Nonlinear magnetic buoyancy instability and turbulent dynamo

Stratified disks with strong horizontal magnetic fields are susceptible to magnetic buoyancy instability (MBI). Modifying the magnetic field and gas distributions, this can play an important role in galactic evolution. The MBI and the Parker instability, in which cosmic rays exacerbate MBI, are often studied using an imposed magnetic field. However, in galaxies and accretion discs, the magnetic field is continuously replenished by a large-scale dynamo action. Using non-ideal MHD equations, we model a section of the galactic disc (we neglect rotation and cosmic rays considered elsewhere), in which the large-scale field is generated by an imposed $\alpha$-effect of variable intensity to explore the interplay between dynamo instability and MBI. The system evolves through three distinct phases: the linear (kinematic) dynamo stage, the onset of linear MBI when the magnetic field becomes sufficiently strong and the nonlinear, statistically steady state. Nonlinear effects associated with the MBI introduce oscillations which do not occur when the field is produced by the dynamo alone. The MBI initially accelerates the magnetic field amplification but the growth is quenched by the vertical motions produced by MBI. We construct a 1D model, which replicates all significant features of 3D simulations to confirm that magnetic buoyancy alone can quench the dynamo and is responsible for the magnetic field oscillations. Unlike with an imposed magnetic field (arXiv:2305.03318,arXiv:2212.03215), the nonlinear interactions do not reduce the gas scale height, so the consequences of the magnetic buoyancy depend on how the magnetic field is maintained.


INTRODUCTION
Magnetic fields are present in many astrophysical systems, often exhibiting highly random structure driven by turbulent flows.It is also very common to find magnetic fields coherently organised on large scales, partially aligned to large-scale shear flows in rotating objects such as spiral galaxies, convecting stars or accretion disks.The large-scale magnetic fields are produced and maintained by dynamo processes which amplify weak seed magnetic fields (Chapter 9 of Shukurov & Subramanian 2021) to generate dynamically significant large-scale magnetic fields (Moffatt 1978;Parker 1979;Krause & Rädler 1980).
The linear (kinematic) stage of the magnetic field amplification, when the Lorentz force is negligible and the magnetic fields still do not affect plasma flows, is well understood but the nonlinear behaviours of magnetic fields and plasma flows remain a subject of active research (Brandenburg & Subramanian 2005;Shukurov & Subramanian 2021).Magnetic fields produced by the mean-field dynamo are subject to further instabilities, and the magnetic buoyancy instability (MBI) is one of the strongest and most ubiquitous of them.The dynamo and magnetic buoyancy instabilities are most often studied separately although it is understood that they can af-★ E-mail: Y.Qazi@newcastle.ac.uk fect each other.In particular, plasma flows produced by the MBI become helical under the action of the Coriolis force, which leads to the mean-field dynamo action (Parker 1992;Hanasz & Lesch 1998;Moss et al. 1999;Thelen 2000;Tharakkal et al. 2023a) in addition to the dynamo associated with the turbulence driven in galaxies by supernova explosions.As we discuss below, the dynamo action also affects the MBI.It is difficult to distinguish the effects of the mean-field dynamo and magnetic buoyancy as they become strongly interdependent in the nonlinear state.
The MBI is expected to develop in any magnetized, strongly stratified astrophysical object and has a relatively short time scale of the order of the sound or Alfvén crossing time based on the magnetic field scale height, of order 10 7 yr in the Solar neighbourhood of the Milky Way (Section 2.8.2 of Shukurov & Subramanian 2021, and references therein).In spiral galaxies, the MBI enhanced by cosmic rays (which produce significant pressure but have negligible weight) is known as the Parker instability (Parker 1966).
The linear stage of the MBI has been extensively explored with and without cosmic rays but its nonlinear states have attracted less attention.The Parker instability has been studied by Hanasz & Lesch (1998), Hanasz et al. (2002Hanasz et al. ( , 2004) ) and Kosiński & Hanasz (2007) with emphasis on the associated mean-field dynamo action, including the nonlinear steady state.In application to accretion discs, the role of the MBI and its effect on the magnetic field structure is discussed by Johansen & Levin (2008) and Gaburov et al. (2012).
More recently, the nonlinear magnetic buoyancy and Parker instabilities of an imposed azimuthal magnetic field has been examined by Tharakkal et al. (2023b,a).They show that magnetic buoyancy in a non-rotating system leads to a substantial redistribution of gas, magnetic field and cosmic rays, resulting in large (a few kiloparsecs) scale heights of the magnetic field and cosmic rays and a correspondingly small (of order a hundred parsec) scale height of the gas, which is supported against gravity mainly by thermal pressure.Rotation leads to somewhat smaller scale heights of the nonthermal constituents (and their larger contribution to the pressure gradient near the midplane).More importantly, the nonlinear MBI in a rotating system can completely change the initial magnetic field even leading to its reversal.
A natural question to ask then is whether a hydromagnetic dynamo can sustain the magnetic field against these effects and how this affects the evolution of the instability.How might this change our understanding of galactic dynamos and the structure of galactic magnetic fields?
To address such questions, we study the nonlinear states of the mean-field dynamo in a disc, paying attention to the effects of the MBI.We explore a range of models with various relations between the dynamo and MBI time scales.The mean-field dynamo action is introduced via an explicit -effect in the induction equation, assumed to be driven by background turbulent flows.Anticipating that rotation would introduce significant qualitative complexity (particularly the enhanced dynamo action arising from the motions induced by the magnetic buoyancy), we commence with a simpler model in which we neglect the overall rotation.Although the -effect emerges physically only in a rotating system, the specific effects of rotation is our subject for a subsequent paper.

