Steady states of the Parker instability: the effects of rotation

We model the Parker instability in vertically stratified isothermal gas using non-ideal MHD three-dimensional simulations. Rotation, especially differential, more strongly and diversely affects the nonlinear state than the linear stage (where we confirm the most important conclusions of analytical models), and stronger than any linear analyses predict. Steady state magnetic fields are stronger and cosmic ray energy density higher than in comparable nonrotating systems. Transient gas outflows induced by the nonlinear instability persist longer, of order 2 Gyr, with rotation. Stratification combined with (differential) rotation drives helical flows, leading to mean-field dynamo. Consequently, the nonlinear state becomes oscillatory (while both the linear instability and the dynamo are non-oscillatory). The horizontal magnetic field near the midplane reverses its direction propagating to higher altitudes as the reversed field spreads buoyantly. The spatial pattern of the large-scale magnetic field may explain the alternating magnetic field directions in the halo of the edge-on galaxy NGC 4631. Our model is unique in producing a large-scale magnetic structure similar to such observation. Furthermore, our simulations show that the mean kinetic helicity of the magnetically driven flows has the sign opposite to that in the conventional non-magnetic flows. This has profound consequences for the nature of the dynamo action and large-scale magnetic field structure in the coronae of spiral galaxies which remain to be systematically explored and understood. We show that the energy density of cosmic rays and magnetic field strength are not correlated at scales of order a kiloparsec.


INTRODUCTION
The Parker instability is a magnetic Rayleigh-Taylor or magnetic buoyancy instability modified by cosmic rays that carry negligible weight but exert significant pressure.The instability is an important element of the large-scale dynamics of the interstellar medium (ISM) as it affects the vertical distributions of the gas, magnetic fields and cosmic rays and can drive gas outflows, thereby affecting the star formation.In our previous work (Tharakkal et al. 2022a), we explored the development of the instability, with a focus on its nonlinear saturation, in a non-rotating disc with imposed unstable distributions of the gas, magnetic field and cosmic rays.Among the essentially nonlinear features of the instability are a transient gas outflow in the weakly nonlinear stage and a strong redistribution of magnetic fields, cosmic rays and thermal gas, resulting in a thinner thermal gas disc and very large scale heights and low energy densities of the magnetic field and cosmic rays.In this paper, we address the effect of rotation on the Parker instability.
Rotation is known to reduce the growth rate of the weak perturbations but it does not suppress the instability completely (Zweibel & Kulsrud 1975;Foglizzo & Tagger 1994, 1995;Matsuzaki et al. 1998;Kowal et al. 2003).However, rotation introduces a fundamentally new feature to the system: under the action of the Coriolis force, the gas flows produced by the instability become helical and can drive mean-field dynamo action that generates a magnetic field at a large scale comparable to that of the initial unstable configuration.Hanasz (1997), Hanasz & Lesch (1997, 1998) and Thelen (2000a) simulate numerically the mean-field dynamo action driven by the magnetic buoyancy with and without cosmic rays, while Moss et al. (1999) present an analytical formulation.A striking feature of the nonlinear evolution of a rotating system, noticed by Machida et al. (2013) in their simulations of the galactic dynamo using ideal magnetohydrodynamics (MHD), is the possibility of quasi-periodic magnetic field reversals at the time scale of 1.5 Gyr, both near the disc midplane and at large altitudes.This appears to be an essentially nonlinear effect that relies on rotation since the linear instability does not develop oscillatory solutions and the nonlinear states are not oscillatory without rotation (Tharakkal et al. 2022a).Foglizzo & Tagger (1994, their Section 7.1) find that the Parker instability can be oscillatory in a certain range of the azimuthal wave numbers.Machida et al. (2013) relate the reversals to the magnetic flux conservation, but we note that the large-scale magnetic flux is not conserved when the mean-field dynamo is active.Our simulations of the nonlinear Parker instability in a rotating system suggest a different, more subtle explanation that relies on the correlations between magnetic and velocity fluctuations not dissimilar to those arising from the -effect that drives the meanfield dynamo action (see below).Large-scale magnetic fields whose horizontal direction alternates with height emerge in the simulations of mean-field dynamo action by Hanasz et al. (2004).This spatial pattern may be related to the field reversals near the midplane.
Table 1.The list of simulation runs discussed: the numerical resolutions along each axis, the angular velocity and rotational shear, and the instability growth rate computed for   and   .We explore the effects of rotation on the Parker instability in a numerical model similar to that of Tharakkal et al. (2022a), quantifying both its linear and nonlinear stages and identifying the roles of the Coriolis force and the velocity shear of the differential rotation.We consider the instability in a local rectangular box with parameters similar to those of the Solar neighbourhood of the Milky Way.The structure of this paper is as follows.Section 2 describes briefly the numerical model, and in Section 3 we consider the linear stage of the instability.Section 4 presents a detailed comparison of the distributions of the thermal and non-thermal components of the system in the nonlinear, saturated stage of the instability and how they change when the rotational speed and shear rate vary. in Section 5, we clarify the mechanism of the magnetic field reversal and Section 8 discusses the effects of rotation on the systematic vertical flows.The mean-field dynamo action of the motions induced by the instability is our subject in Section 6 where we discuss the kinetic and magnetic helicities.

