2 A global disc simulation

2.1 Description of the model

In order to study the possibility of dynamo action of the form envisaged by VC2001 we use the recent global hydromagnetic accretion disc simulation by Arlt & Rüdiger (2001). Unlike earlier global MHD disc simulations by Armitage (1998) and Hawley (2000), the radial extent does here not cover the entire disc, so inflow and outflow boundary conditions have to be applied in the radial direction. A possible advantage of this is that a larger spatial resolution per unit length is now possible. Like Armitage (1998), Arlt & Rüdiger (2001) also use the ZEUS-3D code which, in turn, is very similar to the code used by Hawley (2000). A version of ZEUS that is close to the one used by Arlt & Rüdiger is described by Stone & Norman (1992a,b) and Stone et al. (1992). The equations that are being solved are

$$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho {\vec u})=0,$$ (1)

$$\frac{\partial{\rho \vec u}}{\partial t}+\nabla\cdot(\rho{\ve... ...ec u}) =-\nabla p-\rho \nabla\Phi+{\vec J}\times {\vec B}+...,$$ (2)

$$\frac{\partial {\vec B}}{\partial t}=\nabla\times({\vec u} \times {\vec B}) +\eta\nabla^2\vec B,$$ (3)

where $\rho$, ${\vec u}$, and ${\vec B}$ are density, velocity, and magnetic field, respectively; p is the pressure, $\Phi$ is the gravitational potential (solely from a central mass $M_\ast$), ${\vec J}=\vec{\nabla}\times{\vec B}/\mu_0$ is the current density, $\mu_0$ is the vacuum permeability, and $\eta$ is a constant magnetic diffusivity as implemented by D. Elstner, Potsdam. The dots on the right hand side of Eq. (2) indicate the presence of numerical viscosity for removing energy at small scales. No energy equation appears as we deal with an isothermal model where $p=c_{\rm s}^2\rho$. The sound speed, $c_{\rm s}$, is 7% of the Keplerian orbital velocity in the middle of the simulated ring. Cylindrical polar coordinates, $(r,\phi,z)$, are used and the entire azimuthal range from $\phi =0$ to $2\pi$is considered. The model covers vertically the range from z=-1 to +1which corresponds to 1.5-3 pressure scale heights, depending on the value of r.

In most of the models we apply closed boundary conditions for the flow on the z-boundaries, with the magnetic field penetrating the boundaries at right angles. One of the models (Model IX; see below) uses open boundaries in the z-direction. The inner radial boundary is open in all models; the mass flux is monitored on the inner boundary and fed back into the domain on the outer radial boundary with either a homogeneous or a Gaussian infall pattern for the density. The maximum infall velocity is constant in time and over zand is either 10^-2 or 10^-3 of the sound speed. The $\phi$-direction has periodic boundaries.

The initial configuration contains a relaxed disk with a slow outflow on the inner boundary due to numerical viscosity. The density scale height varies between $H_\rho=0.33$ and 0.66 between the radii r=4 and 6. In the absence of magnetic fields the system is hydrodynamically stable; this was verified numerically for up to 30 orbits without producing any visible changes at the end of the simulation. The magnetic field imposed to this configuration is merely a vertical B_z field with zero net flux through the vertical surfaces.

Table 1 summarizes the global runs used for this analysis. The same numbering of the models as in Arlt & Rüdiger (2001) is used. Two of the simulations are new on this list: Model Va is a repetition of the configuration of Model V, but with an initial magnetic field of mixed parity. The parity of the B_z field was zero at the beginning as was the parity of the emerging $B_\phi$ field. All other models start with an initial field parity of -1 (antisymmetry). Model IX is similar to Model V, but uses outflow boundary conditions for the vertical direction instead of closed boundaries. In Table 1 the magnetic diffusivity $\eta$ is also given for all runs. Since the considerations presented here were made after the actual production runs had been performed, only a limited number of fully three-dimensional snapshots were available, which is the reason for a relatively coarse sampling.

 

 

Table 1: Global simulations used for the analysis. The duration $t_{\rm end}$ of each run is given in orbits. The infall velocities are given in units of the sound speed $c_{\rm s}$. Symmetries refer to the type of initial B_z fields. The kinetic energy is based on the poloidal velocities only and is an average over the last two orbits.