BASIC EQUATIONS AND NUMERICAL METHODS
We model isothermal gas and magnetic field in a section of a within a Cartesian box, with ,  and  representing the radial, azimuthal and vertical directions, respectively.The simulation domain extends 4 kpc in each horizontal direction and 3 kpc vertically, centered at the galactic midplane.We have tested computational boxes of various sizes from 0.5 to 8 kpc to confirm that we capture all essential features of the evolving system.The grid resolution is 256×256×384 mesh points, yielding the grid spacing about 15.6 pc along each dimension.The size of the domain is larger than the expected vertical and horizontal scales of the instability, and the resolution is sufficient to obtain convergent solutions to model the instability (Tharakkal et al. 2023b).

Basic equations
We solve the system of isothermal non-ideal compressible MHD equations using the sixth-order in space and third-order in time finitedifference Pencil Code (Brandenburg & Dobler 2002;Pencil Code where D/D = / + u • ∇,  and u are the gas density and velocity, A is the magnetic vector potential, with B = ∇ × A for the magnetic field (thus, ∇ • B ≡ 0), and  s is the speed of sound assumed to be constant.equation (3) includes the term B responsible for the mean-field dynamo action, where  is specified in Section 2.3.The key parameters used are shown in Table 1.
The traceless rate of strain tensor τ has the form τ   = 1 2 (    +     ) − 1 3        (where   = /  and summation over repeated indices is understood).The shock viscosity   =  shock  shock in equation ( 2), with  shock ∝ |∇ • u − | (where  shock is a constant given in Table 2 and ∇ • u − is the divergence of the velocity field where it is negative and vanishes otherwise) regularises shock fronts propagating perpendicular to steep pressure gradients, and differs from zero only in regions of converging flow (see Gent et al. 2020, for details).Following Gent et al. (2020), we also include the terms with   in equations ( 1) and (2) to ensure the momentum conservation.Similarly to   , the shock-capturing term in the continuity equation has   =  shock  shock , where  shock is a constant given in Table 2.The hyper diffusion with the coefficients  6 and  6 is used to resolve grid-scale numerical instabilities, with τ (5) (Brandenburg & Sarson 2002;Gent et al. 2021).

Boundary conditions
The boundary conditions in both horizontal directions  and  are periodic for all variables.To prevent an artificial inward advection of the magnetic energy through the top and bottom of the domain at  = ±1.5 kpc, we impose the conditions   =   = 0 and   / = 0, i.e.,    / =    / =   = 0 at || = 1.5 kpc . (4) The boundary conditions for the velocity field are stress-free,   / =   / = 0 at || = 1.5 kpc .
(5) To permit vertical gas flow across the boundaries without exciting numerical instabilities, the boundary condition for   imposes the boundary outflow speed across the ghost zones outside the domain whereas an inflow speed at the boundary tends smoothly to zero across the ghost zones (Gent et al. 2013).The density gradient is kept at a constant level at the boundaries, with the scale height intermediate between that of the Lockman layer and the galactic halo, and we note that the value of the scale height adopted has a negligible effect on the results.
In this paper we neglect the effects of rotation and velocity shear and excite an  2 -dynamo by imposing an explicit -dependent effect in the induction equation, The vertical extent of the dynamo-active layer is ℎ  on each side of the midplane; the smaller is ℎ  , the stronger is the vertical gradient of the magnetic field and the more it is buoyant.The dynamo intensity (in particular, the growth rate of the large-scale magnetic field) depends on the dimensionless number In the Solar neighbourhood of the Milky Way,  0 ≃ 0.5 km s −1 and ℎ  ≃ 0.5 kpc (e.g., p. 317 of Shukurov & Subramanian 2021).We consider a range of values for these two parameters, as summarised in Tables 1, 2 and A1.The purpose of the form of () adopted at || > ℎ  /2 is to ensure that |()| decreases smoothly from its maximum || =  0 at || = ℎ  /2 to negligible values at || ≫ ℎ  since the -effect dynamos are sensitive to discontinuities in ∇ which can cause unphysical oscillations of the magnetic field.Figure 1 shows the effective (local) value of R defined as in equation ( 9) but with |()| of equation ( 8) rather than  0 , confirming that the imposed dynamo action is confined to the layer || ≲ ℎ  where   >  ,c , and we have determined that  ,c = 3.8 in our case.
The resulting form of the mean-field dynamo is known as the  2 -dynamo; it is weaker than the -dynamo that involves the differential rotation in galaxies but our goal here is to explore the interaction of the dynamo and MBI rather than to reproduce in detail the galactic conditions.The  2 -dynamo is a relatively simple and wellunderstood mean-field dynamo mechanism; the temporal and spatial scales of the magnetic field produced by  2 -dynamo are known accurately (e.g., Sokoloff et al. 1983).This will help us to distinguish the effects of the dynamo action and the MBI.
The initial conditions represent isothermal hydrostatic equilibrium with the gas density profile given a negligible initial magnetic field.The sound speed  s = 15 km s −1 corresponds to the temperature  = 2 × 10 4 K for 50 per cent of hydrogen ionized, so that the mean molecular weight is 0.74 (assuming 70% of hydrogen and 30% of neutral helium by mass).The resulting initial density scale height | ln /| −1 is about 400 pc at a distance from the midplane.
The seed magnetic field applied as an initial condition represents as Gaussian random noise in the vector potential component   with a mean amplitude proportional to  1/2 () and the maximum strength 10 −6 µG at  = 0.This field has   = 0.A random initial magnetic field leads to shorter transients than in the case of a unidirectional initial field.

MEAN-FIELD DYNAMO TRANSFORMED BY MAGNETIC BUOYANCY
The system that we explore supports both the mean-field dynamo and magnetic buoyancy instabilities.They can be distinguished clearly only during their linear stages when the corresponding perturbations grow exponentially.However, the effects of the dynamo action and magnetic buoyancy are strongly intertwined at the later stages when the Lorentz force becomes dynamically significant and the system evolves into a stationary state.It is, correspondingly, harder to identify the role of each individual instability in the nonlinear behaviour.To facilitate such an identification, we use a range of models with various relations between the time scales of the dynamo action and the MBI.The mean-field dynamo driven by the MBI has been explored earlier (see references in Section 1) but the effects of the MBI on the dynamo action driven in spiral galaxies by the turbulence in the gas disc have received little attention.Our approach is designed to address this aspect of the evolution of magnetic fields in galaxies and accretion discs.The effects we find are rather unexpected: the buoyancy of the magnetic fields produced by the dynamo does not just modify the dynamo but can change entirely the spatial structure and evolution of the large-scale magnetic field.
The parameters of the models explored are presented in Tables 2  and A1.Since the dynamo action is relatively insensitive to the kinematic viscosity  while the MBI responds to both viscosity and magnetic diffusivity , we consider systems with  >  (Table 2)  Table 2. Simulation runs with the kinematic viscosity  = 0.008 kpc km s −1 , magnetic diffusivity  = 0.03 kpc km s −1 and various values of   of equation (9) based on the values of  0 and ℎ  that appear in equation ( 8);  D is the rate of exponential growth of the magnetic field strength during the linear phase of the dynamo,  B is the growth rate of the magnetic field obtained when the magnetic field becomes buoyant and   is the growth rate of the root-mean-square gas speed.A missing entry means that the corresponding variable does not grow.The last column shows the period  of nonlinear oscillations of the azimuthal magnetic field (the -component), a missing entry in this column means that the oscillation period was not measured, in particular, when it is very long.and  <  (Table A1).The results and conclusions presented below apply in both cases.
The growth rate of the most rapidly growing magnetic field mode under the  2 -dynamo action with  depending on  alone is given by (Section 4.iii of Sokoloff et al. 1983) for large and moderate values of   .As shown in Fig. 2, the growth  2.
rates obtained for the early stages of the magnetic field growth fit this dependence with high accuracy.The resultant growth rates from equation ( 11) are positive when   >  ,c = 2.The threedimensional numerical solutions have a somewhat larger critical value  ,c ≈ 3.8 (obtained by interpolating the data points shown in Fig. 2 using the cubic spline for 3 ≤   ≤ 7 ) and smaller  D because of the magnetic diffusion in the -and -directions.The scale over which the magnetic field strength in this mode decreases by the factor e is given by , and the spatial structure is oscillatory with the wavelength  D = 4ℎ  /  .For   ≃ 10, we have  D ≃  D ≃ ℎ  .The dynamo also supports a wide range of modes at scales larger than about  D /2 but they grow at rates lower than that of equation ( 11).Meanwhile, the scale of the MBI along the magnetic field is 1-2 kpc, significantly larger than the values of ℎ  that we use.We consider a wide range of   to ensure that the magnetic fields produced by the dynamo and the MBI can have sufficiently different growth rates and spatial scales to be distinguishable in their linear stages.
Our simulations start with a weak random magnetic field that launches the dynamo action.The growth rate of the MBI is of the order of  A /ℎ  (with  A the Alfvén speed), so there is a relatively long period when the magnetic field is too weak to be unstable because of its buoyancy.The growth rate of the MBI exceeds that of the dynamo when  A ≳  0 /4; this happens when the magnetic field is relatively strong.Depending on the system parameters, the second phase of the system evolution when the magnetic field generated by the dynamo becomes unstable with respect to the MBI may or may not be clearly discernible.
Figure 3 presents Model R10h2, in which the growth rates and characteristic scales of the mean-field dynamo and MBI differ noticeably. D ≃ 0.13 kpc for the scale at which magnetic field strength decreases and  D ≃ 0.25 kpc for the wavelength at which the field direction changes.The evolving spatial structure of the magnetic field is illustrated in Fig. 4. At early times after the short-term transients from the initial condition have dissipated, the magnetic field is weak, has a relatively small scale imposed by the dynamo and is confined to a thin disc || ≲ ℎ  where the -effect is imposed.At a later time,  = 1.4 Gyr, the magnetic field spreads out of the disc because of the magnetic diffusion and the first signs of the magnetic buoyancy emerge, particularly manifested in the enhancement of   (the right-hand column of Fig. 4).At that time, the growth rate of the magnetic field increases (Fig. 3).The effect of the magnetic buoyancy becomes even more significant at  = 1.5 Gyr -at this stage, the exponential field amplification typical of a linear instability continues at the enhanced rate.In the strongly nonlinear stages at  > 1.5 Gyr, the magnetic field structure characteristic of the mean-field dynamo is completely wiped away despite the continued action of the -effect.
To illustrate the three-dimensional magnetic field structure, including the magnetic loops produced by the MBI and the magnetic field generated by the dynamo transformed by the MBI, Fig. 5 shows the evolving three-dimensional structure of magnetic lines at various scales.The total magnetic field B in this figure and in Fig. 3 is separated into the contributions B B of the larger scales (characteristic of the MBI) and of the smaller ones B D (driven by the imposed dynamo action) using the Gaussian smoothing, where the integration extends over the whole domain volume with the smoothing kernel  ℓ (ξ) = (2ℓ 2 ) −3/2 exp[−|ξ| 2 /(2ℓ 2 )] with ℓ = 200 pc.The smoothing scale is chosen to be equal to   ≃   ≃ ℎ  , the scale of the spatial variations of the magnetic field generated by the dynamo as estimated above, and ℎ  = 200 pc in Fig. 5.We note that the small-scale part B D also contains random magnetic fields produced by nonlinear effects at the later stages of the system's evolution.
The restructuring of the magnetic field by the MBI is quantified in Fig. 6, which shows the two-dimensional power spectra of the -component of the magnetic field at the times indicated above each frame.These confirm the evolution pattern visible in Fig. 4.Over time the dominant horizontal scales 2 −1  and 2 −1  of the magnetic field grow larger: at  ≲ 0.5 Gyr, the spectrum centres around   ≈   ≈ 2/ D ≈ 15 kpc −1 , which reduces by  = 2 Gyr to   ≃   ≲ 5 kpc −1 .The latter wavenumbers correspond to scales in excess of 1-2 kpc characteristic of the MBI.As the peak wavenumbers decrease their spread is broader as the MBI excites a wider range of unstable modes.Typical of the  2 -dynamo, the magnetic field components   and   are of about equal strength and scale (velocity shear due to differential rotation would enhance the azimuthal component   in comparison with   and modify the field scales to   <   ).
Figure 7 shows that the distribution in  of the horizontally averaged gas density varies very little with time, only becoming slightly wider in the non-linear phase as the magnetic buoyancy carries matter away from the midplane.This behaviour is quite different than in the simulations of the MBI and Parker instability of an imposed magnetic field.There the scale height of the magnetic field (and cosmic rays) increases in the nonlinear state so strongly that the gas layer loses the magnetic and cosmic-ray pressure support as the system evolves and becomes much thinner, as shown in Fig. 11d of Tharakkal et al. (2023b) and Fig. 4 of Tharakkal et al. (2023a).Remarkably and unexpectedly, the nonlinear states of the MBI (and the Parker instability) are sensitive to the mechanism by which the magnetic field is maintained (imposed versus mean-field dynamo in our case).

NONLINEAR OSCILLATIONS
Figure 8 illustrates an even more fundamental consequence of the MBI: the magnetic field, which grows monotonically at early times, develops oscillations at  ≳ 1.5 Gyr when it becomes strong enough to make the system essentially nonlinear.The figure shows the evolution of the horizontally averaged magnetic field components ⟨  ⟩   and ⟨  ⟩   from Model R10h2, normalised to their maximum values at each time to better expose the field structure at early times when it is still weak.
The magnetic field generated by the linear dynamo is confined to a relatively thin layer || ≲ ℎ  and grows monotonically.However, it spreads to larger altitudes because of the buoyancy (to achieve the scale height of order 1 kpc).When fully nonlinear the magnetic field becomes oscillatory, reversing its direction at intervals of order 0.5 Gyr.There is a phase shift in Model R10h2of around a quarter of a period between the -and -components of the magnetic field.The period of the oscillations is presented in the last column of Table 2.
The field scale height, the timing of the onset of the oscillations and their period vary with ℎ  , but the oscillations are a general property of the models with sufficiently strong magnetic fields.
Figure 9 shows the variation of the oscillation period with ℎ  in both the three-dimensional simulations presented in Table 2 and the 1D model of Section 4.1 for  0 = 1.5 km s −1 .Table 2 shows that  2) and the one-dimensional model of Section 4.1 with  0 = 1.5 km s −1 (crosses).The fit to the three-dimensional simulation results is shown dash-dotted its form is specified in the legend with ℎ  in kiloparsecs.
the period  does not exhibit any significant variation with  0 for a fixed ℎ  , but does change systematically with ℎ  for a fixed   , the dimensionless measure of the kinematic dynamo intensity.This is understandable since the oscillations occur in the nonlinear stage which is not very sensitive to the strength of the kinematic dynamo.The increase of the oscillation period with ℎ  can be attributed to the change in the characteristic time of the system, the sound (or Alfvén) crossing time ℎ  / s if the oscillations are driven by the magnetic buoyancy or the magnetic diffusion time ℎ 2  / if they are a nonlinear dynamo phenomenon.When normalised to the magnetic diffusion time ℎ 2  /, the oscillation period can be approximated for all the models explored (either three-or one-dimensional) as (shown dash-dotted in Fig. 9) The residual variation, ℎ 0.35  , is weak enough to suggest that the oscillations are mainly driven by the dynamo action triggered by the magnetic buoyancy as discussed in Section 5.
Magnetic field reversals, which are likely to be a manifestation of similar nonlinear oscillations, also occur in the simulations of the MBI of an imposed magnetic field (Tharakkal et al. 2023a).Johansen & Levin (2008) and Gaburov et al. (2012) also observe a reversal of the large-scale magnetic field in their simulations of the magnetic buoyancy in accretion discs.As we argue in Section 5, they are produced by the mean-field dynamo action at || ≳ ℎ  associated with the helical gas flows driven by the Lorentz force.

