Steady states of the Parker instability

We study the linear properties, non-linear saturation, and a steady, strongly non-linear state of the Parker instability in galaxies. We consider magnetic buoyancy and its consequences with and without cosmic rays. Cosmic rays are described using the ﬂuid approximation with anisotropic, non-Fickian diffusion. To a v oid unphysical constraints on the instability (such as boundary conditions often used to specify an unstable background state), non-ideal magnetohydrodynamic equations are solved for deviations from a background state representing an unstable magnetohydrostatic equilibrium. We consider isothermal gas and neglect rotation. The linear evolution of the instability is in broad agreement with earlier analytical and numerical models; but we show that most of the simplifying assumptions of the earlier work do not hold, such that they provide only a qualitative rather than quantitative picture. In its non-linear stage the instability has signiﬁcantly altered the background state from its initial state. Vertical distributions of both magnetic ﬁeld and cosmic rays are much wider, the gas layer is thinner, and the energy densities of both magnetic ﬁeld and cosmic rays are much reduced. The spatial structure of the non-linear state differs from that of any linear modes. A transient gas outﬂow is driven by the weakly non-linear instability as it approaches saturation.


INTRODUCTION
The magnetic buoyancy (or magnetic Rayleigh-Taylor) instability (Newcomb 1961), modified and enhanced by cosmic rays, is known as the Parker instability (Parker 1958(Parker , 1966(Parker , 1979)).The horizontal magnetic field in a gas layer confined by gravity can be unstable with respect to undular modes which grow exponentially on a timescale comparable to the sound or Alfvén crossing time over the gas density scale height.For the observed scale height 0.5 kpc for the warm interstellar gas in the Solar neighbourhood, with the sound and Alfvén speeds both at about 10 km s −1 , the time scale is of order 5×10 7 yr.It is much shorter than the galactic life-time, and several effects have been explored that can make the galactic discs less unstable including the roles of cosmic ray diffusion (Kuznetsov & Ptuskin 1983;Kuznetsov 1987;Heintz & Zweibel 2018) and rotation (Zweibel & Kulsrud 1975;Foglizzo & Tagger 1994, 1995).We note in this connection that widely used heuristic models of the stratified interstellar medium, informed by observations of the interstellar medium, are likely to be unstable (e.g., Lachièze-Rey et al. 1980).The Parker instability can hardly be completely suppressed in spiral galaxies, so it is important to explore its nonlinear, steady states in order E-mail: anvar.shukurov@ncl.ac.uk to understand why the gas distributions observed in spiral galaxies are not destroyed by the instability.
Aside from its effect on the vertical distributions of the interstellar gas, magnetic field and cosmic rays, the Parker instability plays a significant role in the evolution of galaxies.It contributes to driving galactic outflows (winds and fountains) and thereby to what is known as the 'star formation feedback' on the evolving galactic disc.
The linear stage of the instability has been thoroughly studied and the dispersion relation has been obtained for a wide range of physical models and parameter regimes (e.g., Giz & Shu 1993;Foglizzo & Tagger 1994, 1995;Kim et al. 1997;Shukurov & Subramanian 2021, and references therein).
However, the nonlinear state of the Parker instability is much less understood, in particular because it can only be studied numerically (Kim et al. 2001).Two-dimensional (2D) simulations of the magnetic buoyancy instability of Matsumoto et al. (1990) and Horiuchi et al. (1988) show how the instability saturate differently across different parameter ranges.We compare our three-dimensional (3D) results with these studies, and find that the solutions may be rather different from any 2D linear modes and their superpositions.Matsumoto et al. (1990) find two types of nonlinear 2D solutions; oscillatory for the azimuthal (along the unstable magnetic field) wave number ky exceeding some critical value, or shock-wave dominated for ky less than critical.When the initial magnetic pressure is much smaller than the thermal pressure (Pm/P th < 0.3), the nonlinear oscillations couple to form long waves which eventually decay due to shocks.3D simulations described here do not reproduce such oscillations or shocks at scales comparable to the instability scales, and the proportions of magnetic and cosmic ray energy density and pressure decrease over time near the midplane.The final nonlinear stages of our models show the magnetic field loops to have less ordered, more small-scale structure in the xyplane than these 2D simulations.Hanasz & Lesch (2000) and Hanasz et al. (2002) study the nonlinear Parker instability with rotation in three dimensions, focusing on the evolution of the magnetic field structure and the mean-field dynamo driven by the instability, adopting a very weak initial magnetic field.Hanasz et al. (2002) observe a decrease in the magnetic field strength at later times (following its growth due to the dynamo action) with an increase in its scale height.Machida et al. (2013) study the 3D nonlinear magnetic buoyancy instability (in the absence of cosmic rays) in toroidal geometry.These authors also observe a gradual reduction in the magnetic field strength accompanied by an increase in its scale, possibly slowed down by the dynamo action.Both Hanasz et al. (2002) and Machida et al. (2013) report the development of vertical current filaments in the late stages of the instability, which could be due to a strongly nonlinear magnetic dissipation in the former case and the ideal magnetohydrodynamic approximation in the latter case.We do not observe any such features.
We explore the instability systematically, starting with the magnetic buoyancy instability and then, to clarify their role, adding cosmic rays in the fluid approximation (based on the advection-diffusion equation).Using our linear stage results where the perturbations remain weak, we test various approximations used in the linear stability analyses (which rely on widely diverse assumptions and simplifications).In order to explore the nonlinear states of the instability, we consider imposed and fixed background gas, magnetic field and cosmic ray distributions and solve fully nonlinear equations for the deviations from that state.
The paper is structured as follows.In Section 2, we set out the basic equations (Section 2.1) and specify the physical and numerical parameters used in the simulations.In Section 2.2, we provide a brief description of the gravitational field models used, and in Section 2.3, the boundary conditions are discussed.The results are presented in Sections 3.1 and 3.2 for the linear and nonlinear stages of the instability, respectively (Sections 3.1 also contains a discussion of the assumptions on which analytical studies rely).Section 4 contains a discussion of the implications of our results with emphasis on the nonlinear, statistically steady state, including the relative distributions of the gas, magnetic fields and cosmic rays (Section 4.1), systematic vertical flows (Section 4.2 and the force balance (Section 4.3).Section 5 summarises our conclusions.Appendix A justifies our approach to the implementation of the unstable background state, while the parameter space is explored in Section 4.4 where we discuss the influence on the results of the scale heights of the gas and non-thermal pressure components in the background state, of the gravity profile and the role of the gas viscosity and magnetic diffusivity.

EXPERIMENTAL DESIGN
As the initial state, we consider a plane-parallel magnetohydrostatic equilibrium in the galactic gravitational field, i.e., a stratified layer of thermal gas, horizontal magnetic field and cosmic rays.The instability is caused by the magnetic buoyancy, which depends on the vertical gradient of the magnetic field strength (Hughes & Cattaneo 1987), and so leads to the removal of the initial magnetic field from the layer's midplane, tending to reduce its gradient.As a result, the initial magnetic field is rapidly lost from the system.In galaxies and accretion discs, the large-scale magnetic field is continuously replenished by the dynamo action near the disc's midplane (Shukurov & Subramanian 2021).However, most simulations of the Parker instability do not include the dynamo action as the source of the unstable magnetic field.Therefore, an equilibrium state introduced as an initial condition is rapidly destroyed by the instability in such simulations.This has prevented Rodrigues et al. (2016) from analysing the steady, strongly nonlinear state of the system.Alternatively, boundary conditions can be used to impose a steady background state.However, this would constrain unphysically the evolution of the system as the fixed boundary conditions would require that the deviations from the background state vanish at the boundaries.
Therefore, our approach is to derive and solve (fully nonlinear) equations for deviations from the background state.In fact, this is the standard approach to explore the linear Parker (or any other) instability analytically, but we extend it to capture a fully nonlinear evolution of the perturbations when their magnitude is no longer small.The boundary conditions for the deviations are not restrictive (we use periodic boundary conditions in the horizontal planes), so that the perturbations can evolve freely.In the nonlinear state of the instability, the magnitude of the deviations from the background state is comparable to that of the background state, altering it fundamentally; so it is important to make the model fully flexible to allow for the possibility of such a strong modification.
In this paper, we consider an isothermal gas and neglect rotation to establish a reference model to allow us to identify the effects of radiative cooling, rotation, dynamo action and supernova activity with the associated multi-phase gas structure, which will be discussed elsewhere.

Basic equations
We solve numerically the non-ideal MHD equations for the gas density ρ, its velocity U , total pressure P (which includes the thermal, magnetic and cosmic-ray contributions), magnetic field B = ∇ × A, its vector potential A (with the advective gauge Φ = η∇ • A) and the energy density of cosmic rays cr: where D/Dt = ∂/∂t + U • ∇ is the Lagrangian derivative, g is the gravitational acceleration, Pcr = (γcr − 1) cr, with γcr = 4/3 for the ultra-relativistic cosmic rays, and F is the diffusive flux of cosmic rays discussed below.The viscous term in equation ( 2) has the form ρ where W is the traceless rate of strain tensor, and we use the cgs system with c the speed of light.The advection-diffusion equation for cosmic rays (4) is derived by, e.g., Skilling (1975) and Drury & Völk (1981) and used in this form by, e.g., Gupta et al. (2021) and Rodrigues et al. (2016).
Each variable is represented as the sum of a background equilibrium value, identified with the subscript zero, and a deviation from it, denoted with the prime or a lower-case symbol, (5) The background state and how it is supported in a steady state is described in detail in Section 2.2.It is convenient to include only the deviations in the gas and cosmic ray pressures into P , whereas the magnetic field deviations contribute to both the magnetic pressure and magnetic tension and thus are included as a part of the Lorentz force.However, the background pressure P0 includes all three pressure components.We stress again that the deviations from the equilibrium state are not assumed to be weak: the nonlinear governing equation for them are obtained by subtracting those for the background state from equations ( 1)-( 4) as discussed in Appendix A, where we also present equations governing the background state which ensure that it is static, ∂/∂t = 0.The background magnetic field can be thought of as being maintained by the dynamo action in the disc; in this case, the assumption of its time-invariance is justified if the dynamo time scale is much shorter than the growth time of the Parker instability.The development of the Parker instability when the dynamo time scale is comparable to or longer than the instability time scale will be discussed elsewhere.
In the absence of a background velocity U0, equations (1)-(3), written for the deviations from the background state, reduce to (Appendix A) where we recall that the pressure perturbation P contains only the gas and cosmic ray pressure perturbations (when the latter is included); the magnetic pressure perturbation is accounted for within the Lorentz force.Equation (4) leads to the following equation for the deviation of the cosmic ray energy density from the background state: where F is the evolving cosmic ray flux perturbation derived from the non-Fickian diffusion model for cosmic rays (Snodin et al. 2006;Rodrigues et al. 2016) as (see Appendix A for the derivation) where κij and κ0 ij are, respectively, the total and background diffusion tensors, where the summation convention is understood, Bi = Bi/|B| and B0 i = B0 i /|B0| are the components of the unit vectors in the direction of B and B0, respectively, and the tensors κij and κ0 ij have the same constant diffusivities κ and κ ⊥ , respectively parallel and perpendicular to the relevant magnetic field.The non-Fickian description of the cosmic ray diffusion ensures their finite propagation speed controlled by the parameter τ (this may be important given the very large cosmic ray diffusivity), and prevents numerical singularities in the cosmic ray propagation near the magnetic null points (Snodin et al. 2006).The parameters that control the diffusion cosmic rays parameters are adopted as τ = 0.01 Gyr, κ ⊥ = 3.16 × 10 25 cm 2 s −1 and κ = 1.58 × 10 28 cm 2 s −1 as given in Rodrigues et al. (2016); Ryu et al. (2003, and references therein).We also include a weak isotropic component of the diffusion tensor κ0 = 0.5κ ⊥ to allow for unresolved random magnetic fields (in fact, this part of the diffusion tensor hardly affects the solutions).Since U0 = 0, equation ( 9) follows as the difference between equations (A2) and (A5) given in Appendix A.

The background state
The background state represents the magnetohydrostatic equilibrium in the galactic gravitational field, where the total pressure consists of the thermal gas pressure P th0 , magnetic pressure Pm0 and cosmic ray pressure Pcr0: and we note that the magnetic tension vanishes in this state because the background magnetic field is unidirectional, B0 = (0, B0(z), 0) in the Cartesian fame (x, y, z) with g aligned with the z-axis.
It is convenient to introduce the pressure parameters βm and βcr as the ratios of magnetic Pm and cosmic ray Pcr pressures to the thermal pressure P th : where We consider isothermal gas.In many of our simulations (as listed in Table 1), the sound speed cs is 18 km s −1 corresponding to the temperature of T = 3.2 × 10 4 K, which is within the range for the warm interstellar gas, but we also consider cs ∈ [9, 29] km s −1 to explore the effect of the gas ionization (through the analogue of sound speed) on the instability.In the background state, βm and βcr are constant as listed in Table 1, but these ratios vary in space and time for the deviations from the background state as the instability develops.We mainly consider two combinations of these parameters; βm = 1 and βcr = 0, or βm = βcr = 0.5.
For a given gravitational profile g(z) we solve equation ( 13) to yield the corresponding density, magnetic field and cosmic ray profiles, given the constants βm and βcr.Following some studies of the Parker instability (e.g., Giz & Shu 1993;Rodrigues et al. 2016) we use the gravitational acceleration profile of the form with H = 500 pc, Σ = 10 2 M pc −2 for the surface mass density and G Newton's gravitational constant, which applies to the self-gravitating stellar disc.This gravity profile is used for the sake of comparison with the earlier results but most of our results are based on the gravity field appropriate for the Solar vicinity of the Milky Way, which also includes the contribution from the dark matter halo (Kuijken & Gilmore 1989), where a1 = 4.4 × 10 −9 cm s −2 , a2 = 1.7 × 10 −9 cm s −2 , z1 = 0.2 kpc and z2 = 1 kpc.Li et al. (2017) propose an alternative form for the gravitational acceleration better applicable at |z| 2 kpc, while Ferrière (1998, equation 35) suggests an extension of equation ( 18) to other galactocentric distances.For the gravity profile (18), the isothermal gas density, the strength of the magnetic field (directed along the y-axis) and cosmic ray distribution in the magnetohydrostatic equilibrium follow as (see Appendix A) where we adopt ρ0(0) = 7 × 10 −25 g cm −3 .Although the density distribution in z deviates from an exponential, it is useful to characterise it using the fitted exponential scale height h given in Table 1, where we summarise the model parameters used in the simulations.Similar expressions for the gravity profile ( 17) can be found in Rodrigues et al. (2016).Although the qualitative picture of the development and saturation of the instability is the same under both gravity profiles ( 17) and ( 18), there are important quantitative differences in both the linear eigenmodes and some features of the nonlinear state: in particular, the outflow speed.
The gravity profile (17) has been adjusted to provide similar values of gz to those of (18) within the computational domain used here, |z| Z0 = 1.75 kpc (Fig. 1), so we do not expect any strong differences in the solutions obtained under the two profiles.However, some rather subtle differences may still occur.
The background magnetohydrostatic equilibrium is unstable.A weak initial velocity perturbation is introduced to launch the instability and the perturbations are then evolved using equations ( 6)-( 9).17) (solid) and Eq. ( 18) (dotted), normalized to their magnitude at z = Z 0 .

Boundary conditions and numerical implementation
We use the Pencil Code (Pencil Code Collaboration et al. 2021) to simulate a local Cartesian box within the gas layer; here we neglect rotation (for effects of rotation see Tharakkal et al. 2022).We simulate magnetised interstellar gas and cosmic rays in a region of the size 4 × 4 × 3.5 kpc 3 (symmetric about the midplane located at z = 0) along the x, y and zaxes, respectively, and use periodic boundary conditions in x and y.The extent of the domain along the background magnetic field (y) is chosen to be large enough to accommodate the dominant modes of the instability (see Section 3.1).
The boundary conditions at the top and bottom of the domain, z = ±Z0 with Z0 = 1.75 kpc, are given in equations ( 22)-( 25) and chosen to minimise their effect on the interior of the domain while remaining physically justifiable.Their effect on the processes within the domain is weak or negligible also because Z0 exceeds by a factor 2-3 the scales of the background gas, magnetic field and cosmic ray distributions.
The magnetic field perturbations satisfy the boundary conditions This allows Bz = 0 at |z| = Z0, and the open magnetic lines can support the gas and cosmic ray flows across the boundary.The boundary condition for ρ at |z| = Z0 imposes an exponential decrease at |z| > Z0 (i.e., within the ghost zones used to impose the boundary condition), at the scale that corresponds to the evolving vertical thermal pressure gradient at |z| = Z0.
The horizontal velocity perturbations vanish at the top and bottom of the domain, The vertical velocity at |z| = Z0 should be handled carefully to allow unrestricted gas outflow, while constraining the gas inflow to within numerically stable limits.To suppress numerical instabilities (associated with spurious strong advection of gas from the outside of the domain), any negative uz at z = Z0 and positive uz at z = −Z0 (an inflow, [uz(x, y, z) sign z] |z|=Z 0 < 0) are gradually quenched to zero within the three ghost points at |z| > Z0.At those horizontal positions on the top and bottom boundaries where the gas flows out of the domain ([uz(x, y, z) sign z] |z|=Z 0 > 0), the vertical derivatives of uz of all relevant orders are assumed to Table 1.The list of simulations with the details of the gravitational field profile and parameters used, including the magnetic field strength B 0 (0) and cosmic ray energy density cr0 (0) in the background state at z = 0.The pressure ratios in the background state βm and βcr are defined in equation ( 15) and h is the effective exponential height of the gas in the background state.The viscosity ν and magnetic diffusivity η are given as the multiples of 1 kpc km s −1 = 3 × 10 26 cm 2 s −1 , and cs is the sound speed.The grid spacing along the x, y and z-axes are given in the penultimate column, with the computational box size of (4 × 4 × 3.5) kpc 3 .
Sim1 Eq. ( 17) 0.5 1.0 0 7 0.1 0.03 18 (15, 7, 13) Sim2 Eq. ( 17) 0.5 0.5 0.5 5 0.1 0.03 18 (15, 7, 13) Sim3 Eq. ( 17 vanish in the ghost zone.Thus, the gas that flows out through |z| = Z0 retains its speed in the ghost zone.These conditions, imposed within the ghost zone rather than on the boundary, can be written as inflow: uz → 0 outflow: where the higher-order derivatives are involved when hyperdiffusion is used. The viscosity and magnetic diffusivity are global constants in our simulations, and we considered a selection of their values shown in Table 1; in all cases, the magnetic Prandtl number is Pm = ν/η ≈ 3, as in Rodrigues et al. (2016).The magnitudes of the transport coefficients are close to (but, in all models except for Sim7, smaller than) their turbulent values in the ISM; this is appropriate since our simulations do not not include the turbulence driven by supernovae.As the instability saturates, small-scale fluctuations get stronger making the system more susceptible to numerical instabilities.The fluctuations at the resolution scale are regularised using sixth-order hyperdiffusion with the diffusivity ν6 = 4 × 10 −11 kpc 6 Gyr −1 (e.g., Brandenburg & Sarson 2002;Haugen & Brandenburg 2004), constraining the maximum mesh Reynolds number to Re∆ = cs∆x 5 /ν6.The magnitude of ν6 was chosen to ensure that Re∆ is close to unity along the lowest-resolution dimension, x.Shocks arising in the system are regularised using the second-order shock diffusion (for implementation, see Gent et al. 2020).To avoid negative gas density values that may arise in intense divergent flows, we apply a density floor and impose the restriction ρ0 + ρ ρmin with ρmin = 0.01 M kpc −3 .The magnitude of ρmin chosen is about one order of magnitude lower than the initial minimum density in the simulation box.
Table 1 presents a summary of the simulations which are based on the two forms of the gravity acceleration, given in equations ( 17) and ( 18), and a selection of the pressure ratios of the magnetic field, βm, and cosmic rays, βcr, to the thermal pressure in the background state.These factors control the intensity of the instability and the contribution of cosmic rays to it; the case βcr = 0 represents the magnetic buoyancy instability.We list the sound speed cs in each model.We compare a simulation with the fiducial values (15, 7, 13) pc for the grid spacing (∆x, ∆y, ∆z) along the corresponding axes, to a simulation with (7, 7, 7) pc.

RESULTS
Our main conclusions are drawn from the numerical models Sim5 and Sim6 in Table 1.They differ in the contribution of cosmic rays described using equations ( 9) and (10): Sim5 is the case of the magnetic buoyancy instability, enhanced by cosmic rays in Sim6.The total pressure in the background state remains unchanged since βm + βcr = 1 in each case.

The linear instability
The exponential growth in the perturbations of magnetic field and gas speed in the linear phase of the instability at t 0.4 Gyr (and after the decay of the initial perturbations down to the leading eigenfunction) is clearly visible in Fig. 2a,b.As shown in Fig. 2c, the total magnetic and cosmic ray energy densities vary little with time as long as the perturbations remain much weaker than the background fields; their subsequent nonlinear evolution is discussed in Section 3.2.The growth rate of the instability in the models of Table 1 is quantified in Fig. 3, which presents both the dimensionless growth rate Γh/v and its unit v/h for various choices of the characteristic speed v.
In agreement with earlier analytical and numerical models, cosmic rays make the system more unstable and the growth rate is smaller in Model Sim5 which has no cosmic rays, than in Sim6 where the cosmic rays and magnetic field have equal pressures in the background state.The growth rate Γ of the root-mean-square (r.m.s.) velocity and magnetic fields in these models increases by approximately 20% due to the cosmic rays from 19 Gyr −1 to 25 Gyr −1 , similar to the increase in Γ in the analytical models of Giz & Shu (1993) and Ryu et al. (2003) and the simulations of Rodrigues et al. ( 2016) (the lat- Figure 2. The evolution of the root-mean-square magnitudes at the midplane z = 0 of (a) the magnetic field perturbation |b|, normalised to B 0 (0), and (b) gas speed in the models Sim5 (solid) and Sim6 (dotted) at z = 0.The linear stage of the instability ends at about t 0.4 Gyr.Panels (c) and (d) show the normalised total magnetic and cosmic ray energy densities at the midplane for the same models, B(0)/B 0 (0) xy (solid for Sim5 and dotted for Sim6) and cr(0)/ cr0 (0) xy (dotted, for Sim6) respectively.
ter authors provide a detailed comparison of the growth rates in several models).
The growth rate of the Parker instability is often assumed to scale with (and, often, to be of the order of) the inverse sound or Alfvén crossing time over the density scale height, cs/h or VA/h, respectively.Our results show that this scaling is rather inaccurate as the growth rate depends significantly on other parameters of the system.Figure 3 shows the dimensionless growth rates for models from Table 1 normalised by various crossing times, Γh/v, where v = cs, v = VA and v = c eff with We note that the horizontal spread of the data points indicates a variation of Γh/v by a factor of two between various models for any choice of v. Secondly, if the scaling of Γ with  26) (squares).The horizontal axis presents the corresponding inverse time scale v/h.The model codes are shown within the symbols as in Table 1, e.g.S5 represents Model Sim5.
v/h were perfect, the data points would be form a horizontal line for the most appropriate choice of v.This does not occur for any choice of v, but the scaling with c eff (shown with squares) is marginally better, especially if the Model Sim1 is disregarded.The models that are not included in Fig. 3 have similar growth rates and do not provide additional information.
Analytical and numerical studies of the Parker instability rely on a wide diversity of assumptions and, perhaps not surprisingly, different models claim different properties of the most rapidly growing mode (we note that many studies are based on ideal magnetohydrodynamics).For example, Parker & Lerche (1969) assume that the perturbations are only weakly dependent on z (kz → 0).The most rapidly growing mode of Lachièze-Rey et al. (1980) has kx = 0, (kyh) 2 = 0.5-0.6.Giz & Shu (1993) assume that the perturbations have finite kx and ky and the amplitudes of both bz and uz vary similarly with z while |uz| ∝ ρ −1/2 0 (z) at large |z|; the most rapidly growing mode has a very large kx.Foglizzo & Tagger (1994) adopt kx = 0 and ky = 0 and the magnitude of the perturbations is assumed to be proportional to ρ −1/2 0 (z).In the model of Parker & Lerche (1969), the system is the most unstable at the largest possible kx, and Giz & Shu (1993) assume that the most unstable mode has a large radial wave number, kxh = 10 3 (their Fig. 2), whereas Heintz & Zweibel (2018) assume that kx → 0. Kim & Ryu (2001) argue that kx = kz = 0 in the most rapidly growing mode.Shukurov & Subramanian (2021, Section 2.8.2) assume that the perturbation amplitudes of the density, thermal pressure and Lagrangian displacement scale with ρ −1/2 0 (z) while the amplitude of b is independent of z, and that the perturbations have kx = 0 but finite ky.The wave numbers of the unstable modes along the background magnetic field, ky, are limited from above (Foglizzo & Tagger 1994), For a finite magnetic diffusivity η, the radial wave number must be smaller than (Foglizzo & Tagger 1994) kx cs hη .
(28) Our simulations are consistent with some of these assumptions and conclusions.The spatial structure of the most unstable linear mode is characterised by Fig. 4 where we present the 2D power spectra of uz and bz in the (kx, ky)-plane, averaged over |z| < 1.75 kpc, for the Models Sim5 (without cosmic rays) and Sim6 (with cosmic rays).The spectra of these variables are nearly identical.The presence of cosmic rays does not change kx much but increases ky significantly (by a factor of two).The energy corresponding to the most unstable mode is also higher by a factor of 4 in Sim6 when compared to Sim5.This confirms that cosmic rays make the system more unstable.The values of kx are well within the upper limit kx ≈ 40 kpc −1 of equation ( 28).
The size of the computational domain along the y-axis, 4 kpc, is larger than that, so it can accommodate the most unstable modes (see also Section 4.4.3where we also discuss the wave numbers of the most unstable mode).The spectrum of the perturbations in this model, shown in Fig. 4b, has the maximum at ky 3.1 kpc −1 , in agreement with equation ( 27).As discussed in Section 4.4.3, the power spectrum of the perturbations has a pronounced maximum at ky ≈ 3 kpc −1 but a broad range of modes at kx 3 kpc −1 carry similar energies for the parameter values of Model Sim6.
Giz & Shu (1993) discuss continuous modes assuming that ρ0u 2 z = const.Figure 5 shows the variation with z of the horizontal average ρ0 u 2 z xy (normalized to its mean value across the computational domain; we note that u 2 z xy grows exponentially with time). 1 The variation is by an order of magnitude across the full range of z and is not monotonic.Shukurov & Subramanian (2021) consider solutions where b is independent of z and the perturbations in the other variables decrease as exp[−z/(2h)] (in particular, ρ0u 2 z = const).While the perturbations in our simulations do not satisfy these assumptions closely, the one-dimensional perturbation power spectra in z show that the most unstable modes of uz and bz have kz ≈ 3.5 kpc −1 in Models Sim5, Sim 6 and Sim9, corresponding to the scale 2π/kz = 1.8 kpc.As shown in Fig. 6, the magnitude of perturbations in different quantities vary differently with z, in contradiction to the assumption of the standard stability analysis.In particular, ρ and cr have oscillations in z superimposed on a large-scale variation.The perturbations in the velocity and magnetic field components ux and bx (which have relatively small amplitudes) and uz and bz (having large amplitudes) are oscillatory in z.Mean- while, uy and by do not have pronounced oscillations in z.
The values of kz presented characterise mainly the oscillatory parts of the variation and should be interpreted with caution.
The stratification of the system affects the spatial structure of the linearly unstable modes; in particular, their wave numbers vary with position and time.For example, Fig. 7 shows that the wave number ky of the dominant mode changes with both z and t, generally increasing with time at a fixed z.During the late stages of the exponential growth phase, t 0.3 Gyr, ky increases from 0 to 3 kpc −1 with distance from the midplane and approaches ky ≈ 3 kpc −1 in the nonlinear phase at all distances (see also Fig. 9).Correspondingly, the growth rate of the instability is no longer a constant but becomes a function of z.As shown in Fig. 8, the growth rate of the vertical velocity increases with |z| from 24 Gyr −1 at z = 0 to 32 Gyr −1 at z = 1.5 kpc.This makes the linear evolution more complicated than the linear analytical models suggests but the properties of the dominant mode are broadly consistent with certain caveats discussed here.
Figure 9 illustrates the development of the spatial structures of the magnetic field, gas density and cosmic rays in the linear stage of the instability in panels (a) and (c), where the perturbations are weak at all |z|, and at the end of the linear stage in panels (b) and (d), where the perturbations are quite strong at |z| h = 0.4 kpc but still are weaker closer to the midplane.The magnetic field lines are aligned with the y-axis initially, and the instability excites the undular modes which deform the magnetic field lines and cause the gas and cosmic rays to be redistributed in z.As the system evolves into the nonlinear state, magnetic fluctuations develop at progressively smaller scales.
Model Sim9 has a somewhat higher spatial resolution than Model Sim6 (with all other parameters being similar).The wave numbers and growth rates obtained from the two models are quite similar confirming that model Sim6 resolves well the spatial structures developing in the system.Model Sim8 is similar to Model Sim5 but has a lower sound speed (thus, a lower thermal pressure) and, consequently, a lower gas scale height.This strengthens the instability and makes the gas flow much more irregular at greater heights in the nonlinear state.

The nonlinear instability
The evolution of the r.m.s.velocity and magnetic fields shown in Fig. 2 marks three distinct stages of the instability development.The initial stage of exponential growth is followed by a short transitional (weakly nonlinear) period, which starts at t = 0.4 Gyr and lasts for around 0.1 Gyr during which the   growth slows down and the instability saturates to the final steady state.The total magnetic and cosmic ray energy densities, which remain nearly constant throughout the linear phase, start decaying during the transitional phase.In the statistical steady state they retain only a few percent of the initial values.Meanwhile, the thermal pressure in the system remains almost constant throughout the simulation, except for a small decrease (by 4%) in the nonlinear state.
The changes in the system as the instability evolves are illustrated in Fig. 10 which shows the horizontal cross-section of the vertical velocity.The pattern of the perturbations at the linear stage is quite regular and, of course, the perturbations average to zero in the horizontal planes (Fig. 10a).However, this pattern undergoes substantial modification during the transitional stage when a systematic transient inflow develops at this altitude (Fig. 10b).An outflow at higher al-titudes (visible in Fig. 11g and discussed in Section 4.2) is followed by a rather chaotic pattern in the nonlinear state (Fig. 10c).The rather regular magnetic loops characteristic of the linear instability (Fig. 9) evolve into an ostensibly random pattern.The model of Rodrigues et al. (2016), which is evolved as far as the very early nonlinear stage, has a similar transition to a chaotic state.
Figure 12 shows the 2D power spectra of the perturbations in the vertical velocity and magnetic field which are useful to compare with Fig. 4. As might be expected, nonlinear effects generally broaden the power spectra, with some notable differences between Models Sim5 and Sim6.In Model Sim5, uz develops modes with small ky but large kx while the linear mode remains prominent.The maximum in the Fourier spectrum of bz shifts to higher values of kx at the same ky as in the linear stage.Meanwhile, the maxima in the power spectra of both uz and bz in Model Sim6 shift to lower ky.However, the spectrum of uz is dominated by very small kx and ky whereas the spectrum of bz is maximum at (kx, ky) = (3.1,1.5) kpc −1 .Altogether, the similarity between the spatial structures of the linear modes and the nonlinear solution is limited and not straightforward.
The vertical redistribution of the gas, magnetic field and cosmic rays is shown in Fig. 11, which presents the horizontally averaged deviations from the background distributions in Panels (a)-(c) and (g) as well as the distributions of the total gas density, magnetic field and cosmic rays in the linear, transitional and nonlinear phases of the instability.In the linear stage the perturbations, periodic in horizontal planes, have vanishing horizontal averages.Non-vanishing horizontal averages arise only due to nonlinear effects.The perturbation in the gas density in the transitional and nonlinear stages is positive near the midplane and negative away from the disc: the nonlinear instability leads to a reduction in the density scale height, so the gas disc becomes thinner.As a result, the mean gas density at the midplane increases from 7 × 10 −25 g cm −3 to 1.08 × 10 −24 g cm −3 as the instability develops.The average energy densities of the total magnetic field and cosmic rays at the midplane are reduced by more than 75%, expanding their vertical profiles.Similar behaviour occurs in the simulations of Heintz et al. (2020, e.g., their Fig. 10) which capture the weakly nonlinear, transitional stage of the instability.
Figure 11g shows the vertical velocity planar averages uz xy at different times.The nonlinear effects drive a systematic inflow at |z| 1 kpc and an outflow at |z| 1 kpc in the transitional stage which, however, is transformed into a weaker inflow at t = 0.9 Gyr while the system still adjusts towards the steady state.These features are further discussed in Section 4.2.At t = 1.6 Gyr, the system of Model Sim6 reaches a statistically steady state with a residual inflow with | uz xy |/cs 0.5 at |z| 1 kpc (which, however, carries little mass).The magnetic field and cosmic ray energy density continue to decrease, saturating at the midplane values B(0) xy /B0(0) ≈ 0.16 and cr(0) xy / cr0(0) ≈ 0.02, while the gas density increases to ρ(0) xy /ρ0(0) ≈ 1.6.
The magnetic buoyancy instability is driven by the vertical gradient of the magnetic field strength, and it might be expected that it would saturate via reducing the gradient to a marginal value.However, the system follows a much more dramatic path, getting rid of the magnetic field altogether.As the instability produces strong vertical magnetic perturbations, the cosmic rays are channelled out from the disc and diffuse along the magnetic field at a high rate.The instability results in wide distributions of both magnetic field and cosmic rays enveloping a relatively thin thermal gas disc.A similar kind of evolution also occurs in the simulations of Heintz et al. (2020) and Girichidis et al. (2022).However, the form of the nonlinear state is sensitive to such features of the system as rotation, and to the mechanism maintaining the unstable magnetic field in the disc.As we show in Tharakkal et al. (2022), the instability in a rotating system can lead even to the reversal of the magnetic field in the disc.

DISCUSSION
In this section we discuss the implications of our results for the dynamics of the interstellar medium and describe the force balance in the system.
Table 2.The Pearson correlation coefficients r between various energy densities in the statistically steady state at t = 1.6 Gyr presented as a, b, where a and b refer to the altitudes z = 0.5 and 1.5 kpc, respectively.

Cross-correlation between energy densities
Understanding the relative spatial distributions of the gas, magnetic field and cosmic rays in the interstellar medium is crucial for the interpretation of the observations in the radio and other wavelength ranges.In particular, the assumption of a tight correlation between cosmic rays and magnetic fields (e.g., the energy equipartition) is routinely used in the interpretations of synchrotron observations (Seta & Beck 2019, and references therein).The correlation or anti-correlation between the thermal electron density and magnetic field can affect significantly the Faraday rotation (Beck et al. 2003).
Table 2 presents the Pearson cross-correlation coefficient r between the fluctuations in the energy densities i = i − i xy , with • • • xy for the horizontal average and i = th, cr, m, k for the thermal, cosmic-ray, magnetic and kinetic ( k = 1 2 ρu 2 ) energy densities, in the nonlinear stage: The relations between these variables are different near the midplane and at a higher altitude, so both are presented (the results at z < 0 are similar).The fluctuations in thermal energy density are significantly anti-correlated with those in the magnetic and cosmic ray energy densities near the midplane and with the cosmic ray energy fluctuations at the higher altitude.This appears to reflect the horizontal pressure balance where the total average pressure is independent of x and y.The kinetic energy density shows no significant correlation with any other variable.The distribution of the cosmic rays is only weakly correlated with that of the magnetic field near the midplane (r = 0.58) but not away from it (r = −0.46).
Figure 13 shows the cross-correlation coefficient of the fluctuations in the magnetic field and cosmic ray energy densities at various heights in the nonlinear stage of Model Sim6.The cross-correlation coefficient of these two quantities increases with |z| reaching a maximum at |z| 0.5 kpc and then decreases to become negative (anti-correlation) at |z| 1 kpc.
To clarify the cause of the correlations, consider the magnetic field and cosmic rays.The magnetic field and cosmic ray energy density can be decomposed into their horizontal where xy and σ 2 = ( ) 2 xy are the variances of the magnetic (the factor 8π can be omitted here) and cosmic ray energy densities.Since Both terms contribute similarly to the correlation coefficient: r(b 2 , cr) = 0.23 and r(B • b, cr) = 0.43 at z = 0.5 kpc and −0.48 and −0.18 at z = 1.5 kpc, respectively.The variance of the magnetic energy density can be represented as  Beck et al. (2003) (see also Sect.13.2 of Shukurov & Subramanian 2021, where typographical errors of the original publication are corrected), the anti-correlation of the total magnetic field strength and thermal gas density can significantly affect the interpretation of the Faraday rotation observations in terms of the magnetic field strength, leading to underestimated field strength when the two variables are assumed to be uncorrelated.The weak correlation between the total magnetic and cosmic ray energy densities near the midplane, r ≈ 0.6, and a similarly weak anti-correlation at a higher altitude, r ≈ −0.5, are inconsistent with the assumption of the local energy equipartition between magnetic fields and cosmic rays at the scales of the fluctuations produced by the Parker instability, of order 1 kpc and less.
The energy equipartition between cosmic rays and magnetic fields is often justified by arguing that cosmic rays can accumulate in the galactic disc until they achieve the energy density comparable with that of the magnetic field.Then their pressure can open up magnetic lines locally and they can be released from the galaxy, after which their pressure decreases down to a level controlled by the magnetic field.Such a self-regulation is argued to lead to their energy equipartition.Our model contains all the elements of these processes in their full nonlinear implementation, and yet no equipartition occurs.In particular, vertical magnetic fields that facilitate the vertical diffusive transport of the cosmic rays are produced by the magnetic buoyancy without any need for the cosmic ray pressure to affect the structure of the magnetic field (even though cosmic rays can enhance the magnetic buoyancy instability).

Vertical flows
The redistribution of the gas, magnetic field and cosmic rays in the unstable system involves systematic vertical flows.As discussed above, the vertical velocity averaged over a horizontal plane necessarily vanishes in the linear stage of the instability but, remarkably, systematic flows emerge in the transitional stage as a nonlinear effect.Figure 14 clarifies the balance between the inflow (zuz < 0) and outflow (zuz > 0) as the instability develops presenting, for 0 z 1.75 kpc (the picture at z < 0 is similar), separately the vertical velocity averaged over those regions in the (x, y)-planes where uz > 0 and uz < 0, denoted uz + and uz −, respectively.In the linear stage, t = 0.3 Gyr, the outflows balance the inflows and the horizontally averaged vertical velocity uz xy = uz − + uz + is negligible.As the nonlinear effects become stronger, at t = 0.6 Gyr, a systematic inflow develops around z = 1 kpc and outflow with the maximum outflow speed of 7 km s −1 is maintained at higher altitudes.At a later time, t = 0.9 Gyr, when the nonlinear effects are still stronger but the system continues evolving, an inflow dominates with −4 km s −1 at z = 1.5 kpc.In the advanced nonlinear stage at t = 1.6 Gyr, the r.m.s.inflow speed saturates at −8 km s −1 at z = 1.5 kpc.
The fact that the Parker instability can drive systematic vertical flows can be of significance for the galactic evolution and regulation of star formation in the disc.Figure 15 shows the evolution of the mass flux at different heights.In the transitional phase (0.4 < t < 0.5 Gyr), the system develops a strong inflow near the midplane.During this period, the gas is redistributed into a thinner layer.As the system evolves into the nonlinear phase (t > 0.5 Gyr), the mass flux reduces, through decaying oscillations to become negligible in the late nonlinear stage.The outflow at higher altitudes discussed above involves dilute gas and does not transport any significant gas mass in this model.

Force balance
Figure 16 shows the horizontally averaged values of the thermal, magnetic, cosmic ray pressure gradients (the averaged magnetic tension force is negligible at all times) as well as the gravitational force at three different stages of the evolution.During the linear stage (Fig. 16a) the horizontal averages of the magnetic field and cosmic rays vanish, and the average force balance is exactly the same as in the background distributions.From the transitional phase to the nonlinear state (Fig. 16b,c, respectively), the contributions of the magnetic and cosmic ray pressures decrease systematically as their scale heights increase, whereas the gradient of the thermal pressure increases because the gas scale height decreases (see Fig. 11).To clarify the nature of the gas outflow prominent during the transitional stage of the instability, we show in Fig. 17 the evolution of the magnitude of the horizontally averaged Lorentz force per unit volume at different heights.During the linear phase of the instability, the average magnetic force remains equal to that in the background (initial) state but it increases sharply during the transitional phase as significant vertical magnetic fields emerge.This increase is stronger at high altitudes (by a factor of 30).Since the gas density decreases with |z|, the increase in the force per unit mass is even stronger, and the magnetic force is clearly the driver of the systematic gas outflow in the transitional state (Section 4.2).As the instability saturates, the Lorentz force decays leaving for the thermal pressure gradient alone to balance the gravitational force.
In the initial state, the magnetic and cosmic ray pressures vary with z exactly as the thermal pressure, so both βm and βcr of Eq. ( 16) are independent of z.Analytical analyses of the linear Parker instability usually rely on the assumption that both βm and βcr remain independent of z as the instability develops.Figure 18 shows the variation with z of the horizontally averaged pressure ratios.In the linear stage, both ratios remain unchanged at |z| 1 kpc, which justifies the assumption of the analytical studies, but increase significantly at |z| 1 kpc, especially βm.As the system evolves into the linear and nonlinear stages, both ratios are reduced to negligible values at z 1 kpc due to the increase in the thermal pressure and the reduction in magnetic and cosmic rays pressures near the midplane caused by the redistribution of the gas towards the midplane and escape of the magnetic field and cosmic rays to greater heights.As a result, the outer layers, |z| 1 kpc, are strongly dominated by the magnetic field and cosmic rays -and yet the magnetic force and the gradient of the cosmic pressure are both negligible in comparison with the thermal pressure gradient in the nonlinear state.

Sensitivity to parameters
The instability is sensitive to a wide range of parameters including the scale height of the unstable magnetic field, the ratios of the magnetic and cosmic ray pressures to the thermal pressure, βm and βcr, and the form of the gravity profile g(z).In this section we discuss the effects of the key parameters and assumptions, with emphasis on the nonlinear states.The models Sim1-Sim4 have the gravity profile (17) while Sim5-Sim9 use g(z) from equation (18).To assess the role of the cosmic ray propagation governed by equation (4), we also consider a model where the only effect of the cosmic rays is to contribute to the total pressure, assuming that this contribution is equal to the magnetic pressure (thus, just doubling the magnetic pressure in such a model).Parameters of the models discussed below can be found in Table 1.

Gas scale height and nonthermal pressures
The purpose of Models Sim1-Sim4 is to explore the role of the gas (and magnetic field) scale heights and nonthermal pressures in the background state.The initial (background) state in Models Sim3 and Sim4 has a relatively large gas scale height h = 1 kpc (and a correspondingly larger scale height of the background magnetic field) and lower magnetic and cosmic ray pressures than in the reference model Sim2 (with the same gravity profile) while Model Sim1 represents the case of the pure magnetic buoyancy instability (βcr = 0).The adopted values of βm, βcr and h require a higher sound speed in these models (see equation ( 9) of Rodrigues et al. 2016).In this sense, these models allow for the presence of the hot interstellar gas.
Figure 19 shows the evolution of the magnetic and cosmic ray energy densities for these two systems -it is useful to compare it with Fig. 2. As expected, the instability is weaker (lower growth rate of the perturbations) in Models Sim3 and Sim4 than in Sim2.The system Sim3 (βm = βcr = 0.5) has a higher growth rate than Sim4 (βm = βcr = 0.25) and, in the transitional stage (1 t 2 Gyr) loses the magnetic field and cosmic rays faster.The amount of the magnetic field and cosmic ray energies lost in the transitional stage is also slightly smaller in the system with weaker instability.

Gravity profile
The difference between the gravity profiles ( 17) and ( 18) is not strong and yet it can be of a physical significance since the gravity fields vary in magnitude and form between different galaxies and between different locations within a specific galaxy.In the linear stage, the development of the instability mostly depends on the system properties at |z| h, where h = 0.3-0.5 kpc in most of the models considered here, so it might by reasonable to expect the modes of the linear instability to be broadly similar under both gravity profiles but it is useful to identify any subtler changes.In Models Sim2 and Sim6, which only differ in the gravity profiles, the growth rate of the instability is about 25 Gyr −1 .
In Fig. 20, we compare the 2D power spectra of the vertical velocity uz and magnetic field bz in the linear instability phase for Models Sim2 and Sim6; the higher value of h in Model Sim2 reflects the somewhat weaker gravity.As also shown in Table 1, both models have similar wave number ky of the most unstable mode although the spectrum in ky around the maximum at ky ≈ 3 kpc −1 is wider in Model Sim2 (suggesting that this mode is less dominant).Furthermore, the mode structures in x and z differ significantly: (17) leads to significantly larger values of kx and a weaker variation of the solution with z.The nonlinear state can be expected to be more sensitive to the system properties at greater distances from the midplane where the two gravity profiles differ stronger (see also Li et al. 2017).However, this difference within our simulation domain is not strong enough to affect, in particular, the vertical flows.In Models Sim2 and Sim6, the outflow speeds during the transitional stage at t = 0.6 Gyr are uz xy = 7 km s −1 and 12 km s −1 at z = 1.5 kpc.Similarly, the maximum magnitude of the downward velocity in the late nonlinear stage at t = 1.25 Gyr is 5 km s −1 and 8 km s −1 , respectively.

The size of the computational domain
The size of the computational domain needs to be large enough to accommodate the most rapidly growing mode.Because of the periodic boundary conditions in x and y, only a discrete set of modes can be excited in the numerical model, with kx,y being multiples of 2π/Lx,y, where Lx and Ly are the domain sizes along the x-and y-axes, respectively.In the models discussed above, the smallest admissible horizontal wave number is kmin = 2π/Lx = 2π/Ly = 1 2 π kpc −1 .The smallest wave number along the z-axis (along which the boundary conditions are non-periodic) is 2π/Lz ≈ 1.8 kpc −1 .In the nonlinear state, the solution is not a perfectly periodic function of (x, y) within the domain, so that the constraint related to the periodic boundary conditions is much less important.
In the reference model Sim6, the most unstable linear mode has (kx, ky, kz) ≈ (1.6, 3.1, 1.8) kpc −1 .This mode has the largest admissible wavelengths in the x-and z-directions and two complete wavelengths fit along the y-direction.In order to confirm that the size of the domain does not affect excessively the parameters of the most unstable mode, we explored Figure 18.From Model Sim6 the horizontally averaged ratios to the thermal pressure, see Eq. ( 16), of (a) magnetic pressure βm and (b) cosmic ray pressure βcr.Vertical distributions are shown at times specified in the legend, which represent the initial state, linear, transitional and nonlinear stages of the instability, respectively.the linear phase of the instability using lower-resolution simulations with the same parameters as Model Sim6 but with bigger domains in x and y.The results are presented in Table 3. Figure 22 shows how the power spectra of the linear fluctuations depend on the horizontal size of the computational domain.The power spectra in x (Fig. 22a) do not vary much between the cases Lx = 8 and 20 kpc and suggest that a wide range of modes with kx 2 kpc −1 are excited by the linear instability.The energy at short wavelengths, kx 2 kpc −1 , is overestimated in the smaller domain with Lx = 4 kpc but a broad maximum at kx ≈ 1.6 kpc −1 occurs for all values of Lx.The power spectra in y (Fig. 22b) have a prominent maximum at ky ≈ 3 kpc −1 in all cases, corresponding to a single most unstable mode.The power spectra in z shown in Fig. 22c have a strong maximum at the smallest available wave number kz ≈ 3.6 kpc −1 for the two smaller domain sizes (representing the large-scale vertical variations shown in Fig. 6) but modes with kz 3 kpc −1 dominate when (Lx, Ly) = (20, 24) kpc: the solution develops a nearly uniform component in z.Remarkably, the horizontal size of the domain affects the vertical structure of the solution.The modes at large wave numbers in all three directions that have an approximately power-law spectrum apparently result from the instability since they are sensitive to the domain size but nonlinear effects, however weak at t 0.4 Gyr, may also contribute.We conclude that the domain size of the models presented in Table 1 is sufficient to capture the properties of the most rapidly growing modes.

The role of dissipation
The only difference of Model Sim7 from the reference Model Sim6 is the higher (by a factor of 10) gas viscosity ν and magnetic diffusivity η (apart form a lower spatial resolution justified by the fact that stronger dissipation suppresses small-scale structures).The magnitudes of the transport coefficients in Model Sim7 are close to their turbulent values in galaxies.The instability growth rates in the two model are compared in Fig. 21; the growth rate is hardly affected by the change.However, stronger dissipation leads to a noticeably larger magnetic field strength and cosmic ray energy density in the steady state; the effects of dissipation in the nonlinear state of the instability perhaps deserve further analysis.

CONCLUSIONS
The nonlinear, saturated state of the Parker instability is different from what might be expected.It is very different from its linear modes.The instability is driven by the gradient of the magnetic field energy density along the gravitational acceleration, and it might be expected that the saturated state would have a reduced gradient of the magnetic field strength.However, the models described here behave differently: not only is the magnetic field gradient reduced, but its strength also reduces throughout the system.The instability produces ubiquitous local vertical magnetic fields which facilitate rapid removal of the cosmic rays from the system by their rapid diffusion along the magnetic field.As a result, the statistically steady state is notable for a small scale height of the gas
The top and bottom boundaries are open for the cosmic ray fluid, allowing it to escape along open magnetic lines at |z| = ±Z0,

Figure 3 .
Figure3.The dimensionless growth rate Γh/v of the r.m.s.magnetic field perturbation is shown for various models, normalized with various crossing times over the density scale height h: v = cs (diamonds), v = V A (circles) and v = c eff of equation (26) (squares).The horizontal axis presents the corresponding inverse time scale v/h.The model codes are shown within the symbols as in Table1, e.g.S5 represents Model Sim5.

Figure 4 .
Figure 4. 2D power spectra in the (kx, ky)-plane at t = 0.3 Gyr (the linear stage of the instability) at |z| 1.75 kpc for the perturbations as indicated in the colour bar labels of vertical velocity uz in kpc 2 km 2 s −2 and the vertical magnetic field bz in kpc 2 µG 2 .Panels (a) and (c) depict Model Sim5 and (b) and (d) present Model Sim6.

Figure 5 .
Figure 5.For Sim6 during the linear and weakly nonlinear stages of the instability at times indicated in the legend the vertical variation of ρ 0 u 2 z xy normalized by its volume averaged energy density ρ 0 u 2 z .

Figure 6 .
Figure 6.The vertical profiles of the perturbation magnitudes normalised by their volume averaged r.m.s.values throughout the domain in Model Sim6 for different variables as denoted in the legends at (x, y) = (−1.2,−1.8) pc and t = 0.3 Gyr.The profiles at fixed (x, y) are shown because the horizontal averages of the perturbations are negligible during the linear phase of the instability.

Figure 7 .Figure 8 .
Figure 7. Time-latitude diagram of the azimuthal wave number ky for y during the linear instability of Model Sim6.

Figure 9 .
Figure 9. Upper row: the gas density distribution (colour) and magnetic lines, with the local field direction indicated, in the cross-section at x = 800 pc in Model Sim6 at (a) t = 0.30 Gyr and (b) t = 0.45 Gyr.Lower row: as above but for the cosmic ray energy density (colour) at (c) t = 0.30 Gyr and (d) t = 0.45 Gyr.The separation of the magnetic lines does not represent the field strength.

Figure 10 .
Figure 10.The vertical velocity perturbation uz from Model Sim6 at z = 0.4 kpc during (a) the linear stage of the instability (t = 0.3 Gyr), (b) the transitional stage (t = 0.6 Gyr) and (c) the nonlinear state (t = 0.9 Gyr).

Figure 11 .
Figure11.From Model Sim6 the horizontally averaged vertical profiles of (left column (a)-(c)) the normalised perturbations and (centre column (d)-(f )) the corresponding total profiles.The solid, dotted, dashed and dash-dotted lines Timing of each profile is as indicated in the legend of (f ).respectively.Panel (g) shows the horizontally averaged vertical velocity normalised to the sound speed uz xy /cs.

Figure 12 .
Figure 12.As Fig. 4 at |z| 1.75 kpc, but t = 0.9 Gyr (the nonlinear stage of the instability) for the perturbations as indicated in the colour bar labels.Panels (a) and (c) depict Model Sim5 and (b) and (d) represent Model Sim6.

Figure 13 .
Figure 13.The vertical variation of the Pearson cross-correlation coefficient between the magnetic ( m) and cosmic ray ( cr) energy densities in Model Sim6 averaged over the time interval 1.2 t 1.7 Gyr with the cadence ∆t = 0.1 Gyr.
Figure13shows the cross-correlation coefficient of the fluctuations in the magnetic field and cosmic ray energy densities at various heights in the nonlinear stage of Model Sim6.The cross-correlation coefficient of these two quantities increases with |z| reaching a maximum at |z| 0.5 kpc and then decreases to become negative (anti-correlation) at |z| 1 kpc.To clarify the cause of the correlations, consider the magnetic field and cosmic rays.The magnetic field and cosmic ray energy density can be decomposed into their horizontal averages and fluctuations, B = B + b and cr = + , where B = B xy and = cr xy , such that B xy = B, b xy = 0, xy = and xy = 0.The correlation coefficient of the total magnetic and cosmic ray energies follows as

since b 2 B
• b xy = 0 for symmetric probability distributions of the Cartesian components of b.To provide an illustration, if the magnetic field fluctuations are isotropic and the Cartesian component bi of b are statistically independent, each having the standard deviation σi, we have b 2 xy = 3σ 2 i and b 4 i xy = 3σ 4 i

Figure 14 .Figure 15 .
Figure 14.The variation of the horizontally averaged vertical velocity component, uz (solid lines) at z > 0 in Model Sim6 at (a) t = 0.3 Gyr, (b) 0.6 Gyr and (c) 1.6 Gyr.The horizontally averaged outflow (inflow) speed is shown dotted (dash-dotted) while the horizontally averaged velocity, their sum, is shown solid.