Run	Grid (z, r, $\phi$)	Radial boundary condition	r-range	$T_{\rm orb}$	$t_{\rm end}$	Init. parity	$\eta$
II	$31\times61\times351$	homogeneous accretion $u_{\rm in}=-0.001c_{\rm s}$	4.0-6.0	0.159	14.7	antisym.	0.001
III	$31\times61\times351$	homogeneous accretion $u_{\rm in}=-0.01c_{\rm s}$	4.0-6.0	0.159	9.7	antisym.	0.001
VIII	$31\times61\times351$	homogeneous accretion $u_{\rm in}=-0.001c_{\rm s}$	3.0-7.0	0.103	22.4	antisym.	0.001
V	$31\times61\times351$	Gaussian accretion $u_{\rm in}=-0.001c_{\rm s}$	4.0-6.0	0.159	18.4	antisym.	0.01
Va	$31\times61\times351$	Gaussian accretion $u_{\rm in}=-0.001c_{\rm s}$	4.0-6.0	0.159	12.4	mixed	0.01
VI	$31\times61\times351$	Gaussian accretion $u_{\rm in}=-0.01c_{\rm s}$	4.0-6.0	0.159	16.1	antisym.	0.01
IX	$31\times61\times351$	Gaussian accretion $u_{\rm in}=-0.01c_{\rm s}$, z open	4.0-6.0	0.159	16.9	antisym.	0.01

2.2 The Vishniac-Cho correlation

In the theory of VC2001 large scale magnetic field generation is possible when there is a finite value of the correlation

$$C_{\rm VC}\equiv\left\langle\omega_\phi\nabla_\phi u_z\right\... ...\langle\omega_\phi^2\rangle\langle(\nabla_\phi u_z)^2\rangle},$$ (4)

where $\nabla_\phi\equiv r^{-1}\partial_\phi$, and $\vec\omega={\rm curl} \, {}\vec u$ is the vorticity. We begin by showing $C_{\rm VC}$ as time series; see Fig. 1. The temporal averages of these time series (excluding the initial five snapshots of each of the runs which correspond to a lapse of 2 to 4 orbits) are listed in Table 2. We found that the simulations can tentatively be divided into two groups: the first group comprises the low-$\eta$ models showing decreasing correlations (left column of Fig. 1), the other group has modestly large values of $\eta$ and is characterized by a non-vanishing $C_{\rm VC}$ that is on the average about 5 to 10 times larger than in the low-$\eta$ cases (right column of Fig. 1).

A scatter plot of $\nabla _\phi u_z$ versus $\omega _\phi$ is shown in Fig. 2 for the data of a snapshot from Model V, which has the largest (negative) correlation. The plot contains points in the $(z,\phi)$ plane at r=5which is in the middle of the computational domain. The plot looks rather noisy, but one sees nevertheless a slight negative correlation.

Figure 1: Correlation of $\nabla _\phi u_z$ versus $\omega _\phi$in global disc simulation; here of Models II and III (left), as well as Models V and Va (right).

Figure 2: Scatter plot of $\nabla _\phi u_z$ with $\omega _\phi$for a snapshot of Model V at $t=14.8\,T_{\rm orb}$. The straight line gives the least square fit, which has a correlation coefficient of $C_{\rm VC}=-0.13$ in this example.

2.3 Resulting magnetic field configurations

We first discuss the overall field structure. Horizontal slices of the field at z=-0.39are shown for Models VIII (less resistive) and V (more resistive) in Figs. 3 and 4, respectively. The former figure exhibits a spiral pattern whilst the latter is rather dominated by intermediate scale structures or eddies, which is probably directly a consequence of the larger magnetic diffusivity in that case.

Figure 3: Projection of magnetic field vectors in a horizontal slice at z=-0.39of Model VIII after 19.9 orbits.

Figure 4: Projection of magnetic field vectors in a horizontal slice at z=-0.39of Model V after 17.6 orbits.

Next we derive a number of averaged quantities from the simulations. Throughout this section we denote azimuthal averages by an overbar, e.g. $\overline{\vec{B}}=\int\vec{B}\,{\rm d} {}\phi/2\pi$. In Fig. 5 are shown the energies contained in the large scale field, $M_{\rm mean}=\int\overline{\vec{B}}{}^2{\rm d} {}V/(2\mu_0)$, and the energies of the remaining fluctuations, $M_{\rm fluct}=\int\vec{b}{}^2{\rm d} {}V/(2\mu_0)$, where $\vec{b}=\vec{B}-\overline{\vec{B}}$. Like in Fig. 1, the temporal evolution behaviours separate into the same two groups: the low-$\eta$ (less resistive) runs which show significant energies in the large scale field, and models with larger $\eta$ that are more resistive, but better able to generate fluctuation energies of at least 50% of the large scale energy. We note that there is one model (Model IX, not shown) where at the end the energy of magnetic fluctuations exceeds the large scale magnetic energy. The fact that the energy of the mean field is typically larger than that of the fluctuating field is somewhat surprising. A possible reason could be that the memory of the initial mean field has not yet been lost. It is also possible, however, that it is because of the global geometry and the shear that a strong large scale field is more easily established when the magnetic Reynolds number is large.

