Analytic Model and Magnetohydrodynamic Simulations of Three-dimensional Magnetic Switchbacks

Parker Solar Probe observations reveal that the near-Sun space is almost filled with magnetic switchbacks (“switchbacks” hereinafter), which may be a major contributor to the heating and acceleration of solar wind. Here, for the first time, we develop an analytic model of an axisymmetric switchback with uniform magnetic field strength. In this model, three parameters control the geometry of the switchback: height (length along the background magnetic field), width (thickness along radial direction perpendicular to the background field), and the radial distance from the center of switchback to the central axis, which is a proxy of the size of the switchback along the third dimension. We carry out 3D magnetohydrodynamic simulations to investigate the dynamic evolution of the switchback. Comparing simulations conducted with compressible and incompressible codes, we verify that compressibility, i.e., parametric decay instability, is necessary for destabilizing the switchback. Our simulations also reveal that the geometry of the switchback significantly affects how fast the switchback destabilizes. The most stable switchbacks are 2D-like (planar) structures with large aspect ratios (length to width), consistent with the observations. We show that when plasma beta (β) is smaller than one, the switchback is more stable as β increases. However, when β is greater than 1, the switchback becomes very unstable as the pattern of the growing compressive fluctuations changes. Our results may explain some of the observational features of switchbacks, including the large aspect ratios and nearly constant occurrence rates in the inner heliosphere.


INTRODUCTION
One of the most striking findings made by Parker Solar Probe (PSP) is that the nascent solar wind is almost filled with magnetic switchbacks ("switchbacks" hereinafter) (Kasper et al. 2019;Bale et al. 2019).Switchbacks are local polarity reversals of the radial magnetic field and are typically Alfvénic structures, with highly correlated velocity and magnetic field fluctuations and nearly constant magnetic field strength (McManus et al. 2020;Woolley et al. 2020;Larosa et al. 2021).For a comprehensive description of switchbacks and their properties, please refer to the review paper (Raouafi et al. 2023a) and references therein.
PSP observations have revealed that switchbacks significantly modify the properties of solar wind turbulence, highlighting their influential role in the solar wind dynamics (Dudok de Wit et al. 2020;Bourouaine et al. 2020;Martinović et al. 2021;Shi et al. 2022b).Besides, as large-amplitude fluctuations, switchbacks are believed to be an important energy source for the heating and acceleration of the solar wind (Halekas et al. 2023).Indeed, observations show that the plasma properties, including proton temperature and temperature anisotropy (Farrell et al. 2020;Larosa et al. 2021;Woolley et al. 2021;Woodham et al. 2021;Luo et al. 2023;Huang et al. 2023;Laker et al. 2024), are different inside and outside the switchbacks, with a possible enhancement of proton temperature inside the switchbacks.

