text:
# Anomalously large critical regions in power-law random matrix ensembles

\[

## Abstract

We investigate numerically the power-law random matrix ensembles. Wavefunctions are fractal up to a characteristic length whose logarithm diverges asymmetrically with different exponents, 1 in the localized phase and 0.5 in the extended phase. The characteristic length is so anomalously large that for macroscopic samples there exists a finite critical region, in which this length is larger than the system size. The Green’s functions decrease with distance as a power law with an exponent related to the correlation dimension.

\]

The power-law random banded matrix (PRBM) model was introduced by Mirlin et al.  to describe a 1D sample with random long-range hopping. The model attracted a lot of attention because it can approximately represent a variety of different physical systems: an integrable billiard with a Coulomb scattering center , two interacting particles in a 1D random potential , quantum chaos in a billiard with a non-analytic boundary , and the Luttinger liquid at finite temperatures . As the model describes a whole family of critical theories, it is also interesting for the study of intermediate spectral statistics .

The model is represented by $`N\times N`$ real symmetric matrices whose entries are randomly drawn from a normal distribution with zero mean and a variance depending on the distance of the matrix element from the diagonal

$$(H_{ij})^2=\frac{1}{1+(|ij|/b)^{2\alpha }}.$$
(1)
From field theoretical considerations , the PRBM model was shown to undergo a sharp transition at $`\alpha =1`$ from localized states for $`\alpha >1`$ to delocalized states for $`\alpha <1`$. This transition is supposed to be similar to an Anderson metal-insulator transition, presenting multifractality of eigenfunctions and non-trivial spectral compressibility at criticality. The parameter $`b`$ plays a role analogous to dimensionality establishing the character of the transition.

From a perturbative treatment of the non-linear $`\sigma `$-model and from renormalization group calculations Mirlin considered different regimes depending on the exponent $`\alpha `$. These models are only justified for $`b1`$, although the other limiting case can also be studied following Levitov . For $`\alpha 3/2`$ the eigenstates are localized, but in contrast to usual tight-binding models they do not decay exponentially, but with a power-law with the same exponent $`\alpha `$ as the hopping matrix elements . For $`1<\alpha <3/2`$ the states are still localized, but with a different typical length. In both cases, the wavefunctions are expected to have integrable power-law tails and the dimensionality to be equal to zero. For very large system sizes, the nearest neighbor normalized level spacings $`s`$ should be distributed according to the Poisson law, $`P_\mathrm{P}(s)=\mathrm{exp}(s)`$.

For $`\alpha <1`$, the previous theories predict that all states should be delocalized, independently of the value of $`b`$. The inverse participation ratio (IPR) should then be proportional to the system size, as for Gaussian orthogonal ensembles (GOE), although higher order cumulants of the IPR should not necessarily behave in the same way as GOE . For extended states, the nearest neighbor level spacing distribution approximately follows Wigner-Dyson law  , $`P_{\mathrm{WD}}(s)=(\pi s/2)\mathrm{exp}(\pi s^2/4).`$

At the critical value $`\alpha =1`$, we expect to have multifractal eigenstates and critical level statistics, intermediate between Poisson and Wigner-Dyson. Recent numerical calculations at criticality and for $`b=1`$ have shown that the nearest level spacing distribution differs from the typical one at a metal-insulator transition .

Our aim in this paper is to explore in detail the characteristics of the PRBM model in the interesting regime around $`\alpha 1`$ and for the case $`b=1`$ where theoretical treatments are not strictly applicable. We do so by studying systematically the nearest level spacing distribution function $`P(s)`$, the behavior of the Green function (GF) as a function of distance and the fractal properties of the states. We consider the PRBM Hamiltonian given by Eq. (1) with the exponent $`\alpha `$ ranging between 0.3 and 2. In order to reduce edge effects, we use periodic boundary conditions, so that our systems are rings of length $`L`$. We have checked that the results do not differ significantly from those obtained for linear chains away from the ends.

We first calculate the eigenvalues of the Hamiltonian matrix by direct numerical diagonalization in order to study $`P(s)`$. We consider a small energy window $`(0.4,0.4)`$ at the center of the band. For this calculation, the system size varies from $`L=200`$ to $`5000`$ and the number of random realizations is equal to $`10^6/L`$. We quantitatively analyze the scaling behavior of $`P(s)`$ through the following function of its variance

$$\eta (L,\alpha )=\frac{\mathrm{var}(s)0.286}{10.286},$$
(2)
which describes the relative deviation of $`\mathrm{var}(s)`$ from the Wigner-Dyson limit due to the finiteness of the system. In Eq. (2) $`\mathrm{var}(s)=s^2s^2`$, and 0.286 and 1 are the variances of Wigner-Dyson and Poisson distributions, respectively. We check that the use of other scaling variables , instead of $`\eta `$, do not alter the results.

Figure 1 shows the $`\alpha `$ dependence of $`\eta `$ for different system sizes: $`L=200`$ ($``$), 500 ($``$), 1000 ($``$) and 3000 () and 5000 ($``$). The value of $`\eta `$ at $`\alpha =1`$ is independent of system size and very close to the GOE value, instead of being intermediate between Wigner-Dyson and Poisson like in a standard metal-insulator transition. This is compatible with the numerical calculations of Ref. , where it was found that the spectral statistics at $`\alpha =1`$ differs from that typical of such a transition. For $`\alpha >1`$ there is a strong increase of $`\eta `$ with both system size and $`\alpha `$ towards the Poisson limit. For $`\alpha <1`$, $`\eta `$ is basically independent of system size and very similar to the GOE value. The small value of $`\eta `$ at criticality makes it difficult to decide whether there is a crossover at $`\alpha =1`$, as in standard metal-insulator transitions. It also indicates that, for $`b=1`$, the present transition is characterized by a very large conductance. This characteristic conductance would decrease with decreasing $`b`$ .

Next, we calculate the real space GF $`G(i,j;E)`$ as the inverse of the Hamiltonian matrix. For each value of $`\alpha `$ and $`L`$, we obtain $`G(i,j;0)`$ at the center of the band $`E=0`$ for different disorder realizations. We first average with respect to the initial position and then we make an ensemble average of the logarithm of this quantity, obtaining

$$F(n)\mathrm{ln}\underset{i}{}|G(i,i+n;0)|\mathrm{ln}L.$$
(3)
In this case, the system size varies from $`L=300`$ to $`4000`$ and the number of random realizations is equal to $`10^6/L`$.

In Fig. 2 we plot $`F(n)F(0)`$ as a function of $`n`$ for $`L=4000`$ and for different values of the exponent $`\alpha =0.7`$, 0.8, 0.85, 0.9, 0.95 and 1. We subtract $`F(0)`$ in order to remove effects associated with variations in the average density of states. The central region (tail) of the results can be fitted very well by a sum of two power-law terms with the same exponent $`\beta `$

$$F(n)F(0)=\mathrm{ln}\left\{\frac{A}{n^\beta }+\frac{A}{(Ln)^\beta }\right\},$$
(4)
where $`A`$ is a constant proportional to the density of states. Each contribution in Eq. (4) corresponds to propagation along one arm of the ring. The exponent $`\beta `$ is closely related to the matrix exponent $`\alpha `$ and weakly dependent on system size. The fit of Eq. (4) to the experimental points is so good that it is indistinguishable from the data on the scale used in Fig. 2. Up to the large scales used states are not truly extended. The results for $`\alpha >1`$ (not shown) are similar to those represented in the figure, for $`\alpha 1`$.

We have tried different extrapolation schemes in order to obtain the large scale behavior of the exponent $`\beta `$. We found that a logarithmic dependence on system size $`L`$ fits fairly well the results, as can be seen in Fig. 3 where we plot $`\beta `$ (solid symbols) versus $`L`$ for $`\alpha =0.8`$, 0.9, 0.95, 1, 1.05 and 1.1 on a semilogarithmic scale. The data can be fitted within errors by straight lines. The value of $`\beta `$ for the critical case, $`\alpha =1`$, is size independent and equal to $`0.25`$. For $`\alpha <1`$, $`\beta `$ decreases with size towards 0, while for $`\alpha >1`$ it increases towards $`\alpha `$. In the limit of macroscopic sizes, let us say, $`L10^8`$, we have a finite regime where the value of the exponent $`\beta `$ is different from zero and so states are not conducting. For infinitely large sizes $`\beta `$ is zero below the transition, $`0.25`$ at the transition and $`\alpha `$ above the transition.

To clarify the nature of the wavefunctions which give rise to this anomalous behavior of the GF, we calculate their fractal dimensions in the region around $`\alpha =1`$. We obtain the eigenfunctions of the Hamiltonian matrix by direct numerical diagonalization for system sizes ranging between $`L=500`$ and 10000. For $`0.75\alpha 1.25`$, we found that the wavefunctions show a complex structure with many sharp peaks of different heights and irregularly spaced which suggests the possibility of multifractal characteristics. We calculate the generalized fractal dimensions $`D_q`$ by the standard box counting procedure through the expression

$$D_q=\frac{1}{q1}\underset{\delta 0}{lim}\frac{\mathrm{ln}\left[\chi _q(l)\right]}{\mathrm{ln}\delta },$$
(5)
where $`\delta =l/L`$ and $`\chi _q(l)`$ is the $`q`$-th moment of the probability density of the wavefunction in the boxes of size of $`l`$ . For the evaluation of $`D_1`$ from the previous equation, one has to expand $`\chi _q(l)`$ around $`q=1`$. In practice we compute Eq. (5) at $`\delta =0.1`$.

