{'fraction_non_alphanumeric': 0.05139703345981373, 'fraction_numerical': 0.02330305469318922, 'mean_word_length': 4.4528735632183905, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
In-situ switchback formation in the expanding solar wind

Recent near-sun solar-wind observations from Parker Solar Probe have found a highly dynamic magnetic environment, permeated by abrupt radial-field reversals, or"switchbacks."We show that many features of the observed turbulence are reproduced by a spectrum of Alfv\'enic fluctuations advected by a radially expanding flow. Starting from simple superpositions of low-amplitude outward-propagating waves, our expanding-box compressible MHD simulations naturally develop switchbacks because (i) the normalized amplitude of waves grows due to expansion and (ii) fluctuations evolve towards spherical polarization (i.e., nearly constant field strength). These results suggest that switchbacks form in-situ in the expanding solar wind and are not indicative of impulsive processes in the chromosphere or corona.


INTRODUCTION
The recent perihelion passes of Parker Solar Probe (PSP) have revealed a highly dynamic near-sun solar wind Kasper et al. 2019). A particularly extreme feature compared to solar-wind plasma at greater distances is the abundance of "switchbacks": sudden reversals of the radial magnetic field associated with sharp increases in the radial plasma flow (Neugebauer & Goldstein 2013;Horbury et al. 2018;Horbury et al. 2020). Such structures generally maintain a nearly constant field strength |B|, despite large changes to B. It remains unclear how switchbacks originate and whether they are caused by sudden or impulsive events in the chromosphere or corona (e.g., Roberts et al. 2018;Tenerani et al. 2020).
In this letter, our goal is to illustrate that turbulence with strong similarities to that observed by PSP develops from simple, random initial conditions within the magnetohydrodynamic (MHD) model. Using numerical simulations, we show that constant-|B| radial-field reversals arise naturally when Alfvénic fluctuations grow to amplitudes that are comparable to the mean field. We hypothesize that the effect is driven by magnetic-pressure forces, which, by forcing |B| to be nearly constant across the domain (Cohen & Kulsrud 1974;Vasquez & Hollweg 1998), cause large-δB ⊥ /B 0 fluctuations to become discontinuous (Roberts 2012;Valentini et al. 2019). The effect is most pronounced around β ∼ 1 (where β is the ratio of thermal to magnetic pressure). Due to the radial expansion of the solar wind, initially small-amplitude waves propagating outwards from the chromosphere naturally evolve into such conditions (van Ballegooijen et al. 2011;Perez & Chandran 2013;Montagud-Camps et al. 2018). We solve the MHD "expanding box" equations (Grappin et al. 1993), which approximate a small patch of wind moving outwards in the inner heliosphere. We deliberately keep the simulations simple, Electronic address: jonathan.squire@otago.ac.nz initializing with random superpositions of outwards propagating waves, assuming a constant wind velocity, and using an isothermal equation of state. That this generic setup reproduces many features of the turbulence and field statistics seen by PSP strongly suggests that switchbacks form in-situ in the solar wind and are not remnants of impulsive or bursty events at the Sun.