Figure 5: Comparison of energies in the large scale ($\phi$-averaged) magnetic field (solid lines) and the small-scale magnetic field (dotted lines) for Models II and III (left), as well as Models V and Va (right). Note that the fluctuations are larger in the cases where $\eta$ is larger.

Another distinction between the two groups of runs is given by the magnetic Taylor microscale, $\lambda _{\rm M}$, which we define here via $\lambda_{\rm M}^2=\langle\vec{B}^2\rangle/\langle\mu_0^2\vec{J}^2\rangle$. In Fig. 6 we show the value of $\lambda _{\rm M}$for the same four models as in Fig. 1. The quantity $2\pi\lambda_{\rm M}$ characterized the typical thickness of flux structures. Its significance is that in runs with dynamo action $\lambda _{\rm M}$ tends to increase with time until it reaches saturation (e.g. Brandenburg et al. 1996). Conversely, when the field is amplified just by field compression the value of $\lambda _{\rm M}$ decreases somewhat with time. This is what happened in the more resistive runs ($\eta=0.01$). These were actually the runs that did show evidence for a finite Vishniac-Cho correlation. In contrast, the less resistive runs ( $\eta=0.001$) do not show any such trend.

A useful quantity for assessing the importance of helicity in the large scale field is to look at the nondimensional quantity

$$\varepsilon_{\rm C}={C_{\rm mean}/k_1\over M_{\rm mean}},$$ (5)

where $C_{\rm mean}=\int\overline{\vec{J}}\cdot\overline{\vec{B}}\,{\rm d} {}V$is the current helicity of the mean field. In BD2001 the value of $\varepsilon\rm _C$ was found to be of order unity (between 1-2) for the models with a halo. In Fig. 7 we show the value of $\varepsilon\rm _C$for the global accretion disc runs. The less resistive models show negative current helicity whilst the more resistive ones show vanishing or positive current helicity toward the end of the run. In any case, the values of $\vert\varepsilon\rm _C\vert$ are smaller compared with the models studied in BD2001 where $\vert\varepsilon_{\rm C}\vert={\cal O}(1)$. This is not surprising, because in BD2001 the kinetic helicity was close to the maximum value, which would be an unrealistic assumption for any astrophysical body.

Figure 6: Comparison of the magnetic Taylor microscale $\lambda _{\rm M}$ for Models II and III (left), as well as V and Va (right).

Figure 7: Evolution of the nondimensional large-scale current helicity parameter (upper disc plane) for Models II and III (left), as well as Models V and Va (right).

We now discuss the sign of $\varepsilon_{\rm C}$. In BD2001 the sign of the kinetic helicity was positive, and so was the sign of the small scale current helicity. The sign of the large scale current helicity is typically opposite, i.e. negative in that case. In the present case, where we consider the upper disc plane, the kinetic helicity is negative, so we would expect a positive value of $\varepsilon_{\rm C}$, which is not what we find (except for one of the more resistive cases, Model Va). However, the unusual sign of $\overline{\vec{J}}\cdot\overline{\vec{B}}$ is primarily related to an unusual sign of the effective $\alpha$ found in the present simulations. This is because in the steady state the energy-generating effect, $\alpha\overline{\vec{B}}^2$, must balance turbulent diffusion, $\eta_{\rm T}\overline{\vec{J}}\cdot\overline{\vec{B}}$. Therefore the sign of $\alpha$ must coincide with the sign of $\overline{\vec{J}}\cdot\overline{\vec{B}}$. In discs, however, the sign of $\alpha$ is negative (in the upper disc plane), so $\overline{\vec{J}}\cdot\overline{\vec{B}}$ must also be negative, and hence $\varepsilon_{\rm C}$ is negative, as observed. The perhaps most convincing explanation for the negative $\alpha$ is that intense parts of a flux tube contract (to maintain pressure balance along field lines), but are also most buoyant. If this contraction is stronger than the expansion associated with the rise into a less dense medium, $\alpha$ will be negative (Brandenburg 1998, see also Rüdiger & Pipin 2000).