Shi et al.
Although there is no doubt that understanding switchbacks is crucial for a complete understanding of the solar wind, how are the switchbacks generated and how do they propagate in the solar wind is still under debate.There are three widely accepted theories for the generation mechanisms of switchbacks, namely the interchange reconnection happening in the solar corona (Yamauchi et al. 2002;Fisk & Kasper 2020;Zank et al. 2020;Drake et al. 2021;He et al. 2021;Telloni et al. 2022;Upendran & Tripathi 2022;Hou et al. 2023), the velocity shear (Landi et al. 2006;Toth et al. 2023) that may arise due to coronal jets (Sterling & Moore 2020;Magyar et al. 2021;Raouafi et al. 2023b) or motion of the magnetic footpoint between source regions of fast and slow streams (Schwadron & McComas 2021), and the natural evolution of Alfvén waves in the expanding solar wind (Belcher 1971;Hollweg 1974;Squire et al. 2020;Shoda et al. 2021;Mallet et al. 2021).However, which mechanism is dominant is still under debate.As different generation mechanisms take place at different radial distances to the Sun, it is thus important to study the radial evolution of the occurrence rate of switchbacks that may suggest which generation mechanism is most effective.Unfortunately, no solid conclusion has been made on this topic either.Mozer et al. (2020) suggest that the switchback occurrence rate does not change much with radial distance, while other studies indicate a dependence of the occurrence rate on the duration of the switchback, i.e. the occurrence rate of larger switchbacks increases with the radial distance and the occurrence rate of shorter-duration switchbacks shows a decreasing trend (Tenerani et al. 2021;Jagarlamudi et al. 2023).Moreover, there is evidence showing that the occurrence rate of switchbacks is smaller inside the Alfvén radius than outside (Pecora et al. 2022;Bandyopadhyay et al. 2022).Hence, it is likely that a subset of the switchbacks observed in-situ are generated in the lower corona while the others are generated locally (Tenerani et al. 2021).Another important question is how does a switchback evolve in the solar wind, because the dynamic evolution of switchbacks inevitably affects the characteristics of switchbacks observed in-situ, including their occurrence rates as well as their topology (Laker et al. 2021).Magnetohydrodynamic (MHD) simulations (Tenerani et al. 2020) show that a 2D Alfvénic switchback maintains stable for quite a long time (tens of Alfvén crossing time) and is eventually destabilized by the growth of density fluctuations, i.e. parametric decay instability (PDI).In addition, the Parker spiral (Johnston et al. 2022;Squire et al. 2022) and Hall effect (Tenerani et al. 2023) also modify the evolution and dissipation of the switchbacks.
The objective of this study is to numerically investigate the dynamic evolution of 3D switchbacks.The 3D effect has not been fully investigated in previous studies because of the difficulty to establish a 3D analytic model that can satisfy both the divergence-free condition and a uniform magnetic field strength.Squire & Mallet (2022) developed a timeintegration method that grows the amplitude of the initially weak spherically-polarized magnetic field fluctuations (with normalized magnetic field fluctuation δB/|B| ≈ 0.2) while keeping the uniform-magnitude constraint and divergencefree condition.Here, we construct an analytic solution of a single 3D switchback, so that we are able control the geometry of the switchback and investigate the effect of the switchback's geometry on its evolution.In addition, we analyze how the plasma beta affects the destabilization process of the switchbacks.
The paper is organized as follows: In Section 2, we describe in detail the analytic model of the 3D switchback.In Section 3, we show the results of a series of 3D MHD simulations conducted based on the analytic model.In Section 4, we summarize the study and discuss possible future investigation.