One-dimensional model of nonlinear oscillations
The kinematic stage of the  2 -dynamo at  ≲ 1.3 Gyr in Model R10h2 is non-oscillatory.Indeed, the kinematic  2 -dynamo is known to develop an oscillatory solution only of a long period (Sokoloff et al. 1983;Shukurov et al. 1985;Rädler & Bräuer 1987) or because of the boundary effects (Baryshnikova & Shukurov 1987).In this section, we propose a nonlinear one-dimensional model of the  2 -dynamo with advection due to the magnetic buoyancy and demonstrate that it not only admits oscillatory solutions but reproduces quantitatively and in fine detail the oscillatory behavior of the magnetic field discussed above.
The -and -components of the mean-field dynamo equation for the mean magnetic field ⟨B⟩ and the mean velocity ⟨u⟩ ≡ U = (0, 0,   ) can be written as where we assume that all variables only depend on  and  (the infinite slab approximation) and we have omitted brackets denoting averaging to simplify the notation in this section.As follows from ∇ • B = 0,   = const in this approximation.The advection velocity   satisfies the Navier-Stokes equation where  is the vertical acceleration due to gravity and the second term on the right-hand side is the Archimedes force resulting from magnetic buoyancy.We neglect the time variation of the gas density, adopting  =  0 at all times but, in the spirit of the Boussinesq approximation, include density variation  ′ in the Archimedes force.Consider a region of density  =  0 +  ′ containing a magnetic field of a strength  +  surrounded by the gas of the density  0 with magnetic field  (here  is the mean field strength and  is its local perturbation).The pressure balance in an isothermal gas then leads to Equations ( 15)-( 17) are solved numerically in 0 <  <  0 with  0 = 1.5 kpc.  ,   and   are symmetric about  = 0 (a quadrupolar magnetic structure, evident in Fig. 8 and known to dominate in a thin layer -e.g., Section 11.3.1 of Shukurov & Subramanian 2021) with At  =  0 we apply an impenetrable boundary condition for   , and vacuum boundary conditions for the magnetic field with justified by the fact that the turbulent magnetic diffusivity increases with || (see Section 11.3 of Shukurov & Subramanian 2021, for details).Larger vertical sizes were used to confirm the domain was large enough to prevent any spurious boundary effects over the simulation period.We use the form of  given in equation ( 8) and an initial magnetic field of about 10 −6 µG in strength.Equations ( 15)-( 17) are nonlinear and the dynamo action can saturate due to the magnetic buoyancy alone.Figure 10 compares the root-mean-square strengths of the magnetic field in the Model R10h2 (solid, the reference model discussed above), and in the one-dimensional model with the same values of   and ℎ  (dash-dotted): the agreement in the growth rate of the magnetic field and the development of its steady state is quite satisfactory.
However, magnetic buoyancy is not the only nonlinearity and, as discussed in Section 5, the growing magnetic field also suppresses the -effect via the magnetic current helicity  m .Therefore, we also  considered a modification of the one-dimensional model that includes the quenching of the -effect by the magnetic field, a widely used heuristic form of the dynamo nonlinearity where  in equations ( 15)-( 16) is replaced with where  = ( 2  +  2  ) 1/2 and  0 = 3 µG.Apart from modifying the transition to the steady state and the magnetic field strength at which the dynamo action saturates, both evident in Fig. 10, the additional nonlinearity affects the period of the magnetic field oscillations, confirming once again the nonlinear nature of the oscillatory behaviour.
Figures 11 and 12 compare the magnetic field evolution, including its oscillations obtained in Model R10h2 with those from the onedimensional model.Solution of equations ( 15)-( 17) with  a function of  alone and independent of  does oscillate but the period of the  oscillations is shorter than in the three-dimensional simulations.The additional -quenching of equation ( 21) makes the oscillation period longer and dependent on  0 .In the lower panels, we have adjusted  0 to obtain the oscillations at about the same period as in the threedimensional simulation although this affects the saturation level of the magnetic field.
We do not attempt to achieve a precise match between the threedimensional and one-dimensional results being content with the fact that the one-dimensional model justifies further our conclusion that the magnetic oscillations are an essentially nonlinear phenomenon that relies on an interaction of the mean-field dynamo and magnetic buoyancy not envisaged earlier.