We shall now return to the question of whether there is any evidence for the presence of a dynamo effect as envisaged by VC2001. We therefore need to look at the possibility of magnetic helicity fluxes through the domain.

 

 

Table 2: Results after the magnetic field was switched on. All values are temporal averages excluding the first five snapshots (about 2 to 4 orbital periods). The first group of models (II, III, and VIII) comprises less resistive runs, whilst the second group (Models V, Va, VI, and IX) refers to the more resistive ones.

Run	$u_{\rm rms}$	$E_{\rm kin}$	$2\eta C_{\rm mean}/(\mu_0 u_{\rm rms}E_{\rm kin})$		$2\eta C_{\rm fluc}/(\mu_0 u_{\rm rms} E_{\rm kin})$		$Q_{\rm mean}/$	slope	$C_{\rm VC}$
			North	South	North	South	$(\mu_0 u_{\rm rms} E_{\rm kin})$
II	3.7	$3.6\times 10^5$	$-7.8\times10^{-6}$	$+1.5\times10^{-6}$	$-5.0\times10^{-6}$	$+5.6\times10^{-6}$	+0.00006	-0.0024	-0.012
III	7.5	$6.9\times 10^5$	$-4.0\times10^{-6}$	$+9.3\times10^{-7}$	$-2.9\times10^{-6}$	$+2.3\times10^{-6}$	+0.00007	-0.0005	+0.003
VIII	7.6	$4.9\times 10^5$	$-5.4\times10^{-7}$	$+3.5\times10^{-7}$	$-2.8\times10^{-6}$	$+3.5\times10^{-6}$	+0.00029	-0.0018	-0.005
V	3.8	$1.4\times 10^6$	$-2.2\times10^{-5}$	$+3.9\times10^{-5}$	$-2.0\times10^{-7}$	$+1.1\times10^{-5}$	-0.00013	-0.0086	-0.050
Va	5.3	$3.9\times 10^5$	$+6.9\times10^{-5}$	$+6.6\times10^{-5}$	$+2.4\times10^{-4}$	$-2.6\times10^{-5}$	-0.00038	-0.0096	-0.086
VI	4.0	$1.7\times 10^6$	$-5.0\times10^{-5}$	$+3.3\times10^{-5}$	$-5.0\times10^{-6}$	$+2.1\times10^{-6}$	-0.00157	$-0.012\phantom{0}$	-0.060
IX	2.9	$2.9\times 10^5$	$-1.9\times10^{-4}$	$+1.6\times10^{-4}$	$-1.4\times10^{-4}$	$+1.9\times10^{-4}$	+0.00007	-0.0053	-0.047

2.4 Magnetic helicity flux

In VC2001 it was argued that, although the overall magnetic helicity is small, there could still be a significant (spatially constant) flux of magnetic helicity vertically through the domain. The numerical procedures for evaluating gauge invariant magnetic helicity and magnetic helicity flux in cylindrical geometry with open boundaries in the r and zdirections are not yet available. However, for the present purpose most important is the contribution from the large scales. If we adopt horizontal averages (over r and $\phi$), the mean fields are one-dimensional and the mean magnetic vector potential can be obtained simply by integration. The corresponding magnetic helicity and magnetic helicity fluxes of the mean field can then be calculated quite easily (see the appendix of BD2001). In Fig. 8 we plot, for the four models, the magnetic helicity flux, $Q_{\rm mean}=Q_{\rm mean}^{(2)}-Q_{\rm mean}^{(1)}$, out of the domain through the two boundaries at z=z₁ and z₂. Here, $Q_{\rm mean}^{(1)}$ and $Q_{\rm mean}^{(2)}$ denote the upward helicity fluxes at z=z₁and z₂, respectively; see Appendix A.

Figure 8: Estimates of the magnetic helicity flux $Q_{\rm mean}$ out through the vertical boundaries. The dotted lines give the zero value.

