Magnetic Helicities and Dynamo Action in Magneto-rotationally Driven Turbulence

TL;DR

This study investigates the relationship between magnetic helicity and flux in magneto-rotational turbulence-driven dynamos, finding no direct correlation and highlighting the nonlinear nature of the dynamo process.

Contribution

It introduces a gauge-invariant method to measure magnetic helicity and flux in open domains within MRI-driven dynamo simulations, revealing their complex relationship.

Findings

Total helicity remains small despite magnetic flux generation.

No direct link between magnetic flux and helicity flux.

Helicity measurements are clarified using a new gauge approach.

Abstract

We examine the relationship between magnetic flux generation, taken as an indicator of large-scale dynamo action, and magnetic helicity, computed as an integral over the dynamo volume, in a simple dynamo. We consider dynamo action driven by Magneto-Rotational Turbulence (MRT) within the shearing-box approximation. We consider magnetically open boundary conditions that allow a flux of helicity in or out of the computational domain. We circumvent the problem of the lack of gauge invariance in open domains by choosing a particular gauge -- the winding gauge -- that provides a natural interpretation in terms of average winding number of pairwise field lines. We use this gauge precisely to define and measure the helicity and helicity flux for several realizations of dynamo action. We find in these cases, that the system as a whole does not break reflectional symmetry and the total helicity remains small even in cases when substantial magnetic flux is generated. We find no particular connection between the generation of magnetic flux and the helicity or the helicity flux through the boundaries. We suggest that this result may be due to the essentially nonlinear nature of the dynamo processes in MRT.

Equations

$$v_z=0,\frac{\partial v_x}{\partial z}=\frac{\partial v_y}{\partial z}=0,$$

$$\frac{\partial B_z}{\partial z}=0,B_x=B_y=0.$$

$$-0.5H<x<0.5H,0<y<\pi H,-3H<z<3H,$$

$$\rho ={\rho }_0\exp (-z^2/H^2).$$

$$-8H<x<8H,-6H<y<6H,-2H<z<2H,$$

$$\frac{\gamma }{\gamma -1}\kappa \rho \mathcal{R}\frac{dT}{dz}=\mp \sigma T^4$$

$$\mathbf{B}=\nabla \times \mathbf{A}\mathrm{with}\nabla \cdot \mathbf{B}=0.$$

$$\mathbf{A}\to \mathbf{A}+\nabla \Psi$$

$${\partial }_t\mathbf{B}=\nabla \times (\mathbf{u}\times \mathbf{B}),$$

$${\partial }_t\mathbf{A}=(\mathbf{u}\times \mathbf{B})+\nabla \Phi ,$$

$$H={\int }_{V}dV(\mathbf{A}\cdot \mathbf{B}),$$

$$\frac{dH}{dt}=-{\int }_S(F+F_a)dS\mathrm{with}F=(\mathbf{B}\cdot \mathbf{n})(\mathbf{A}\cdot \mathbf{u}+\Phi )\mathrm{and}{F}_a=(\mathbf{A}\cdot \mathbf{B})(\mathbf{u}\cdot \mathbf{n}),$$

$$H^R={\int }_V(\mathbf{A}+{\mathbf{A}}_{\mathbf{p}})(\mathbf{B}-{\mathbf{B}}_{\mathbf{p}})dV$$

$${\nabla }_H\cdot {\mathbf{A}}^w=0,\mathrm{where}{\nabla }_H\equiv ({\partial }_x,{\partial }_y,0).$$

$$\mathbf{J}=\nabla \times \mathbf{B}=\nabla \times (\nabla \times \mathbf{A})=-{\nabla }^2\mathbf{A}+\nabla (\nabla \cdot \mathbf{A}).$$

$$\mathbf{J}=-{\nabla }^2{\mathbf{A}}^w+\nabla ({\partial }_z{A}_z^w).$$

$$J_z=-{\nabla }_H^2{A}_z^w$$

$${\nabla }_H^2{\Phi }^w=-{\nabla }_H\cdot (\mathbf{u}\times \mathbf{B}).$$

$${\mathbf{A}}^w\to {\mathbf{A}}^w+{\mathbf{e}}_{z}f(z),{\Phi }^w\to {\Phi }^w+g(z),$$

$${\int }_{V}dVf(z)B_z=\int dzf(z)\int dxdy{B}_z=0,$$

$$f(z){\int }_{S}dxdy{B}_z{u}_z+g(z){\int }_{S}dxdy{B}_z=0,$$

$$\mathbf{B}=(b_1(z),b_2(z),0).$$

$$A_z=C_{1}x{b}_2(z)+C_{2}y{b}_1(z),$$