HELICAL FLOWS AND DYNAMO DRIVEN BY A BUOYANT MAGNETIC FIELD
Figure 13 shows the mean kinetic helicity coefficient of the gas flow (this does not include the imposed -effect) defined as computed using the gas velocity u in the simulations and the flow correlation time  obtained using the autocorrelation function of |u| as explained Appendix B. Figure 13 shows  k additionally averaged over  and  (the horizontal average), ⟨ k ⟩   .Apart from the -effect, large-scale magnetic fields are subject to the turbulent magnetic diffusivity with the coefficient Since Model R10h2 does not include rotation, the mean helicity of the flow ⟨u•(∇×u)⟩ can only be driven by the Lorentz force.The Coriolis force in a stratified, rotating system is the cause of the conventional -effect with  k > 0 for  > 0 and  k (−) = − k () (e.g., Section 7.1 of Shukurov & Subramanian 2021).While the antisymmetry of  k in  is evident in Fig. 13, the sign of  k is opposite to that of the -effect produced by the Coriolis force.The mean-field dynamo action imposed near the midplane, with  as given in equation ( 8), has the conventional sign,  > 0 at  > 0. However, the magnetic field generated by the dynamo drives helical motions that have the opposite sign of the mean helicity, which leads to the saturation of the dynamo action in the layer near the midplane where the -effect is imposed (e.g., Section 7.11 of Shukurov & Subramanian 2021).
Because of the magnetic buoyancy, the magnetic field spreads out from the layer || ≤ ℎ  and drives motions with the 'anomalous' mean helicity,  k < 0 at  > 0 and  k > 0 at  < 0 producing the picture shown in Fig. 13.Indeed, the gas flow becomes helical only at later times  > ∼ 1.5 Gyr when the magnetic field has become strong enough and has spread to large ||.Since the system is highly nonlinear, significant parts of the velocity and magnetic fields are randomised as shown in Figs. 4 and 5. Tharakkal et al. (2023a, their Fig. 13) similarly find that the sign of the mean kinetic helicity is anomalous, with  k < 0 at  > 0, in their simulations of the magnetic buoyancy of an imposed magnetic field in a rotating system.Thus, the mean magnetic field generated near the midplane spreads to larger altitudes where it drives helical flows and an ensuing meanfield dynamo, but the sign of the -effect is opposite to the sign of the imposed -effect near the midplane.As we show in Section 4.1, the resulting complex dynamo system generates oscillatory magnetic fields.
As with any other dynamo action, the dynamo driven by the Lorentz force away from the midplane saturates by reducing the magnitude of the -effect, so that the total -coefficient is modified as  =  k +  m , where (Section 7.11 of Shukurov & Subramanian 2021) where b = B − ⟨B⟩ is the deviation of the total magnetic field B from its mean ⟨B⟩.To compute b and then  m , we derived ⟨B⟩ by smoothing the total field with a Gaussian kernel as in equation ( 12) with a smoothing length of ℓ = 200 pc.As expected,  m , shown in the lower panel of Fig. 13, has the sign opposite to that of  k and a comparable magnitude.The blue dash-dotted line in the upper row of panels is the sum of the kinetic and magnetic helicities  k () +  m (), obtained using equations ( 22) and ( 24), and thus not including the imposed  of equation ( 8).
The green dash-dotted line in the lower row of panels represents the total magnetic diffusivity  t +  obtained using equation ( 23).