METHODS
We solve the isothermal MHD equations in the "expanding box" frame (Grappin et al. 1993), which moves outwards in the radial (x) direction at the mean solar-wind velocity, while expanding in the perpendicular (y and z) directions due to the spherical geometry. We impose a mean anti-radial (sunward) field B 0 = −B x0x , with initial Alfvén speed v A = B x0 / 4πρ. The mass density ρ, flow velocity u, and magnetic field B evolve according to Here T = (0, 1, 1) and L = (2, 1, 1), a(t) = 1 +ȧt is the current perpendicular expansion withȧ/a the expansion rate, and ∇ ≡ (∂ x , a −1 ∂ y , a −1 ∂ z ) is the gradient operator in the frame co-moving with the expanding flow. The equation of state is isothermal, but with a sound speed c s (t) that changes in time to mimic the cooling that occurs as the plasma expands (for adiabatic expansion, c 2 s ∝ a(t) −4/3 ). A detailed derivation of equations (1)-(3) is given in, e.g., Grappin et al. (1993); Dong et al. (2014). The anisotropic expansion causes plasma motions and magnetic fields to decay in time. Those that vary slowly compared toȧ/a scale as B  Chandran et al. (2011) adapted to match the slow wind seen during PSP's first perihelion (see Chen et al. 2020), which is shown with the thick dotted line. The Athena++ simulations retain nearly constant β as they evolve (the inset zooms out to show β0.02-GS, β0.1-GS, and β35-GS). Panel b illustrates the evolution of the normalized fluctuation amplitude, B 2 ⊥ / B x 2 (solid lines) and ρu 2 ⊥ / B x 2 (dashed lines). The dot-dashed line (almost directly behind the solid line for SW-Gaus) shows the WKB expectation for waves without dissipation: The dotted lines, which are read with the right-hand axis, show the normalized cross-helicity σ c = √ ρu · B / ρu 2 + B 2 . u x ∝ a 0 , u ⊥ ∝ a −1 , with ρ ∝ a −2 (Grappin & Velli 1996). In contrast, for small-amplitude Alfvén waves with frequencies ȧ/a, B ⊥ /B x ∝ a 1/2 and u ⊥ /v A ∝ a 1/2 . By appropriate parameter choices, we ensure that it is primarily this latter WKB effect that leads to large δB ⊥ /B 0 waves in our simulations. The former non-WKB effect preferentially grows B ⊥ / √ ρ compared to u ⊥ , thus acting like a reflection term for largest-scale waves and aiding in the generation of turbulence.
2.1. General considerations In order to simulate the propagation of Alfvénic fluctuations from the transition region to PSP's first perihelion, our simulations would ideally have the following two properties: (i) Large expansion factor. As a parcel of solar-wind plasma flows from r i 1R to r f 35R , it expands by a factor exceeding 35 because of super-radial expansion in the corona, leading (at least in the absence of dissipation) to a large increase in normalized fluctuation amplitudes.
(ii) Minimum dissipation. The damping of inertial-range Alfvénic fluctuations is negligible, aside from nonlinear damping processes (turbulence, parametric decay, etc.). Unfortunately, these properties are difficult to achieve: (i) causes the numerical grid to become highly anisotropic, which can cause numerical instabilities, while (ii) necessitates large numerical grids and careful choice of dissipation properties. Further, probing smaller scales deeper in the inertial Table 1 Properties of the simulations studied in this work (β is shown in figure 1a).
L i ⊥ and L f ⊥ are the initial and final perpendicular box size.
2.2. Numerical methods For the reasons discussed above, we use two separate numerical methods to solve equations (1)-(3), which have different strengths and provide complementary results. The first -a modified version of the pseudospectral code Snoopy (Lesur & Longaretti 2007) -allows for larger expansion factors, making it possible to probe the growth of turbulence and switchbacks from low-amplitude waves. The second -a modified version of the finite-volume code Athena++ (Stone et al. 2008;White et al. 2016) -is better suited for capturing shocks and sharp discontinuities, but develops numerical instabilities for a 4. In Snoopy, we solve equations (1)-(3) using ln ρ and p ≡ ρu as variables, applying the∇ operator by making the k y and k z grids shrink in time. We use sixth-order hyper-dissipation to regularize p, B, and ln ρ, with separate time-varying parallel and perpendicular diffusion coefficients chosen to ensure that the Reynolds number remains approximately constant as the box expands. In Athena++, we use variables B ⊥ = B ⊥ /a and u ⊥ = u ⊥ /a, which leads to MHDlike equations with a modified magnetic pressure (Shi et al. 2019). We use a simple modification of the HLLD Reimann solver of Mignone (2007) with second-order spatial reconstruction, without explicit resistivity or viscosity.