$$A_x=(C_1+1)G_2(z),A_y=(C_2-1)G_1(z),$$

$$\mathbf{B}(z)=(\sin z,\cos z,0).$$

$$\frac{dH}{dt}=-{\int }_S{F}_{T}dS+{\int }_{V}WdV,$$

$$F_T=F+\eta \mathbf{n}\cdot (\mathbf{A}\times \mathbf{J}),\mathrm{and}W=-2\eta (\mathbf{B}\cdot \mathbf{J}).$$

$$L={\int }_S\sigma T^{4}dS,$$

$$Q=\eta {\int }_{V}dV{\mathbf{J}}^2,$$

$$\mathcal{F}(t)={\int }_S{F}_T(t)dS.$$

Full text

accretion disc - MRI - MHD - dynamos - turbulence

1 Introduction

The magnetic helicity is a measure of the topological complexity of a magnetic field; in an ideal fluid bounded by a co-moving flux surface, it is a conserved quantity. This results makes precise the idea that in the absence of diffusion, magnetic field lines cannot reconnect and therefore, all their knots and linkages must be preserved (Moffatt, 1978). By contrast, dynamo action describes the sustained generation of magnetic fields by inductive processes within the bulk of an electrically conducting fluid; it often involves efficient reconnection or enhanced diffusivity, and more often than not, involves substantial changes in field topology. In a plasma with high electrical conductivity these two processes must find some way to compromise. How this compromise is reached is the subject of considerable debate. If the magnetic diffusivity is small, or more precisely, if the magnetic Reynolds number $Rm$ is large, one would expect that the field topology be nearly preserved, or, alternatively, that substantial topological changes might require timescales of the order of the (long) diffusion time. This is incompatible with the requirement that a dynamo operating at high $Rm$ should be able to change the field topology on the (short) dynamical time-scale. The problem is particularly acute in the case of large-scale dynamo action. The latter describes the generation of substantial amounts of magnetic flux. One can imagine that if the flux bundles making up the magnetic field are at all linked, increasing the flux in the bundles by dynamo action will lead to a corresponding increase in magnetic helicity in blatant violation of the general idea that in a highly conducting fluid the magnetic helicity should be (nearly) conserved. This has lead to a sustained theoretical effort to find ways to reconcile helicity conservation on the one hand, and efficient flux generation on the other (Kulsrud & Anderson, 1992; Vainshtein & Cattaneo, 1992; Gruzinov & Diamond, 1994; Cattaneo & Hughes, 1996; Field & Blackman, 2002; Blackman & Field, 2002; Shukurov et al., 2006; Sur et al., 2007).

One popular choice is based on the observation that if the region in which the dynamo operates is not enclosed by flux surfaces then helicity need not be conserved. If the magnetic field “pokes through” some of the boundaries, under the right conditions, there could be a non-zero flux of helicity in or out of the dynamo region that could alleviate the constraint posed by helicity conservation and allow efficient flux generation (Blackman & Field, 2000; Vishniac & Cho, 2001; Subramanian & Brandenburg, 2004; Shukurov et al., 2006; Ebrahimi & Bhattacharjee, 2014). At some basic level the idea is that the magnetic field could keep itself reasonably disentangled by expelling all the unwanted tangles though the boundaries. The particularly attractive feature of this idea is that the process of expulsion does not require fast reconnection or turbulent diffusivity, in principle, it could operate even at infinite $Rm$. It should be noted that if this general idea is correct one would expect the presence of open boundaries to be a necessary condition for flux generation, and some evidence of a connection between large-scale dynamo action and helicity fluxes through the boundaries. In theory, such a connection could be sought in numerical models of large-scale dynamos with open boundaries in an effort better to understand the relationship between dynamo action and helicity expulsion. One problem, however, immediately arises that makes a straightforward application of these considerations problematic. In a domain with open magnetic boundaries the helicity is not gauge invariant, and consequently nor are the fluxes. Different choices of gauge could lead to widely different results thereby undermining the whole analysis. Fortunately, there appear to be at least two natural ways to circumvent this problem. One is to replace the helicity with some other quantity that is related to the helicity, in other words it is still a measure of topological complexity, that is gauge invariant in an open field configuration. This avenue was explored by Berger & Field (1984) and led to the introduction of the relative helicity. One would hope that even though the relative helicity is not quite the helicity it may still be informative about dynamo action. The other is to choose a particular gauge and stick to it. This approach is very simple and very straightforward, but it requires some justification of why a particular gauge is selected, and more importantly some understanding of what is the geometric or topological interpretation of the helicity in that specific gauge (Prior & Yeates, 2014). In particular one would want to choose a gauge that is relevant to dynamo processes.