Electromotive force and transport coefficients
To verify, refine and further justify our interpretation of the results, particularly the dynamo action driven by the MBI, we have computed the components of the (pseudo-)tensor    and tensor    using the method of singular value decomposition (SVD) introduced by Bendre et al. ( 2020) and Bendre & Subramanian (2022).In this approach, the time series for the Cartesian components of the electromotive force E  = ⟨u × b⟩  are by E  =    ⟨  ⟩ −    (∇ × ⟨B⟩)  . Explicitly, are solved to determine the elements of the tensors    and    , which are assumed to be independent of time.This assumption is valid in either the early stages of the exponential growth of the magnetic field or in the later, stationary state of the system.In these calculations horizontal averaging is used, ⟨B⟩ = ⟨B⟩   , as displayed in Fig. 8, such that the tensor elements are functions of  alone.The horizontal average of the vertical component of the magnetic field vanishes due to the horizontal periodic boundary conditions.Hence, the analysis is applied only to the horizontal components of the magnetic field.
We also estimated the components of    and    using the alternative IROS method described by Bendre et al. (2023).The results obtained with the two independent methods are quite consistent with each other, and we present the SVD results below.
The governing equation for the mean magnetic field (when the mean gas velocity vanishes, as in our case) has the form The diagonal elements of the -tensor represent the -effect, with  k +  m ≈ 1 2 (   +    ).If the flow is isotropic in the (, )-plane    is antisymmetric (   = −   ), and the off-diagonal elements of    represent the transfer of the mean magnetic field along the -axis at the effective speed   = −   due to the increase in the turbulent magnetic diffusivity with || resulting mainly from the increase of the random flow speed (turbulent diamagnetism -e.g., Section 7.9 of Shukurov & Subramanian 2021).The diagonal components of the tensor    represent the turbulent magnetic diffusion.Since the imposed -effect of equation ( 8) and the magnetic diffusivity  in equation ( 3) are not associated with any explicit velocity and magnetic field components, the elements of    and    obtained with the SVD method are not sensitive to them.
Figure 14 presents the resulting components of the tensors    and    for the nonlinear stage of the evolution.The yellow stripes represent one standard deviation of the variables obtained from five estimates, each resulting from the sampling of every fifth entry in the time series of E containing 2500 data points at each , measured with the time interval 1 Myr.
Confirming the results and arguments of Section 5,    +   is significant in magnitude, antisymmetric with respect to the midplane  = 0 and mostly negative at  > 0. The magnitudes of    and   are close to  k +  m obtained using equations ( 22) and ( 24) and increase with ||.
The off-diagonal components of    are quite close to the expected antisymmetry,    = −   and reflect the transfer of the mean magnetic field towards  = 0 associated with the increase of the turbulent magnetic diffusivity with ||, which opposes the buoyant migration of the magnetic field away from the midplane.The estimated values of    and   shown in the lower row of Fig. 14 are significantly larger than  = 0.03 kpc km s −1 , which shows that the velocity field develops a rather strong random component by  = 1.5 Gyr.The components of    obtained with the SVD agree less with equation ( 23) although    is reasonably consistent with the analytical estimate at || ≲ 1 kpc.The accuracy of the SVD results is expected to be lower for the magnetic diffusion tensor because they involve the noisier spatial derivatives of the mean magnetic field.
To confirm our conclusion that the dynamo action and the associated complex behaviour of the mean magnetic field are essentially nonlinear phenomena, we have verified that the components of both    and    fluctuate around the zero level during the linear stage of the evolution without any significant effect on the evolution.