We concentrate on the correlation dimension $`D_2`$, which reflects the scaling of the density-density correlation function. We found that, as for the exponent $`\beta `$, a logarithmic dependence of $`D_2`$ with $`L`$ fits fairly well the numerical data. We also found the following empirical relation between $`\beta `$ and $`D_2`$

$$\beta +D_2=\{\begin{array}{cc}1\hfill & \text{for }\alpha <1,\hfill \\ \alpha \hfill & \text{for }\alpha 1.\hfill \end{array}$$
(6)
To show graphically this relationship, we plot in Fig. 3, together with $`\beta `$ (solid symbols), the value of $`1D_2`$ (empty symbols) as a function of $`L`$ for $`\alpha =0.8`$, 0.9, 0.95, 1, 1.05 and 1.1, from bottom to top. In agreement with Eq. (6), the data for $`\beta `$ and $`1D_2`$ overlap for $`\alpha <1`$, while they are approximately shifted by an amount $`\alpha 1`$ for $`\alpha >1`$. This shift is consistent with the fact that in this regime $`D_2`$ tends towards zero, while $`\beta `$ goes to $`\alpha `$.

At $`\alpha =1`$, $`D_2`$ is practically size independent and equal to $`0.8`$. It very slowly tends to 0 for $`\alpha >1`$ and to 1 for $`\alpha <1`$, but even for (finite) macroscopic sizes remains between 0 and 1 for a finite width interval of the parameter $`\alpha `$. The value of $`D_2/d`$ at criticality is larger than for the standard Anderson metal-insulator transition in 3D , $`1.68/3=0.56`$, which indicates that the characteristic critical adimensional conductance for the PRBM model with $`b=1`$ is larger than for the Anderson model. This is in agreement with the small value found for the variance of the nearest level spacing at criticality, Fig. 1.

As for $`D_2`$, the other fractal dimensions are size independent at criticality. We obtained that the information (or entropy) dimension $`D_1`$ is equal to 0.88. Other fractal dimensions are $`D_2=1.27`$, $`D_1=1.14`$, $`D_3=0.71`$ and $`D_4=0.65`$. Trivially the similarity dimension $`D_0=d`$. In the range of dimensions studied, the deviations with respect to the weak multifractality limit, where $`(D_q/d)1q`$, are of the order of 10 %.

Away from the critical point, the wavefunctions are fractal up to a given length, which depends drastically on $`\alpha `$ and diverges as $`\alpha `$ goes to unity. This length $`\xi _{\mathrm{frac}}`$ is defined as the system size for which the correlation dimension becomes either 0 or 1 in a plot like Fig. 3. In Fig. 4 we show $`\mathrm{log}\xi _{\mathrm{frac}}`$ as a function of $`\alpha `$. We see that $`\mathrm{log}\xi _{\mathrm{frac}}`$ diverges at $`\alpha =1`$ in an asymmetric way, unlike in standard metal-insulator transitions. In the regime $`\alpha >1`$, a perturbative treatment of the non-linear $`\sigma `$-model predicts the existence of a length scale $`\xi `$

$$\mathrm{log}\xi c\frac{2\alpha 1}{2\alpha 2},$$
(7)
which plays the role of the localization length for the PRBM model . Note that $`\xi `$ is not a length characterizing an exponential decay of the wavefunction, since the GF, for system sizes larger than $`\xi _{\mathrm{frac}}`$, decreases with distance as a power law with an exponent equal to $`\alpha `$. For this reason it is the logarithm of $`\xi _{\mathrm{frac}}`$ rather than $`\xi _{\mathrm{frac}}`$ itself what diverges in Eq. (7) with a critical exponent and so plays the role of a localization length in a standard metal-insulator transition. Eq. (7) was obtained for the regime $`b1`$ and the constant $`c`$ was estimated to be of the order of $`\mathrm{log}b`$, although one can expect it to be qualitatively valid for any value of $`b`$. Our data in the localized regime can be fitted very well by a curve of the form given by Eq. (7), as can be seen in Fig. 4, where the right continuous curve corresponds to this equation with $`c=1.67`$. This value of $`c`$ is larger than predicted by the theory if extended to our case $`b=1`$.

For $`\alpha <1`$, the data of Fig. 4 cannot be fitted by Eq. (7), which corresponds to a critical exponent equal to unity, and there is no other theoretical prediction available. In this regime the data can be fitted to a curve of the form

$$\mathrm{log}\xi \frac{c}{|\alpha 1|^\nu },$$
(8)
with the critical exponent $`\nu `$ equal to 0.5. The left solid line in Fig. 4 corresponds to Eq. (8) with $`\nu =0.5`$ and the fitting parameter $`c=1.92`$. It is the first time to our knowledge that a metal-insulator transition is approach by two different critical exponents, which is in contradiction with a single parameter scaling theory. The possibility of the existence of two characteristic lengths has been already pointed out in a different one-dimensional model . We note that we are able to obtain two exponents numerically because the transition point is known exactly ($`\alpha =1`$), which constitutes a key advantage of this model.

To summarize, we have explored the critical region of the PRBM model in the case $`b=1`$ inaccessible to field theoretical methods. Our results, although compatible with theoretical predictions, show new light on the specific form of the GF and on the multifractal structure of the wavefunctions near the critical point. The wavefunctions are fractal up to a given length which diverges asymmetrically with two different critical exponents at $`\alpha =1`$. It is the first time that such an asymmetry has been obtained in part due to the exact knowledge of the position of the critical point.The possibility of a similar asymmetric behavior in other metal-insulator transitions should be reconsidered. We note that a fitting procedure of the raw data without an exact knowledge of the critical value of the transition parameter is compatible with equal exponents in the diverging lengths, above and below the transition. The fractal behavior is correlated with a power law decay of the GF with small exponents. The sum of the exponent characterizing the GF decay and the correlation dimension is approximately equal to unity in the localized phase and to $`\alpha `$ in the extended phase. All this should have important consequences in the problem of two interacting particles and in the other problems related to the PRBM model although their study is clearly outside the scope of this work.

###### Acknowledgements.

We would like to thank the Spanish DGESIC, project numbers BFM2000-1059 and 1FD97-1358, for financial support.  Present address: Department of Physics, California State University at Bakersfield, USA.


format: text

text:
# A simple solution to color confinement

## Abstract

We show that color confinement is a direct result of the nonabelian, i.e. nonlinear, nature of the color interaction in quantum chromodynamics. This makes it in general impossible to describe the color field as a collection of elementary quanta (gluons). A quark cannot be an elementary quanta of the quark field, as the color field of which it is the source is itself a source hence making isolated (noninteracting) quarks impossible. In geometrical language, the impossibility of quarks and gluons as physical particles arises due to the fact that the color Yang-Mills space does not have a constant trivial curvature.

One major problem in contemporary particle physics is to explain why quarks and gluons are never seen as isolated particles. A lot of effort has gone into trying to resolve this puzzle over a period of years, including lattice QCD, dual Meissner effect, instantons, etc. For a review, see .

We will take a different route than normally used, to eliminate the problem before it arises.

Usually, most particle physicists use “fields” and “particles” interchangeably, i.e. as denoting the same things. That is because the almost universal use of Feynman diagrams gives the false impression that particles (quanta) are always exchanged, even when they do not exist. The use of Feynman diagrams can be justified in mildly nonlinear theories (weak coupling limit) but breaks down for strongly coupled nonabelian theories. (And also for strongly coupled abelian theories with sources.) In quantum chromodynamics (QCD) it is at first sight a puzzle why the color force should be short-range, and especially why gluons are not seen as free particles, as the nonbroken SU(3) color symmetry seems to demand massless quanta, which naively would have infinite reach. However, as we shall see, there are no quanta.

In quantum field theory a particle , i.e. a quantum of a field, is defined through the creation and annihilation operators, $`a^{}`$ and $`a`$. For instance, in quantum electrodynamics (QED), the entire electromagnetic field can be seen as a collection of superposed quanta, each with an energy $`\omega _k`$. The hamiltonian of the electromagnetic field (omitting the zero-point energy) can be written

$$H=\underset{k}{}N_k\omega _k,$$
(1)
where

$$N_k=a_k^{}a_k,$$
(2)
is the “number operator”, i.e. giving the number of quanta with a specific four-momentum $`k`$ when operating on a (free) state,

$$N_k|\mathrm{}n_k\mathrm{}>=n_k|\mathrm{}n_k\mathrm{}>,$$
(3)
and, because all oscillators are independent,

$$|\mathrm{}n_k\mathrm{}>=\underset{k}{}|n_k>,$$
(4)
where $`n_k`$ is a positive integer, the number of quanta with that particular momentum. The energy in the electromagnetic field is thus the eigenvalue of the hamiltonian (1). The reasoning for fermion fields is the same, but then the number of quanta in the same state can be only 0 or 1.

Assuming that QCD is the true theory of quark interactions, a problem arises, as it is generally not possible to write the color fields in terms of superposed harmonic oscillators. It is not possible to make a Fourier expansion and then interpret the Fourier coefficients as creation/annihilation operators, as the color vector potentials $`A_\mu ^b`$ ($`b1,\mathrm{},8`$), even without a quark current, are governed by nonlinear evolution equations,

$$D^\mu F_{\mu \nu }=j_\nu ,$$
(5)
with quark current $`j_\nu g_s\overline{\psi }\gamma _\nu \psi =0`$, we get, in component form

$$(\delta _{ab}^\mu +g_sf_{abc}A_c^\mu )(_\mu A_\nu ^b_\nu A_\mu ^b+g_sf_{bde}A_\mu ^dA_\nu ^e)=0,$$
(6)
where $`g_s`$ is the color coupling constant (summation over repeated indices implied).

When we have an abelian dynamical group, as in QED, all the structure constants $`f_{abc}`$ are zero, and a general solution to Eq.(6) can be obtained by making the Fourier expansion

$$A_\mu ^{QED}=d^3k\underset{\lambda =0}{\overset{3}{}}a_k(\lambda )ϵ_\mu (k,\lambda )e^{ikx}+a_k^{}(\lambda )ϵ_\mu ^{}(k,\lambda )e^{ikx},$$
(7)
where $`ϵ_\mu `$ is the polarization vector.

However, for a theory based on a nonabelian group, like QCD, this is no longer true, due to the nonlinear nature of Eq.(6) when $`f_{abc},f_{bde}0`$,

$$A_\mu ^{bQCD}d^3k\underset{\lambda =0}{\overset{3}{}}a_k^b(\lambda )ϵ_\mu (k,\lambda )e^{ikx}+a_k^b(\lambda )ϵ_\mu ^{}(k,\lambda )e^{ikx}.$$
(8)
The color fields can be represented by harmonic oscillators only in the trivial, and physically empty, limit when the strong interaction coupling “constant” tends to zero, $`g_s0`$ (or equivalently when $`Q^2\mathrm{}`$ because of asymptotic freedom). Hence, no elementary quanta of the color interaction can exist. This means that no gluon particles are possible, and that Eq.(1) does not hold for color fields. (As another elementary example, there is nothing within QCD which resembles the photo-electric effect in QED, i.e. no “gluo-electric” effect!). The fields are always there, but their quanta, gluons and quarks, are relevant only when probed at sufficiently (infinitely) short distances. Generally, quarks do not exchange gluons, but the fermion fields $`q`$ react to the color fields given by $`A_\mu `$. Fields are primary to particles.

So far we have only banished gluons. To also banish quarks as physical particles we note that a quark field is the source of a color field, but this color field is itself a source of a color field. Hence, a quark field is never removed from other sources, is always interacting, and can never be considered to be freely propagating, resulting in that it can never be represented by harmonic oscillator modes. This means that no quark field quanta (quarks) can ever exist.

In QED things are very different. An electric charge gives rise to an abelian field, which is not the source of another field. Hence an electrically charged field can be removed from other sources and thus exist as a physical particle. Thus, the observability of, e.g. an electron is ultimately due to the fact that that electromagnetic quanta (photons) can exist as real particles.

A more mathematical treatment of the physical picture given above is provided by geometry. A case analogous to the one we are studying appears in quantum field theory on a curved spacetime , where it is well known that fields are more fundamental than particles. Indeed, there it can be shown that the very concept of a particle is, in general, useless . Actually nonabelian gauge fields and quantum field theory on a curved background have a lot in common. The total curvature, and also the dynamical coupling to “matter fields” through the covariant derivative, is given by one part coming from the Yang-Mills connection (i.e. gauge potential) and one part coming from the Riemannian (Levi-Civita) connection . Only when both the gauge field curvature and the spacetime curvature  are zero, or at most constant, can a particle be unambiguously defined. The former is constant for abelian quantum field theory, the latter is zero on a Minkowski background with inertial observers, and constant for some special (and static) spacetimes. The curvature in gauge space is given by the field strength tensor,

$$F_{\mu \nu }^b=_\mu A_\nu ^b_\nu A_\mu ^b+g_sf_{bcd}A_\mu ^cA_\nu ^d.$$
(9)
This is the analog in gauge space to the Riemann curvature tensor ($`R_{\mu \nu \sigma \rho }`$) for spacetime. The properties of $`F_{\mu \nu }`$ under a gauge transformation, $`U`$, is

$$F_{\mu \nu }F_{\mu \nu }^{}=UF_{\mu \nu }U^1,$$
(10)
i.e. the gauge curvature generally transforms as a tensor in gauge space. (It is also a tensor in flat spacetime, but not in a general Riemannian spacetime.)

However, we see directly that for an abelian gauge theory, like QED,

$$F_{\mu \nu }F_{\mu \nu }^{}=F_{\mu \nu },$$
(11)
as $`U`$ commutes with $`F_{\mu \nu }`$. This means that the curvature is constant (invariant) in gauge space for an abelian field, i.e. that $`F_{\mu \nu }`$ transforms as a scalar in gauge space. It also means that the field is a gauge singlet, which only reflects that it has no “charge” and that the fields have no self-interactions (abelian gauge fields $``$ sources of fields).

For nonabelian fields, like QCD, the gauge curvature, $`F_{\mu \nu }`$ transforms as a tensor, i.e. is covariant, not invariant, and is thus generally different at different points in gauge space. The color-electric and color-magnetic fields, $`𝐄^b`$ and $`𝐁^b`$, which are the components of $`F_{\mu \nu }^b`$ defined by $`F_{0i}^b=E_i^b`$ and $`F_{ij}^b=ϵ_{ijk}B_k^b`$ ($`i,j,k1,2,3`$) , are thus not gauge independent and cannot be observable physical fields, which is another way of seeing that gluons cannot exist, regardless of coupling strength. (In contrast to usual electric and magnetic fields which are both gauge singlets and observed.) Thus, color confinement is just a special case of the more general requirement that observables be gauge invariant, i.e. independent of the local choice of gauge “coordinates”. (Physically, $`F_{\mu \nu }^b`$ color singlet, implies that color gauge fields are sources of color gauge fields.)

We see that the unbroken nonabelian gauge theories of gravity and QCD are strictly incompatible with the concept of (“charged”) particles. In practice, however, this only rules out gluons and quarks as physical particles, as spacetime curvature (or, equivalently, observer accelerations) is normally completely negligible in experimental settings in particle physics. The difference can be traced to the fact that the dynamical curvature is directly related to the (nonlinear) coupling strength, which is enormously much larger for QCD than for gravity. Leptons can exist as physical particles as QED has abelian gauge dynamics and weak (nonabelian) SU(2) is broken, i.e. absent from the point of view of particle detectors.

It also follows, as a direct corollary of the argument above, that hadrons must be color singlets (i.e. color neutral) as they otherwise could not exist as physical particles. It would be interesting to continue the analogy with gravity and speculate that the hadrons are “grey holes” (as the color stays inside). The curvature induced by the color fields would then give the hadron (or confinement) radius. This would require solutions to the coupled $`q`$-$`A`$ system (with fully dynamical quark fields), which is a very hard, and unsolved problem. (Strictly, also gravity should be included, perhaps in a Kaluza-Klein fashion, the lagrangian then containing both $`F_{\mu \nu }F^{\mu \nu }`$, now with covariant spacetime derivatives, and $`R=R_\mu ^\mu `$, the Ricci-scalar.) Although this is a nice picture, which may/may not be true, it is not necessary for the purpose of excluding quarks and gluons as physical particles, for that the argument given in this article is sufficient.

In conclusion, what we have done is to provide a “Gordian knot”-type of solution to color/quark confinement.
We assume only that:

1) QCD is the correct theory of quark-field interactions