We address some of these issues in the present paper. Our approach is straightforward. We choose for our analysis the results of two numerical studies of dynamo action driven by the Magneto-Rotational-Instability (MRI) in shearing-box simulations (Bodo et al., 2014, 2015, Balbus, 2003, for a general discussion on MRI see e.g.). One of the reasons for choosing this particular setup is that with appropriate boundary conditions, the geometry of shearing-boxes is such that the radial and azimuthal magnetic fluxes are not conserved. Thus large-scale dynamo action can be easily demonstrated, and indeed it is, by monitoring the generation of these fluxes. We then choose a specific gauge that leads to a natural topological interpretation for the resulting helicity, compute the helicity and the helicity fluxes in this gauge and compare them with dynamo activity. We would like to note that other quantities that involve the magnetic vector potential and that may have topological interpretations have been proposed as relevant to dynamo action in various geometries. A notable example is the Vishniac and Cho flux (Vishniac & Cho, 2001) that may be important in the accretion disk geometry. In this paper, however, we choose to focus on the magnetic helicity alone. We believe this is an appropriate starting point since it is the magnetic helicity (properly defined as an integral quantity) that is both an ideal invariant and carries topological information about the magnetic field. Finally we close the paper with a discussion of our findings.

2 The Model

We will perform the analysis of magnetic helicity behavior for two sets of stratified shearing box simulations of MRI driven turbulence with zero net vertical magnetic flux, already presented elsewhere (Bodo et al., 2014, 2015). The two sets differ essentially for the equation of state: In the first set we employed an isothermal equation of state (Bodo et al., 2014), while in the second set we used the perfect gas equation of state with a finite thermal diffusivity (Bodo et al., 2012, 2013, 2015). The magnetic and mechanical boundary conditions are the same for both: periodic in the azimuthal direction $y$, shear periodic in the radial direction $x$, while in the vertical direction we assume that the boundaries are in hydrostatic balance, impenetrable and stress free, giving

[TABLE]

and also that the magnetic field is purely vertical, giving

[TABLE]

We should note that these conditions allow a net flux of magnetic helicity through the boundaries. It is instructive to compare these two cases because they give rise to different kinds of dynamo action. In the perfect gas cases, there are situations in which the dynamo is of the large-scale type capable of generating substantial amounts of toroidal flux, whereas in the isothermal cases the dynamo always appears to be small-scale. Although it is clear that the differences in dynamo activity are associated with the development of vigorous convection in some of the perfect gas cases, why the presence of convection in MRI driven turbulence should enable the generation of large-scale fields is not currently understood.

Both sets of simulations were performed with the PLUTO code (Mignone et al., 2007). As is often the case in many astrophysical simulations there is no explicit viscosity or magnetic diffusivity, this is done to extract the largest possible dynamic range for the given computational resources. However, it does not mean that the evolution of the system is ideal, amongst other things that would preclude the possibility of dynamo action, rather it implies that dissipative processes are controlled by the resolution. Thus, for the simulations to be presented here, the characteristic size of dissipative structures is of the order of a few gridpoints. This property notwithstanding, numerical dissipations acts very much like physical dissipation from the point of view of mediating reconnection and allowing changes in magnetic field topology. Furthermore, since PLUTO is written in conservative form for the total energy, whatever kinetic or magnetic energy is lost to numerical dissipation it is correctly accounted for in terms of a corresponding increase in internal energy.

Finally, we note here that in the stratified shearing-box geometry, gravity is anti-symmetric about the mid-plane. As a result, quantities like the horizontally averaged kinetic helicity density is likewise antisymmetric about the mid-plane, so in a sense, the regions above and below the mid-plane are individually not reflectionally symmetric, even though the layer as a whole is. However, in the cases we consider here, the normalized horizontally averaged kinetic helicity density is very small, smaller than half of a percent throughout most of the interior, increasing to approximately 2 percent in thin boundary layers near the horizontal boundaries. So the systems we consider are almost reflectionally symmetric.

2.1 The isothermal case

The computational domain covers the region

[TABLE]

where $H$ is the scale height. The grid size is $128\times 384\times 768$ with 128 grid-points per scale height. As mentioned above, there is no net vertical magnetic flux threading the box. The equilibrium density distribution is given by

[TABLE]

