{'fraction_non_alphanumeric': 0.038873141667148725, 'fraction_numerical': 0.023042363240199766, 'mean_word_length': 4.584365241735759, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Turbulent generation of magnetic switchbacks in the Alfv\'enic solar wind

One of the most important early results from the Parker Solar Probe (PSP) is the ubiquitous presence of magnetic switchbacks, whose origin is under debate. Using a three-dimensional direct numerical simulation of the equations of compressible magnetohydrodynamics from the corona to 40 solar radii, we investigate whether magnetic switchbacks emerge from granulation-driven Alfv\'en waves and turbulence in the solar wind. The simulated solar wind is an Alfv\'enic slow-solar-wind stream with a radial profile consistent with various observations, including observations from PSP. As a natural consequence of Alfv\'en-wave turbulence, the simulation reproduced magnetic switchbacks with many of the same properties as observed switchbacks, including Alfv\'enic v-b correlation, spherical polarization (low magnetic compressibility), and a volume filling fraction that increases with radial distance. The analysis of propagation speed and scale length shows that the magnetic switchbacks are large-amplitude (nonlinear) Alfv\'en waves with discontinuities in the magnetic field direction. We directly compare our simulation with observations using a virtual flyby of PSP in our simulation domain. We conclude that at least some of the switchbacks observed by PSP are a natural consequence of the growth in amplitude of spherically polarized Alfv\'en waves as they propagate away from the Sun.

The classical idea of the solar wind is based on pressure-driven acceleration that leads to a transonic outflow (Parker 1958(Parker , 1965Velli 1994). Early observation of supersonic plasma velocities in interplanetary space (Neugebauer & Snyder 1962, 1966 supported this idea. This thermally-driven wind model is, however, incomplete in that it cannot reproduce the wellestablished anti-correlation between solar-wind velocity and coronal (freezing-in) temperature (Geiss et al. 1995;von Steiger et al. 2010) nor the large wind velocities measured in fast-solar-wind streams near earth (Durney 1972). In addition, the solar-wind mass flux remains nearly constant regardless of wind speed or solar activity (Goldstein et al. 1996;Wang 1998;Cohen 2011;Cranmer 2017) in contrast to the thermally-driven model that predicts a sensitive dependence of the mass flux to the coronal temperature (Parker 1965;Lamers & Cassinelli 1999;O'Fionnagáin & Vidotto 2018). To explain these observations, as well as the mass and energy budget across the transition region, a self-consistent description of coronal heating and wind acceleration via magnetic field needs to be considered (Hammer 1982;Withbroe 1988;Hansteen & Leer 1995;Hansteen & Velli 2012). Two different types of magnetically driven solar-wind models have been proposed: wave/turbulence-driven (WTD) models and reconnection/loop-opening (RLO) models (see e.g. Cranmer 2009).
Compositional analysis of the solar wind indicates that (a part of) the slow solar wind may have a different origin than the fast solar wind. The first-ionizationpotential (FIP) bias, the degree of relative enhancement of low first-ionization-potential elements, is observed to be large and variable in the slow solar wind (von Steiger et al. 2000;Stansby et al. 2020). A large FIP bias is also observed in closed-loop regions such as helmet streamers (Raymond et al. 1997;Feldman et al. 1998) and active regions (Widing & Feldman 2001;Brooks & Warren 2011;Brooks et al. 2015;Baker et al. 2018;Doschek & Warren 2019), which suggests that slow-solar-wind streams with large FIP bias are formed by the leakage of closed-loop plasma. Active-region outflows (Sakao et al. 2007;Harra et al. 2008;Brooks & Warren 2011;Brooks et al. 2015) or streamer blobs (Sheeley et al. 1997;Wang et al. 1998;Viall & Vourlidas 2015) are possible observed signatures of the leakage of the closedfield plasma. Therefore, magnetic reconnection and the resultant opening of loops are possibly central to origin of the slow solar wind. Models of this category are called reconnection/loop-opening (RLO) models (Fisk et al. 1999;Fisk 2003;Antiochos et al. 2011;Higginson et al. 2017;Réville et al. 2020a;Wang 2020). RLO models may be particularly relevant to active stars, for which reconnection and flares occur more frequently (Candelaresi et al. 2014).
Winds from coronal holes are believed to be driven by waves and turbulence, because the energy released by magnetic reconnection in coronal holes is likely insufficient to accelerate the solar wind (Cranmer & van Ballegooijen 2010;Lionello et al. 2016). This conclusion, however, is challenged by the Parker Solar Probe's recent observations of large numbers of "magnetic switchbacks", which are defined as abrupt, large-angle rotations of the magnetic field that yield local reversals in the radial component of the magnetic field, in solar wind emanating from a low-latitude coronal hole Kasper et al. 2019).
Parker Solar Probe (hereafter PSP) (Fox et al. 2016) is a mission to observe the near-Sun solar wind by measuring electromagnetic fields (FIELDS, Bale et al. 2016), the distribution functions of thermal electrons, alpha particles, and protons (SWEAP, Kasper et al. 2016), high-energy particles (IS IS, McComas et al. 2016) and scattered white light (WISPR, Vourlidas et al. 2016). One of the most important early results from PSP is the ubiquitous presence of sudden local magnetic polarity reversal events called magnetic switchbacks . The presence of switchbacks was first reported in the high-latitude wind (Balogh et al. 1999;Matteini et al. 2014), and later was also found in equatorial wind at 1 au (Gosling et al. 2009) and 0.3 au (Horbury et al. 2018). The detailed observational properties are summarized in Section 5.
Possible scenarios for the origin of switchbacks include magnetic transient events (reconnection and/or jets) in the solar atmosphere (Yamauchi et al. 2004a;Roberts et al. 2018;He et al. 2020;Tenerani et al. 2020;Sterling & Moore 2020;Zank et al. 2020), shear-wave interaction (Landi et al. 2006;Ruffolo et al. 2020;Shi et al. 2020) and large-amplitude turbulence (Squire et al. 2020). If magnetic transient events are the principal source of switchbacks, then this would imply that RLOlike events play a significant role even in the coronal-hole wind, which would pose a challenge to the conventional understanding that waves and turbulence dominate the heating and acceleration in coronal-hole wind.
Motivated by the above background, in this work we aim to investigate the theoretical properties of magnetic switchbacks in a wave/turbulence-driven solarwind model.
For this purpose, we perform a three-dimensional direct numerical simulation of the turbulence-driven solar wind from the coronal base to 40R , thus directly connecting coronal dynamics and solar-wind fluctuations.

Basic equations
We solve the three-dimensional magnetohydrodynamic equations with gravity and heating (e.g. Priest 2014): where ρ is the mass density, p is the pressure, v is the velocity field, g is the gravitational acceleration, B is the magnetic field, e int is the internal energy per unit volume and L is the heat injection per unit volume and time, respectively. The internal energy and pressure are related by where γ = 5/3 is the specific heat ratio of adiabatic gas. Mathematically, Eq.
(3) ensures that ∇ · B = 0 as long as the magnetic field is divergence-free at any given time. In numerically solving the Eq.s (1)-(4), however, numerical error gives rise to non-zero ∇ · B that could result in an unphysical solution. To overcome this issue, a scalar variable ζ is added to the basic equations (Dedner et al. 2002).
Eq.s (1), (2), (4), (6), (7) are the basic equations used in this work. We solve the basic equations in the spherical coordinate system (r, θ, φ), in which the nabla operator is given by where e r,θ,φ denotes the unit vector in r, θ, φ direction. The extension of the simulation domain in the θ direction is only 4.32 × 10 −2 rad, and thus, we can assume, without loss of generality, that θ ≈ π/2 or equivalently sin θ ≈ 1. The ∇ operator in this local spherical coordinate system is given by Using Eq. (9) and after some algebra, the basic equations (1), (2), (4), (6), (7) are given in the form of a conservation law, where U , F r,θ,φ and S are the conserved variables, corresponding fluxes in each direction, and source terms, respectively, which are given by where δ α,β is the Kronecker delta andB = B/ √ 4π. e and p T are the total energy per unit volume and total pressure given by Assuming that hydrogen is ionized throughout the simulation domain, the equation of state for a fully-ionized hydrogen plasma is used to obtain the temperature: where n p and m p are the proton number density and mass, respectively. g is the absolute value of the gravitational acceleration g at radial distance r.
where G is the gravitational constant and M is the solar mass. The energy source term L is inlcuded to mimic the thermal conduction that works to homogenize the temperature profile. To speed up the simulation, we implement exponential cooling to a given reference temperature profile by setting where e int,ref is the internal energy corresponding to the reference temperature T ref , where T i and T o represents the asymptotic behaviors of the reference temperature in the inner and outer sides of the simulation domain, which are given by We set the relaxation time τ cnd = 1 s to ensure that the thermal relaxation is faster than the turbulent evolution. Note that this is comparable to the time scale of thermal conduction in the solar wind (∼ 4 s, assuming a number density of 10 5 cm −3 , temperature of 10 6 K, and temperature scale length of 10 5 km).

Simulation domain and boundary conditions
The simulation domain is a quadrangular pyramid rooted in the coronal base. The range of the simulation domain is defined by 1.02 ≤ r/R ≤ 40, −θ max ≤ θ ≤ θ max , −φ max ≤ φ ≤ φ max , where θ max = φ max = 2.16 × 10 −2 rad. The horizontal size of the simulation domain is 30, 000 km at the coronal base, which is comparable to the scale of super granulation (Leighton et al. 1962). We employ (16000, 216, 216) uniform grid points in the (r, θ, φ) directions.
Periodic boundary conditions are imposed in the θ and φ directions. Because the solar wind is supersonic and super-Alfvénic at r/R = 40, free boundary conditions are imposed at the outer boundary. At the inner boundary, fixed boundary conditions are imposed on the density, radial velocity, temperature and radial magnetic field: ρ inner = 4.0 × 10 −16 g cm −3 , T inner = 1.0 × 10 6 K, v r,inner = 0.0 km s −1 and B r,inner = 1.1 G, where the subscript "inner" denotes the value at the inner boundary. The transverse velocity and magnetic field are given in terms of Elsässer variables z ± θ,φ defined by respectively.
To ensure that the velocity and magnetic field are divergence-free at the inner boundary, the Elsässer variables are given in terms of stream functions.
where z ± θ,φ,inner denote the θ-and φ-components of Elsässer variables at the inner boundary and ψ ± inner are the corresponding stream functions.
The stream function of the upward component ψ + inner is given such that it has broadband spectra in space and time: where k θ and k φ denote the normalized wavenumbers in the θ and φ directions, f is the frequency, and C (k θ , k φ , f ) stands for the amplitude of the mode. Reduced MHD simulations show that the spectral index of upward Elsässer (Alfvén-wave) energy with respect to perpendicular wave number k ⊥ = k 2 θ + k 2 φ is approximately −0.8 at the coronal base, when the perpendicular outer scale of the turbulence in the chromosphere is sufficiently small (see Figure 7 of Chandran & Perez 2019). Considering this result, the injected Alfvén waves are assumed to exhibit a spectral index of −2/3 with respect to k ⊥ and a broken power-law spectrum with respect to f . This is achieved by setting which gives a flat frequency spectrum at f < f inj and a 1/f spectrum at f > f inj . We set n θ = n φ = 27 for the following reason. To represent a broadband spectrum, a large number needs to be chosen for n θ and n φ . Although the maximum value for n θ and n φ (the Nyquist wavenumber) is 108 in our case, care needs to be taken so that the maximum wavenumber mode imposed at the inner boundary does not suffer from immediate numerical dissipation. To avoid the numerical dissipation, any imposed wave modes need to be resolved by at least a few grid points. For this reason, we set the maximum wavenumber of the imposed mode n θ and n φ to be 27, which is resolved by 8 grid points and thus is numerically non-diffusive. f inj corresponds to the energy-injection scale in terms of frequency. We set f inj = 10 −3 Hz, which corresponds to the timescale of the solar granulation (Hirzberger et al. 1999), which is believed to be the primary source of the Alfvén waves launched by the Sun (Steiner et al. 1998).
The frequency bins in Eq. (21) are uniformly distributed in 10 −4 Hz ≤ f ≤ 10 −2 Hz. The amplitude of ψ + is tuned so that the rms amplitude of the velocity fluctuations of the upward Alfvén waves is 30 km s −1 at the lower boundary. At the lower boundary, the stream function of the downward component is assumed to vanish (ψ − inner = 0), which yields a discontinuity in the amplitude of the downward Alfvén waves. However, since the amplitude of downward Alfvén waves is small near the bottom boundary ( 10 km s −1 ), an unphysical discontinuity is unlikely to affect the dynamics of the solar wind.

ALFVÉNIC SLOW SOLAR WIND
In approximately 10 5 s, the system reaches a quasisteady state (QSS) that is not influenced by the initial conditions. Since our interest is the dynamics of turbulence in the quasi-steady solar wind, we analyze the simulation data in the QSS with time duration of τ dur = 1.2 × 10 4 s and time cadence of ∆t = 6 s. The actual time step in the simulation is in the order of 10 6 .
In Table 1, we show the averaged properties of the simulated solar wind at several radial distances. According to the the values shown in Table 1, the sim-ulated solar wind exhibits an averaged mass-loss rate of 2.4-2.7 × 10 −14 M yr −1 (which is in the observed range of 2-3 × 10 −14 M yr −1 , see Figure 1 of Cranmer et al. 2017) and is characterized by its relatively low speed v r 400 km s −1 and high normalized cross helicity σ c 0.95, and thus is categorized as Alfvénic slow solar wind (Marsch et al. 1981;D'Amicis & Bruno 2015;D'Amicis et al. 2019). Since the first perihelion of PSP is dominated by Alfvénic slow solar wind , it is worthwhile to compare our simulation results with the encounter-1 observations.

Overview on rθ-plane
To see the radial evolution of the solar wind and turbulence therein, visualization of the vertical (rθ) plane is useful. Figure 1 shows a snapshot of the QSS in the rθ plane. Left and right panels visualize the whole simulation domain (1 ≤ r/R ≤ 40) and the near-Sun region (1 ≤ r/R ≤ 10), respectively. The panels that display the outward and inward Elsässer variables show the complex nature of turbulent transport in the solar wind. Comparing the spatial structure of z + φ and z − φ at 1 < r/R < 10, one finds finer structures in z − φ . This structure difference was previously found in a numerical simulation of fast solar wind (Shoda et al. 2019). Farther from the sun where the solar wind is super Alfvénic, the chaotic nature of z − φ disappears, and instead large-scale structure with respect to θ is evident.
Local reversals of the radial magnetic field are highlighted by arrows in the map of B r /B. In the near-Sun solar wind, the radial magnetic field is found to exhibit small fluctuations. Beyond r ∼ 25R , B r occasionally locally reverses to form magnetic switchbacks ). Figure 1 also shows that large fluctuations in B r are associated with local enhancements in v r , which are also highlighted by arrows. These switchbacks are discussed in more detail in Section 5.

Radial profiles of the mean field
To analyze the large-scale properties of the simulated solar wind, we average the plasma properties and magnetic field over the θφ plane and time. For any given variable X, its horizontal average is defined by while its time average is by where t = 0 is the time when the system reached QSS. The averaging time τ dur is set to 1.2×10 4 s. The average with respect to time and horizontal space is denoted by The root-mean-square fluctuation amplitude of X, δX rms , is defined by The black dashed lines in Figure 2 show the radial profiles of the averaged density ρ ave , radial velocity v r,ave , and radial magnetic field B r,ave . Also shown in the grey solid lines are the corresponding non-averaged (snapshot) profiles at the center of the simulation domain (θ = φ = 0). Observations from SOHO/LASCO (Brueckner et al. 1995), SOHO/UVCS (Kohl et al. 1995, MLSO/MK-III, STEREO-A/SECCHI/COR2 (Howard et al. 2008), Akatsuki (Nakamura et al. 2011), PSP/FIELDS , PSP/SWEAP (Kasper et al. 2016) are also plotted (see caption for corresponding papers). All PSP data are retrieved from Kasper et al. (2019) as the most probable values at given radial distance in PSP's encounter 1. In converting the proton/electron number density to the mass density, we assumed that the helium-to-proton ratio is 5% in number. We note that, although the typical slow solar wind is found to exhibit smaller helium abundance (∼ 2.5%, Kasper et al. 2007), highly Alfvénic slow solar wind can be helium-rich with abundance as large as 5% (Huang et al. 2020).
Although there are broad similarities between our simulation and the observations, there are also several differences. In the density plot, the measurements from DeForest et al. (2018) systematically lie above our model possibly because they estimated the density of blobs while our model shows the density of the ambient solar wind. A systematic gap also exists between our model and the velocity measurements from Miralles et al. (2001). A possible reason for this discrepancy is that they reported the outflow velocity of O VI which is often found to be systematically larger than the bulk speed of the solar wind (Kohl et al. 1998).
The radial magnetic field on average scales as B r ∝ r −2 . This is because the mean field is forced to be aligned with the simulation domain that expands radially in the spherical coordinate system. In the snapshot data (grey line), one finds multiple sudden decreases of B r (magnetic switchbacks) that appear at r/R 10.

Radial profiles of the fluctuation amplitudes
Since the energy for the wind heating and acceleration lies in the turbulent fluctuations, how the amplitude of the turbulence varies in r is worth investigating. Figure  3 shows the radial profiles of the amplitudes of the fluctuations in the Elsässer fields (z ± = v ∓ B/ √ 4πρ ) and magnetic field. In the top panel, solid lines show the rms amplitudes of outward and inward Elsässer fields, and the dashed lines show the rms amplitudes of the transverse components of the Elsässer fields: z ± ⊥ = z ± θ e θ +z ± φ e φ . Since the fluctuations are nearly transverse, the solid and dashed lines nearly overlap. Crosses are PSP encounter-1 observations from Chen et al. (2020). Also shown by the dotted lines are asymptotic radial scalings of z ± from Chen et al. (2020) (see also Bavassano et al. 2000). High Alfvénicity is maintained throughout the simulation domain, possibly because the weaker Elsässer field decays more rapidly (Dobrowolny et al. 1980;Pouquet et al. 1986). The bottom panel shows the profiles of the radial and transverse magnetic-field fluctuations. The radial-field fluctuation is approximately one order of magnitude smaller than the transverse fluctuation, which is con- sistent with the transverse nature of the Alfvén wave. According to Vasquez & Hollweg (1998), the interaction between two oblique Alfvén waves with the same group velocity drives a second-order parallel-field fluctuation given by consistent with earlier work by Barnes & Hollweg (1974) on spherically polarized (or constant-|B|) Alfvén waves. The red dashed line in the bottom panel of Figure 3 shows a root-mean-square version of Eq. (27), which nicely overlaps the black dashed line, supporting the idea that the radial-field fluctuation is driven nonlinearly by oblique Alfvén waves. Eq. (27) predicts that the radial-field fluctuation can be as large as the mean radial field when the Alfvén-wave fluctuation is comparable to the mean field. In contrast to most reduced-MHD simulations (e.g.  Shoda et al. (2019) showed that the PDI plays an essential role in driving turbulence in the fast solar wind. However, it is unclear whether the PDI is important also in the Alfvénic slow solar wind. Here, we discuss the origin and role of density fluctuations, following an analysis similar to that of Shoda et al. (2019). Figure 4 shows the radial profile of the fractional density fluctuation δρ rms /ρ ave in our model (black solid lines). The growth rate of the parametric decay instability γ PDI is also shown in the top panel (dashed line). In calculating γ PDI we numerically solved the Goldstein-Derby dispersion relation (Goldstein 1978;Derby 1978) with frequency f = f inj = 10 −3 Hz and added the suppression terms from wind acceleration and expansion (Tenerani & Velli 2013;Shoda et al. 2018a). Comparing the growth rate of the PDI and the magnitude of density fluctuations, one finds that the density fluctuations grow rapidly where the growth rate of the PDI is large. This spatial correlation supports the idea of densityfluctuation generation by the PDI, consistent with what is found in fast-wind simulations (Shoda et al. 2018a(Shoda et al. , 2019. The bottom panel of Figure 4 compares the rootmean-square amplitudes of the simulated and observed density fluctuations. The observed values are from Akatuski (Nakamura et al. 2011), PROBA2/SWAP (Berghmans et al. 2006;Seaton et al. 2013) and PSP/FIELDS ) (see caption for the corresponding papers). Radio scintillation observations from Miyamoto et al. (2014) exhibit a large scatter. To take into account the underestimation of the local density-fluctuation amplitude by the positive-negative cancellation along the line of sight, we focus on the upper envelope of the measurements (the blue dashed line) rather than each value (blue circles). The upper envelope is globally consistent with our model. EUV observations (green diamonds) are systematically higher possibly because we ignore the slow magnetosonic waves from below the transition region (DeForest & Gurman 1998;Ofman et al. 1999;Kiddie et al. 2012) and also the large cross-field variation in the background density along different magnetic flux tubes that has been inferred from Comet-Lovejoy observations (Raymond et al. 2014).
Type-III radio burst observations (orange squares) show a similar magnitude but a somewhat different trend at r/R < 5. In spite of these differences, overall, the observations are broadly consistent with our model. We note, however, that MHD simulations may overestimate the density fluctuations because they neglect collisionless damping (Barnes 1966; but see Schekhochihin et al 2016, Meyrand et al 2019), and thus care needs to be taken in the interpretation. We also note that fastmode waves (Cranmer & van Ballegooijen 2012) and corotating interaction regions (Cranmer et al. 2013) may also contribute to the observed density fluctuations.

DIRECT COMPARISON WITH PSP DATA
Since our simulation domain extends beyond the perihelions of Parker Solar Probe, a direct comparison between our simulation and PSP observations is possible. Here we compare our simulation results with PSP encounter-1 data.
Around the first perihelion, PSP was nearly corotating with the Sun. Although the effect of rotation is not included in our simulation, we simply assume that the simulation domain is co-rotating with the Sun. The location of PSP in the simulation domain is therefore assumed to be fixed in time. For any variable X, the simulated time series of encounter-1 X sim (t) is calculated using the equation X sim (t) = X (r, θ, φ, t)| r=35.5R , θ=−θmax, φ=−φmax .
Observed encounter-1 data of PSP are retrieved from the Coordinated Data Analysis Web (CDAWeb). Magnetic fields in RTN coordinates are taken from the Level-2 data of the fluxgate magnetometer (MAG), part of the FIELDS instrument suite. The proton density, bulk velocity and thermal velocity are taken from the Level-3 data of the Solar Probe Cup (SPC), part of the SWEAP instrument suite, as 0th-, 1st-and 2nd-moments of the reduced distribution function.
Two sets of two-hundred minutes of data are compared with our simulation: the interval from 2018/11/05 12:40:00 to 2018/11/05 16:00:00 and the interval from 2018/11/06 06:20:00 to 2018/11/06 09:40:00. The former period corresponds to a relatively active phase that exhibits large fluctuations in B r , and the latter to a relatively quiet phase that exhibits small fluctuations in B r . Table 2 shows the quantitative difference in the magnetic-field properties in these two periods and the simulated data. The data from Nov 5th are characterized by large δB r,rms / |δB rms | with little 'variance anisotropy': δB r,rms / |δB rms | ≈ 0.577 and δB r ≈ δB t ≈ δB n . On the other hand, both the data from Nov 6th and the data from our simulation exhibit much smaller values of δB r,rms / |δB rms |. Since our simulation captures neither high-frequency fluctuations originating from kinetic physics nor sub-grid-scale fluctuations, we take a 24-second time average of the PSP data before comparing with our simulation results. Figure 5 shows a direct comparison between observed (left and middle column) and simulated data (right column) of PSP. In this comparison, the sign of the simulated magnetic field is changed (B sim → −B sim ) so that mean magnetic-field polarity is consistent with observations: during this period PSP observed negativepolarity-dominated wind while our simulation assumes a positive-polarity-dominated wind. A change in the sign of B does not affect the comparison because the structure of field line is invariant. In converting bulk velocity v to proton bulk velocity v p , we assume no differential flow between different species, thus v = v p . The proton number density is obtained by a simple assumption that the solar wind is composed of fully ionized hydrogen plasma.
Considering the lack of artificial fine tuning in the simulation, the similarity between the Nov-6th observations and simulation results is surprising. Various observational properties are also seen in the simulated time series: approximately constant density and field strength, the one-sided nature of the fluctuations in v r and B r , turbulent fluctuations in v and B, and the Alfvénic correlation between B and v shown in the bottom two panels of Figure 5.  Table 2. This table shows magnetic-field properties in the two observational periods in the encounter 1 and the corresponding simulated data displayed in Figure 5. Note that the averaged and root-mean-squared values are defined in terms of time average over 200 min. Meanwhile, the Nov-5th observation exhibits several substantial differences from the simulation data. Mag-netic switchbacks (where B r > 0, see Section 5 for definition) that last more than a few minutes are observed multiple times, which are not seen in the simulation data. As a result of the emergence of switchbacks, the magnetic field experiences more discontinuous jumps than the simulation predicts. However, as noted above and shown in Table 2, δB r,rms /|δB rms | is smaller in the simulation than in the Nov-5 data by a factor of more than 2. Given that the fractional magnetic-field fluctuation is significantly larger in the Nov-5 data, it is perhaps not surprising that switchback features are more prominent than in the simulation. However, it remains to be seen whether a numerical simulation with larger δB r,rms /|δB rms | (and perhaps higher numerical resolution, as discussed above) could explain the data, or whether additional physical ingredients (such as impulsive forcing at the coronal base) are needed. It remains to be seen whether a numerical simulation with higher resolution and/or amplitude could could explain the data, or whether additional physical ingredients (such as impulsive forcing at the coronal base) are needed.
In Figure 6, we compare a single switchback observed in our simulation at r/R = 35 with an individual switchback observed by PSP on November 6, 2018 (see Zank et al. 2020, for a similar comparison). For the PSP measurements, we show both the full data and the data averaged over 24 s, the latter of which is to remove the high-frequency fluctuations beyond the capability of our model. The simulated and observed switchbacks exhibit similar characteristics, including similar durations and amplitudes, consistent with the hypothesis that the switchbacks seen by PSP emerge naturally from the dynamics of solar-wind turbulence, even without impulsive forcing near the Sun. The detailed statistical properties of switchbacks are discussed in the next section.

MAGNETIC SWITCHBACKS
Observationally, magnetic switchbacks are defined as regions in which the radial field is reversed against the ambient region. Since the mean radial magnetic field is positive in our simulation, switchbacks are defined as regions in which B r < 0, and switchback boundaries are defined as surfaces on which B r = 0. In the previous section, we have already shown the presence of switchbacks in our simulation.
In this section, after reviewing the observational properties of magnetic switchbacks, we present several detailed analyses of switchbacks, focusing on their physical properties and the comparison with observations.

Observational properties
Before presenting our analysis of the switchbacks in our simulations, we first summarize the main observa- tional characteristics of the magnetic switchbacks observed by PSP: 1. Switchbacks are Alfvénic in the sense that the jump in the magnetic field ∆B is associated with a jump in the velocity ∆v given by ±∆B/ √ 4πρ, where ρ is the plasma density, and the +/− sign corresponds to Alfvén waves propagating antiparallel/parallel to the background magnetic field in the frame co-moving with the bulk solar wind. The particular value of this +/-sign for each switchback is opposite to the sign of B r , so that switchbacks always propagate away from the Sun in the plasma frame. As a consequence, the switchbacks are associated with enhanced values of the radial component of the plasma velocity (Kasper et al. 2019;Bale et al. 2019;Horbury et al. 2020). We note that switchbacks are sometimes 2. Switchbacks are approximately "spherically polarized," in the sense that the strength of the magnetic field is almost constant as the field abruptly rotates Farrell et al. 2020).
3. Switchbacks appear to be field-aligned, elongated structures with a typical aspect ratio of 10 near the first perihelion Laker et al. 2020), where the aspect ratio is defined as the ratio of the length scales parallel and perpendicular to the magnetic field. Wit et al. 2020). As a result, although the filling factor of B r -reversed regions is 6% during PSP's first perihelion encounter , the "active phase" of the solar wind that contains switchbacks occupies at least 75% of the total observation time near the first perihelion .

A waiting-time analysis indicates that switchbacks appear in clusters (Dudok de
5. Switchbacks are S-shaped, locally folded magneticfield lines. The propagation directions of the electron strahl (Kasper et al. 2019;Whittlesey et al. 2020), small-scale Alfvén waves (McManus et al. 2020), and proton-alpha differential velocity (Yamauchi et al. 2004b) support this picture.
6. The proton core parallel temperature, the fieldaligned temperature obtained from the core of the proton distribution function, is found to be constant across switchbacks , possibly because switchbacks move through the plasma at the local Alfvén speed, which prevents them from strongly heating any particular parcel of plasma before propagating past it. We note, however, that the active phase mentioned above is characterized by an enhanced proton parallel temperature (Woodham et al. 2020).
7. The volume filling factor of switchbacks is observed to increase as the heliospheric distance increases from ∼ 35R to ∼ 50R (Mozer et al. 2020). In light of the comparative rarity of switchbacks near 1 au, the volume filling factor should begin to decrease somewhere between 50R and 1 au.
The theory of magnetic switchbacks should explain these observational properties.

Structure in the θφ plane
To see the localized nature of magnetic switchbacks, we show in Figure 7 the simulation data in the θφ plane at r/R = 35. The four panels show the normalized radial magnetic field B r /B, density fluctuation δρ/ρ = ρ/ρ − 1, radial Poynting-flux fluctuation δS r /S r = S r /S r − 1, and radial velocity v r . The white solid lines in the top-right and lower-left panels and the black solid lines in the lower-right panel correspond to the boundaries of magnetic switchbacks (where B r = 0). Note that the overline denotes an average in the θφ plane.
A careful observation of the θφ-plane structure reveals several physical properties of the switchbacks: 1. As seen in the top-left panel of Figure 7, in which the switchback (blue patch) is much smaller than the rθ plane, the switchbacks are highly-localized structures in that their horizontal spatial extent is much smaller than the energy-containing scale of the turbulence, which is comparable to the size of simulation domain in the θ and φ directions. any correlation between the emergence of switchbacks and density fluctuations, since no clear enhancement in the density fluctuation is seen in the switchback region of the bottom-right panel.
These characteristics are more quantitatively discussed in the following sections. The structure of a magnetic switchback can be clearly seen in a plot of its one-dimensional spatial profile. Figure 8 shows an example of the one-dimensional structure of a switchback along the white dashed and dotted lines in the left-top panel of Figure 7. An interesting point is that the total magnetic field is nearly constant across the switchback in both the θ and φ directions, which is inferred from the B profile (black solid line) in Figure 8. To investigate whether switchbacks are closer to rotational discontinuities (RD) or tangential discontinuities (TD) (see, e.g., Horbury et al. 2001), we define the field-aligned and field-normal length scales (Squire et al. 2020) whereb = B/B andn is a unit vector normal tob (n·b = 0). Takingn to be in the θφ plane, at each given radial distance, lb and ln are calculated at the boundaries of switchbacks (where B r = 0) and then averaged over at each radial distance. To be specific, the total number of detected boundaries is 192 at r/R = 10 and 155002 at r/R = 40, and lb and ln values calculated at these boundary grids are averaged for each θφ plane. The quantities lb and ln represent the length scales of the field deflection in the parallel and perpendicular directions. If the boundary of a switchback is TD-like, the maximum-variation direction is nearly perpendicular, which yields lb > ln. The radial variations of averaged lb and ln are shown in Figure 9. The field-aligned length scale lb is found to be statistically larger than the field-normal length scale ln. We thus conclude that switchbacks are on average more TD-like than RD-like, which is consistent with the observational result that the discontinuities in the solar wind are mostly TD-like (Horbury et al. 2001).

Alfvénic nature
Important physical properties of magnetic switchbacks can be inferred from the correlations between their magnetic fluctuations and other quantities. One of the important correlations is found between the radial-field fluctuation δB r and the radial-velocity fluctuation δv r , which are connected by the following Alfvénic relation: where B r,ave and v Ar,ave are the averaged radial magnetic field and radial Alfvén speed v Ar = B r / √ 4πρ. In addition to the radial velocity enhancement, some of the observed magnetic switchbacks are associated with density and Poynting-flux enhancements ).
Here we investigate whether these correlations are found in our simulation. Figure 10 shows scatter plots of normalized B r fluctuations (−δB r /B r,ave ) and normalized radial-velocity enhancements δv r /v Ar,ave (top panels), normalized density fluctuations δρ/ρ (middle panels), and normalized Poynting-flux fluctuations δS r /S r (bottom panels) at different radial distances, with the color scale indicating the probability of finding each combination of variables in the simulation domain. The radial component of the Poynting flux S r is defined by The grey lines in the top panels show the Alfvénic relation, Eq. (30). Note that the data shown in Figure 10 is restricted to switchback regions in which δB r + B r,ave < 0. Several features are found in these scatter plots: 1. Fluctuations in the radial magnetic field are always associated with fluctuations in the radial velocity. All the switchback events approximately satisfy the Alfvénic relation, Eq. (30), indicating the Alfvénic nature of magnetic switchbacks.
2. No clear correlations are found between density fluctuations and switchbacks, in that the majority of switchback events lie near δρ = 0, regardless of the switchback amplitude −δB r /B r,ave .
3. Magnetic switchbacks always exhibit larger-thanaverage Poynting flux. At the same time, within the switchback population, the switchback amplitude is not correlated with the magnitude of Figure 10. Scatter plots of the normalized radial-field fluctuation −δBr/Br,ave and radial-velocity enhancement (δvr/vAr,ave, top panels), density fluctuation (δρ/ρave, middle panels), and Poynting-flux fluctuation (δSr/Sr,ave, bottom panels). Left, center and right panels correspond to r/R = 10, 20, and 40, respectively. The color represents the normalized kernel density estimate of the points. Grey lines in the top panels show the Alfvénic-correlation relation δBr/Br,ave = −δvr/vAr,ave.
the Poynting-flux fluctuation. This means that switchbacks with arbitrary amplitude can emerge once the radial Poynting flux exceeds a critical threshold value that depends on heliocentric distance.
These results indicate that the switchbacks in our simulation are associated with large-amplitude, unidirectional Alfvén waves that exhibit strong v-B correlation, weak density fluctuations, and larger-thanaverage radial Poynting flux.

Magnetic compressibility
The constant-B nature of switchbacks found in Figure  8 is more quantitatively seen from the magnetic com-pressibility, which is given by where B = |B| and the overline denotes a spatial average. Black and red symbols in Figure 11 show the radial evolution of magnetic compressibility in the whole θφ plane and near switchbacks (in the square region the maximum and minimum values of C B 2 , respectively, and the points represent the mean values. The magnetic compressibility defined in the whole θφ-plane is much smaller than unity which is consistent with observations in Alfvénic solar wind. Interestingly, the magnetic compressibility near switchbacks is similar to or even smaller than that defined in the whole θφ plane. In terms of B, this analysis shows that the switchbacks are approximately "spherically polarized," in the sense of having nearly constant B, consistent with many of the switchbacks observed by PSP. Some of the observed switchbacks, on the other hand, have a compressional nature. If we define "compressional events" as having C B 2 > 0.02, then ∼ 0.1-1% of switchbacks are compressional, which is much smaller than the observed fraction (∼ 27%, see Larosa et al. 2020). Our results suggest that the explanation of such compressional switchbacks lies beyond the mechanisms seen in the simulations presented here and requires additional physics.

Aspect ratio
PSP measurements suggest that switchbacks are elongated along the background magnetic field, with a typical aspect ratio of 10 Laker et al. 2020). It is worth investigating whether this large aspect ratio is reproduced in our simulation. We used the following method to measure the aspect ratio of mag- netic switchbacks. First, at each time step, we find the minimum-B r grid point in a given rθ plane. When the minimum B r is negative, we define the minimum-B r point as a switchback region. For each point in a switchback region, the adjacent grid points that satisfy B r < 0 are iteratively found and collectively defined as the switchback region. In this way, we obtain a certain connected domain that satisfies B r < 0 everywhere. The ratio between the radial extent l SB and transverse extent l SB ⊥ of the domain is defined as the aspect ratio: l SB /l SB ⊥ . We note that multiple switchbacks can be found at given φ and t, among which we only focus on the one with minimum B r . Figure 12 shows the measured aspect ratio as a function of radial distance. Although the data are highly scattered, on average the aspect ratio is larger than unity, indicating elongated structure along the meanfield direction (r axis). The aspect ratio also seems to decrease with r. However, we need to note that the value of aspect ratio is highly influenced by numerical resolution. Most of the switchbacks are resolved by a few grid points in the perpendicular (θ and φ) directions and thus are broadened by numerical dissipation. That the aspect ratio is smaller than the typical observed value (∼ 10, see Horbury et al. 2020;Laker et al. 2020) is possibly due to insufficient resolution. The actual value of the aspect ratio should be investigated in future highresolution simulations. 5.6. Filling factor and amplitude histogram Figure 13. Filling factor of field deflection events versus radial distance. Blue, green and red symbols correspond to the threshold angles of θBR = 60 • , θBR = 75 • , and θBR = 90 • , respectively.
The number of field-deflection events, the events with large fluctuations in B r including switchbacks, is found to increase with distance from the Sun regardless of the deflection angle (see Figure 3 of Mozer et al. 2020) This trend is observed in our simulation data. To show this, we define the filling factor of field-deflection events via the following procedure. We first count the number of grid points that satisfy the criteria θ BR > 60 • , θ BR > 75 • , and θ BR > 90 • , where θ BR is the angle between the radial direction and the local magnetic field, in each θφ plane over the 2000 time steps we analyze. The filling factor of each field-deflection event is then defined as a ratio of the cumulative counts to the total number of θφ-plane grid points over 2000 steps (216 × 216 × 2000 = 93, 312, 000). Figure 13 shows the radial evolution of the filling factors of field-deflection events with θ BR > 60 • (blue), θ BR > 75 • (green), and θ BR > 90 • (red). The deflection-angle thresholds (60 • , 75 • and 90 • ) are chosen arbitrarily to investigate if the behavior depends on the threshold angle. Regardless of the deflection-angle threshold, the filling factor of field-deflection events increases with radial distance. We note, however, that the filling factor of magnetic switchbacks (θ BR > 90 • ) is much smaller than the observed value at r ≈ 35R (∼ 6%, see Bale et al. 2019). Specifically, the filling factor in our simulation is 0.063% at r/R = 35.5 (close to the first perihelion of PSP), and thus, there exists a two-orders-of-magnitude gap between the observation and the current simulation. The smaller filling factor of switchbacks in our simulation may be due to some combination of a smaller turbulence amplitude and insufficient resolution. In our simulation, |δB rms /B r,ave | = 0.695, whereas |δB rms /B r,ave | ranged from 0.670 to 1.280 in the two intervals from PSP's first perihelion encounter listed in Table 2. As shown in Figure 3 of Squire et al. (2020), the switchback volume filling fraction is highly sensitive to the value of |δB rms /B r,ave |. The presence of intervals during PSP's first perihelion encounter with larger turbulence amplitudes than our simulation may thus at least partially account for the difference in switchback filling fractions. In addition, Squire et al. (2020) showed that the number of switchbacks is sensitive to the numerical resolution. In their expanding-box simulations, which are designed to emulate the evolution of Alfvénic turbulence out to r = 35.7R , the switchback filling fraction reaches over 3% in their highest-resolution run, which has 540 × 1120 2 grid points. In our simulations, if we reduce the resolution of the θφ plane from (216, 216) to (96, 96), the filling factor of magnetic switchbacks is reduced by a factor of 10. Thus, it is possible that the filling factor of magnetic switchbacks would increase substantially in a higher-resolution version of the simulation we have presented.
The distribution of switchback amplitudes is also of interest. Figure 14 shows histograms of detected switchback events as functions of normalized switchback amplitude δB r /B r . Here a grid point with negative B r is counted as a switchback event. As shown by this histogram, as well as the number of events, the maximum  Figure 15. Detected propagation speed of magnetic switchbacks (vSB, red crosses). Black lines are averaged radial velocity (vr,ave, solid line) and radial velocity plus Alfvén speed (vr,ave + vA,ave, dashed line).
amplitudes of switchbacks are found to increase with radial distance.

Propagation speed of magnetic switchbacks
Since the simulation box extends globally from the coronal base to 40R , once a magnetic switchback emerges, we are able to measure its propagation velocity by tracing it. Specifically, we measure how the grid position with minimum B r moves in time. Although this rough estimation yields only a limited number of switchback detections, the number of detected events is sufficiently large to discuss their statistical properties. Figure 15 compares the detected switchback velocity v SB with the averaged radial velocity v r,ave (black solid line) and radial velocity plus Alfvén speed v r,ave + v A,ave (black dashed line). Clearly the radial profile of v SB matches v r,ave + v A,ave , supporting the idea that magnetic switchbacks are locally bent field lines that propagate through the plasma in the anti-Sunward direction at the local Alfvén speed.

SUMMARY AND DISCUSSION
In this work, we have performed a direct numerical simulation of the wave/turbulence-driven solar wind. Our simulation yields Alfvénic slow solar wind, the same type of solar wind observed in the first perihelion encounter of PSP. The radial profiles of the density, velocity, and magnetic field are in agreement with remote-sensing observations and in-situ measurements ( Figure  2).
The turbulence in our simulations is generated by launching outward-propagating Alfvén waves into our simulation domain through the inner (coronal-base) boundary. These waves become turbulent as they propagate away from the sun and ultimately form magnetic switchbacks. A possible scenario for switchback generation inferred from our simulation is summarized as follows. 2. Outward-propagating Alfvén waves nonlinearly evolve towards a constant-B or "spherically polarized" state ( Figure 11, see Cohen & Kulsrud 1974;Barnes & Hollweg 1974;Vasquez & Hollweg 1998).
3. Once δB/B increases to values ∼ 1, the ongoing nonlinear drive towards spherical polarization generates discontinuities in the field, because it is (or at least appears to be) impossible to form a 3D state with constant-B and δB/B ∼ 1 (Figure 8, see Barnes 1976;Valentini et al. 2019).
The switchbacks in the simulation have several of the properties exhibited by observed switchbacks, including spherical polarization (Figures 8, 11), Alfvénic v-b correlation (Figure 10), and a volume filling factor that increases with r ( Figure 13). Switchbacks are also found to exhibit a rapid change in field direction (Figure 8), to be more 'TD-like' than 'RD-like' (Figure 9), and to propagate outward at the Alfvén speed in the frame of background flow (Figure 15). From the simulation results and the discussion above, we conclude that at least some of the magnetic switchbacks seen in the solar wind arise as a natural consequence of Alfvén waves and turbulence with sufficiently large amplitude. However, the volume filling fraction of switchbacks in our simulation is almost two orders of magnitude smaller than in PSP's first perihelion encounter with the Sun, possibly because of insufficient numerical resolution and/or the fact that |δB rms /B r,ave | is smaller in our simulation than in the PSP data (see Table 2). Determining the volume filling fraction in simulations with higher-resolution and different values of |δB rms /B r,ave | is an important goal for future research in this area.
In addition to its importance for determining the switchback volume filling factor, running simulations with larger numerical resolution will be important for describing the fine-scale structure of switchbacks, including the local reversal of the cross helicity (McManus et al. 2020). Besides increasing the numerical resolution, there are a number of other ways that the simulation presented here could be further improved. Since an artificial energy source term is used, our model does not accurately capture the thermodynamic properties of switchbacks Woodham et al. 2020). Field-aligned thermal conduction needs to be considered in the future, ideally incorporating the deviation from the Spitzer-Härm heat flux in the weakly collisional, large-Knudsen-number regime (Salem et al. 2003;Bale et al. 2013;Verscharen et al. 2019). Also, kinetic effects are ignored in our MHD treatment. Considering the scale gap between the typical switchback size and the proton gyroradius, the MHD approximation is likely sufficient for describing the dynamics of switchbacks. However, the temperature is anisotropic in the solar wind , which can affect Alfvén-wave dynamics (Squire et al. 2016;Tenerani & Velli 2018). A fluid approximation with temperature anisotropy (Chandran et al. 2011;Hirabayashi et al. 2016;Tenerani & Velli 2018) is a possible future generalization of the current simulation.
The authors thank Marco Velli, Kosuke Namekata, Shinsuke Takasao, Jono Squire, Romain Meyrand, Al-fred Mallet, Stuart Bale, Justin Kasper, Tim Horbury, and Gary Zank for valuable discussions and comments. The authors are also grateful to the anonymous referee for valuable comments. We acknowledge the NASA Parker Solar Probe Mission and the FIELDS team led by S. Bale and SWEAP team led by J. Kasper for use of data. Parker Solar Probe was designed, built, and is now operated by the Johns Hopkins Applied Physics Laboratory as part of NASA's Living with a Star (LWS) program (contract NNN06AA01C). Support from the LWS management and technical team has played a critical role in the success of the Parker Solar Probe mission. Numerical computations were carried out on the Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. MS is supported by a Grant-in-Aid for Japan Society for the Promotion of Science (JSPS) Fellows and by the NINS program for cross-disciplinary study (grant Nos. 01321802 and 01311904) on Turbulence, Transport, and Heating Dynamics in Laboratory and Solar/ Astrophysical Plasmas: "SoLaBo-X." BC is supported in part by NASA grant NNN06AA01C to the Parker Solar Probe FIELDS Experiment and by NASA grants NNX17AI18G and 80NSSC19K0829.