SUMMARY AND IMPLICATIONS
The nonlinear interaction of the mean-field dynamo and magnetic buoyancy leads to profound changes in the evolution of the largescale magnetic field transforming its properties far beyond what can be envisaged from a study of the early stages of the magnetic field amplification (when the Lorentz force is still negligible).A buoyant magnetic field not only spreads naturally to larger altitudes from a relatively thin layer where it is generated by the dynamo, but the buoyancy can change dramatically the magnetic field that has caused it, reversing its direction and leading to nonlinear magnetic oscillations.Neither the dynamo system in a plane layer nor the magnetic buoyancy instability exhibit any oscillatory behaviours on their own.
The nonlinear development of the magnetic buoyancy instability of an imposed uni-directional magnetic field has been studied by Tharakkal et al. (2023b,a).We have extended those results to the case of a buoyant magnetic field supported by the mean-field dynamo action.Both systems behave similarly as they reverse the direction of the magnetic field, and we attribute this to the secondary meanfield dynamo action driven by the helical flows associated with the magnetic buoyancy.
The ability of magnetic buoyancy to produce the -effect and thus a mean-field dynamo action is well known but its synergy with the dynamo action of a turbulent flow driven independently has not been explored earlier.
The following picture emerges from the studies of Tharakkal et al. (2023b,a) and our work.
When the buoyant magnetic field is not helical (e.g., unidirectional) and there is no rotation, magnetic buoyancy just redistributes the large-scale magnetic field to larger altitudes reducing very strongly its pressure gradient and leaving the support of the gas layer against the gravity to the thermal pressure gradient and contributions from turbulence and random magnetic fields if present (Tharakkal et al. 2023b).Rotation changes the picture because the gas flows that accompany magnetic buoyancy become helical driving a mean-field dynamo that overwhelms the imposed magnetic field leading to its reversal, apparently signifying the onset of nonlinear oscillations (Tharakkal et al. 2023a).As we show here, similar oscillations occur independently of the rotation if only the buoyant magnetic field is produced by the mean-field dynamo.In this case, the magnetic field is helical, so the associated Lorentz force powers helical flows independently of the rotation and its Coriolis force.Because of the nonlinearity, flows over a wide range of scales develop and are randomized (Rodrigues et al. 2016).
Notably, the sign of the -effect caused by the Lorentz force is opposite to that of the conventional -effect resulting from the action of the Coriolis force.
Our simulations are performed in a relatively large but still limited part of a gas layer (2×2 kpc horizontally) using Cartesian coordinates.The computational domain is large enough to accommodate the most rapidly growing MBI mode, and it is likely that the results would not be much different in cylindrical coordinates where the unstable magnetic field would no longer be unidirectional.Therefore, it is reasonable to expect that our main conclusions apply to galactic and accretion discs as a whole, at least at some distance from the disc axis where the curvature is not very strong.
These results can change our understanding of the nonlinear evolution of large-scale magnetic fields in spiral galaxies and accretion discs.In the case of the mean-field dynamos in disc galaxies, the conventional picture (Shukurov & Subramanian 2021) is that of a non-oscillatory, monotonically growing large-scale magnetic field whose direction is controlled by the initial conditions.The strength of the magnetic field can decrease with time if the interstellar gas is depleted (Rodrigues et al. 2015) but this does not affect its direction.
Magnetic buoyancy and its transformation of the dynamo action enrich this picture with the possibility that the magnetic field in a spiral galaxy or accretion disc can become oscillatory at a later stage when its energy density becomes comparable to the turbulent and thermal energy densities.Both the strength and direction of the large-scale magnetic field can vary in time in this state.Such an oscillatory behaviour has been reported in the simulation of the meanfield dynamo action by supernova-driven turbulence under largescale shear by Gent et al. (2023).In their case, a single sign reversal of the mean magnetic field occurs close to dynamo saturation.This corresponds to a change in the sign of , which is due primarily to a late increase in the strength of  m , so could be a signature of MBI in their system.
The intensity and consequences of the interaction of the magnetic dynamo and buoyancy depend on the relative intensities of the two processes and may vary with location within a galaxy and between galaxies.The dynamo efficiency is controlled by the scale height of the gas and rotation frequency.The consequences of magnetic buoyancy also depend on these parameters but in a different manner.The implications of our findings for the large-scale magnetic fields in astrophysical discs require a dedicated exploration.