BASIC EQUATIONS AND THE NUMERICAL MODEL
We use a model very similar to that of Tharakkal et al. (2022a), with the only difference being that we now consider rotating systems, with either a solid-body or differential rotation.We consider the frame rotating at the angular velocity of the centre of the domain with the -axis aligned with the gravitational acceleration and the angular velocity , the -axis directed along the azimuth and the -axis parallel to the radial direction of the local cylindrical frame.Vector -components are occasionally referred to as radial, while -components are called azimuthal.
The non-ideal MHD equations are formulated for the gas density , its velocity , total pressure  (which includes the thermal, magnetic and cosmic-ray contributions), magnetic field  and its vector potential , and the energy density of cosmic rays  cr .The initial conditions represent an unstable magneto-hydrostatic equilibrium, and the corresponding distributions  0 ,  0 and  cr,0 in  are maintained throughout the simulation as a background state.We solve for the deviations from them, denoted  for the density,  for the velocity,  for the total pressure,  for the magnetic field and  for its vector potential, and  cr and  for the cosmic-ray energy density and flux.Cosmic rays are described in the fluid approximation with non-Fickian diffusion, so we have separate equations for their energy density and flux.The governing equations are solved numerically in a rectangular shearing box of the size 4×4×3.5 kpc 3 along the ,  and  axes, respectively, with the mid-plane at  = 0 and || ≤ 1.75 kpc.The boundary conditions are periodic in , sliding-periodic in  and allow for a free exchange of matter through the top and bottom of the domain as specified in detail by Tharakkal et al. (2022a).
The total velocity is given by  =  0 + , where  0 =  ŷ is the mean rotation velocity in the rotating frame with the shear rate  =  dΩ/d, and  is the deviation from this, associated with the instability.For a solid-body rotation,  = 0, we have  0 = 0.Both  and Ω are assumed to be independent of  and  < 0 for realistic galactic rotation profiles.We neglect the vertical gradient of Ω and ; for its observed magnitude of order −15-25 km kpc −1 (Section 10.2.3 of Shukurov & Subramanian 2021, and references therein), Ω and  only vary by about 10-15 per cent within || 1.5 kpc.
The presence of rotation only affects the momentum and induction equations, so equations ( 1), ( 4)-( 6), ( 9) and ( 10) for the mass conservation and cosmic rays of Tharakkal et al. (2022a) still apply and only the momentum and induction equations are augmented with terms containing Ω and : where D/D = / + ( 0 + ) • ∇ is the Lagrangian derivative,  is the gravitational acceleration and  is the viscous stress tensor.The Kepler gauge for the vector potential, as described by Oishi & Mac Low (2011) (see also Brandenburg et al. 1995), is appropriate for this shearing box framework.We use the gravity field  = −() ẑ obtained by Kuĳken & Gilmore (1989) for the Solar vicinity of the Milky Way and consider an isothermal gas with the sound speed  s = 18 km s −1 and temperature  = 3.2 × 10 4 K.In the background state (identified with the subscript zero, this is also the initial state), both the magnetic and cosmic ray pressures are adopted to be half the thermal pressure,  m,0 / th,0 =  cr,0 / th,0 = 0.5, where  th,0 =  2 s  0 (0),  m,0 =  2 0 (0)/(8) and  cr,0 =  cr0 (0)/3 are the thermal, magnetic and cosmic ray pressures, respectively, and  0 (0) = 5 µG.The gas viscosity  (included in ) and magnetic diffusivity  are chosen as  = 0.1 kpc km s −1 and  = 0.03 kpc km s −1 , respectively, to be somewhat smaller than the turbulent values in the ISM (see Tharakkal et al. 2022a, for further details and justification).
Table 1 presents the simulation runs discussed in this paper.The value of Ω near the Sun is close to 30 km s −1 kpc −1 (referred to as the nominal value hereafter), and  = −Ω when the rotational speed is independent of the galactocentric distance (a flat rotation curve), |×| = const.Model Ω00N is identical to Model Sim6 of Tharakkal et al. (2022a), Model Ω30N only differs by the solid-body rotation at the nominal angular velocity, Model Ω30S adds the large-scale velocity shear (differential rotation), whereas Model Ω60S has both the angular velocity and its shear doubled.The averages at  = const (horizontal averages) are denoted • • • h .
Figure 1 presents a pictorial summary of the changes in the magnetic field and gas density as the instability develops through its linear stage and then saturates in Model Ω30S.During the linear phase, at  = 0.3 Gyr, the magnetic field and gas density retain the structure of the imposed fields with weak perturbations in .By the weakly nonlinear stage at  = 0.6 Gyr, both the gas density and magnetic field are strongly perturbed to the extent that the mean azimuthal magnetic field   h starts reversing.The reversal is complete in the late nonlinear stage at  = 1.6 Gyr and magnetic loops are prominent.We explain and detail these processes below.