2) particles (quanta) are represented by $`a`$ and $`a^{}`$
which leads directly to the result that QCD can have no elementary quanta (gluons). If a specific fundamental quantum does not exist within a certain (supposedly correct) theory, it neither can be detected in experiments. As the quark fields always generate color fields, which in turn act as sources of other color fields, the quark fields can never be considered to be noninteracting. Hence no expansion in harmonic oscillator modes is possible, which means that no quark field quanta (quarks) can exist. Only if I) QCD is wrong, or II) quanta are not necessarily described by harmonic oscillator modes (which would mean that Einstein’s relation $`E=\mathrm{}\omega `$ does not hold), or both, can gluons and quarks exist as physical particles. In geometrical terms the curvature of Yang-Mills (color) space makes quarks (as particles) impossible. This proves that the theory of QCD automatically forbids particles with color charge, hence implying color (gluon/quark) “confinement”. In a way, there is nothing to confine in terms of particles.


text:
# The effects of vertical outflows on disk dynamos

## 1 Introduction

Magnetic fields play an important role in accretion disks. Inside the disk, for example, they allow the magneto-rotational instability to develop (Balbus & Hawley 1998). This instability can lead to the development of turbulence which in turn transports angular momentum radially outwards, thus allowing the matter to accrete. On the other hand, large scale magnetic fields are often considered necessary for launching winds or jets from accretion disks. The origin of such large scale magnetic fields is still unclear: they may be advected from the surrounding medium or be generated by a dynamo inside the disk. Advection appears unlikely as turbulence leads to enhanced viscosity and magnetic diffusivity, so that the two are of the same order (Prandtl number of order unity; Pouquet et al. 1976). Thus, turbulent magnetic diffusivity can compensate the dragging of the field by viscously induced accretion flow (van Ballegooijen 1989; Lubow et al. 1994; Heyvaerts et al. 1996). Dynamo action is a plausible mechanism for producing magnetic fields in accretion disks (Pudritz 1981; Stepinski & Levy 1988). However, dynamo magnetic fields can affect winds from accretion disks (Campbell et al. 1998; Campbell 1999; Brandenburg et al. 2000). On the other hand, the wind can also affect the dynamo. In particular, an enhancement by wind was suggested by Brandenburg et al. (1993) in the context of galactic dynamos. In the same context, Elstner et al. (1995) also consider the effects of shear in the wind.