Simulation parameters and initial conditions
Our simulations are listed in Table 1. Each starts with randomly phased, "outward-propagating" (z + ) Alfvén waves with u ⊥ = B ⊥ and u x = 0. The initial normalized fluctu- in each simulation varies between 0.2 and 1. To explore the range of power spectra that might be present in the low corona (see, e.g., van Ballegooijen et al. 2011; Chandran & Perez 2019), we consider three types of initial magnetic power spectra E ⊥0 (k) in our simulations: (i) isotropic spectra peaked at large The magnetic-field structure in final snapshots, concentrating on β1-Gaus-HR (panels a-c). Panel a illustrates the field-line structure of a strong switchback that appears in the simulation, with the color indicatingb x . This switchback appears to be a tangential discontinuity with clear 3D structure to the field lines, while a number of rotational discontinuities are visible in other regions as sharp bends in the field lines. Right-hand panels illustrate the perpendicular turbulence structure (normalized B z and u z ), showing its sharp discontinuities and Alfvénic correlation. Panels c-f show simulated "flybys" through various simulations, scaled to solar-wind timescales using equations (5) and (6). The flybys sample the direction (1, 1/ √ 32, −1/2π) (grey line in panel a) to approximate slow azimuthal motion through a fast radial-wind outflow. Black lines show |B|, and blue, yellow, and red show B x , B y , and B z respectively, each normalized to B x . The dotted-blue line shows B x + u x 4πρ 1/2 (normalized by B x ), illustrating the radial jets associated with switchbacks (Kasper et al. 2019). Table 1, corresponding to the one-dimensional perpendicular and parallel spectra in the critical-balance model (Goldreich & Sridhar 1995;Cho et al. 2002). We avoid imposing any specific features on the flow or field (e.g., large-scale streams; Shi et al. 2019). Two additional simulation parameters are c 2 s (t), which determines β through β = 8πc 2 s (t) ρ/B 2 , and the expansion rate which is the ratio of outer-scale Alfvén time L x /v A to expansion time (ȧ/a) −1 . In our simulations, Γ sim is constant because v A ∝ a −1 , whileȧ and L x are constant. 1 Since most of the fluctuations propagate away from the Sun at speed v A in the plasma frame, we map position x in the simulation to time t in the perihelion (zero-radial-velocity) PSP frame by setting We then convert from simulation units to physical units by equatingȧ/a in the simulations with U/r and setting r = 35R , U = 300 km s −1 , and v A = 115 km s −1 , values that approximate the conditions during PSP's first perihelion pass (Kasper et al. 2019;Bale et al. 2019). The parallel box length L x then corresponds via equation (4) to the physical length scale 4.8 × 10 6 Γ sim km. An outward-propagating Alfvén wave , illustrating how C B 2 is minimized for β ∼ 1 and when the expansion rate is slower. Panel b shows how abrupt changes in field direction, as measured by P(|δ b x | > 1), become more common as δB ⊥ /B 0 grows. The inset shows the good correlation between the two switchback-fraction measures, fb x >0 and P(|δ b x | > 1). Simulation β1-ȧ0 has fb x >0 = P(|δ b x | > 1) = 0 (indicated with a downwards arrow).
with parallel wavelength L x has a frequency Hz (6) in the spacecraft frame. The smallest resolvable parallel wavelength, approximately 4 times the grid scale, corresponds to an Alfvén-wave frequency ( where N is the number of grid points in the x direction. The Snoopy simulations, labeled SW-#, are designed to optimize property (i) of Section 2.1. They start from smallamplitude waves, B 2 ⊥0 /B 2 x0 = 0.2, and expand by a factor of a f /a i = 10 to reach large normalized amplitudes. We prescribe non-adiabatic temperature evolution c 2 s (t) = 0.35t −1 for these runs, causing β to increase from β ≈ 0.2 to β ≈ 1 (figure 1a), in approximate agreement with near-sun predictions from the model of Chandran et al. (2011) with parameters chosen to match conditions seen by PSP (see Chen et al. 2020). The amplitude increases broadly as expected, with some dissipation, as shown in figure 1b. In contrast, the Athena++ simulations (labeled β#-#) are designed to explore the basic physics of expanding-box turbulence at a range of β values. They are limited to a f /a i 4 by numerical instabilities, which necessitates starting from larger wave amplitudes, B 2 ⊥0 /B 2 x0 = 1.0. The fluctuations rapidly (within less than an Alfvén time) develop broadband turbulence and near-constant |B|, resembling that at later times in the Snoopy runs. The temperature evolution in these runs is adiabatic (c 2 s (t) ∝ a −4/3 ), so β remains nearly constant (figure 1a). We choose Γ sim = 0.5 for most simulations so that wave growth is mostly in the WKB regime, although it proved necessary to use an elongated box (L x = 4,ȧ = 0.5; Γ sim = 2) in the Snoopy simulations because of the enhanced dissipation over longer simulation times. We also explore a slower expansion rate Γ sim = 0.2 (β1-ȧ0.2) and the same initial conditions with no expansion at all (β1-ȧ0).

RESULTS
The evolved state from a selection of runs is shown in figure 2, with a focus on β1-Gaus-HR. The field-line visu-alization shows a strong switchback with a complex, threedimensional structure, as well as sharp, angular bends along other field lines. In the spacecraft-like trace shown below (panel c), we see that the switchback (around t = 0) results in a full reversal of the magnetic field to the antisunward directionb x = B x /|B| ≈ +1. The spacecraft remains in the B x > 0 region for approximately 1000s. The jumps at the boundaries of the switchback, with a scale of ∼ 100 s, are resolved by only 8 to 10 grid cells. Because of dissipation, this ∼ 100 s scale is likely around the simulation's minimum resolvable scale. Figure 2c also shows that |B| remains fairly constant despite large changes to B, as ubiquitously observed in solarwind data (Belcher & Davis 1971;Barnes 1981;Bruno & Carbone 2013). 2 Time traces in the lower panels illustrate several other relevant simulations. SW-Gaus resembles β1-Gaus-HR, but with larger time scales and smoother variations due to the numerical method and lower resolution. In contrast, at high-β (β35-GS), the tendency to maintain constant |B| is almost eliminated, although there remains significant variation in B x . With no expansion (β1-ȧ0) the lack of growth of δB ⊥ /B 0 means there are only small fluctuations in B x and no reversals, while |B| is nearly perfectly constant.
More quantitative measures from all simulations are illustrated in figure 3. We measure "switchback fraction" with two methods. First, we simply count the fraction of cells witĥ b x > 0, denoting this as fˆb x >0 . Second, to quantify the prevalence of small-scale sharp changes in field direction, we calculate the PDF ofb x increments, |δ b x | ≡ |b x (x+ )−b x (x)| across scale = | | = L ⊥ /68 (∼ 8 grid cells), and define P(|δ b x | > 1) as the proportion of increments with |δ b x | > 1. Two other useful statistics are the "magnetic compression," which measures the degree of spherical polarization (|B| = constant), 3 and the normalized fluctuation amplitude, which we define as B 2 ⊥ / B x 2 = B 2 y + B 2 z / B x 2 . We choose B x to normalize the fluctuation amplitudes, as opposed to B 2 x or |B| , because B x ∝ a −2 is predetermined by the expansion and B 2 ⊥ / B x 2 is not bounded from above . Figure 3a plots C B 2 and fˆb x >0 for all simulations and shows only modest changes to C B 2 as the fluctuations evolve. Comparing our primary set of β ≈ 1 runs with β0.02, β0.1, and β35, we see -surprisingly -that C B 2 is minimized around β ∼ 1. We also observe that C B 2 increases with expansion rate (explaining the larger C B 2 in Snoopy runs), probably due to non-WKB growth forcing B ⊥ > u ⊥ at large Γ sim , and that C B 2 increases as the initial spectrum becomes steeper. Figure 3b plots normalized fluctuation amplitude and P(|δ b x | > 1), illustrating that the increasing B 2 ⊥ caused by expansion causes the field to become more discontinuous. The good correlation with fˆb x >0 (inset) shows that similar processes lead to radial-field reversals. Comparing P(|δ b x | > 1) between cases β35-GS and β1-GS, which have identical initial wave spectra but very different C B 2 , tentatively supports the hypothesis that magnetic-pressure forces, in striving to keep |B| constant (Vasquez & Hollweg 1998;Roberts 2012), lead to a discontinuous field with switchbacks. 4 We also see that: (i) at low β or with shallower initial spectra, the solutions are similarly discontinuous but the lower fˆb x >0 and P(|δ b x | > 1) is caused by these simulations not reaching such large amplitudes due to increased numerical dissipation (see also figure 1b); (ii) lower resolution simulations have lower switchback fraction, especially in the Snoopy set (SW-Gaus-LR); and (iii) slower expansion (β1-ȧ0.2) causes larger fˆb x >0 and P(|δ b x | > 1) for similar amplitudes.
Magnetic-field and velocity spectra are shown in figure 4. Figure 4a, which pictures SW-Gaus, shows how the spectrum flattens to around E(k ⊥ ) ∼ k −1.5 ⊥ . This is expected and in line with previous results because non-WKB growth from expansion acts like a reflection term for large-scale waves, aiding in the development of turbulence (Grappin & Velli 1996). Figure  4b shows the evolution of the spectral slope from a number of simulations, illustrating the general trend towards slopes around ∼ k −1.5 ⊥ . Although there remain some differences depending on initial conditions, and a larger difference between kinetic and magnetic spectra than seen in the solar wind (Chen et al. 2020), the spectra in runs with nonzeroȧ are quite similar to solar-wind spectra, despite the limited resolution of the simulations. In contrast, simulation β1-ȧ0 exhibits continual steepening of the spectrum due to grid dissipation and shows no tendency to develop solar-wind-like turbulence (this is also true withȧ = 0 and large-amplitude B 2 ⊥0 /B 2 x0 = 4 initial conditions, which have fˆb x >0 0; not shown). Finally, we note that the turbulence exhibits scale-dependent anisotropy similar to standard Alfvénic turbulence (Goldreich & Sridhar 1995;not shown). Figure 5 illustrates the shape of the most intense struc- There is a clear evolution towards spectra around k −1.5 ⊥ , with shallower velocity spectra (dotted lines and dashed lines for Snoopy and Athena++ cases, respectively). The green curve with the top time axis shows β1-ȧ0, which -although it exhibits very low magnetic compressibility -does not evolve towards solar-wind spectra, and does not develop an excess of magnetic energy.
tures in β1-Gaus-HR (other cases are similar). We compute n ⊥ =x×b and n 2 =n ⊥ ×b, and use these to calculate the PDF of structures' length scales along field lines (|b · ∇b| −1 ) and in the perpendicular directions (|n ⊥ · ∇b| −1 and |n 2 · ∇b| −1 ). A rotational discontinuity is characterized by the scale of variation alongb being comparable to variation in the other directions, while a tangential discontinuity varies more rapidly in the directions perpendicular tob. By filtering to include only regions of very large gradient as measured by the Frobenius norm (|∇b| −1 ≤ |∇b| −1 cut = 4L ⊥ /N ⊥ ; solid lines), we see that many of the most intense structures do not exhibit sharp changes alongb, as evidenced by the lack of a sharp peak in the PDF of |b · ∇b| −1 at |∇b| −1 cut (unlike the PDFs of |n · ∇b| −1 ). This shows that the most intense structures are generally tangential discontinuities, consistent with the simulated fly-by and switchback illustrated in figure 2a.
35R , even if the initial wave field is smooth. The answer to this question is crucial for assessing what can be inferred about chromospheric and coronal processes from PSP measurements. Although we cannot rule out switchback formation in the low solar atmosphere, the simplicity of the in-situ formation hypothesis is compelling, so long as it is consistent with observations. Our calculations clearly reproduce important basic features, including sudden Alfvénic radialfield reversals and jets, constant-|B| fields, and spectra around E(k) ∼ k −1.5 . On the other hand, our simulations do not quite reach fˆb x >0 ≈ 6%, as observed in PSP encounter 1 . This discrepancy, however, may simply result from lack of expansion or scale separation: fˆb x >0 depends on resolution (figure 3), particularly for the Snoopy runs, which were designed to approximate the evolution of a solar-wind plasma parcel. Similarly, simulation β1-ȧ0.2, which probes smaller scales, 5 exhibits larger fˆb x >0 and P(|δ b x | > 1) for similar δB ⊥ /B 0 . More detailed comparison to PSP fluctuation statistics will be left to future work (Dudok de Wit et al. 2020;Chhiber et al. 2020;Horbury et al. 2020).
The range of magnetic compressibility C B 2 in our simulations broadly matches PSP observations on the relevant scales. A more intriguing prediction is the dependence of C B 2 on β (figure 3a), which may be a fundamental property of MHD turbulence. While the tendency for magnetic pressure to reduce δ|B| has been studied in a variety of previous works (e.g., Cohen & Kulsrud 1974;Barnes 1981;Vasquez & Hollweg 1996, 1998Matteini et al. 2015), we believe that our observation that this tendency is strongest at β ∼ 1 is new. The increase in C B 2 at high β is unsurprising: thermal pressure forces dominate and interfere with the magnetic-pressure driven motions that would act to reduce δ|B| (absent extra kinetic effects; Squire et al. 2019). At low β, however, the larger C B 2 is more puzzling; we speculate on two possible reasons. First, the theory of Cohen & Kulsrud (1974) for Alfvénic discontinuities predicts a nonlinear time that approaches zero at β ∼ 1 because magnetic pressure resonantly drives parallel compressions that reduce δ|B|. It is, however, unclear how this mechanism operates in a spectrum of oblique waves. Second, predicted parametric-decay growth rates increase at low β (Sagdeev & Galeev 1969;Goldstein 1978;Malara et al. 2000;Del Zanna 2001), which may disrupt the constant-|B| state and increase C B 2 (although this remains unclear; see Cohen & Dewar 1974;Primavera et al. 2019). Parametric decay processes may also aid in the development of turbulence (Shoda et al. 2019). This physics, along with the role of magnetic pressure in producing discontinuous fields and switchbacks, will be investigated in future work.

CONCLUSION
We show, using expanding-box compressible MHD simulations, that a spectrum of low-amplitude Alfvén waves propagating outwards from the sun naturally develops into turbulence with many similarities to that observed in the recent perihelion passes of Parker Solar Probe. Some features that are reproduced by our simulations include the presence of abrupt radial-magnetic-field reversals (switchbacks) associated with jumps in radial velocity, nearly constant magneticfield strength, and energy spectra with slopes around k −1.5 ⊥ . We present two sets of simulations with complementary numerical methods: the first (SW-# in table 1) follows a parcel of plasma through a factor-of-10 expansion, with waves growing from linear to nonlinear amplitudes; the second (β#-#) explores the dependence of expansion-driven turbulence on initial conditions, β, and expansion rate.
Our key results can be summarized as follows: (i) Switchbacks can form "in situ" in the expanding solar wind from low-amplitude outward-propagating Afvénic fluctuations that grow in normalized amplitude due to radial expansion. This indicates that PSP observations are broadly consistent with a natural state of large-amplitude imbalanced Alfvénic turbulence. The strongest switchbacks are tangential discontinuities.
(ii) The magnetic field develops correlations between its components in order to maintain constant |B| (Vasquez & Hollweg 1998). We show that this effect, which has been extensively documented in observations, is strongest at β ∼ 1, absent at high β, and reduced for β 1. (iii) Imbalanced Alfvénic turbulence, with energy spectra ∼ k −1.5 ⊥ and scale-dependent anisotropy, develops from homogenous, random, outward-propagating Alfvén waves in the compressible expanding box model. We hypothesize that (i) results from (ii) -i.e., that switchbacks result from the combination of spherical polarization and expansion-induced growth in δB ⊥ /B 0 -because discontinuities are a natural result of large-δB ⊥ /B 0 fluctuations with constant |B| (Roberts 2012;Valentini et al. 2019).
To more thoroughly assess the in-situ switchback-formation hypothesis and compare detailed statistical measures with observations, two important numerical goals for future work are to run to larger expansion factors without excessively fast expansion or high dissipation, and to reach the steady state where turbulent dissipation balances expansion-driven wave growth. This will require long simulation times at high res-olution. Including mean-azimuthal-field growth (the Parker spiral) and/or improved expansion models (e.g., Tenerani & Velli 2017) may also be necessary. Theoretically, significant questions remain regarding the mechanisms for reducing δ|B| (Cohen & Kulsrud 1974;Vasquez & Hollweg 1998).
We thank Stuart Bale, Chris Chen, Tim Horbury, Justin Kasper, Lorenzo Matteini, Anna Tenerani, and Marco Velli for helpful discussions. Support for JS and RM was provided by Rutherford Discovery Fellowship RDF-U001804 and Marsden Fund grant UOO1727, which are managed through the Royal Society Te Apārangi. VC was supported in part by NASA grants NNX17AI18G and 80NSSC19K0829 and NASA grant NNN06AA01C to the Parker Solar Probe FIELDS Experiment.
High-performance computing resources were provided by the New Zealand eScience Infrastructure (NeSI) under project grant uoo02637 and the PICSciE-OIT TIGRESS High Performance Computing Center and Visualization Laboratory at Princeton University. This work was started at the Kavli Institute for Theoretical Physics in Santa Barbara, which is supported by the National Science Foundation under Grant No. NSF PHY-1748958.