THE LINEAR INSTABILITY
The linear phase of the Parker instability in the absence of rotation is discussed in detail in our previous work (Tharakkal et al. 2022a), where we compare the growth rate and the spatial structure The evolution of the gas density and magnetic field in Model Ω30S is illustrated for its three significant epochs: (a) the linear stage, (b) beginning of the magnetic field reversal in the early nonlinear stage and (c) the advanced nonlinear state (the specific simulation times are indicated for each frame).Selections of magnetic lines are shown (with colour representing the local magnetic field strength in µG) in the ( , , )-space at the time indicated to the left of each frame.The horizontal average of the azimuthal magnetic field   h in µG is shown with colour on the vertical (, )-plane as it evolves continuously (rather than at discrete times used for the magnetic lines).The gas density distribution is shown with colour on the vertical ( , )-planes (in g cm −3 ) for each time.
of the most rapidly growing mode with those obtained in a range of analytical and numerical models.In this section, we focus on the modifications of the exponentially growing perturbations caused by the rotation and velocity shear.
Figures 2a,b show the evolution (in both the linear and nonlinear stages) of the root-mean-square (r.m.s.) magnitudes of the perturbations in the magnetic field and velocity, while Panels (c) and (d) show how the total magnetic field strength  r.m.s. and the mean cosmic ray energy density  cr at  = 0, respectively, evolve in the models of Table 1.As expected (Shu 1974;Zweibel & Kulsrud 1975;Foglizzo & Tagger 1994, 1995;Hanasz & Lesch 1997), the instability growth rate Γ (given in Table 1) decreases systematically with the angular velocity.The stretching of the magnetic lines along the radial () direction by the Coriolis force enhances the magnetic tension thus opposing the instability while the differential rotation shears the perturbations to reduce the radial wavelength also suppressing the instability (Foglizzo & Tagger 1994).
The spatial structure of the unstable modes is illustrated in Fig. 3, which presents the two-dimensional power spectra of the perturbations affected by the solid-body (c-d) and differential (e-f) rotation and compares them with the non-rotating case (a-b).The spectra of the velocity and magnetic field perturbations are identical when Ω = 0 but noticeable differences develop in rotating systems.In agreement with the analysis of Shu (1974), the dominant azimuthal wave number   decreases under the influence of rotation.The solidbody rotation leads to wider spectra in the radial and azimuthal wave numbers, consistent with the weaker variation of the instability growth rate with   in a rotating system (Fig. 1 of Foglizzo & Tagger 1994).Since the Coriolis force couples the radial and azimuthal motions, the spectra in   and   are more similar to each other than in the case Ω = 0.However, the velocity shear strongly reduces the range of   while the perturbations have significantly larger radial wave numbers   than in the cases Ω = 0 and  = 0.

THE SATURATED STATE
Figure 2 also shows that the nonlinear development of the instability and its statistically steady state are strongly affected by the rotation and velocity shear.Solid-body rotation does not affect much the magnitude of the magnetic field perturbations at  1 Gyr, presented with the solid and dash-dotted curves in Panel (a), but reduces the velocity perturbations shown in Panel (b).Understandably, the velocity shear enhances both (the dotted curves) by stretching the radial magnetic fields which, in turn, affect the motions.The case of faster rotation and correspondingly stronger shear confirms this tendency (dashed curves).
Panels (c) and (d) of Fig. 2, which show the total magnetic field strength and cosmic ray energy density at  = 0, suggest that the structure of the magnetic field is changed profoundly by rotation and, especially, by the velocity shear.For example, the magnitude of the magnetic field perturbations in Model Ω30S shown with the dotted curve in Panel (a) is less than twice larger than at Ω = 0 (solid curve), but the total magnetic field at  = 0 shown in Panel (c) is almost an order of magnitude stronger since the perturbation is better localised near  = 0 (see below).The instability still removes both the magnetic field and cosmic rays from the system as in the case Ω = 0, but at a much lower efficiency that depends on both the angular velocity and the rotational shear.
As compared to the case Ω = 0, the system retains stronger magnetic field under the solid-body rotation but less cosmic rays, as shown with the solid and dash-dotted curves in Fig. 2(c,d).Figure 4 clarifies the details of the changes effected by rotation and velocity The evolution of the root-mean-square magnitudes at the midplane  = 0 of (a) the magnetic field perturbation | |, normalised to  0 (0) (the strength of the background magnetic field at  = 0), and (b) gas speed in the Models Ω00N (solid, no rotation), Ω30N (dash-dotted, solid-body rotation at the nominal Ω), Ω30S (dotted, differential rotation at the nominal Ω and ) and Ω60S (dashed, doubled Ω and ).Similarly, panels (c) and (d) show the horizontally averaged total magnetic and cosmic ray energy densities at  = 0 for those models, normalized to the respective midplane values in the background state,    (0)/ 0 (0) and  cr   (0)/ cr0 (0), respectively.Moreover, as we discuss below, the gas flow becomes helical in a rotating system (see Section 6), supporting the mean-field dynamo action.As a result, a large-scale radial magnetic field   , clearly visible in Fig. 5(d,f), emerges in a rotating system.The velocity shear changes the nonlinear state qualitatively.Firstly, the scale heights of  and  cr near the midplane are even smaller at  = 0.6-0.9Gyr in Panels (h) and (i) than at the comparable times in Panels (e) and (f).Secondly, and more importantly, the vertical profile of the magnetic field strength evolves to become more complicated at  = 1.6 Gyr in Panel (h), and the cosmic ray distribution reflects this change.The energy density of cosmic rays in Model Ω30S,  cr h (0) = 0.2 cr0 at  = 1.6 Gyr (Fig. 4i) is ten time larger than in Model Ω00N.Differential rotation helps to confine cosmic rays because it drives dynamo action generating strong horizontal magnetic field, and this slows down the escape of cosmic rays as they spread along larger distances guided by the magnetic field.
The change in the vertical profile of  h in Model Ω30S at  = 1.6 Gyr reflects the reversal of the horizontal magnetic field near the midplane discussed and explained in Section 5.