In this paper, we consider the generation of magnetic field by a dynamo acting in a thin accretion disk. We provide a study of the effects of vertical velocities. We show how vertical velocities can enhance the dynamo, allowing a larger growth rate and leading to super-exponential growth of the magnetic field in the nonlinear regime.

The vertical velocities considered can have several origins. They can, for example, be due to a wind emanating from the disk or to magnetic buoyancy inside the disk. Note that we are interested here in magnetic buoyancy as generating vertical outflows. We neglect the associated dynamo effect (Moss et al. 1999).

The paper is organized as follows. In Sect. 2, we describe the two complementary approaches used and specify our choice of parameters. In Sect. 3, we discuss the linear results without vertical outflows. Results with vertical outflows are discussed in Sect. 4 for the linear regime and in Sect. 5 for the nonlinear regime. Sect. 6 presents our conclusions.

## 2 The models

Simulations of dynamo-generated turbulence in accretion disks have shown that, if stratification is included, a large scale field can be generated whose structure and dependence on boundary conditions closely resemble what is known from mean-field $`\alpha \omega `$-dynamo theory (Brandenburg 1998). It has been suggested that nonlinear effects can suppress the $`\alpha `$-effect to insignificant values (Vainshtein & Cattaneo 1992), but as recent simulations have shown, this does not affect the very possibility of mean-field dynamo action, but merely the time-scale on which the large scale field is established (Brandenburg 2001). With these justifications in mind, we adopt the mean-field approximation. The evolution of the magnetic field is given by

$$\frac{B}{t}=\times (V\times B+\alpha B\eta \times B),$$
(1)
where $`\eta `$ is the turbulent magnetic diffusivity and $`\alpha `$ the coefficient quantifying the $`\alpha `$-effect. We ignore anisotropies of $`\alpha `$ because only the toroidal diagonal component of the $`\alpha `$-tensor matters in the so-called $`\alpha \mathrm{\Omega }`$-regime. The antisymmetric part of the $`\alpha `$-tensor corresponds to turbulent diamagnetism which tends to push magnetic fields into the disk, and therefore just lowers the effective wind speed. Furthermore, as shown by Gabov et al. (1996, Fig. 6), turbulent diamagnetism is relatively unimportant for dynamo numbers less than 200, and our results mostly refer to this range. Given the substantial uncertainties in the mean-field transport coefficients, we decided to adopt the simplest possible approach. In order to minimize the number of parameters in our model, we also neglect the radial component of the velocity and thus our velocity field is $`V=(0,r\mathrm{\Omega },V_z)`$ in cylindrical coordinates ($`r,\phi ,z`$) (see comments in Sect. 2.2 on the consequences of this assumption).

We use two approaches to solve this equation. Firstly, a 2D numerical simulation based on a finite-difference scheme is applied (Sect. 2.1). Secondly, we develop a 1D free-decay mode expansion (Sect. 2.2). Details of these two complementary approaches are given below. Each of these two methods has its own advantages. The finite-difference simulation allows investigation of the nonlinear behaviour and applies to 2D geometry. On the other hand, the free-decay mode expansion gives the full eigenstructure and provides a tool for qualitative analysis. Furthermore, boundary conditions corresponding to an ideal vacuum can be implemented easily in the 1D model, whereas our 2D model where the disk is embedded in a halo can provide only an approximation to a vacuum outside the disk. Our two approaches are therefore complementary.

The problem is governed by three dimensionless numbers which appear in the dimensionless form of Eq. (1). Let $`h`$, the half-thickness of the disk, be a characteristic length-scale, $`\alpha _0`$ a characteristic value of $`\alpha `$, $`V_{z0}`$ a characteristic vertical velocity and $`\eta _0`$ a characteristic value of $`\eta `$. We also define the rotational shear $`S(r)=r\mathrm{d}\mathrm{\Omega }/\mathrm{d}r`$. The dimensionless numbers are the azimuthal magnetic Reynolds number

$$R_\omega =\frac{Sh^2}{\eta _0},$$
(2)
the magnetic Reynolds number based on the $`\alpha `$-effect

$$R_\alpha =\frac{\alpha _0h}{\eta _0},$$
(3)
and the vertical magnetic Reynolds number

$$R_v=\frac{V_{z0}h}{\eta _0}.$$
(4)
Note that $`\alpha _0`$ and $`S`$ depend on $`r`$ and thus the magnetic Reynolds numbers $`R_\omega `$ and $`R_\alpha `$ also depend on $`r`$. However, we assume later that $`\alpha _0S`$ is approximately constant. In this work, we adopt the $`\alpha \omega `$-approximation for the dynamo. In that case, for thin disks, the growth rates depend on the dynamo number $`𝒟_\alpha _\omega `$ only and not on the individual values of $`R_\alpha `$ and $`R_\omega `$. However, note that the ratio of the radial to the azimuthal magnetic field $`|B_r|/|B_\phi |`$ is about $`\sqrt{|R_\alpha |/|R_\omega |}`$ (Ruzmaikin et al. 1988). Moreover, Pudritz (1981) showed that the dynamo number $`𝒟`$ is approximately constant in a thin accretion disk.

### 2.1 The finite-difference model

To ensure that the magnetic field $`B`$ is solenoidal, we evolve the induction equation in terms of the vector potential $`A`$, where $`B=\times A`$,

$$\frac{A}{t}=V\times B+\alpha B\eta \mu _0j,$$
(5)
where $`j=\times B/\mu _0`$ is the current density due to the mean magnetic field and $`\mu _0`$ is the magnetic permeability. The $`\alpha `$-term $`\alpha B`$ together with the velocity term $`V\times B`$ drives the mean-field dynamo. We restrict ourselves to axisymmetric solutions of Eq. (5). We do not make the thin-disk approximation in this model but we consider thin disks by choosing an appropriate aspect ratio (see below). We do not make the $`\alpha \omega `$-dynamo assumption explicitly i.e. we retain the $`\alpha `$-effect in the equations for $`A_r`$ and $`A_z`$, but we restrict ourselves to the $`\alpha \omega `$-dynamo regime by our choice of magnetic Reynolds numbers $`|R_\omega |/|R_\alpha |=10^4\text{}10^6`$.

Our computational domain contains the disk embedded in a halo: $`0rL_r=2`$, $`|z|L_z=1`$, where the disk occupies the region $`0r1.5`$, $`|z|0.15`$, $`R=1.5`$ being the outer radius of the disk and $`h=0.15`$ its half-thickness. This results in a disk aspect ratio of $`h/R=0.1`$ and the upper boundary of the domain at $`z7h`$. Using a resolution of $`101\times 101`$ grid points and a uniform grid, the resulting mesh spacing is $`\delta r=\delta z=0.02`$. Equation (5) is solved using a sixth order finite-difference scheme in space and a third order Runge–Kutta time advance scheme.

We control the symmetry of our solutions (dipolar or quadrupolar) by using appropriate initial conditions. On the boundaries of the computational domain (i.e. far away from the disk) we impose the conditions that the normal component of $`A`$, $`A_{}`$, vanishes together with the normal derivative of the tangential components, $`A_{}`$,

$$A_{}=0,\frac{}{n}A_{}=0.$$
(6)
This is similar to the “normal-field” condition, where the tangential components of the magnetic field and the normal component of the electric field $`E`$ vanish, $`B_{}=0`$, $`E_{}=0`$. To verify that the boundary conditions do not affect our results, we have also tried perfectly conducting boundary conditions and found differences in the magnetic field structure only in the close vicinity of the boundaries.

The $`\alpha `$-coefficient $`\alpha (r,z)`$ is antisymmetric about the midplane and vanishes outside the disk. We adopt