The mean outward flux, expressed in "dynamical'' units, $\mu_0 u_{\rm rms} E_{\rm kin}$, is small ( ${\sim}10^{-4}$). A small outward flux was also found in the case of helical turbulence (BD2001). Furthermore, the signs tend to be different in the cases where $\eta$ is small ( $Q_{\rm mean}>0$ for $\eta=10^{-3}$) and where it is larger ( $Q_{\rm mean}<0$ for $\eta=10^{-2}$). It is doubtful that this result indicates any significant departure from zero, because the magnetic helicity on the two sides about the midplane of the disc are expected to be different. Thus, equal losses or gains on the two surfaces (z=z₁ and z₂) should result in zero net magnetic helicity flux. It is therefore now necessary to determine the mean upward fluxes of magnetic helicity on the two sides, $Q_{\rm mean}^{(2)}$ and $Q_{\rm mean}^{(1)}$. Its average is denoted by $Q_{\rm mean}^{\rm (up)}={\textstyle{1\over2}}[Q_{\rm mean}^{(2)}+Q_{\rm mean}^{(1)}]$. If there is indeed a systematic upwards flux through the two boundaries, this quantity should be finite and positive. As expected, this quantity turns out to be small. In order to check whether this is the result of some cancellation, we need to consider the magnetic helicity fluxes in smaller sub-volumes.

A difficulty associated with calculating magnetic helicity and magnetic helicity fluxes separately in two sub-volumes (e.g., above and below the midplane) is that we want to make sure that the sum of the two is equal to the total magnetic helicity calculated earlier. This will be the case provided the magnetic helicity in each sub-domain is calculated using the same gauge that also made the helicity of the full domain gauge invariant. This then also allows one to calculate the integrated magnetic helicity fluxes out of each sub-domain. The corresponding formulae are given in Appendix B.

It turns out that the helicity fluxes out of each sub-volume are actually large but of opposite sign. This means that there is actually a large magnetic helicity flux through the midplane, but not through the upper and lower boundaries. Having fixed the gauge such that $\int\overline{\vec{A}}\cdot\overline{\vec{B}}\,{\rm d} {}z$ is equal to the helicity of BD2001 for the full domain, we can also calculate the local magnetic helicity fluxes. We denote these by $Q_{\rm mean}^{\rm (mid)}$ (if evaluated at the midplane), or by $Q_{\rm mean}^{(z)}$ (if calculated for all values of z). In Fig. 9 we plot $Q_{\rm mean}^{\rm (mid)}$ and compare with the averaged boundary fluxes $Q_{\rm mean}^{\rm (up)}$. It turns out that $Q_{\rm mean}^{\rm (mid)}$ is indeed mostly positive, as predicted by VC2001, but this flux is not sustained all the way to the boundaries: $Q_{\rm mean}^{\rm (up)}$ is virtually zero by comparison. An exception is Run Va, where $Q_{\rm mean}^{\rm (mid)}$shows large variations about zero and $Q_{\rm mean}^{\rm (up)}$ begins to deviate systematically from zero. (We recall that this is the run where the initial field had mixed parity.)

Figure 9: Estimates of the integrated magnetic helicity fluxes at the midplane, $Q_{\rm mean}^{\rm (mid)}$, (solid line) and on the boundaries, $Q_{\rm mean}^{\rm (up)}$ (dotted line).

In order to see whether the magnetic helicity flux at the midplane is typical of the entire interior of the simulation domain, we plot in Fig. 10 the vertical distribution of the magnetic helicity flux, which was derived from horizontal averages of field and flow and then averaged in time (again excluding the first 2 to 4 orbits). Figure 10 shows that a positive (i.e. upwards) flux of magnetic helicity is indeed typical of the interior of the entire domain, and that it vanishes only near the boundaries. Thus the boundaries seem to play an important role, which may of course be unrealistic. We note, however, even Model IX with open boundaries does show a rapid drop of magnetic helicity flux near the z-boundaries. Clearly, an abrupt change of this flux implies generation and destruction of magnetic helicity near the vertical boundaries.

To summarize the global disc simulations, the correlation anticipated by VC2001 is present, provided the magnetic diffusivity is not too small. If the magnetic diffusivity is smaller, the magnetic field tends to be stronger and can become more important and may hence be able to suppress this correlation. Nevertheless, these are also the cases which show most clearly a systematic magnetic helicity flux within the simulation domain, even though it is unable to leave or enter it through the boundaries. It is difficult to say whether the effect of VC2001 was really responsible for the field generation found in the disc simulations. We recall that in the present simulations there is also some evidence for a $\alpha$-effect, although it is based on a rather noisy correlation between the turbulent electromotive force and the mean field; see Arlt & Rüdiger (2001). In the following section we isolate the VC-effect by studying a more idealized model with no net helicity. This model also allows for longer runs and therefore clearer statistics.