The solutions are characterized by cyclic patterns of mean azimuthal field propagating away from the equatorial plane and becoming more evident at high altitudes, while the equatorial region is dominated by small scale fluctuations. This behavior has been observed in all isothermal simulations (see e.g. Davis et al., 2010; Shi et al., 2010; Gressel, 2010; Flaig et al., 2010; Guan & Gammie, 2011) and is best described in term of two coupled dynamos operating in different regions: a small-scale turbulent dynamo operating in the mid-plane region, where most of the mass is concentrated, and a mean-field type dynamo that operates in the tenuous overlying regions and generates magnetic structures in the form of upward propagating dynamo waves (Gressel, 2010). Our data set covers an interval of time of $300\;\Omega^{-1}$, where $\Omega$ is the shearing box rotational frequency. In Fig. 1 we plot the resulting volume averaged Maxwell stresses (in units of $\rho_{0}H^{2}\Omega^{2}$) as a function of time.

2.2 The perfect gas case

In the second series of simulations we adopted the perfect gas equation of state, and also introduced thermal conduction so that the heat generated by dissipation of the MRI turbulence can be carried to the top and bottom boundaries, where it is radiated away. The simulations are initialized with an isothermal equilibrium configuration, with a density distribution given by Eq. 4, and are then evolved until they reach a stationary configuration in which all the heat dissipated by turbulence is radiated away at the top and bottom boundaries. The computational domain covers the region

[TABLE]

where $H$ is the scale height of the initial isothermal atmosphere. The final quasi-steady state is completely different from the initial one and bears no memory of it, however we keep a reference to the scale height of the initial isothermal configuration in order to make possible a comparison between these simulations and the isothermal ones. The computational domain is covered by a grid of $512\times 384\times 128$ grid-points. In this case, because of the introduction of thermal conduction, we need an additional boundary condition, in addition to those described above, for specifying either the temperature or the heat flux at the boundaries. Our choice has been to put the heat flux equal to the black body energy flux at the surface temperature, namely at $z=\pm 2H$ we impose

[TABLE]

where $T$ is the temperature, $\kappa$ is the thermal diffusivity, $\cal R$ is the gas constant and $\sigma$ is the Stefan-Boltzmann constant. We found two distinctive regimes: conductive and convective, corresponding respectively to large and small values of the thermal diffusivity. In the conductive regime, the heat generated by dissipation is transported by thermal conduction, the temperature and density profiles are not very different from the isothermal cases and also dynamo action has similar properties. At low thermal diffusivities the layer becomes unstable to convective motions that carry all the heat, the density profiles become flat and we observe an extremely effective large-scale dynamo capable of generating substantial amounts of toroidal flux and, correspondingly, much larger values of the Maxwell stresses. This behavior has been also confirmed by Gressel (2013); Hirose et al. (2014) and Hirose (2015). Two broad classes of dynamo solutions can be identified: the ones in which the toroidal flux is symmetric about the mid-plane and the others in which it is anti-symmetric. The symmetric solutions have higher Maxwell stresses. The system appears to swap between these states randomly but, typically, spends more time in the symmetric, higher-efficiency states. The data set covers a time interval of $2000\;\Omega^{-1}$. In Fig. 2 we show the behavior of the volume averaged Maxwell stresses (in units of $\rho_{0}H^{2}\Omega^{2}$) as a function of time; we have marked in red the time intervals in which we observe the antisymmetric solution.

3 Magnetic Helicities

The primary objective of the present work is to explore the relationship between magnetic helicity and large-scale dynamo activity. Before this task can be undertaken we need to provide ourselves with a workable definition of the magnetic helicity in a domain with open magnetic boundaries. Let ${\bf B}$ and ${\bf A}$ be the magnetic field intensity and the associated vector potential satisfying

[TABLE]

Clearly the above expression does not define ${\bf A}$ uniquely but only up to the addition of the gradient of an arbitrary scalar function. Thus the gauge transformation

[TABLE]

leaves ${\bf B}$ unchanged. In an ideal fluid the evolution of ${\bf B}$ is governed by the induction equation

[TABLE]

where u is the fluid velocity. Using (7) we can uncurl the induction equation to obtain a corresponding evolution equation for the vector potential, namely

[TABLE]

where $\Phi$ is an arbitrary function related to the the choice of gauge. The magnetic helicity is defined by

[TABLE]

where the integral is over the volume containing the magnetofluid. The evolution equation for $H$ can readily be obtained from (9) and (10),

[TABLE]