$$\alpha (r,z)=\{\begin{array}{cc}\alpha _0(r)\mathrm{sin}(\pi z/h)& \text{for }|z|h,\hfill \\ 0& \text{for }|z|>h.\hfill \end{array}$$
(7)
The radial profile of $`\alpha (r,z)`$ is smoothly cut off at $`r>R`$ and at $`r_i=0.15`$ where the rotational shear is very strong.

As appropriate for accretion disks, we adopt a softened Keplerian angular velocity profile in $`r`$ in the whole computational domain (including the halo as well as the disk), assuming $`\mathrm{\Omega }`$ to be $`z`$-independent:

$$\mathrm{\Omega }(r)=\sqrt{\frac{GM}{r^3}}\left[1+\left(\frac{r_0}{r}\right)^n\right]^{(n+1)/2n},$$
(8)
where $`G`$ is Newton’s gravitational constant, $`M`$ the mass of the central object, $`r_0=0.05`$ the softening radius, and $`n=5`$. At $`r=0`$, $`\mathrm{\Omega }`$ vanishes as $`r^{3/2}`$.

The turbulent magnetic diffusivity is given by an interpolation between $`\eta _{\mathrm{disk}}`$ in the disk midplane and $`\eta _{\mathrm{halo}}`$ in the halo,

$$\eta (r,z)=\eta _{\mathrm{halo}}+(\eta _{\mathrm{disk}}\eta _{\mathrm{halo}})\xi (r,z).$$
(9)
Here, the profile $`\xi (r,z)`$ defines the disk,

$$\xi (r,z)=g(rr_0,d)g(|z|z_0,d),$$
(10)
where $`g(x,d)`$ is a smoothed Heaviside step function with a smoothing half-width $`d`$ set to 8 grid zones. So $`\xi (r,z)`$ is equal to unity deep in the disk and vanishes in the halo.

To obtain the magnetic Reynolds numbers of Eqs. (2)–(4) we choose $`\eta _0=\eta _{\mathrm{disk}}`$.

In accretion disks, the dynamo number is approximately constant with radius. In this model this is achieved by setting

$$\alpha _0(r)=\frac{𝒟\eta _{\mathrm{disk}}^2}{S(r)h^3}.$$
(11)
### 2.2 The free-decay mode expansion

The lowest order approximation for Eq. (1) in a thin disk results from retaining derivatives with respect to $`z`$ alone. In the $`\alpha \omega `$-approximation and with a constant $`\eta `$, Eq. (1) then becomes (see Ruzmaikin et al. (1988) for a form with $`V_z=0`$)

$`{\displaystyle \frac{B_r}{t}}`$ $`=`$ $`{\displaystyle \frac{}{z}}(\alpha B_\phi )+\eta {\displaystyle \frac{^2B_r}{z^2}}{\displaystyle \frac{}{z}}(V_zB_r),`$ (12)
$`{\displaystyle \frac{B_\phi }{t}}`$ $`=`$ $`SB_r+\eta {\displaystyle \frac{^2B_\phi }{z^2}}{\displaystyle \frac{}{z}}(V_zB_\phi ),`$ (13)
where once again $`S(r)=rd\mathrm{\Omega }/dr`$ and where $`B_r`$ and $`B_\phi `$ are respectively the radial and azimuthal components of the magnetic field. Note that we have neglected $`V_r`$, the radial component of the velocity in the disk, but the inclusion of this extra term would be straightforward (Moss et al. 2000).

We consider solutions with exponential time dependence

$`B_r(z,t)`$ $`=`$ $`\stackrel{~}{B_r}(z)\mathrm{exp}(\gamma t),`$ (14)
$`B_\phi (z,t)`$ $`=`$ $`\stackrel{~}{B_\phi }(z)\mathrm{exp}(\gamma t),`$ (15)
where $`\gamma `$ is the complex growth rate. Measuring $`z`$ in units of the half-thickness $`h`$ of the disk, $`t`$ in units of the diffusion time $`h^2/\eta `$, the coefficient $`\alpha `$ in units of a characteristic value $`\alpha _0`$ and $`V_z`$ in units of a characteristic value $`V_{z0}`$, the dimensionless form of Eq. (12) and Eq. (13) is given by

$`\gamma \stackrel{~}{B_r}`$ $`=`$ $`R_\alpha {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}(\alpha \stackrel{~}{B_\phi })+{\displaystyle \frac{\mathrm{d}^2\stackrel{~}{B_r}}{\mathrm{d}z^2}}R_v{\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}(V_z\stackrel{~}{B_r}),`$ (16)
$`\gamma \stackrel{~}{B_\phi }`$ $`=`$ $`R_\omega \stackrel{~}{B_r}+{\displaystyle \frac{\mathrm{d}^2\stackrel{~}{B_\phi }}{\mathrm{d}z^2}}R_v{\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}(V_z\stackrel{~}{B_\phi }),`$ (17)
where we have kept the same notation for the dimensionless variables and where $`R_\omega `$, $`R_\alpha `$ and $`R_v`$ are given by Eqs. (2), (3) and (4) with $`\eta _0=\eta `$.

We expand $`\stackrel{~}{B_r}`$ and $`\stackrel{~}{B_\phi }`$ with respect to the orthonormal base formed by the free-decay modes (Galerkin expansion),

$`\stackrel{~}{B_r}={\displaystyle \underset{n=0}{\overset{N}{}}}a_nB_r^{0,n},\stackrel{~}{B_\phi }={\displaystyle \underset{n=0}{\overset{N}{}}}b_nB_\phi ^{0,n},`$ (18)
where $`a_n`$ and $`b_n`$ are complex coefficients. The free-decay modes are the solutions of the equation

$$\gamma _n^0B^{0,n}=\frac{\mathrm{d}^2}{\mathrm{d}z^2}B^{0,n}.$$
(19)
Vacuum boundary conditions are imposed at the disk surface. We impose them for each free-decay mode so that $`\stackrel{~}{B_r}`$ and $`\stackrel{~}{B_\phi }`$, and thus $`B_r`$ and $`B_\phi `$, satisfy these conditions

$$B_r^{0,n}|_{z=1}=B_\phi ^{0,n}|_{z=1}=0.$$
(20)
At the midplane, the boundary conditions are connected to the symmetry of the mode considered. For quadrupolar modes we have

$$\frac{\mathrm{d}B_r^{0,n}}{\mathrm{d}z}|_{z=0}=\frac{\mathrm{d}B_\phi ^{0,n}}{\mathrm{d}z}|_{z=0}=0,$$
(21)
and for dipolar modes we have

$$B_r^{0,n}|_{z=0}=B_\phi ^{0,n}|_{z=0}=0.$$
(22)
Thus, the normalized eigenfunctions of free decay and the corresponding growth rates are, for quadrupolar symmetry

$`B_{r,\phi }^{0,n}`$ $`=`$ $`\sqrt{2}\mathrm{cos}\left(\left(\frac{1}{2}+n\right)\pi z\right)\mathrm{with}n=0,1\mathrm{}N,`$ (23)
$`\gamma _n^0`$ $`=`$ $`\left(\frac{1}{2}+n\right)^2\pi ^2,`$ (24)
and for dipolar symmetry

$`B_{r,\phi }^{0,n}`$ $`=`$ $`\sqrt{2}\mathrm{sin}(n\pi z)\mathrm{with}n=0,1\mathrm{}N,`$ (25)
$`\gamma _n^0`$ $`=`$ $`n^2\pi ^2.`$ (26)
Note that the free-decay dipolar mode $`n=0`$ would correspond to a uniform steady vertical magnetic field. The free-decay modes are normalized such that

$$_0^1\left(B_{r,\phi }^{0,n}\right)^2dz=1$$
(27)
(except for the dipolar mode with $`n=0`$). The free-decay modes are degenerate, forming doublets $`(B_r^{0,n},0),(0,B_\phi ^{0,n})`$ decaying at the same rate $`\gamma _n^0`$.

After making this expansion, the differential equations become an algebraic set of equations forming an eigenvalue problem of $`2(N+1)`$ equations

$`(\gamma \gamma _m^0)a_m`$ $`=`$ $`R_\alpha {\displaystyle \underset{n=0}{\overset{N}{}}}b_nY_{nm}R_v{\displaystyle \underset{n=0}{\overset{N}{}}}a_nW_{nm}^r,`$ (28)
$`(\gamma \gamma _m^0)b_m`$ $`=`$ $`R_\omega a_mR_v{\displaystyle \underset{n=0}{\overset{N}{}}}b_nW_{nm}^\phi ,`$ (29)
with $`m=0,1\mathrm{}N`$.

The coefficients $`W_{nm}^r`$, $`W_{nm}^\phi `$ and $`Y_{nm}`$ are defined by

$`W_{nm}^r`$ $`=`$ $`{\displaystyle _0^1}B_r^{0,m}(V_zB_r^{0,n})^{}dz,`$ (30)
$`W_{nm}^\phi `$ $`=`$ $`{\displaystyle _0^1}B_\phi ^{0,m}(V_zB_\phi ^{0,n})^{}dz,`$ (31)
$`Y_{nm}`$ $`=`$ $`{\displaystyle _0^1}B_r^{0,m}(\alpha B_\phi ^{0,n})^{}dz,`$ (32)
where the prime denotes the derivative with respect to $`z`$. We normalize the eigenfunctions (18) such that

$$\underset{0}{\overset{1}{}}\left[\left(\frac{\stackrel{~}{B_r}}{R_\alpha }\right)^2+\stackrel{~}{B_\phi }^2\right]dz=\underset{n=0}{\overset{N}{}}\left[\left(\frac{|a_n|}{R_\alpha }\right)^2+|b_n|^2\right]=1.$$
(33)
In this model, the truncation level needed depends on the values of $`R_v`$ and $`𝒟`$. Typically, $`N=4`$ is enough to compute the eigenvalues (growth rates) of the problem whereas $`N=9`$ is usually required to compute the eigenfunctions. The precision of all the results presented here was checked.

### 2.3 Choice of parameters

In this section, we discuss properties that are common to the two models.

The rotation of the thin disk is assumed to be close to Keplerian and independent of $`z`$ for the finite-difference model (see Sect. 2.1). Only the assumption of the independence of $`\mathrm{\Omega }`$ on $`z`$ is needed in the free-decay mode expansion.

We use two configurations: one where the magnetic diffusivity is the same in the disk and in the halo around the disk ($`\eta _{\mathrm{halo}}=\eta _{\mathrm{disk}}`$) and one with a low conductivity in the halo or a vacuum around the disk. Note that assuming a vacuum or low conductivity in the halo above the disk is a rather good approximation for galaxies (Sokoloff & Shukurov 1990; Poezd et al. 1993). It also allows for a connection with established theory. Note finally that the low-conductivity approximation ($`j0`$) is a good approximation for potential fields ($`j=0`$) and hence for a particular case of force-free fields ($`j\times B=0`$).

The boundary conditions of the free-decay mode expansion imply a vacuum above the disk. Technically, to achieve $`\eta _{\mathrm{halo}}=\eta _{\mathrm{disk}}`$ in this model, we restrict the dynamo action to an inner layer of the disk. The outer part of the disk is thus playing the role of the halo. The results simply need to be rescaled for a proper comparison.

In the finite-difference model, we study the case $`\eta _{\mathrm{halo}}=\eta _{\mathrm{disk}}`$ and the case of a low conductivity halo with $`\eta _{\mathrm{halo}}/\eta _{\mathrm{disk}}=20`$. Note that, with this model, to approximate a vacuum in the halo we have to choose $`\eta _{\mathrm{halo}}/\eta _{\mathrm{disk}}`$ as high as possible. But, this choice is restricted because the higher this ratio, the smaller the time step.

In this paper, the model with $`\eta _{\mathrm{halo}}=\eta _{\mathrm{disk}}`$ is referred to as the “disk+halo” model whereas the model with a vacuum or low conductivity above the disk is referred to as the “disk+vacuum” model.

In both models, the profile for $`\alpha `$ is taken to be

$$\alpha =\alpha _0\mathrm{sin}\left(\pi \frac{z}{h}\right)\text{for}|z|h$$
(34)
(see Eq. (7) for details concerning the finite-difference model), except in Fig. 1 where we take a linear dependence of $`\alpha `$ with $`z`$ to compare with previous work. The profile for $`V_z`$ is linear in $`z`$,

$$V_z=V_{z0}\frac{z}{h}.$$
(35)
The accuracy of the models depends on the number of grid points in $`r`$ and $`z`$ for the finite-difference model and on the level of truncation for the free-decay mode expansion.

In this work, the two control parameters are the dynamo number $`𝒟`$ and the vertical magnetic Reynolds number $`R_v`$. For reference, we give here possible values for these two parameters. While dynamo numbers are rather low in galactic disks, $`|𝒟|10`$ (Ruzmaikin et al. 1988), in accretion disks, they can be as high as $`|𝒟|10^210^4`$, depending on the Mach number of the turbulence, see e.g. Pudritz (1981). The value of $`R_v`$ obviously depends on the origin of the vertical outflow. With magnetic buoyancy, a rough estimate can be obtained assuming that the vertical motions occur at the Alfvén speed. Estimating the Alfvén speed from magnetic equipartition, one gets a vertical magnetic Reynolds number $`R_v`$ of order unity.

In the following, the figures obtained with the finite-difference model are labelled with “FD” and the ones obtained with the free-decay mode (Galerkin) expansion with “GE”.

## 3 Linear results for $`R_v=0`$

In this section, we present results with no vertical velocity in the disk ($`R_v=0`$). Such a study is an interesting test for the quality of the results of the free-decay mode expansion. It also gives a point of comparison between the numerical simulation and the free-decay mode expansion.

We first compare the results of the free-decay mode expansion (with $`R_v=0`$) with the results published by Stepinski & Levy (1991). We thus take a linear dependence of $`\alpha `$ on $`z`$, $`\alpha =\alpha _0z/h`$. The calculated growth rates are plotted in Fig. 1. The agreement is good but a discrepancy is observed at large $`|𝒟|`$ where the increase of the growth rate with $`|𝒟|`$ is faster according to Stepinski & Levy than in our free-decay mode expansion for both quadrupolar and dipolar symmetry ($`𝒟<0`$). This discrepancy seems to be due to an insufficient number of grid points in Stepinski & Levy’s model. Indeed, if we increase the truncation level of the free-decay mode expansion, no change is observed in the results whereas if we decrease the truncation level, we recover the results of Stepinski & Levy. Note that the results discussed in this paper are for $`|𝒟|<400`$ where the agreement is excellent.

We now want to compare the finite-difference simulation with the free-decay mode expansion. In particular, we want to test our method to mimic, in the case of the free-decay mode expansion, a model for $`\eta _{\mathrm{halo}}=\eta _{\mathrm{disk}}`$ starting with the “disk+vacuum” model (see Sect. 2). For an infinitely extended halo with $`\eta _{\mathrm{halo}}=\eta _{\mathrm{disk}}`$, the growth rates of dipolar and quadrupolar modes should be exchanged when $`𝒟`$ reverses sign, i.e. the graph of $`Re(\gamma )`$ as a function of $`𝒟`$ should be symmetric with respect to the vertical axis (Proctor 1977). The finite-difference simulation shows this symmetry with a good precision (Fig. 2a). The results of the free-decay mode expansion (Fig. 2b) are less symmetric. This is not surprising since the ratio between the height of the computational domain and the height of the disk is about 7 in the finite-difference simulation whereas it is only 2 in the free-decay mode expansion.

We now discuss our results for the “disk+halo” and “disk+vacuum” configuration. Let us first consider the “disk+halo” configuration (Figs. 2a and b). For moderate dynamo number $`|𝒟|300`$, the leading mode is quadrupolar and non-oscillatory if $`𝒟<0`$ ($`R_\alpha >0`$) and it is dipolar and non-oscillatory if $`𝒟>0`$ ($`R_\alpha <0`$). In the following, we thus study the effect of vertical motions for $`𝒟<0`$ in quadrupolar symmetry and for $`𝒟>0`$ in dipolar symmetry. Above $`|𝒟|300`$, the symmetry of the leading mode changes and it becomes oscillatory. The critical dynamo number is $`|𝒟_{\mathrm{crit}}|5`$.

Let us now consider the “disk+vacuum” configuration (Figs. 2c and d). The results obtained with the two approaches are rather different. With the free-decay mode expansion, the dipolar mode is always oscillatory and subcritical for $`|𝒟|400`$. On the other hand, with the finite-difference simulation, the dipolar mode becomes supercritical for $`𝒟50`$ and is non-oscillatory for all $`𝒟>0`$. This is because the dominant dipolar mode in the finite-difference simulation is the “forgotten mode” i.e. a dipolar mode corresponding to $`n=0`$ in Eqs. (25) and (26). The “forgotten mode” cannot be obtained with our 1D calculation (free-decay mode expansion). However, we reproduced the corresponding profile $`\gamma (𝒟)`$ with a modified 1D model which takes into account the finite aspect ratio of the disk. The model is based on the assumption of separability of $`r`$\- and $`z`$-dependence and yields equations similar to the thin-disk equations (12), (13), but with additional correction terms of order $`(h/r)^2`$. Thus, the dipolar modes of Fig. 2c and Fig. 2d are of a different nature. As for the quadrupolar mode, the results obtained with the two approaches are not so different. The branch for $`𝒟<0`$ is rather similar in both cases with a critical dynamo number $`50`$. The branch for $`𝒟>0`$ is all subcritical with the finite-difference simulation whereas it becomes supercritical for $`𝒟170`$ with the free-decay mode expansion. However, we have checked that, in the finite-difference model, if we take a steeper transition between $`\eta _{\mathrm{disk}}`$ and $`\eta _{\mathrm{halo}}`$, the branch for $`𝒟>0`$ moves upwards, resembling the one for the free-decay mode expansion.

The qualitative behaviour for the case of low conductivity in the halo is in many respects similar to that of the “disk+halo” configuration (Figs. 2a and c). But there are quantitative differences. There is no intersection point where the symmetry of the leading mode would change for $`|𝒟|700`$. The diagram does not look symmetric anymore.

The results of the finite-difference model presented here are in agreement with the accretion disk models of Rüdiger et al. (1995) and von Rekowski et al. (2000). Note, however, that these models have an inner radius and that at this radius a perfectly conducting boundary condition is used, which leads to a vertical current layer on the boundary and potentially lowers the excitation threshold for magnetic fields.

## 4 Linear behaviour with vertical velocities

In this section, we study the linear behaviour of the $`\alpha \omega `$-dynamo in the presence of vertical velocity ($`R_v0`$). As discussed in Sect. 3, we present results for quadrupolar symmetry with $`𝒟<0`$ (Sect. 4.1) and results for dipolar symmetry with $`𝒟>0`$ (Sect. 4.2.). An explanation of these results is given in Sect. 4.3.

### 4.1 Quadrupolar modes

We consider quadrupolar modes with $`𝒟<0`$. For sufficiently low dynamo numbers ($`|𝒟|50`$), the growth rate decreases with $`R_v`$ whereas at larger dynamo numbers, the growth rate increases for small $`R_v`$ and then decreases at large $`R_v`$, forming a maximum (Fig. 3). This maximum is a robust feature, the existence of which does not seem to depend on the boundary conditions or on the model. Indeed, we find this maximum with the finite-difference simulation (Figs. 3a and c) and with the free-decay mode expansion (Figs. 3b and d) in both configurations, “disk+halo” and “disk+vacuum” (see Sect. 2.3 for the definition of these configurations). The boundary conditions at the surfaces of the disk affect the position of the maximum in $`R_v`$ and the dynamo number above which it occurs. In particular, the smaller the magnetic diffusivity of the halo, the smaller this dynamo number. Nevertheless, the required values for $`𝒟`$ to get a maximum are larger than the typical values in galaxies. When $`𝒟`$ is increased, the maximum is more pronounced and its position is shifted towards larger $`R_v`$. Note that the values of $`R_v`$ at which the maximum occurs are between 1 and 10. They thus fit within the range of values estimated for magnetic buoyancy. Finally, in the range of dynamo numbers studied here, $`|𝒟|300`$, the vertical velocity does not change the symmetry or temporal behaviour of the dominant mode.

We now discuss the magnetic field structure. Figure 4 shows the eigenfunctions for $`B_r`$ and $`B_\phi `$ for $`𝒟=200`$ and for different $`R_v`$ in the “disk+halo” configuration. With this choice of the dynamo number, a maximum in the curve $`\gamma (R_v)`$ occurs around $`R_v=7.5`$. An obvious feature is that the vertical velocity pushes the zeros of $`B_r`$ and $`B_\phi `$ outwards. For the dynamo to work, it is essential that $`B_r`$ changes sign in the disk since only then can magnetic flux of the opposite sign to that in the midplane leave through the disk surface (Ruzmaikin et al. 1988). When $`R_v`$ takes large values, the radial component of the magnetic field $`B_r`$ has no zero inside the disk and thus the dynamo cannot work ($`\gamma <0`$). This is the technical reason for the decrease of $`\gamma `$ with $`R_v`$. The physical reason for it is investigated in Sect. 4.3. The increase of the growth rate at small $`R_v`$ will be explained in Sect. 4.3. We note however that when we decrease the truncation level of the free-decay mode expansion down to $`N=0`$ (one doublet), the maximum disappears.

The magnetic field configuration in 2D geometry is shown in Fig. 5. It shows both the effect of increasing the dynamo number and the effect of increasing the vertical magnetic Reynolds number. When one decreases the dynamo number from $`𝒟=50`$ to $`𝒟=200`$, in the absence of any vertical velocity, the magnetic field becomes more concentrated at the outer edge of the disk (Figs. 5a and 5b). If the vertical magnetic Reynolds number is then increased from $`R_v=0`$, the magnetic field is radially redistributed over the whole disk (Fig. 5c) and then advected towards the halo (Fig. 5d). From panel (b) to (d), the growth rate increases as shown by the curve for $`𝒟=200`$ in Fig. 3a. Thus, the magnetic field diffusion becomes less efficient as the vertical scale of the magnetic field increases with $`R_v`$. This leads to an increase of the growth rate up towards its maximum at $`R_v=7.5`$ (Fig. 5d). Note that in the halo the poloidal magnetic field lines become more and more aligned with lines of constant angular velocity as $`R_v`$ increases. Note also that the mode structure of Fig. 5d (where Re($`\gamma `$) is maximum for $`𝒟=200`$) is very similar to the mode structure for $`𝒟=50,R_v=3`$, where that mode has attained its maximum growth rate.

### 4.2 Dipolar modes

We now study the linear behaviour with vertical velocity for positive $`𝒟`$ where the dipolar non-oscillatory mode is dominant. Figure 6 shows that there is no maximum in the growth rates; at $`R_v`$ of order $`1`$ or even less the dynamo is switched off.

The magnetic field structure in 2D geometry is shown in Fig 7. The magnetic field structures for $`𝒟=50`$ and $`R_v=0`$ and $`0.75`$ show that the field is advected vertically outwards and the poloidal field lines are aligned with the $`\mathrm{\Omega }`$ contours lines very quickly.

### 4.3 The origin of the dynamo enhancement

In this section, we explain why a vertical velocity can increase the growth rate of the dynamo. Therefore, we plot the different terms of Eq. (16) and Eq. (17) obtained with the free-decay mode expansion. To avoid any ambiguity with the choice of normalization, in particular when comparing results at different $`R_v`$, we rewrite the equations in the following way. We multiply Eq. (16) and Eq. (17) by $`B_r`$ and $`B_\phi `$ respectively, and we integrate over the disk thickness. We thus obtain

$`\gamma =\gamma _\alpha ^r+\gamma _{\mathrm{diff}}^r+\gamma _{\mathrm{stret}}^r,`$ (36)
$`\gamma =\gamma _\omega ^\phi +\gamma _{\mathrm{diff}}^\phi +\gamma _{\mathrm{stret}}^\phi ,`$ (37)
with

$`\gamma _\alpha ^r`$ $`=`$ $`R_\alpha {\displaystyle \frac{_0^1B_r(\alpha B_\phi )^{}dz}{_0^1B_r^2dz}},`$ (38)
$`\gamma _\omega ^\phi `$ $`=`$ $`R_\omega {\displaystyle \frac{_0^1B_rB_\phi dz}{_0^1B_\phi ^2dz}},`$ (39)
$`\gamma _{\mathrm{diff}}^{r,\phi }`$ $`=`$ $`{\displaystyle \frac{_0^1B_{r,\phi }B_{r,\phi }^{\prime \prime }dz}{_0^1B_{r,\phi }^2dz}}={\displaystyle \frac{_0^1\left(B_{r,\phi }^{}\right)^2dz}{_0^1B_{r,\phi }^2dz}},`$ (40)
$`\gamma _{\mathrm{stret}}^{r,\phi }`$ $`=`$ $`R_v{\displaystyle \frac{_0^1B_{r,\phi }(V_zB_{r,\phi })^{}dz}{_0^1B_{r,\phi }^2dz}}{\displaystyle \frac{R_v}{2}}{\displaystyle \frac{_0^1V_z^{}B_{r,\phi }^2dz}{_0^1B_{r,\phi }^2dz}},`$ (41)
where we have used a vacuum boundary condition at $`z=h`$ to rewrite $`\gamma _{\mathrm{stret}}`$, $`\gamma _{\mathrm{diff}}`$. This form of the equations allows us to investigate the different physical contributions to the growth rate $`\gamma `$. Note that Eqs. (36) and (37) are not explicitly time-dependent. The terms $`\gamma _\alpha ^r`$ and $`\gamma _\omega ^\phi `$ are the contributions to the growth rate $`\gamma `$ from the $`\alpha `$-effect and $`\omega `$-effect. The terms $`\gamma _{\mathrm{diff}}^{r,\phi }`$ are the contributions from dissipation which are always negative. Finally, the terms $`\gamma _{\mathrm{stret}}^{r,\phi }`$ are the contributions to $`\gamma `$ from stretching the magnetic field. The profile for $`V_z`$ is always antisymmetric with respect to the disk midplane and $`V_z^{}>0`$ for $`z>0`$. Thus, these stretching terms $`\gamma _{\mathrm{stret}}^{r,\phi }`$ are always negative. In particular, with a linear profile for $`V_z`$, they are both equal to $`R_v/2`$.

The general picture is as follows. The positive contribution to $`\gamma `$ comes from $`\gamma _\alpha ^r`$ and $`\gamma _\omega ^\phi `$ whereas dissipation has a negative contribution. In the presence of a vertical velocity, an additional negative contribution comes from the stretching of the magnetic field $`\gamma _{\mathrm{stret}}^{r,\phi }`$. But, at the same time, the stretching reduces dissipation by increasing the scale of the magnetic field (see Eq. (40)). Thus, depending on the relative values of the reduction of dissipation and of the negative stretching term, the growth rate is enhanced or not (a maximum occurs or not). Obviously, when the magnetic field is completely stretched, dissipation cannot be reduced anymore. However, this does not necessarily correspond to the maximum of the growth rate as the terms $`\gamma _\alpha ^r`$ and $`\gamma _\omega ^\phi `$ evolve themselves with $`R_v`$ (which quantifies the vertical velocity). These different quantities are plotted in Fig. 8 in a case where a maximum occurs and in a case with no maximum.

We note that dissipation increases again after being reduced by the stretching of the magnetic field in the vertical direction. This increased dissipation is due to the formation of boundary layers near $`z=h`$. We also note that the term $`\gamma _\alpha ^r`$ decreases whereas $`\gamma _\omega ^\phi `$ increases. The evolution of these terms is linked to the disappearance of the positive area of $`B_r`$ (with the convention of sign of Fig. 4). Finally, note that at large $`R_v`$, the growth rate $`\gamma `$ always decreases as the stretching terms $`\gamma _{\mathrm{stret}}^{r,\phi }`$ become dominant.

In short, the growth rate can increase with $`R_v`$, thus forming a maximum, because dissipation is reduced by the stretching of the magnetic field.

## 5 Nonlinear behaviour with vertical velocities

The nonlinear behaviour has only been investigated with the finite-difference model. We restrict ourselves here to the model “disk+halo” ($`\eta _{\mathrm{halo}}=\eta _{\mathrm{disk}}`$) and to negative dynamo numbers considering only quadrupolar symmetry which, for $`R_v=0`$, is the dominant symmetry for moderate values of $`𝒟<0`$. Solutions are all saturated and non-oscillatory, as in the corresponding linear modes. We study three different cases of nonlinearity: magnetic buoyancy, $`\alpha `$-quenching and a combination of both. For a comparison of nonlinear solutions with a different parametrization of magnetic buoyancy see earlier works by Torkelsson & Brandenburg (1994a, b).

### 5.1 Magnetic buoyancy

We parameterize the effect of magnetic buoyancy by assuming that the vertical velocity depends on the horizontal magnetic field,

$$V_z(z)=V_{z0}\frac{z}{h}\sqrt{B_r^2+B_\phi ^2}/B_{\mathrm{eq}},$$
(42)
where $`B_{\mathrm{eq}}`$ is the equipartition magnetic field $`B_{\mathrm{eq}}=\sqrt{\rho \mu _0}c_\mathrm{s}\mathrm{Ma}_{\mathrm{turb}}`$, $`\rho `$ is the density, $`c_\mathrm{s}`$ the sound speed and $`\mathrm{Ma}_{\mathrm{turb}}`$ a turbulent Mach number, treated here as a free parameter. The profiles of density $`\rho `$ and sound speed $`c_\mathrm{s}`$ are taken from a simple model of a cool accretion disk embedded in a hot corona (Brandenburg et al. 2000; see also Dobler et al. 1999).

We define an effective vertical magnetic Reynolds number

$$R_v^{(\mathrm{eff})}=R_v\sqrt{B_r^2+B_\phi ^2}/B_{\mathrm{eq}}_{\mathrm{disk}},$$
(43)
where $`_{\mathrm{disk}}`$ denotes averaging over the disk. Fig. 9 shows the magnetic energy and $`R_v^{(\mathrm{eff})}`$ as a function of time for $`𝒟=200`$ and $`R_v=7.5`$. This value of $`R_v`$ corresponds to the location of the maximum of the growth rate $`\gamma (R_v)`$ in the linear model (cf. Fig. 3a, upper curve). Starting from an initial magnetic field strength $`B_{\mathrm{init}}=10^5`$, $`B`$ grows and therefore $`R_v^{(\mathrm{eff})}`$ increases; so the growth rate roughly follows the upper curve in Fig. 3a. Thus, as long as $`R_v^{(\mathrm{eff})}`$ is approximately less than the position of the maximum of $`\gamma `$, i.e. the horizontal magnetic field is smaller than the equipartition field, $`\gamma `$ increases with time, which results in the super-exponential growth of the magnetic energy seen for $`10t28`$ (Fig. 9). After $`t28`$, when the horizontal field is comparable with the equipartition field, $`\gamma `$ decreases and eventually the magnetic field approaches its saturation level. This level is reached when $`R_v^{(\mathrm{eff})}R_v21.5`$, where $`R_v`$ is the zero of $`\gamma (R_v)`$ in Fig. 3a, upper curve.

In this subsection we assume that the only nonlinearity is due to magnetic buoyancy and thus appears in the vertical velocity term only. In particular, we do not consider here any dependence of $`\alpha `$ or $`\eta `$ on $`B`$. In this case the dimensionless form of the induction equation in the steady state reads

$$\widehat{}B+\times \left[\left(R_v\sqrt{B_r^2+B_\phi ^2}/B_{\mathrm{eq}}\right)u\times B\right]=0,$$
(44)
where $`\widehat{}`$ is a linear operator and $`u`$ is a vector function of position, $`u=(0,0,z/h)`$ in the case considered here. Let $`b`$ be the solution of Eq. (44) for $`R_v=1`$. Then, for arbitrary $`R_v`$, the solution is $`b/R_v`$. From this, it follows that the saturation value of magnetic energy is

$$E_{\mathrm{mag}}^{\mathrm{sat}}=\frac{1}{R_v^2}E_{\mathrm{mag}}^{\mathrm{sat}}|_{R_v=1}\frac{1}{R_v^2},$$
(45)
and the structure of the saturated magnetic field is independent of $`R_v`$. This is confirmed with our numerical model.

The magnetic field is extended into the whole halo and almost completely stretched. Therefore, adding a wind term such that

$$V_z(z)=\left(V_{z0}^{\mathrm{buoy}}\sqrt{B_r^2+B_\phi ^2}/B_{\mathrm{eq}}+V_{z0}^{\mathrm{wind}}\right)\frac{z}{h}$$
(46)
and fixing $`R_v^{\mathrm{buoy}}=R_v^{(\mathrm{eff})}`$ while increasing $`R_v^{\mathrm{wind}}=R_v`$, cannot lead to enhancement of magnetic energy. The wind is carrying magnetic field away and $`E_{\mathrm{mag}}`$ is decreasing.

### 5.2 $`\alpha `$-quenching

We now parameterize the back-reaction of the magnetic field on the turbulent $`\alpha `$-effect by making $`\alpha `$ depend on $`B`$,

$$\alpha (r,z)=\alpha _0(r)\mathrm{sin}\left(\pi \frac{z}{h}\right)\frac{\xi _\alpha (r)}{1+B^2/B_{\mathrm{eq}}^2},$$
(47)
where $`B_{\mathrm{eq}}`$ is again the equipartition magnetic field.

Figure 10 shows the magnetic energy as a function of $`R_v`$ for the three dynamo numbers $`𝒟=50`$, $`𝒟=100`$ and $`𝒟=200`$. From this we can see the following features:

* $`E_{\mathrm{mag}}`$ has a maximum as a function of $`R_v`$ in the nonlinear regime;
* but, this maximum is shifted towards much smaller $`R_v`$ (of order $`0.1`$), compared to that of the growth rate in the linear regime, and it is also less pronounced;
* the value of $`R_v`$ where the maximum in $`E_{\mathrm{mag}}`$ occurs, is almost independent of $`𝒟`$;
* the magnetic energy scales roughly with the dynamo number, $`E_{\mathrm{mag}}𝒟`$;
* all modes remain non-oscillatory when increasing $`R_v`$.

The effect of vertical velocity on the structure of the magnetic field for the model with $`\alpha `$-quenching can be seen in Fig. 11 where it is shown for $`𝒟=50`$. But looking at the magnetic field configuration for $`𝒟=100`$ and $`𝒟=200`$, one finds that the field structure changes only weakly with the dynamo number.

### 5.3 $`\alpha `$-quenching together with magnetic buoyancy

In this section, we combine $`\alpha `$-quenching, Eq. (47), with magnetic buoyancy, Eq. (42), and compare the results with the case of $`\alpha `$-quenching with wind. In Fig. 12, magnetic energy is plotted against $`R_v^{(\mathrm{eff})}`$ for these two models and for $`𝒟=50`$. Evidently, the two curves behave similarly, supporting the notion that wind and magnetic buoyancy have similar effects on the dynamo. For $`R_v^{(\mathrm{eff})}1`$, the field structure is solely determined by $`\alpha `$-quenching and thus $`E_{\mathrm{mag}}`$ is independent of $`R_v^{(\mathrm{eff})}`$. For $`R_v^{(\mathrm{eff})}1`$, on the other hand, $`\alpha `$-quenching becomes unimportant and $`E_{\mathrm{mag}}(R_v^{(\mathrm{eff})})^2`$, as discussed in Sect. 5.1. The largest disagreement between the two curves occurs for intermediate values of $`R_v^{(\mathrm{eff})}`$, i.e. where both $`\alpha `$-quenching and magnetic buoyancy are important, and where the field has not yet reached its asymptotic shape (Fig. 11c).

Perfect agreement between the two curves cannot be expected, because the averaging in Eq. (43) was somewhat arbitrary; the agreement could perhaps be improved by choosing a different definition of $`R_v^{(\mathrm{eff})}`$. The fact that with magnetic buoyancy the magnetic energies are larger for a given value of $`R_v^{(\mathrm{eff})}`$ can therefore not be interpreted as a relative enhancement of dynamo action.

## 6 Conclusions

Disk dynamos have been studied intensively over the past 20 years, but only recently the effects of outflows and winds have been included (e.g. Campbell 1999; Dobler et al. 1999; Brandenburg et al. 2000). It was immediately clear that these changes in the boundary conditions affect the structure of the solutions. The present work is an attempt to provide a more thorough survey of solutions with wind for different values and sign of the dynamo number. One of the effects that we have highlighted is the enhancement of the kinematic dynamo growth rate for intermediate strengths of the wind. In the nonlinear case this corresponds to a phase shortly before saturation where the growth becomes super-exponential.

In the present work we have deliberately isolated the dynamo problem from the wind problem. In reality the two aspects together constitute a more complicated self-consistent dynamo problem that needs to be looked at in more detail in future. One step in that direction was taken by Brandenburg et al. (2000) where the feedback on the flow from the Lorentz force was taken fully into account. A more formidable task would be to obtain self-consistent dynamo action in three dimensions without assuming an $`\alpha `$-effect, thus obtaining dynamo action and outflows in a self-consistent manner from first principles. In the context of accretion disks, this would require simulating turbulent dynamo action, where the turbulence itself is driven by dynamo-generated magnetic fields (Brandenburg et al. 1995). This would be of great interest not only from the astrophysical point of view, but also from that of dynamo theory in general, because, as was recently realized, the catastrophic quenching of the standard $`\alpha \omega `$-dynamo could be alleviated in the presence of open boundary conditions which may substantially improve the working of the dynamo (Blackman & Field 2000). This is simply because the $`\alpha `$-effect generates helical fields, but at the same time magnetic helicity is conserved in a closed domain. Thus, being able to get rid of helical fields through open boundaries may be vital for the dynamo (Brandenburg & Dobler 2001). Furthermore, the production of fields at large scales could be facilitated by the production – and subsequent loss through boundaries – of small scale fields with opposite sign of magnetic helicity. This latter aspect could not be addressed with mean-field models and would require a direct simulation of accretion disk turbulence in a global model.

###### Acknowledgements.

We acknowledge financial support from PPARC (Grant PPA/G/S/1997/00284) and from the Leverhulme Trust (Grant F/125/AL). We thank D. Sokoloff for helpful discussions. The use of the PPARC supported parallel computer GRAND in Leicester and the PPARC funded Compaq MHD Cluster in St Andrews is acknowledged.