ANALYTIC MODEL OF A 3D AXISYMMETRIC MAGNETIC SWITCHBACK
We consider a cylindrical coordinate system (r, ϕ, z) and assume axisymmetric (∂/∂ϕ ≡ 0).In this case, the magnetic field can be expressed by a scalar flux function ψ(r, z) plus a ϕ-component: which automatically satisfies the ∇ • B = 0 condition.We further decompose ψ into a uniform background field ψ 0 = 1 2 B 0 r 2 , which corresponds to B 0 = B 0 êz , and a fluctuation part ψ 1 (r, z).In order to have a switchback, we assume that ψ 1 has the form where ∆B is the strength of the switchback, H is the scale length of the switchback along z direction, i.e. its "height".
To confine the switchback within a region in radial direction, we further require where R 0 is the radial location of the switchback's center, and R 1 is the radial size of the switchback, i.e. its "width".The above model of magnetic field is 2.5D instead of fully 3D as we are assuming symmetry in ϕ.Nonetheless, the parameter R 0 serves as a proxy of the size of the switchback along the third dimension: The larger R 0 is, the more 2D-like, i.e. planar, the structure is.In Appendix A, we present an alternative model which consists of trigonometric functions instead of Gaussian functions.The solution of f (r) to Equation ( 3) is where erf(x) is the error function and the integral constant C is chosen such that f (0) = 0, i.e.
Thus, the analytic form of B r and B z are Equation (6b) indicates that the switchback is centered at (r = R 0 , z = 0) and has a size (R 1 , H) in the r − z plane.Ideally, B ϕ can be used to maintain the constant-|B| condition.However, Equation (6) does not ensure B 2 r + B 2 z ≤ B 2 0 , and there are certain regions where the constant-|B| condition cannot be preserved.The reason is that in the above model, B r → 1/r for large r, which is not fast enough to maintain B 2 r + B 2 z ≤ B 2 0 .This also leads to a difficulty in setting up the boundary condition because we cannot directly implement periodic boundary conditions if B r does not decay to zero at the boundaries.To overcome this difficulty, we modify the flux function with a predefined function p(r) such that B r decays exponentially outside the switchback region and B z does not change in the z=0 plane: Shi et al.
For p(r), we choose the following function where R m is the radial location at which the modulation of the magnetic field starts, and R d is the scale length of the modulation.With the modified flux function, one can calculate the modified magnetic field which results in the modified B r and B z : . We note that this modulation indeed does not change the magnetic field at the z = 0 plane.A different but equivalent approach to do the modulation is to assume B r is modified by a given function p(r) (Equation (9a)) and B z is modified by an unknown function q(r, z), and then solve q(r, z) such that the modified magnetic field preserves the divergence-free condition.This approach is described in Appendix B. We neglect the tildes above B r and B z hereinafter for convenience.
In Figure 1, we visualize the above analytic model with parameters ∆B/B 0 = 2, R 0 = 2, R 1 = 0.5, H = 3, R m = 2, and R d = 1.5.Panel (a) shows the radial profiles of the three components (colored) and magnitude (black) of the magnetic field.Solid lines are z = 0, i.e. through the center of the structure, and dashed lines are z = H/2.In panel (b), the three orthogonal surfaces are color-coded with B z , and the green and white lines are two magnetic field lines.Here, the green line is connected to the background magnetic field, while the white line is confined within the switchback region.That is to say, in our model, not all of the field lines are connected to the Sun.In contrast, in-situ observations show that the magnetic field inside switchbacks is connected to the Sun (Balogh et al. 1999;Kasper et al. 2019).Thus, this model is a combination of switchbacks and flux ropes.However, we note that magnetic flux ropes are occasionally observed near the switchbacks by PSP (Choi et al. 2024).This phenomenon implies the possibility that reconnection takes effect either during the generation of the switchbacks (Drake et al. 2021;He et al. 2021) or during the propagation of the switchbacks in the solar wind.In Appendix C, we briefly discuss how the reconnection between two switchbacks may generate a flux rope.
One limitation of the model is that, in the region |z| > H and R m ≤ r ≲ R m + R d , B z is smaller than B 0 due to the introduced modification (see panel (b) of Figure 1 and Figure 2).Consequently, there is a finite B ϕ in this region, i.e. the background field is twisted azimuthally.This low-B z region extends infinitely in z.Actually, the ∇ • B = 0 condition means it is impossible to find a switchback model that has finite extents in both r (the x − y plane) and z.The proof is simple: Suppose that the fluctuation part of the field totally vanishes outside a finite radial distance R b in the x − y plane and consider a cylinder with r = R b and 0 ≤ z ≤ H b .Clearly there is no magnetic flux across the side surface r = R b , hence the divergence-free condition writes as As the perturbation in B z at the z = 0 surface is negative definite, we immediately see that B z (z = H b ) ≡ B 0 cannot satisfy the above equation, i.e. for any horizontal plane z = H b , there must be some regions where B z < B 0 .