where $S$ is the surface bounding $V$ and ${\bf n}$ is the outward normal. Clearly, if $S$ is a comoving, or impenetrable boundary $F_{a}$ vanishes, furthemore, if $S$ is a flux surface, i.e. $({\bf B}\cdot{\bf n})=0$ on $S$ then $F=0$ and the helicity is conserved and gauge invariant. On the other hand, if $({\bf B}\cdot{\bf n})\neq 0$ on some subset of $S$ there will be a corresponding contribution to $F$ and the helicity will, in general, not be a conserved quantity. Notice that in this case, the flux of helicity through the boundaries defined by $F$, and therefore the helicity itself depend on the choice of gauge. This is a well known problem associated with domains with open magnetic boundaries. As we mentioned in the introduction there are at least two natural ways to proceed. One is to give up the helicity altogether and to replace it with some other quantity that, still has a useful topological interpretation, and is gauge invariant. This approach was pursued by Berger & Field (1984) who noted that the difference between the helicities associated with two distinct magnetic fields could be made gauge independent. This opens the possibility to define in a gauge invariant way the helicity in an open field domain relative to some reference field. Although in principle the reference field need only satisfy some appropriate boundary conditions, a popular choice is a potential field in which case the resulting quantity is called the relative helicity. Because a potential field is current free, a possible topological interpretation for the relative helicity is that it is the helicity associated with the current distribution in the fluid. For a domain bounded by plane surfaces, the relative helicity is defined by

[TABLE]

where ${\bf B}_{p}$ is the potential field with the same normal component as ${\bf B}$ on the boundaries, and ${\bf A}_{p}$ is the associated vector potential.

Another possibility is to pick a particular gage. Given our geometry, a good choice is given by the winding gauge, discussed by Prior & Yeates (2014) and defined by the requirement that the vector potential satisfies

[TABLE]

It is possible to show that in cylindrical domains, i.e. domains in which the vertical boundaries are flux surfaces while the horizontal ones are open, the helicity computed with the winding gauge gives the average pairwise winding of the field lines (Prior & Yeates, 2014). This gives a particularly intuitive picture of the field topology. To see how the helicity and helicity flux can be computed in the winding gauge we begin by noting that in general

[TABLE]

which upon use of (14) and restricting to our geometry, gives

[TABLE]

The $z$ component of equation (16) reads

[TABLE]

which can be solved to find $A_{z}^{w}$ up to an additive function of $z$ , $f(z)$, say. $A_{z}^{w}$ can then be substituted back into (16) to obtain the horizontal components of ${\bf A}^{w}$. We note that, unlike $A_{z}^{w}$, $A_{x}^{w}$ and $A_{y}^{w}$ are uniquely defined. The scalar potential $\Phi^{w}$ in the winding gauge can likewise be computed by taking the horizontal divergence of (10), which gives

[TABLE]

The last equation defines $\Phi^{w}$ up to an additive function of $z$, $g(z)$, say. We now show that for our geometry, the winding helicity $H^{w}$ and the associated boundary fluxes are independent of the choice of $f$ and $g$. We consider a transformation of the form

[TABLE]

its effect on $H^{w}$ is to give rise to an extra term of the form

[TABLE]

since there is no net vertical flux. Similarly from (12) the contributions to the helicity flux have the form

[TABLE]

where the first integral vanishes because $u_{z}$ vanishes on the horizontal boundaries, and the second vanishes because, as before, there is no net vertical flux.

Inspection of equation (17) shows that some care is required when the vertical component of the current vanishes identically. This is, for instance, the case for magnetic fields of the form

[TABLE]

From which it follows that the $z$ component of the vector potential in the winding gauge can be written as

[TABLE]

for any constants $C_{1}$ and $C_{2}$. Substituting (23) into (16) immediately gives the horizontal components of the vector potential in the form

[TABLE]

where $G_{1}^{\prime}(z)=b_{1}(z)$ and $G_{2}^{\prime}(z)=b_{2}(z)$. Clearly the resulting vector potential satisfies the winding gauge requirement (14) and is compatible with the prescribed magnetic field for any choice of $C_{1}$ and $C_{2}$. This ambiguity can easily be removed by considering the specific case of a horizontal field whose amplitude is constant and whose direction rotates uniformly with $z$. This can be written as

[TABLE]

It is reasonable to demand that the helicity density for such a field should be unity which corresponds to the specific choice $C_{1}=C_{2}=0$, which is what we adopt here. Finally, we note that adding a constant (vector) to the vector potential again satisfies all the requirements of the winding gauge but leads to changes in the helicity density. Again the ambiguity can be removed by requiring that uniform horizontal magnetic fields have vanishing helicity density which gives that the (vector) constant must itself be zero.

The previous considerations apply to an ideal fluid. If non ideal effects are included, like a finite magnetic diffusivity, $\eta$, say, the helicity is no longer conserved even in a magnetically closed system. The helicity flux must be modified and a new production term appears. The evolution equation for the helicity then becomes

[TABLE]

where

[TABLE]

In the numerical simulations being considered here, $\eta$ is not explicitly known since the dissipation is entirely numerical. It can, however be estimated for any simulation by noting that in a steady state the energy production rate, $Q$, must balance the luminosity, $L$. The latter can be computed directly from the simulation from the expression

[TABLE]