Figure 1 .
Figure 1.The local strength of the dynamo action, R = | () | ℎ  / for () of equation (8) and various values of   of equation (9) specified in the legend.The critical value of   is shown dotted: the mean-field dynamo is supercritical at those |  | where R () >  ,c .

Figure 2 .
Figure 2. The dimensionless growth rates  = ℎ 2 / of the magnetic field at early times of the evolution when the dynamo action is predominant (symbols, as specified in the legend) and the similarly normalised dynamo growth rate from equation (11) for those values of   where it applies (solid).The critical value of   , corresponding to  = 0, is obtained by interpolating the data points as shown with the cross.

Figure 3 .
Figure 3.The evolving strengths of the large-(solid) and small-scale (dashed) magnetic fields (with the separation scale of ℓ = 200 pc) in Model R10h2 averaged over |  | < ℎ  (blue) and |  | > ℎ  (red).The dashed lines represent the exponential growth with the growth rates in Gyr −1 from Table2.

Figure 4 .
Figure 4.The strength of the horizontally averaged magnetic field components ⟨  ⟩   , ⟨  ⟩   and ⟨  ⟩   (columns from left to right) in the ( , )-plane at various evolutionary stages:  = 0.5 Gyr and 1 Gyr (two upper rows, the linear dynamo phase: the magnetic field strength grows while its spatial structure remains unchanged) and  = 1.5 and 2 Gyr (the two bottom rows, an early development of the MBI and the nonlinear stage) in Model R10h2.