Simulation setup
Based on the analytic model described in Section 2, we conduct a series of 3D MHD simulations using a pseudospectral MHD code (Shi et al. 2019(Shi et al. , 2020(Shi et al. , 2022a;;Tenerani et al. 2023;Shi et al. 2023).The code evolves the MHD equation set in conservation form on a Cartesian grid.The initial condition comprises a uniform density ρ 0 ≡ 1, a uniform thermal pressure P 0 , magnetic field as described in Section 2 with B 0 = 1, and velocity that is anti-correlated with the magnetic field (V = −B/ √ ρ 0 ).We note that the code uses normalized quantities such that permeability µ 0 is eliminated.The parameters for the conducted simulations are summarized in Table 1.In all the runs listed in Table 1, we use the same spatial resolution ∆ = 0.047 along all dimensions.In addition, several fixed parameters include the switchback strength ∆B/B 0 ≡ 2, the width of the switchback R 1 ≡ 0.5, and the parameters for modification of the magnetic field outside the switchback region: R m ≡ R 0 , R d ≡ 1.5.We do not implement explicit resistivity and viscosity.Instead, we apply a numerical filter in the wave-vector space, which is similar to the filters usually adopted in compact-finite-difference schemes (Lele 1992), to all quantities after each time step to maintain numerical stability.
For a detailed description of the filter, please refer to equations (16,17) of Shi et al. (2020).To verify grid convergence of the simulations, several low-resolution runs (not listed in Table 1) based on Run 0 are conducted, including one run with n z = 256, one run with n x = n y = 256, one run with n x = n y = n z = 256, and one run with n x = n y = n z = 256 and weaker numerical filter.We have verified that the linear growth rate of the instability is not affected by the resolution and the filter strength.
Compared with Run 0, Run 1 has doubled R 0 , Run 2 has doubled H, and Run 3 has one half of H, so that we can investigate the effect of the geometry on the switchback evolution.We also carry out one run, labeled as Run 0IC, with exactly the same initial condition with Run 0, but using incompressible version of the simulation code.For Runs 0-3, P 0 = 1.0 so that the plasma beta is β = 2.0.To investigate how β influences the evolution, we conduct Runs 4-8 which have the same initial configuration with Run 0 but different values of P 0 so that β varies in the range 0.1-2.0.

Results
In Figure 2, we show the snapshots of Runs 0 (top), 0IC (middle), & 8 (bottom) color-coded with B z .The plots are 2D slices in x − z plane, across the central axis (y = 0).The black lines are magnetic field lines projected on the 2D plane.In Run 0IC, the switchback remains stable until the end of the simulation (t = 300), while in Run 0 and Run 8, the switchback remains stable for a certain time but eventually destabilizes.These results indicate that compressibility is a necessary condition for the destabilization of the switchbacks, verifying the conclusion of Tenerani et al. (2020), i.e. parametric decay instability is the major mechanism in the destabilization of the switchbacks. .
In Figure 3, we plot the time evolution of the normalized density fluctuation (δρ/ ⟨ρ⟩) (top) and minimum B z in the simulation domain (bottom) in different runs.Here δρ is the root-mean-square and ⟨ρ⟩ is the average of density.Let us first focus on the left column, which shows results of Runs 0 (blue), 1 (orange), 2 (green), and 3 (red).Apparently, the density fluctuations grow exponentially over certain time intervals in all the Runs before they saturate.As the density fluctuations grow, minimum B z gradually decreases below −1.After the density fluctuations saturate, minimum B z starts to rise rapidly and eventually becomes positive.That is to say, the magnetic field is first compressed as the PDI grows and then destabilizes rapidly after the PDI saturates.We note that, due to the numerical diffusivity, min(B z ) gradually increases with time if there is no compressibility or the growth rate of PDI is small.Comparing Runs 0, 2, and 3, which have different values of H (length/height of the switchback), we can see that the growth rate of the PDI is largest in Run 3 and smallest in Run 2, indicating that longer switchbacks are more stable than the shorter ones, consistent with the PSP observations which show that the switchbacks in the solar wind are typically very long structures (Horbury et al. 2020;Laker et al. 2021).Comparing Run 0 and Run 1, we see that the switchback is much more stable in Run 1, in which R 0 , the size of the structure along the third dimension, is larger.That is to say, the switchback is more stable if it is more 2D-like or planar.In conclusion, Runs 0-3 show that the most stable switchbacks have 2D-like structures and large aspect ratios (length to width).Then we investigate the effect of plasma beta by inspecting the right column of Figure 3. From the time profiles of δρ/ ⟨ρ⟩, we see that the evolution in runs with β > 1 is quite different from the evolution in runs with β ≤ 1.The growth rates of the instability are almost identical in Run 0 (β = 2.0) and Run 4 (β = 1.5).When β drops to unity (Run 5), the growth rate decreases significantly.As we further decrease β to 0.4 (Run 6), the growth rate does not change much compared with Run 5. Eventually, as β continues to drop (Run 8, β = 0.1), the growth rate starts to increase and becomes comparable to β > 1.Although the growth rates of the instability are comparable in Runs 0, 4, and 8, the disturbance of the magnetic field is much slower in Run 8 than in Runs 0 & 4 (bottom right panel), implying that the characteristics of the density fluctuation may be very different depending on plasma beta.
where ω and k are frequency and wavenumber of the density fluctuation normalized to the frequency and wavenumber of the mother wave, A = δB/B 0 is the amplitude of the mother wave normalized to the background magnetic field, and β = C 2 s /V 2 A where C s and V A are the sound speed and Alfvén speed of the background plasma.Note that we use the symbol β instead of β since β = 5/6 × β assuming an adiabatic index of 5/3.In panels (d)-(f) of Figure 4, we plot the solved dispersion relation based on Equation ( 11) for (d) A = 0.1, (e) A = 1.0, and (f) A = 2.0 respectively.In each panel, the solid lines are growth rate and the dashed lines are frequency.Colors of the lines correspond to different values of β.Clearly, regardless of the wave amplitude A, the growth rate decreases with β and the wavenumber corresponding to the maximum growth rate decreases toward k 0 .As pointed out by (Jayanti & Hollweg 1993), for β > 1, the instability is more like a beat instability as the wavenumber and frequency of the generated sound wave are close to those of the mother wave.In this case, the secondary transverse fluctuation is dominated by the forward propagating Alfvén wave with wavenumber close to 2k 0 .In panels (a1)-(c2), we show the snapshots color-coded with ρ 1 in Runs 0 (a, β = 2.0), 5 (b, β = 1.0), and 8 (c, β = 0.1) with the magnetic field lines.Here ρ 1 = ρ − ⟨ρ⟩ is the density fluctuation.We can see that, for β = 2.0, during the linear stage (panel (a1)), the density fluctuation is concentrated within the switchback region.For small β(= 0.1) (panels (c1,c2)), the density fluctuation has shorter wavelength and is distributed along the z-direction, similar to that observed in 2D simulations conducted by Tenerani et al. (2020) (see their Figure 6).For intermediate β (Run 5, β = 1.0), we see growth of density fluctuations both inside and outside the switchback, and the wavelength of the growing mode is larger than Run 8.These results are qualitatively consistent with the linear theory.However, a major discrepancy is that, in our simulations, the growth rate of the instability for β > 1 is larger than β < 1 except for very small β (β = 0.1), while in the linear theory, the maximum growth rate monotonically decreases with β.We point out that, the simulated switchback is not a simple monochromatic, circularly-polarized, one-dimensional Alfvén wave, so the stability of this structure is likely to deviate from linear theory described by Equation ( 11) (Malara & Velli 1996).

SUMMARY
In this study, we develop an analytic model of an axisymmetric 3D magnetic switchback with uniform magnetic field strength.A series of 3D MHD simulations are conducted to investigate the evolution of the switchback.We focus on how the geometry and plasma beta modify the switchback evolution.The major conclusion of the simulations is • Compressibility is a necessary condition for the destabilization of the switchbacks (see also Tenerani et al. 2020).
The switchback is stable in the incompressible-MHD simulation.
• The most stable switchbacks are 2D-like structures with large aspect ratios, i.e. planar structures elongated along the background magnetic field direction.
• The plasma beta plays an important role in the destabilization of switchbacks.For β ≪ 1, the structure is unstable due to the classic parametric decay instability.β ≲ 1 is the most stable regime for the structure.For β > 1, the structure becomes unstable again due to the beat instability (Jayanti & Hollweg 1993).Destabilization of the magnetic field in this case is much faster than β < 1.
Our results may explain why the switchbacks observed by PSP are mostly long structures (Horbury et al. 2020;Laker et al. 2021).Moreover, the β-dependence of the stability shown by the simulations suggests a scenario where the occurrence rate of switchbacks in the inner heliosphere (β ≲ 1) does not change significantly with radial distance to the Sun (Mozer et al. 2020;Pecora et al. 2022), but drops significantly when β grows to greater than unity during the propagation of the solar wind.This may explain why much less switchbacks are observed beyond 1 AU, e.g.Ulysses observation (e.g.Balogh et al. 1999), than inside 1 AU.An investigation of the relation between the occurrence rate of switchbacks and the plasma beta will be necessary to verify this point.Last, we note that our simulations start with a switchback that is already generated.An important question is whether and how different generation mechanisms lead to different structures of the switchbacks.However, this is beyond the scope of the current study.Future numerical and observational studies are necessary to answer this question.2020)), closed field lines, i.e. flux ropes, are inevitably included as there is no straightforward way to construct a model with open magnetic field lines only.The existence of flux ropes around switchbacks indicates prior reconnection events.However, we emphasize that switchbacks accompanied by flux ropes are not necessarily generated by interchange reconnection (Drake et al. 2021;He et al. 2021).One possible mechanism is the interaction between two switchbacks during the propagation of the solar wind, as illustrated by Figure 5.When the two switchbacks collide, reconnection happens at two different locations (a).At each reconnection site, the inverse part of one field line reconnects with the normal part of the other filed line.After the reconnection (b), two new switchbacks and a flux rope are generated.We note that, in practice, it may not be easy to trigger reconnection in switchbacks as they are highly Alfvénic.Nonetheless, reconnection is occasionally observed near the switchback boundaries (Froment et al. 2021).Thus, the mechanism shown by Figure 5 should not be excluded.