where the integrals are computed at the top and bottom boundaries and $\sigma$ is Stefan-Boltzmann constant. The former can be approximated by the expression

[TABLE]

where we have assumed that the viscous heating is small compared to Ohmic heating. This is motivated by noting that the numerical Prandtl number for our codes is close to unity, and that in our solutions the Maxwell stresses always substantially exceed the Reynolds stresses. As we shall see presently, even though this is only a working assumption, it gives reasonable estimates.

Having defined the helicity in two possible ways, relative and winding, it is natural to ask which, if any, is more naturally related to dynamo processes. Since there is no completely convincing argument that one choice is far better than the other, the sensible thing to do is to try them both. However the helicity densities associated with the relative helicity or with the winding helicity are practically identical for our set of data. Thus from here on we discuss all relevant quantities in the winding gauge.

4 Results and Analysis

We begin our analysis by specifying the units we adopt to describe the results. In order to be able to compare the isothermal and perfect-gas results we use the same units in both cases. Unfortunately there is no single set of units that is entirely natural to both systems, nevertheless, and as a reasonable compromise we adopt $1/\Omega$ as the unit of time, the isothermal scale height $H$ as the unit of length, and the mid-plane density $\rho_{0}$ as the unit of density. In the perfect-gas case $H$ and $\rho_{0}$ are computed relative to an isothermal atmosphere in hydrostatic balance with the same total mass and vertical extent.

We consider the perfect-gas case first. All quantities are computed over the same time interval as in Figure 2. Figure 3 shows toroidal flux (upper panel) and total helicity (lower panel) as functions of time. In order to get a feeling for the values of the helicity we recall that the helicity in the winding gauge measures the average twist between pairs of field lines in units of the rms square flux. Accordingly, on the right of the figure (lower panel) we have included a scale in these units.

Some comments are in order. At a typical time, the average amount of twist is modest; less than unity actually. Assuming that the magnetic field is concentrated in bundles, and as we shall see presently this is not a bad assumption, the small average twist arises not because the individual bundles are untwisted, rather because they are twisted with either sign roughly equally. Also, the total helicity changes sign with time and there is no evidence of symmetry breaking with a corresponding establishment of a long-lasting helical state of either sign. Occasionally, during the anti-symmetric phases the helicity drops to very small values. Interestingly the toroidal flux remains one signed for much longer periods than the helicity. Thus, if we use the toroidal flux as a measure of large-scale dynamo action, it seems to be uncorrelated, or at any rate unconstrained by the total helicity. For instance, after $t\approx 1100$ the toroidal flux remains substantial and positive while the helicity changes sign over ten times during the same period. Of course, one should entertain the possibility that the helicity is small and uncorrelated to dynamo action because most of it is efficiently expelled through the boundaries. This being the case, there should be some correlation between the time history of the toroidal flux and the helicity flux in or out of the domain. This possibility can easily be eliminated by inspection of the upper panel in figure 4, which shows the corresponding time history of the total helicity flux (i.e. the difference of the fluxes though the two horizontal surfaces), defined by

[TABLE]

Figure 4 also shows the separate contributions to the total flux arising from the ideal and non-ideal parts (lower panel). Very much like the helicity itself, the helicity flux shows rapid fluctuations and no particular preference for either sign. Also, since the area under a typical peak is comparable with a typical value of the helicity, we conclude that the helicity flowing in and out of the volume is comparable to the amount of helicity inside the volume. Again one could argue that, since what is being shown in figure 4 is the total helicity flux, the fluxes though the individual boundaries could be much larger but with a small difference. Actually that is not the case, and indeed, the individual fluxes are comparable to the total. The other contribution to the helicity is the helicity production rate, related to non-ideal effects and given by last term in (26), namely

[TABLE]

This measures how efficiently the turbulence itself can generate or destroy helicity. Figure 5 shows the time history of $S_{H}(t)$ where we have used an effective diffusivity $\eta$ computed by equating the time averages of (28) and (29). Once again we note that the $S_{H}(t)$ fluctuates rapidly and has no particular sign. Furthermore, much as in the case of ${\cal F}(t)$, the helicity produced during a typical fluctuation is comparable to the total helicity in the domain. We should emphasize that this quantity is only a rough estimate, based on the assumptions that the unknown numerical dissipation can be approximated by regular diffusion with an effective diffusivity estimated by the method above. To check the validity of this assumption we note that the time integral of (26) gives

[TABLE]