MAGNETIC FIELD REVERSAL
The reversal of the magnetic field in the nonlinear stage of the instability has been noticed earlier by a few authors (see Section 1) but our simulations identify it as a generic feature of the Parker and magnetic buoyancy instabilities in rotating systems.This process is illustrated in Fig. 5 which shows how the evolution of the large-scale (f) 0.3 Gyr 0.6 Gyr 1.6 Gyr 2.6 Gyr The evolution of the vertical profiles of the horizontally averaged and normalised gas density  h / 0 (0) (left-hand column), magnetic field strength  h / 0 (0) (middle) and cosmic ray energy density  cr h / cr0 (0) (right-hand column).First row: Model Ω00N (no rotation), second row: Model Ω30N (nominal solid-body rotation), third row: Model Ω30S (nominal rotation and shear).The times corresponding to the line styles are given in the legend of each row.Note that the direction of the mean azimuthal magnetic field   h has reversed within a certain distance of the midplane at the later times,  = 1.6 and 2.6 Gyr.
horizontal magnetic field components   h and   h depends on rotation and the velocity shear.
Figure 5a shows again (see also Tharakkal et al. 2022a, for details) that, in a non-rotating system, the azimuthal magnetic field   h decreases with time in strength and its scale height increases, while the radial field   h shown in Fig. 5b is much weaker and varies along  without any systematic pattern.Solid-body rotation causes two major changes: the azimuthal field strength (Fig. 5c) first decreases faster than without rotation but then starts growing and, at late times, is stronger than for Ω = 0.The field direction remains the same as of the imposed field,   h > 0. Meanwhile, the radial field (Fig. 5d) is, at late times, comparable in strength to   h , well-ordered and is predominantly negative,   h < 0. This change is a result of the mean-field  2 -dynamo action driven by the mean helicity of the gas flow as discussed in Section 6.
The differential rotation of Model Ω30S (Fig. 5e,f) changes the evolution even more dramatically: it drives the more efficient dynamo with stronger   h and, remarkably, exhibits a reversal of the large-scale horizontal magnetic field.The reversal starts in the weakly nonlinear phase at  = 0.5 Gyr with a rather abrupt emergence of a relatively strong positive radial magnetic field near the midplane,   h > 0. The velocity shear with  < 0 stretches the positive radial field into a negative azimuthal magnetic field, so that   h starts decreasing and reverses at  = 1.6 Gyr (Fig. 5e).The total horizontal magnetic field strength (   2 h +   2 h ) 1/2 decreases to a minimum before increasing again, as   h decreases to zero and then reemerges with the opposite direction.These changes in the large-scale magnetic field structure start near the midplane and spread to larger altitudes because of the magnetic buoyancy.

The mechanism of the reversal
To understand the process that leads to the reversal of the large-scale azimuthal magnetic field, we consider individual terms in the induction equation written for the deviation from the imposed magnetic field, Figure 6 shows, for Model Ω30S, the evolution of the mean radial and azimuthal components of the first three terms on the right-hand side of this equation, which represent the advection, stretching and compression of the corresponding magnetic field components near the midplane.The stretching terms ( • ∇)  and ( • ∇)  clearly dominate, producing a mean radial field   h > 0 during the weakly nonlinear stage, 0.6  0.8 Gyr, which decreases only slowly at later times (because of diffusion and buoyancy) while being gradually stretched by the differential rotation  < 0 into a negative azimuthal field   h , eventually leading to the reversal of the initially positive   h .This picture is very different from that for Model Ω00N, where the stretching terms in both components rapidly vanish after a negative excursion during the early nonlinear phase (see Figs 5a,b and  7).Under the solid-body rotation, a positive radial field does emerge near  = 0 in the early nonlinear stage but, without the velocity shear, this does not lead to the reversal of the azimuthal field (Fig. 5c,d).
We have analyzed various parts of the averaged stretching term ( • ∇)  h in the -component of equation ( 3) to understand which of them produces a positive radial component of the mean field.We  note that   h = 0 and then ( • ∇)  h = ( • ∇)  h .Thus, Figure 8 shows that the first two terms on the right-hand side of this equation are less significant than the third term, and that     / h > 0 at || 0.2 kpc.The term     / h also contributes to the generation of a positive   h at all .
The positive correlation between   and   /, the main driver in the generation of the positive   h , arises because of: (i) the Coriolis force; and (ii) the emergence of a local minimum of   h at the midplane produced by the buoyancy.To demonstrate this, we express   using the -component of the momentum equation ( 1) with  = −Ω, differentiate the result with respect to , multiply it by   and average to obtain where we have neglected the fluctuations in  when averaging on the left-hand side (which is justifiable since the random gas speed is subsonic) and Ψ combines all other terms: where we neglect the viscosity (represented by the viscous stress tensor ) and  2 =  2  +  2  +  2  .Figures 9a,b show vertical profiles of   h in Models Ω30N (where no reversal occurs) and Ω30S, while Fig. 9c clarifies the form of various terms in equation ( 5).The positive correlation     / h emerges because of the first term on the right-hand side as soon as magnetic buoyancy produces a local minimum of   h at  = 0 (see Fig. 9b), so that  2   / 2 is systematically positive at  = 0.Such a minimum does not develop in the case of solid-body rotation (Fig. 9a) where no reversal of   h happens.As shown in Fig. 9c, the second and third terms in equation ( 5) are smaller in magnitude than the first term near  = 0 and partially compensate each other.The correlation     / h is dominant and positive near  = 0, driving a reversal in the largescale magnetic field near the midplane which then spreads to larger || as shown in Fig. 6e,f because of the magnetic buoyancy.We stress that the minimum of   h at  = 0 can only arise at the nonlinear stage of the instability, because only then do the fluctuations   not average to zero.
4 Gyr in duration (see Fig. 5).This is already a significant fraction of the galactic lifetime; therefore, we did not extend them further to find out if a further reversals would occur at later times.However, periodic reversals occur in a similar model where the unstable magnetic field is generated by an imposed mean-field dynamo action (Y.Qazi et al. 2022, in preparation).It appears that the emergence of the local minimum of   h at  = 0 and its ensuing reversal is related to the mean-field dynamo action (which our imposed field emulates).The dynamo is driven by the mean helicity of the gas flow, and both Models Ω30N and Ω30S support this mechanism (as discussed below).However, the dynamo in Model Ω30N, which has solid-body rotation (so is an  2 -dynamo), is too weak, whereas the differential rotation of Model Ω30S enhances the dynamo enough (making it an -dynamo) to produce the reversal.In the next section, we compute and discuss the mean helicity of the gas flow and other evidence for the mean-field dynamo action in Model Ω30S.