Figure 5 .Figure 6 .
Figure5.The magnetic lines in Model R10h2 of the total field B (left-hand column), separated using equation (12) with ℓ = 200 pc into contributions at the larger scales characteristic of the magnetic buoyancy B B (middle) and the smaller scales of the imposed dynamo B D (right-hand column).The red isosurface corresponds to a gas number density of 0.7 cm −3 .

Figure 7 .
Figure 7.The horizontal average of the gas density in the Model R10h2, in the initial condition ( = 0), in the early MBI phase ( = 1.5 Gyr) and in the nonlinear state ( = 2 Gyr).

Figure 8 .
Figure 8.The evolution of the horizontally averaged magnetic field components ⟨   ⟩   (upper panel) and ⟨   ⟩   (lower panel) in Model R10h2 normalised to their maximum values at each time.

Figure 9 .
Figure 9. Variation of the period of the nonlinear oscillations of the component of the magnetic field ℎ  : from three-dimensional simulations (pluses, also presented in Table2) and the one-dimensional model of Section 4.1 with  0 = 1.5 km s −1 (crosses).The fit to the three-dimensional simulation results is shown dash-dotted its form is specified in the legend with ℎ  in kiloparsecs.

Figure 10 .
Figure10.The evolution of the magnetic field strength in the Model R10h2 (solid/black) and the one-dimensional mean-field model of Section 4.1 without (dash-dotted/blue) and with (dashed/red) -quenching of equation (21) with  0 = 3 µG obtained using the same values of all relevant parameters.

Figure 11 .
Figure 11.The evolution of the horizontally averaged -component of the magnetic field in Model R10h2 (top panel, identical to the upper panel of Fig. 8) and in the one-dimensional model of Section 4.1 without (middle) and with (bottom) -quenching of equation (21) obtained using the common parameters   = 10 and ℎ  = 0.2 kpc.Magnetic field strengths are normalised to unity at each time.

Figure 12 .
Figure 12.As in Fig. 11 but for the -component of the magnetic field.

Figure 13 .
Figure 13.The evolution of the horizontally averaged mean kinetic helicity (upper panel) and mean magnetic helicity (lower panel) coefficients of the random velocity and magnetic fields, equations (22) and (24) respectively, in Model R10h2.The horizontal dotted lines are shown at |  | = ℎ  .

Figure 14 .
Figure 14.The elements of the turbulent transport tensors introduced in equation (25) for Model R10h2 in the nonlinear state at 1.5 ≤  ≤ 3.0 Gyr.The yellow stripes indicate one standard deviation of the variables based on bootstrap resampling of the time series of E. The blue dash-dotted line in the upper row of panels is the sum of the kinetic and magnetic helicities  k () +  m (), obtained using equations (22) and (24), and thus not including the imposed  of equation (8).The green dash-dotted line in the lower row of panels represents the total magnetic diffusivity  t +  obtained using equation (23).