Thus the helicity produced within the volume can also be estimated as the difference between the helicity present in the volume and that expelled (or injected) through the boundaries. Again we note that the above expression is only exact analytically, and not numerically. There are two main sources of possible errors: the first is that the code is not written in conservative form for the helicity, the second is that we only sample the results every unit of time, thus the (numerical) time integral is carried out on a coarsely sampled signal. The first error is controlled by the resolution the second by the cadence of the data collection. For our simulations the former is much, much smaller that the latter. These issues notwithstanding, the time histories of the helicity production computed from the two different methods are shown in figure 6. Given the uncertainties in estimating the values of $\eta$, it is quite remarkable that the two curves track each other as well as they do. It is clear that the curves slowly drift apart, but over a remarkably long time–of the order of several hundred orbits.

It is instructive to examine how the contributions to the total helicity are distributed spatially. To this end, figure 7 shows, in the two top panels, a plot of the normalized helicity density in the mid-plane at two times during which the system is in the symmetric (left panel) and anti-symmetric state (right panel). In the same figure, in the two bottom panels, we show also the corresponding distributions of the azimuthal magnetic field component, a comparison between the top and bottom panels show no apparent relationship between the two distributions. Inspection of the top panels indicates that the helicity density is highly inhomogeneous and is concentrated in oblique sheet-like structures with both right and left handed twist. This left-right symmetry is further evidenced in figure 8 which shows the probability distribution function (pdf) of the helicity density over the entire volume normalized by $<{\bf B}^{2}>$ at the corresponding times. The black curve corresponds to the symmetric state, the red curve to the anti-symmetric state, and we shall return to the blue curve presently. Clearly the modest value of the total helicity arises mostly from the high degree of symmetry of the pdf and not by the fact that the pdf is narrowly distributed around zero. We can make this more precise by noting that for the symmetric case the rms value is approximately 25 times larger than the average value. If we apply this value of 25 to a typical fluctuation in figure 3 we conclude that, if all the flux bundles had the same sign, the average pairwise twist would be close to 2.5 which is also an estimate of the typical twist in a bundle. If we repeat a similar argument for the antisymmetric case, the average value is close to 12 times smaller than in the symmetric case, while the ratio of rms to average is 3.4 times larger which corresponds to a typical twist per bundle of 0.7 or approximately 3.5 times smaller than in the symmetric case. Inspection of figure 7 shows that the most striking difference between the symmetric and anti-symmetric cases is the the smaller size of the typical structure in the latter. This size difference can be used to argue that, in a simple-minded way, the anti-symmetric case can be thought of as two symmetric cases stuck on top of each other and squished back into the same original volume. In the horizontal, that would correspond to 4 copies of the symmetric layer next to each other and shrunk so as to cover the same original area.

We note here, that the high degree of symmetry of the pdf’s can be used to argue against the possibility that regions of positive and negative twist have different filling factors. One could conceive of a situation in which the positive twist, say, is gentle but occupies a large fraction of the volume, while the negative twist is high but concentrated over a smaller region so that the overall twist is small. This would manifest itself in a skewed distribution with nearly zero mean–which is definitely not the case here. We also note that the high degree of symmetry of the pdf’s is preserved even when the distributions are computed separately for the regions above and below the mid-plane. This rules out possibility that the average twist is significant but distributed anti-symmetrically about the mid-plane.

We now turn to the isothermal case. Figure 9 and 10 show respectively the total helicity and helicity flux as functions of time. A few things are immediately obvious, the average twist for the isothermal case is smaller than for the perfect-gas case, also the characteristic time-scale for fluctuations in the fluxes are much shorter than the corresponding time scale for the total helicity. In fact a simple comparison of the values shows that the contributions to the total helicity from boundary fluxes are all but negligible. Whatever twist is present is generated internally. This result is not surprising if we recall that in the isothermal case the average density stratification is close to a Gaussian with most of the matter concentrated in a narrow region about the mid-plane. This is the region in which most of the magnetic energy and most of the Maxwell stresses are concentrated. Accordingly, the contributions to the helicity integral come from the same region with hardly anything flowing in or out of the boundaries. Inspection of the pdf for the isothermal case (blue line in Figure 8) shows that in analogy with the perfect gas case the value of the total helicity is mostly controlled by the symmetry of the pdf.

5 Discussion

We have initiated a numerical study of the relationship between large-scale dynamo action and magnetic helicity in domains with open magnetic boundaries. The purpose of this study is to provide some quantitative evidence to support the general idea that fluxes of magnetic helicity are an important ingredient in the large-scale dynamo process. The problem of the lack of gauge invariance in the definition of helicity and helicity fluxes in open domains has been circumvented by choosing a particular gauge, the winding gauge that, at least in domains with slab geometry, lends itself to a natural physical interpretation in terms of average pairwise twist of field lines (Prior & Yeates, 2014). We have applied these ideas to two related systems consisting of dynamos driven by MRI turbulence in shearing boxes with zero net flux. The difference between the two systems was in the equation of state, isothermal in one case, perfect-gas in the other. The reason to pick these two particular system was because despite their obvious similarities, they are basically the same system except for the equations of state, their dynamo properties are radically different. In the isothermal case there is little or no evidence for large-scale dynamo action, in terms of generation of substantial flux, whereas the evidence in the perfect-gas case is very convincing.