HELICITY AND DYNAMO ACTION
In Models Ω30N, Ω30S and Ω60S, the Coriolis force causes the gas motions to become helical, and the resulting -effect produces a large-scale radial magnetic field   h (e.g., Sect.Here overbar denotes a suitable averaging, and we use the horizontal averages in our discussion, so ũ and b are understood as the deviations from the horizontal averages  h and  h , such that Figure 10 shows the evolution of the kinetic and current helicities and their variation with  obtained using the horizontal averages.As expected, both quantities have odd symmetry in  (e.g., Sect. 11.3.1 of Shukurov & Subramanian 2021).Both are weak throughout the linear phase when the instability-driven perturbations are still weak, but increase significantly in magnitude during the early nonlinear phase at about  = 0.5 Gyr.The kinetic helicity reaches its maximum magnitude near the upper and lower boundaries,  = ±1.6 kpc, during the transitional phase at  = 0.6 Gyr.At a later time,  = 1.9 Gyr, the kinetic helicity reduces to a maximum of |  k | = 340 km 2 s −2 kpc −1 at || = 1.6 kpc.At early stages of the evolution, the current helicity has local extrema Figure 11.The spatial distribution of the mean kinetic helicity  k at  = 0.7 Gyr for four imposed (initial) magnetic field strengths specified by the parameters  m,0 and  cr,0 defined in equation ( 12) and given in the legend.Among the models shown in this figure, cosmic rays are present only in Model Ω30S where ( m,0 ,  cr,0 ) = (0.5, 0.5) (dash-dotted: this is a vertical cross-section of the distribution in Fig. 10a).The vertical profiles of both kinetic and current helicities evolve in a rather complicated manner, with  k < 0 at  > 0 close to the midplane (although the magnitude is small), and  k > 0 at larger  in the case of pure magnetic buoyancy (dotted curve in Fig. 11 representing  = 0.7 Gyr).In Model Ω30S,  k < 0 at  > 0 close to the midplane just before  = 0.7 Gyr.Negative  k at  > 0 is expected from the action of the Coriolis force on the ascending and descending volume elements (Sect.7.1 of Shukurov & Subramanian 2021).However,  k > 0, as it occurs at larger  for all models presented in Fig. 11, is unexpected (see below for a discussion).
The -coefficient of the nonlinear mean-field dynamo is related to the kinetic and current helicities as (Sect.7.11.2 of Shukurov & Subramanian 2021) where, in terms of the horizontal averages, and  0 is the characteristic (correlation) time of the random flow.The relevant time scale  0 differs from the time scale of the linear instability 2/( 0   ) where  0 and   are the characteristic speed and azimuthal wave number of the most unstable mode shown in Figs 2b and Fig. 3e-f, respectively.Instead,  0 is determined by nonlinear effects and has to be measured separately.We calculate the correlation time using the time autocorrelation function  () of   (the vertical velocity   is a representative component since it is directly related to the instability), with the normalized autocorrelation function calculated as where  is the duration of the time series used to compute  ().For a given , the integral in equation ( 11) is calculated for each (, ) and the result is averaged over (, ).Thus defined, the autocorrelation function and the corresponding correlation time depend on .
Figure 12 shows the time autocorrelation of   at three values of , and the form  () = exp(−/ 0 ) provides a good fit, with the fitted values of  0 given in the legend: they vary between 18 Myr at  = 0 and 40 Myr at  = 1.5 kpc.We use the fitted  () to estimate  0 as this provides a more accurate result than the direct integration as in the definition (10).
We use  0 = 30 Myr in equations ( 9), and the results are shown in Fig. 13.The largest in magnitude values | k | ≈ 7 km s −1 are reached during the transition phase around  = 0.6 Gyr near || = 1.5 kpc, whereas | m | is at its maximum around 3 km s −1 during the nonlinear phase at  = 3.6 Gyr.The spatial structure of  k is relatively simple during the early nonlinear phase but becomes more complicated later.Closer to the midplane and at later stages of the evolution,  k > 0 at  > 0 (and  k < 0 at  < 0) as expected, and the region where  k is predominantly positive (albeit small in magnitude) extends to larger || with time (see Fig. 14 representing vertical sections of Fig. 13a).
As expected, the sign of the current helicity is opposite to that of  k at almost all  and , so that the back-reaction of the magnetic field on the flow weakens the dynamo action leading to a (statistically) steady state at  3 Gyr.
The negative sign of  k at  > 0 (corresponding to the positive kinetic helicity  k ) appears to be a specific feature of a system driven by magnetic buoyancy or another magnetically driven instability such as the magneto-rotational instability (MRI).Hanasz & Lesch (1998) argue, using a model of reconnecting magnetic flux ropes, that negative  k at  > 0 can occur in magnetic buoyancy-driven mean-field dynamos.In his analysis of the mean electromotive force produced by the magnetic buoyancy instability in its linear stage, Thelen (2000a, his Fig. 4) finds  < 0 in the unstable region of the northern hemisphere in spherical geometry (corresponding to  > 0 in our case), although the 'anomalous' sign of  k remained unnoticed (Thelen 2000b).However, Brandenburg & Schmitt (1998) find  k > 0 at  > 0 in their analysis of the -effect due to magnetic buoyancy.Brandenburg & Sokoloff (2002) find  k < 0 at  > 0 in simulations of the MRI-driven dynamos (their Section 2 and   in Figs 5, 7, 9 and 11).Kinetic helicity (and the corresponding  k ) of this 'anomalous' sign is also found in the simulations of MRI-driven dynamos of P. Dhang et al. (2023, in preparation) (K.Subramanian 2022, private communication).The origin and properties of the kinetic helicity of random flows driven by magnetic buoyancy and MRI deserves further attention.Our results indicate not only that the kinetic helicity has the anomalous sign but also that it can change in space and time.
The current helicity (Fig. 10b) and the corresponding contribution to the -effect (Fig. 13b) have the opposite signs to, and closely follow both the spatial distribution and evolution of,  k and  k respectively (although the magnetic quantities have smoother spatial distributions than the corresponding kinetic ones).This confirms that the action of the Lorentz force on the flow weakens the dynamo action as expressed by equation ( 8).Together with the removal of the large-scale magnetic field by the Parker instability, this leads to the eventual evolution of the system to the statistically steady state.
Although the gas flows that become helical are driven by the instability, no simple and obvious relation of the mean helicity to the parameters that control the strength of the instability is apparent.Figure 11 shows how the vertical profile of the kinetic helicity  k changes with the magnetic and cosmic ray pressures in the initial (imposed) state, specified in terms of their ratios to the thermal pressure at  = 0, and  cr0 = ( cr − 1) cr0 (0) where  cr = 4/3.To avoid complications associated with the cosmic rays in the system behaviour, only one model of the four illustrated in Fig. 11 contains cosmic rays (Model Ω30S discussed elsewhere in the text).The midplane strengths of the imposed magnetic field  0 (0) corresponding to  m0 = 0.5, 1 and 1.5 are 5, 7 and 9 µG, respectively.When ( m0 ,  cr0 ) = (0.5, 0), the magnetic field is too weak to be unstable and the system remains in the state of magneto-hydrostatic equilibrium,  k = 0. Adding cosmic rays, ( m0 ,  cr0 ) = (0.5, 0.5) (Model Ω30S) destabilises the system producing helical flows discussed above.Adding magnetic rather than cosmic ray pressure, ( m0 ,  cr0 ) = (1, 0), also makes the system unstable, and the resulting mean helicity at larger || is greater than for ( m0 ,  cr0 ) = (0.5, 0.5).A still stronger magnetic field, ( m0 ,  cr0 ) = (1.5, 0) leads to comparable  k the previous two cases in || 1 kpc, except near the midplane.Altogether, it is difficult to identify a clear pattern in the dependence of the magnitude and spatial distribution of the mean helicity of the gas flow driven by the Parker instability; this invites further analysis, both analytical and numerical.
The dimensionless measure of the mean-field dynamo activity in a differentially rotating gas layer is provided by the dynamo number (Section 11.2 of Shukurov & Subramanian 2021) where ℎ is the layer scale height,  is the velocity shear rate ( = −Ω in our case),  is given in equation ( 8) and is the magnetic diffusivity.The first term in this expression is the turbulent diffusivity and  is the explicit magnetic diffusivity from equation ( 2) or (3).As we use the horizontal averages in these relations,  is a function of  and varies with time together with ℎ,  and ; thus defined,  might be better called the local dynamo number, a measure of the dynamo efficiency at a given  and .In Model Ω30S,  = 0.03 kpc km s −1 while the turbulent diffusivity varies, at  = 1 Gyr, from 0.03 kpc km s −1 at  = 0 to 0.5 kpc km s −1 at  = 1 kpc (a nominal turbulent diffusivity in the ISM, where turbulence is mainly driven by supernovae, is 1 kpc km s −1 ).The dynamo amplifies a large-scale magnetic field provided || >  c , where  c is a certain critical dynamo number (see below).
Figure 15 shows how the dynamo number varies with  and .

During the transient phase, ũ2
h is relatively low while || is at its maximum.The resulting dynamo number is as large as | | 10 4 .As the system evolves into the nonlinear state, the turbulent diffusivity increases and the dynamo number reduces in magnitude.At  = 0.6 Gyr,  varies from 4 near the midplane to 6 × 10 3 at  = 1 kpc.At later times,  is larger near the midplane and reduces further in magnitude: at  = 0.9 Gyr,  = 300 near the midplane and 9 at  = 1 kpc.
As shown by Ruzmaikin et al. (1980), the -dynamo in flat geometry generates oscillatory magnetic fields for  > 0, quadrupolar for  180 and dipolar for  550.The behaviour of the largescale magnetic field in Model Ω30S is consistent with these results: it is quadrupolar and oscillatory.

RELATIVE DISTRIBUTIONS OF COSMIC RAYS AND MAGNETIC FIELD
Similar to our analysis in Tharakkal et al. (2022a), we present in Table 2 the Pearson cross-correlation coefficient between the fluctuations in energy densities for different components in model Ω30S at  = 0.5 and 1 kpc for the late nonlinear stage at  = 2.6 Gyr, derived as The only significant entry in the table is the anti-correlation between the magnetic and cosmic ray energy fluctuations at  = 1 kpc where their contribution to the total pressure is noticeable (see Section 8).
There are no signs of energy equipartition between cosmic rays and magnetic fields at kiloparsec scales; nor are there indications of equipartition at the turbulent scales, for either cosmic ray protons (Seta et al. 2018) or electrons (Tharakkal et al. 2022b).

VERTICAL FLOWS AND FORCE BALANCE
Rotation affects significantly the vertical gas flow driven by the instability.As discussed by Tharakkal et al. (2022a) (and also in Model Ω00N), a systematic gas outflow is transient without rotation and only occurs during the early nonlinear stage.Figure 16 shows the horizontally averaged vertical velocity   h in Models Ω30N (solid-body rotation) and Ω30S (differential rotation).In both cases, systematic vertical flows occur at || 1 kpc.The solid-body rotation (Fig. 16a) does not change much the structure of the flow in comparison with the non-rotating system, with a transient outflow during the early nonlinear stage and a weak inflow at later times.In Model Ω30N, the maximum outflow speed is |   h | = 9 km s −1 at  = 0.7 Gyr, followed by the inflow at the speed |   h | = 7 km s −1 at  > 1.4 Gyr.However, differential rotation not only changes dramatically the magnetic field structure and evolution (Fig. 5), but also supports a prolonged period of a systematic gas outflow at 0.6  3 Gyr, which eventually evolves into a weak gas inflow at large || (Fig. 16b).The maximum outflow speed in Model Ω30S is |   h | = 7 km s −1 at  = 0.6 Gyr at large ||, while the later inflow speed is |   h | = 1 km s −1 at  3 Gyr.The pattern of the vertical flows shown in Fig. 16b is not dissimilar to the structure of the magnetic field shown in Fig. 5e-f and the dynamo number (Fig. 15) -especially at later stages,  3 Gyr -suggesting that the magnetic field contributes noticeably to the vertical flow in Model Ω30S.
To understand what drives the vertical flows, we present in Fig. 17 the vertical forces acting during various evolutionary stages of Model Ω30S.It is instructive to compare them with those in nonrotating systems discussed by Tharakkal et al. (2022a).Without rotation, as in Model Ω00N (see also Fig. 12 of Tharakkal et al. 2022a), both magnetic and cosmic ray pressures are reduced significantly as the system evolves into the nonlinear state, and the vertical gas flows are driven by the thermal pressure gradient.This changes in Model Ω30S, where magnetic field, and to a lesser extent cosmic rays, make a stronger contribution to the force balance.Moreover, the gravity force and the thermal pressure gradient balance each other almost completely in the nonlinear state, so that the weaker magnetic and cosmic ray pressures appear to be capable of controlling the vertical velocity pattern, especially at || 0.5 kpc.This is is illustrated in Fig. 18, which shows that the vertical variations of the net vertical force per unit mass are indeed similar in detail to those of the magnetic pressure gradient.
The magnetic and cosmic ray pressure gradients are weak because both non-thermal components of the simulated ISM are much less stratified than the thermal gas.However, their energy densities are large and they dominate over the thermal gas at || 0.5-1 kpc. Figure 19 shows the vertical profiles of the horizontally averaged ratios of the magnetic and cosmic ray pressures to the thermal pressure,  m and  cr respectively, defined as in equation ( 12) but for the evolving quantities.Although each non-thermal pressure component is subdominant near the midplane at all stages of the evolution, each of them exceeds the thermal pressure at larger altitudes as soon as the instability becomes nonlinear,  0.6 Gyr.It is useful to compare Fig. 19 with Fig. 18 of Tharakkal et al. (2022a): rotation somewhat reduces the magnitudes of  m and  cr at large || but leads to the dominance of the non-thermal pressure components at smaller values of || than in a non-rotating system, and leads to a larger contribution from cosmic rays.

DISCUSSION AND CONCLUSIONS
Differential rotation affects the nonlinear state of the Parker instability more strongly than its linear properties.Without rotation, the system loses most of its magnetic field and cosmic rays as it evolves towards the steady state.A solid-body rotation does not change the nonlinear state significantly.However, differential rotation allows the system to retain better both the magnetic field and cosmic rays.The reason for that is the dynamo action (present also under the solid-body rotation but significantly enhanced by the differential rotation) which produces strong (about 2-3 µG) large-scale magnetic field both near the midplane and at large altitudes.As a result, cosmic rays (governed by anisotropic diffusion) spend longer times within the system.
The systematic vertical gas flows are also affected by the rotation, which prolongs the transient outflow at a speed |   h | = 7 km s −1 to the time interval 0.6  3 Gyr.It appears that the magnetic field contributes significantly to driving the outflow.Meanwhile, cosmic rays do not play any significant role in driving the outflow at the scales explored here, || ≤ 1.5 kpc: because of the large diffusivity of cosmic rays, the vertical gradient of their pressure is very small.
Another dramatic effect of the dynamo action is that it leads to a reversal of the large-scale magnetic field, in what appears to be a sign of nonlinear oscillations of the large-scale magnetic field.Neither the Parker instability nor the dynamo are oscillatory by themselves.We have identified the rather subtle mechanism of the reversal and argue that it is an essentially nonlinear phenomenon.
The reversal of the large-scale magnetic field is also reflected in its spatial distribution.The reversal starts near the midplane and then the reversed magnetic field spreads to larger altitudes (see Fig. 5e-f).As a result, the direction of the large-scale magnetic field reverses along  at any given time.An arguably similar pattern of regions with the sign of the Faraday depth alternating along the direction perpendicular to the disc plane is observed in the edge-on galaxy NGC 4631 (Mora-Partiarroyo et al. 2019).The comparison of Figs 5e-f and 5c-d shows that the Parker instability in a dynamo active system can produce rather complicated magnetic field structures.Our use of horizontal averages in Fig. 5 and elsewhere in the text conceals strong localised vertical magnetic fields typical of the magnetic buoyancy (see, e.g., and  cr , respectively, at various times specified in the legend: the linear state,  = 0.3 Gyr (solid), transitional period,  = 0.6 Gyr (dotted), nonlinear state at  = 1.6 Gyr when the magnetic field reversal occurs (dashed) and a late nonlinear state,  = 2.6 Gyr (dash-dotted).
Figure1.The evolution of the gas density and magnetic field in Model Ω30S is illustrated for its three significant epochs: (a) the linear stage, (b) beginning of the magnetic field reversal in the early nonlinear stage and (c) the advanced nonlinear state (the specific simulation times are indicated for each frame).Selections of magnetic lines are shown (with colour representing the local magnetic field strength in µG) in the ( , , )-space at the time indicated to the left of each frame.The horizontal average of the azimuthal magnetic field   h in µG is shown with colour on the vertical (, )-plane as it evolves continuously (rather than at discrete times used for the magnetic lines).The gas density distribution is shown with colour on the vertical ( , )-planes (in g cm −3 ) for each time.

Figure 5 .
Figure 5.The evolution of the horizontally averaged magnetic field components,   h (left-hand column) and   h (right-hand column) in Models Ω00N (a-b),Ω30N (c-d) and Ω30S (e-f).For Ω30S the mean azimuthal field   h decreases after  = 0.6 Gyr, and undergoes a reversal in sign at  ≈ 1.6 Gyr, with the reversal then spreading to higher altitudes.Meanwhile, the mean radial field   h becomes positive and relatively strong near  = 0 rather abruptly at  ≈ 0.5 Gyr and then also spreads away from the midplane.
7.1 of Shukurov & Subramanian 2021).Differential rotation (in Models Ω30S and Ω60S) enhances the dynamo significantly, and we have discovered that this leads to a reversal in the azimuthal magnetic field direction discussed in Section 5.Both types of the turbulent dynamo ( 2 dynamo in Ω30N and  in Ω30S and Ω60S) are driven by the mean kinetic helicity of the gas flow  k = ũ • (∇ × ũ), and the current helicity of the magnetic fluctuations  m = b • (∇ × b) opposes the dynamo instability leading to a reduction of the -coefficient until a steady state is achieved (e.g., Sect.7.11 of Shukurov & Subramanian 2021).

Figure 12 .
Figure 12.The time autocorrelation function of the vertical velocity component, equation (11), for 0 ≤  ≤ 2 Gyr (with the minimum time lag of 10 Myr) at  = 0 (solid), 0.6 (dashed) and 1 kpc (dotted) in Model Ω30S.The correlation time  0 at each  is given in the legend, obtained from the fits of the form  ( ) = exp(−/ 0 ), shown with dotted curves.close to the midplane, where the magnetic field is stronger, |  m | = | b • (∇ × b) h | = 89 µG 2 kpc −1 at  = 0.6 Gyr, || = 0.1 kpc.The extrema move away from the midplane in the nonlinear stage, to reach |  m | = 7 µG 2 kpc −1 at  = 1.2 Gyr, || = 0.5 kpc and |  m | = 5 µG 2 kpc −1 at  = 3 Gyr, || = 1 kpc.The vertical profiles of both kinetic and current helicities evolve in a rather complicated manner, with  k < 0 at  > 0 close to the midplane (although the magnitude is small), and  k > 0 at larger  in the case of pure magnetic buoyancy (dotted curve in Fig.11representing  = 0.7 Gyr).In Model Ω30S,  k < 0 at  > 0 close to the midplane just before  = 0.7 Gyr.Negative  k at  > 0 is expected from the action of the Coriolis force on the ascending and descending volume elements (Sect.7.1 of Shukurov & Subramanian 2021).However,  k > 0, as it occurs at larger  for all models presented in Fig.11, is unexpected (see below for a discussion).The -coefficient of the nonlinear mean-field dynamo is related to the kinetic and current helicities as(Sect.7.11.2 of Shukurov &  Subramanian 2021)

Figure 13 .
Figure 13.The evolution of (a)  k and (b)  m , given in equations (9), in Model Ω30S.

Figure 15 .
Figure 15.The evolution and vertical variation of the dynamo number of equation (13) in Model Ω30S.

Figure 16 .
Figure 16.The evolution and variation with  of the horizontally averaged vertical velocity   h in Models (a) Ω30N and (b) Ω30S.

Figure 19 .
Figure19.The distribution in  of the horizontally averaged ratios of (a) magnetic and (b) cosmic ray pressures to the thermal pressure in Model Ω30S,  m and  cr , respectively, at various times specified in the legend: the linear state,  = 0.3 Gyr (solid), transitional period,  = 0.6 Gyr (dotted), nonlinear state at  = 1.6 Gyr when the magnetic field reversal occurs (dashed) and a late nonlinear state,  = 2.6 Gyr (dash-dotted).

Table 2 .
The cross-correlation coefficient  of the fluctuations in various energy densities in the statistically steady state of Model Ω30S at  = 2.6 Gyr presented as , , where  and  refer to  = 0.5 and 1 kpc, respectively.