Figure 1 .
Figure 1.Visualization of the analytic model of the magnetic field described in Section 2. The model parameters are: ∆B/B0 = −2, R0 = 2, R1 = 0.5, H = 3, Rm = 2, and R d = 1.5.(a) Radial profiles of the r (blue), ϕ (orange), z (green) components, and magnitude (black) of the magnetic field.(b) Two magnetic field lines (green and white) over-plotted on three orthogonal surfaces which are color-coded with Bz.The background field is along z direction.

Figure 2 .
Figure 2. Snapshots of the x − z plane across the central axis (y = 0) color-coded with Bz for Runs 0 (a,top), 0IC (b,middle), and 8 (c,bottom).Black curves are magnetic field lines projected on the 2D plane.

Figure 3 .
Figure 3.Time evolution of (top) the normalized density fluctuations δρ/ ⟨ρ⟩ where δρ is the root-mean-square of the density and ⟨ρ⟩ is the average density and (bottom) the minimum value of Bz in the simulation domain.In the left column, we show the results for Runs 0, 1, 2, 3.In the right column we show the results for Runs 0, 4, 5, 6, 7, 8.

Figure 4 .
Figure 4. Panels (a1)-(c2): Snapshots of the x − z plane across the central axis (y = 0) color-coded with ρ1 for Runs 0 (a,top), 5 (b,middle), and 8 (c,bottom) where ρ1 = ρ − ⟨ρ⟩ is the density fluctuation.Black curves are magnetic field lines projected on the 2D plane.Panels (d)-(f): Linear growth rate (solid) and frequency (dashed) of the parametric decay instability for a monochromatic circularly-polarized Alfvén wave with amplitude (d) A = 0.1, (e) A = 1.0, and (f) A = 2.0.In each panel, different curves correspond to different values of β = C 2 s /V 2 A , where Cs and VA are the sound speed and Alfvén speed of the background field.ω0, k0 are the frequency and wavenumber of the mother Alfvén wave.A = δB/B0 is the amplitude of the mother wave normalized to the background magnetic field.

Figure 5 .
Figure 5. Illustration of two interacting switchbacks.(a) Reconnection happens at two different locations marked by the red dots.(b) Two new switchbacks and a flux rope are generated.

Table 1 .
Simulation parameters.Additional low-resolution runs to verify the grid-convergence are not listed here (see the text).
We thank Dr. Jaye L Verniero for very inspiring discussions.