We find that there is no evidence for symmetry breaking with the formation of a helical state. Rather the helicity fluctuates, changing sign on a time-scale much shorter than the corresponding timescale for sign changes in the magnetic flux. Most of the helicity density is concentrated in isolated flux bundles whose average twist is of order unity. The presence of nearly equal amounts of the right and left handed bundles makes the average twist rather modest. In the perfect-gas case the helicity flux and helicity production rate are comparable, in the isothermal case the helicity flux is all but negligible. These considerations lead us to conclude that we find no compelling evidence of a simple relationship between large-scale dynamo action, measured, say, by the production of toroidal flux, and the helicity or the helicity flux. Even when substantial toroidal flux is present the overall state remains largely untwisted, and the helicity fluxes modest. We would argue that it would be very difficult to use the time traces in figures 3 or 9 to place any useful constrain on the dynamo efficiency. We should note that this result does not preclude that some quantity related to the magnetic vector potential other than the helicity may be a better indicator of large-scale dynamo action. We leave these speculations for a later study.

We close by suggesting that one possible explanation for tour findings is that the system under consideration belongs to the class of essentially nonlinear dynamos in the sense that the velocity field responsible for the dynamo action is entirely magnetically driven–indeed it is the result of the growth of the magneto-rotational-instability to finite amplitude. The velocity in the unmagnetized state is purely a Keplerian shear and as such incapable of sustaining a dynamo. It is only through the intervention of magnetic stresses that turbulence develops that can regenerate the magnetic field. This should be contrasted with the more conventional type of large-scale dynamos in which the basic flow is helical from the start (i.e. it lacks reflectional symmetry) and is capable of dynamo action. In this case one would expect the emergence of a non-reflectionally symmetric magnetic field in which the magnetic helicity and its flux may play a more prominent role. It would be interesting to repeat the kind of analysis presented here on such a system.

6 Acknowledgment

This work was supported in part by the National Science Foundation sponsored Center for Magnetic Self Organization at the University of Chicago. We acknowledge that the results in this paper have been achieved using the PRACE Research Infrastructure resource FERMI based in Italy at the Cineca Supercomputing Center.

Appendix A Magnetic helicity calculation

In order to compute magnetic helicity and its flux, equations (16) and (18) must be solved to determine the vector and scalar potentials in the winding gauge. If all quantities were periodic this task could most easily be achieved by Fourier transform methods. However the data as it stands is not periodic in the $z$ direction and it is only shear periodic in $x$. Furthermore, even if ${\bf B}$ were periodic, A, need not be. We show here how the data can easily be rearranged so that everything is exactly periodic. For the $x$ direction the problem can easily be overcome by restricting our analysis to the times when the shear-periodic domain re-aligns and the data becomes exactly periodic in $x$. Furthermore, the magnetic field can be split into two parts, namely

[TABLE]

where $\overline{{\bf B}}$ is the horizontally averaged magnetic field. The vector potential associated with ${\bf b}$ is itself periodic and poses no problems, the vector potential associated with the average can be calculated separately and is given by (24). In the $z$ direction we can exploit some useful symmetries implied by the boundary conditions in $z$. We recall that these are

[TABLE]

and

[TABLE]

on both boundaries. We can construct periodic data in the $z$ direction by the following procedure. Given the original domain ${\cal D}=\{0<x<L_{x},0<y<L_{y},0<z<L_{z}\}$, we define a new domain $\cal D^{\prime}$, with double size in the vertical direction, ${\cal D^{\prime}}=\{0<x<L_{x},0<y<L_{y},0<z<2L_{z}\}$. For each function $f$ defined on $\cal D$, that satisfies the condition $f=0$ on the boundary, we define a function $f^{\prime}$ on $\cal D^{\prime}$ by copying the original function $f$ in the lower half of $\cal D^{\prime}$, we then reverse $f$ along the $z$ direction, change sign and copy to the upper half of $\cal D^{\prime}$. More precisely we define $f^{\prime}$ as

[TABLE]

For each function $f$ defined on $\cal D$, that satisfies the condition $\partial f/\partial z=0$ on the boundary,we proceed in the same way, but we don’t change sign. We have then

[TABLE]

The functions defined on the domain $\cal D^{\prime}$ are fully periodic and can be treated by fast Fourier transforms techniques.