format: text

text:
# Topological and confining properties of Abelian-projected SU(3)-QCD

## 1 INTRODUCTION

In the present talk, we shall mostly discuss various nonperturbative properties of the effective low-energy theory of the $`SU(3)`$-QCD , which models confinement of quarks as the dual Meissner effect . The partition function of this $`[U(1)]^2`$ magnetically gauge-invariant theory reads

$$𝒵=\left|\mathrm{\Phi }_a\right|𝒟\left|\mathrm{\Phi }_a\right|𝒟\theta _a𝒟𝐁_\mu \delta \left(\underset{a=1}{\overset{3}{}}\theta _a\right)\times $$

$$\times \mathrm{exp}\{d^4x[\frac{1}{4}(𝐅_{\mu \nu }+𝐅_{\mu \nu }^{(c)})^2+$$

$$+\underset{a=1}{\overset{3}{}}\left(\frac{1}{2}\right|(_\mu 2ig_m𝐞_a𝐁_\mu )\mathrm{\Phi }_a|^2+$$

$$+\lambda \left(\right|\mathrm{\Phi }_a|^2\eta ^2)^2)]\}.$$
(1)
Here, $`g_m`$ is the magnetic coupling constant, related to the QCD coupling constant $`g`$ as $`g_m=\frac{4\pi }{g}`$, and $`𝐞_a`$’s are the root vectors of $`SU(3)`$, whose explicit form is $`𝐞_1=(1,0)`$, $`𝐞_2=(\frac{1}{2},\frac{\sqrt{3}}{2})`$, $`𝐞_3=(\frac{1}{2},\frac{\sqrt{3}}{2})`$. Next, in Eq. (1), $`𝐅_{\mu \nu }`$ stands for the field strength tensor of the field $`𝐁_\mu `$ dual to the field $`𝐀_\mu (A_\mu ^3,A_\mu ^8)`$, and $`\mathrm{\Phi }_a=\left|\mathrm{\Phi }_a\right|\mathrm{e}^{i\theta _a}`$ are the dual Higgs fields describing the condensation of Cooper pairs of Abelian-projected monopoles.

Note that the phases $`\theta _a`$’s are related to each other by the constraint $`\underset{a=1}{\overset{3}{}}\theta _a=0`$, imposed by the respective $`\delta `$-function on the right-hand side of Eq. (1). This constraint reflects the dependence of the monopoles of three kinds of each other. Such a dependence is inspired by the fact that the monopole magnetic charges are distributed over the lattice defined by the root vectors, whose sum vanishes.

Finally, in Eq. (1) we have introduced the field strength tensor $`𝐅_{\mu \nu }^{(c)}`$ of an external quark of the colour $`c=R,B,G`$ (red, blue, green, respectively), which obeys the equation $`_\mu \stackrel{~}{𝐅}_{\mu \nu }^{(c)}=g𝐐^{(c)}j_\nu `$. Here, $`\stackrel{~}{𝒪}_{\mu \nu }\frac{1}{2}\epsilon _{\mu \nu \lambda \rho }𝒪_{\lambda \rho }`$, $`j_\nu (x)\underset{C}{}𝑑x_\nu (\tau )\delta (xx(\tau ))`$, and $`𝐐^{(c)}`$’s are the charges of a quark of the colour $`c`$ with respect to the Cartan subgroup of $`SU(3)`$: $`𝐐^{(R)}=(\frac{1}{2},\frac{1}{2\sqrt{3}})`$, $`𝐐^{(B)}=(\frac{1}{2},\frac{1}{2\sqrt{3}})`$, $`𝐐^{(C)}=(0,\frac{1}{\sqrt{3}})`$. The present lattice data  indicate that in the regime of the model (1) corresponding to the real QCD, the coupling constant $`\lambda `$ is much larger than one, namely $`\lambda 65`$. This makes it reasonable to consider this model in the London limit, $`\lambda \mathrm{}`$. In this limit, the model (1) allows for an exact string representation, and the resulting string effective action reads

$$S_c=\pi ^2d^4xd^4yD_m^{(4)}(xy)\times $$

$$\times \left[\eta ^2\overline{\mathrm{\Sigma }}_{\mu \nu }^a(x)\overline{\mathrm{\Sigma }}_{\mu \nu }^a(y)+\frac{8}{3g_m^2}j_\mu (x)j_\mu (y)\right].$$
(2)
Here, $`m=g_H\eta `$ is the mass of the dual vector bosons with $`g_H\sqrt{6}g_m`$ standing for their magnetic charge, which they acquire due to the Higgs mechanism. Next, $`D_m^{(4)}(x)=\frac{m}{4\pi ^2|x|}K_1(m|x|)`$ is the propagator of these bosons, where from now on $`K_\nu `$’s stand for the modified Bessel functions. In Eq. (2), we have also introduced the notation $`\overline{\mathrm{\Sigma }}_{\mu \nu }^a=\mathrm{\Sigma }_{\mu \nu }^a2s_a^{(c)}\mathrm{\Sigma }_{\mu \nu }`$. In this expression, $`s_a^{(c)}`$’s stand for certain numbers equal to $`0`$ and $`\pm 1`$, which obey the relation $`𝐐^{(c)}=\frac{1}{3}s_a^{(c)}𝐞_a`$, and $`\mathrm{\Sigma }_{\mu \nu }^a(x)\underset{\mathrm{\Sigma }_a}{}𝑑\sigma _{\mu \nu }(x_a(\xi ))\delta (xx_a(\xi ))`$ is the vorticity tensor current defined on the closed dual string world sheet $`\mathrm{\Sigma }_a`$ with $`\xi (\xi ^1,\xi ^2)`$. Note that owing to the one-to-one correspondence existing between $`\mathrm{\Sigma }_{\mu \nu }^a`$’s and the multivalued parts of $`\theta _a`$’s, the three vorticity tensor currents are subject to the constraint $`\underset{a=1}{\overset{3}{}}\mathrm{\Sigma }_{\mu \nu }^a=0`$, which stems from the analogous constraint imposed on $`\theta _a`$’s. Finally, in the definition of $`\overline{\mathrm{\Sigma }}_{\mu \nu }^a`$, we have denoted by $`\mathrm{\Sigma }_{\mu \nu }`$ the vorticity tensor current defined on an arbitrary surface $`\mathrm{\Sigma }`$ bounded by the contour $`C`$, which is the world sheet of the open dual string, ending up at a quark and an antiquark.

The $`\mathrm{\Sigma }_{\mu \nu }\times \mathrm{\Sigma }_{\mu \nu }`$-interaction in the action (2) can be shown  to describe confinement of quarks, whereas the $`j_\mu \times j_\mu `$-interaction clearly describes their Yukawa interaction at small distances. Notice also that in what follows we shall be interested in the string effective actions, rather than the measure of integration over world-sheet coordinates $`x_a(\xi )`$’s. The Jacobian emerging during the change of integration variables $`\theta _ax_a(\xi )`$, which should be accounted for in this measure, has been evaluated in Ref. .

## 2 INCLUDING THE $`\mathrm{\Theta }`$-TERM

Let us now add to the Lagrangian of the model (1) the following $`\mathrm{\Theta }`$-term:

$$\mathrm{\Delta }=\frac{i\mathrm{\Theta }g_m^2}{4\pi ^2}\left(𝐅_{\mu \nu }+𝐅_{\mu \nu }^{(c)}\right)\left(\stackrel{~}{𝐅}_{\mu \nu }+\stackrel{~}{𝐅}_{\mu \nu }^{(c)}\right).$$

In the London limit, the string representation of such an extended partition function then reads

$$𝒵_\mathrm{\Theta }^c=\mathrm{exp}[\frac{2}{3}(\frac{(2\pi )^2}{g_m^2}+\frac{(\mathrm{\Theta }g_m)^2}{\pi ^2})\times $$

$$\times \underset{C}{}dx_\mu \underset{C}{}dy_\mu D_m^{(4)}(xy)8(\pi \eta )^2\times $$

$$\times d^4xd^4y\mathrm{\Sigma }_{\mu \nu }(x)D_m^{(4)}(xy)\mathrm{\Sigma }_{\mu \nu }(y)]\times $$

$$\times \mathrm{exp}\{2s_a^{(c)}[d^4xd^4y(\frac{i\mathrm{\Theta }}{3}\stackrel{~}{\overline{\mathrm{\Sigma }}}_{\mu \nu }^a(x)j_\nu (y)\times $$

$$\times _\mu D_m^{(4)}(xy)+2(\pi \eta )^2\mathrm{\Sigma }_{\mu \nu }^a(x)D_m^{(4)}(xy)\times $$

$$\times \mathrm{\Sigma }_{\mu \nu }(y))\frac{i\mathrm{\Theta }}{3}\widehat{L}(\mathrm{\Sigma }_a,C)]\}_{\mathrm{\Sigma }_a}.$$
(3)
Here,

$$\widehat{L}(\mathrm{\Sigma }_a,C)$$

$$d^4xd^4y\stackrel{~}{\mathrm{\Sigma }}_{\mu \nu }^a(x)j_\nu (y)_\mu D_0^{(4)}(xy)$$

is the 4D Gauss linking number of the contour $`C`$ with the closed world sheet $`\mathrm{\Sigma }_a`$, where $`D_0^{(4)}(x)=\frac{1}{4\pi ^2x^2}`$. The average $`\mathrm{}_{\mathrm{\Sigma }_a}`$ is formally defined with respect to $`S_c[\mathrm{\Sigma }=0]`$, but its exact meaning will be discussed below.

The first argument in the first exponential factor on the right-hand side of Eq. (3) shows that due to the $`\mathrm{\Theta }`$-term quarks acquire a nonvanishing magnetic charge , i.e. become dyons. We also see that there appear short- and long-ranged interactions of such dyons with closed dual strings, as well as the short-ranged dyon–open-string interaction. The above-mentioned long-ranged interaction, described by the linking number, can be viewed as a scattering of dyons by the closed dual strings, which thus play the rôle of solenoids carrying electric flux. According to Eq. (3), such a scattering, which is nothing else but the four-dimensional analogue of the Aharonov-Bohm effect , takes place at $`\mathrm{\Theta }3\pi \times (\mathrm{integer})`$. Note that in the $`SU(2)`$-case, the respective critical values of $`\mathrm{\Theta }`$ have been found  to be equal to $`2\pi \times (\mathrm{integer})`$.

Let us now carry out the average $`\mathrm{}_{\mathrm{\Sigma }_a}`$, taking into account that at zero temperature closed dual strings with opposite winding numbers are known to form virtual bound states, called vortex loops . The summation over the grand canonical ensemble of these objects has been performed in Ref.  and yields an effective sine-Gordon theory of two antisymmetric spin-1 tensor fields. In the dilute gas approximation, which is relevant to the reality since vortex loops are only virtual (and therefore small-sized) objects, such a theory enables one to calculate correlators of loops. Averaging the $`\mathrm{\Sigma }_a`$-dependent exponential factor on the right-hand side of Eq. (3) in the sense of this effective theory by making use of the cumulant expansion and accounting in this average only for the contribution of the dominant, bilocal, irreducible correlator of vortex loops, we get

$$𝒵_\mathrm{\Theta }^c=\mathrm{exp}\{8d^4xd^4y[(\pi \eta )^2\mathrm{\Sigma }_{\mu \nu }(x)\times $$

$$\times D_{M_1}^{(4)}(xy)\mathrm{\Sigma }_{\mu \nu }(y)+\frac{i\mathrm{\Theta }}{3}\stackrel{~}{\mathrm{\Sigma }}_{\mu \nu }(x)j_\nu (y)\times $$

$$\times _\mu D_m^{(4)}(xy)]\frac{2}{3}(\frac{(2\pi )^2}{g_m^2}+\frac{(\mathrm{\Theta }g_m)^2}{\pi ^2})\times $$

$$\times \frac{1}{g_H^2+g_D^2}\underset{C}{}dx_\mu \underset{C}{}dy_\mu \times $$

$$\times (g_H^2D_{M_1}^{(4)}(xy)+g_D^2D_0^{(4)}(xy))\}.$$
(4)
Here, $`g_D=\frac{2\pi \sqrt{\zeta }}{\mathrm{\Lambda }^2}`$ describes the contribution to the magnetic charge of the dual vector bosons, stemming from the Debye screening of those in the gas of electric vortex loops. In this formula, $`\mathrm{\Lambda }=\sqrt{\frac{L}{a^3}}`$ is the ultraviolet momentum cutoff with $`L`$ standing for the typical distance between the neighbours in the gas of vortex loops and $`a`$ denoting the typical size of the loop, $`aL`$. Next, $`\zeta \mathrm{e}^{S_0}`$ is the Boltzmann factor of a single vortex loop with the action $`S_0`$ equal to the string tension times the characteristic area of the loop.

Due to the Debye screening, the mass $`m`$ of the dual vector bosons increases. Since owing to the constraint $`\underset{a=1}{\overset{3}{}}\mathrm{\Sigma }_{\mu \nu }^a=0`$ there exists the second independent type of vortex loops, there respectively appears also the second value of the Debye charge, equal to $`g_D\sqrt{3}`$. The two full masses then read $`M_1=\eta \sqrt{g_H^2+g_D^2}`$ and $`M_2=\eta \sqrt{g_H^2+3g_D^2}`$.

As it follows from Eq. (4), the above-mentioned screening changes drastically the interaction of quarks (dyons). Namely, besides the modification of the massive propagator due to the additional contribution to the mass of the dual vector bosons, there also appears a novel massless interaction. It is also worth noting that, as it should be, when the effect of screening is disregarded, i.e. $`g_Dg_H`$, one recovers the classical result, which is nothing else, but the $`\mathrm{\Sigma }_a`$-independent part of Eq. (3).

## 3 APPLICATIONS TO THE STOCHASTIC VACUUM MODEL

Let us now discuss in more details the confining properties of the model (1) in the London limit. Namely, let us consider the bilocal correlator of electric field strengths, $`𝐟_{\mu \nu }=_\mu 𝐀_\nu _\nu 𝐀_\mu `$, which plays the major rôle in the so-called stochastic vacuum model  (see Ref.  for reviews). Within this model, such a correlator is parametrized by the two coefficient functions as follows:

$$f_{\mu \nu }^i(x)f_{\lambda \rho }^j(0)_{𝐀_\mu ,𝐣_\mu ^\mathrm{m}}=\delta ^{ij}\{(\delta _{\mu \lambda }\delta _{\nu \rho }\delta _{\mu \rho }\delta _{\nu \lambda })\times $$

$$\times D\left(x^2\right)+\frac{1}{2}[_\mu (x_\lambda \delta _{\nu \rho }x_\rho \delta _{\nu \lambda })+$$

$$+_\nu (x_\rho \delta _{\mu \lambda }x_\lambda \delta _{\mu \rho })]D_1\left(x^2\right)\}.$$
(5)
In this formula, $`\mathrm{}_{𝐀_\mu }`$ stands for the average over free diagonal gluons, and $`\mathrm{}_{𝐣_\mu ^\mathrm{m}}`$ denotes a certain average over monopoles, which provides the condensation of their Cooper pairs. It is the coupling of the $`𝐁_\mu `$-field, dual to the field $`𝐀_\mu `$, to $`𝐣_\mu ^\mathrm{m}`$’s, which makes both functions $`D`$ and $`D_1`$ nontrivial and finally adequate to those of the real QCD. Namely, these functions turn out to have the following form :

$$D=\frac{m^2M_1}{4\pi ^2}\frac{K_1(M_1|x|)}{|x|},$$

$$D_1=\frac{g_D^2}{\pi ^2\left(g_H^2+g_D^2\right)|x|^4}+\frac{m^2}{2\pi ^2M_1x^2}\times $$

$$\times \left[\frac{K_1(M_1|x|)}{|x|}+\frac{M_1}{2}\left(K_0(M_1|x|)+K_2(M_1|x|)\right)\right].$$

In the infrared limit, $`|x|M_1^1`$, the respective asymptotic behaviours read

$$D\frac{(mM_1)^2}{4\sqrt{2}\pi ^{\frac{3}{2}}}\frac{\mathrm{e}^{M_1|x|}}{(M_1|x|)^{\frac{3}{2}}},$$
(6)
$$D_1\frac{g_D^2}{\pi ^2\left(g_H^2+g_D^2\right)|x|^4}+\frac{(mM_1)^2}{2\sqrt{2}\pi ^{\frac{3}{2}}}\frac{\mathrm{e}^{M_1|x|}}{(M_1|x|)^{\frac{5}{2}}}.$$
(7)
One can see that according to the above-presented asymptotics, the inverse correlation length of the vacuum in the model under study is equal to the smallest of the two full masses, $`M_1`$. Besides that, owing to the screening, the function $`D_1`$ contains the novel nonperturbative $`\frac{1}{|x|^4}`$-term, which might be important for modelling the Lüscher term  in the quark-antiquark potential . Subtracting from the function $`D_1`$ this contribution, which has the same functional form as the pure perturbative contribution in the real QCD, we see that the remained part of the asymptotics (7) together with the asymptotics (6) are in a good agreement with the QCD lattice measurements of the functions $`D`$ and $`D_1`$  (see Ref.  for recent developments and Ref.  for reviews). Namely, due to the preexponential behaviour, $`DD_1`$ in the infrared limit, $`|x|M_1^1`$, whereas in the ultraviolet limit, $`DD_1`$.

It is also worth noting that when the screening is disregarded, i.e. $`g_Dg_H`$, the inverse correlation length of the vacuum goes over to $`m`$, and the classical expressions  (see Ref.  for the $`SU(2)`$-case) for the functions $`D`$ and $`D_1`$ recover. In particular, the nonperturbative screening-motivated $`\frac{1}{|x|^4}`$-contribution to the function $`D_1`$ vanishes in this limit.

Besides the above-discussed modifications of the classical expressions for the correlators of electric field strengths inspired by screening, this effect changes also the classical expression for the propagator of the dual vector bosons. Indeed, the latter one reads

$$B_\mu ^a(x)B_\nu ^b(0)=\delta ^{ab}\delta _{\mu \nu }D_m^{(4)}(x)$$
(8)
with $`B_\mu ^a𝐞_a𝐁_\mu `$, whereas with the account for screening it changes to

$$\delta _{\mu \nu }\{\frac{1}{3}D_m^{(4)}(x)+\frac{g_D^2}{2}[2\frac{g_H^2+2g_D^2}{\left(g_H^2+3g_D^2\right)\left(g_H^2+g_D^2\right)}\times $$

$$\times D_0^{(4)}(x)+\frac{g_H^2}{g_D^2}(\frac{1}{3\left(g_H^2+3g_D^2\right)}D_{M_2}^{(4)}(x)+$$

$$+\frac{1}{g_H^2+g_D^2}D_{M_1}^{(4)}(x))]\}$$
(9)
for $`a=b`$ and

$$\delta _{\mu \nu }\{\frac{1}{3}D_m^{(4)}(x)+\frac{g_D^2}{2}[\frac{2g_D^2}{\left(g_H^2+3g_D^2\right)\left(g_H^2+g_D^2\right)}\times $$

$$\times D_0^{(4)}(x)+\frac{g_H^2}{g_D^2}(\frac{1}{3\left(g_H^2+3g_D^2\right)}D_{M_2}^{(4)}(x)$$

$$\frac{1}{g_H^2+g_D^2}D_{M_1}^{(4)}(x))]\}$$
(10)
for $`ab`$. We see that the screening makes the propagator of the dual vector bosons nonvanishing even for $`ab`$. As it should be, in the limit when the screening is disregarded, $`g_Dg_H`$, these off-diagonal components of the propagator given by Eq. (10) vanish, whereas the diagonal ones given by Eq. (9) go over to the classical expression, so that Eq. (8) recovers.

Contrary to the four-dimensional case, in three dimensions the present lattice data allow one to assume that Abelian-projected monopoles form a gas. Such a gas of $`SU(3)`$-monopoles has for the first time been considered in Ref. , and its partition function reads

$$𝒵=1+\underset{N=1}{\overset{\mathrm{}}{}}\frac{\zeta ^N}{N!}(\underset{n=1}{\overset{N}{}}d^3z_n\underset{a_n=\pm 1,\pm 2,\pm 3}{})\times $$

$$\times \mathrm{exp}\left[\frac{g_m^2}{4\pi }\underset{n<k}{}\frac{𝐞_{a_n}𝐞_{a_k}}{|𝐳_n𝐳_k|}\right].$$
(11)
In this formula, $`\zeta \mathrm{exp}\left(\frac{\mathrm{const}}{g^2}\right)`$ stands for the Boltzmann factor of a single monopole, and $`𝐞_a=𝐞_a`$. The string representation of the Wilson loop in this gas, constructed in Ref. , turned out to be alternative to the one of the $`SU(2)`$-case, found in Ref. . (See Ref.  for the discussion of Polyakov loops and their correlators in the $`SU(2)`$ monopole gas.) By virtue of this representation in the approximation when the monopole gas is so dilute that its density is much less than $`\zeta `$, one can deduce the respective expressions for the functions $`D`$ and $`D_1`$. In the model under study, those are defined by Eq. (5) with the average $`\mathrm{}_{𝐣_\mu ^\mathrm{m}}`$ replaced by the average with respect to the partition function (11). Together with the contribution of the free diagonal gluons to the function $`D_1`$, these two functions read

$$D=12\pi \zeta \frac{\mathrm{e}^{m|𝐱|}}{|𝐱|},$$
(12)
$$D_1=\frac{24\pi \zeta }{(m|𝐱|)^2}\left(m+\frac{1}{|𝐱|}\right)\mathrm{e}^{m|𝐱|}.$$
(13)
Here, $`m=g_m\sqrt{3\zeta }`$ is the Debye mass of the two scalar bosons, dual to the diagonal gluons. Similarly to the above-considered four-dimensional case, we see that Eqs. (12) and (13) well agree with the lattice calculations in the real QCD . In particular, the inverse correlation length of the vacuum is now equal to $`m`$, and at the distances larger than this length, $`DD_1`$ due to the preexponential behaviour.

## 4 CONCLUSIONS

In the present talk, we have discussed various nonperturbative properties of Abelian-projected $`SU(3)`$-QCD in four and three dimensions. In the four-dimensional case, we have firstly considered the respective dual Abelian Higgs type model extended by the introduction of the $`\mathrm{\Theta }`$-term. In this way, the string representation of such an extended model has been derived, which has in particular demonstrated how the $`\mathrm{\Theta }`$-term leads to the appearance of the magnetic charge of external quarks, making out of them dyons. The critical values of $`\mathrm{\Theta }`$, at which the Aharonov-Bohm scattering of dyons over the closed dual strings disappears, have been found. Next, the effect of the Debye screening of the dual vector bosons by virtual electric vortex loops, built out of the closed dual strings, has been taken into account.

Then, the confining properties of the four- and three-dimensional $`SU(3)`$ Abelian-projected theories within the stochastic vacuum model have been addressed. In particular, in the four-dimensional case the rôle of the above-mentioned screening has been discussed. The influence of this effect to the propagators of the dual vector bosons has also been considered. Finally, the bilocal correlator of electric field strengths in the dilute three-dimensional gas of $`SU(3)`$ Abelian-projected monopoles has been evaluated.

In conclusion, the performed investigations have shown that Abelian-projected theories are not only adequate to the description of confinement of quarks in QCD, but possess also a lot of interesting nonperturbative properties themselves.

## 5 ACKNOWLEDGMENTS

The author is indebted to Prof. A. Di Giacomo for useful discussions and hospitality, and INFN for the financial support. He is also greatful to the organizers of the Euroconference ”QCD 00” (Montpellier, France, 6-13th July 2000) for an opportunity to present these results in a very stimulating atmosphere.

## 6 DISCUSSION

Dr. N. Brambilla (University of Heidelberg): You have mentioned that the London limit is in the agreement with the lattice results, but I would like to point out that a lot of the most recent lattice data show that the QCD vacuum looks like a dual superconductor at the border of type-I and type-II. Could you please comment on this issue.

D. Antonov: My statement was based on the lattice data of Ref. , which demonstrated that in the QCD-relevant regime of the effective dual theory considered, the coupling constant of the dual Higgs potential should be much larger than unity.


text:
# Mixing through shear instabilities

## 1 Introduction

Mixing is a fundamental process which profoundly affects the evolution of stars. Nonetheless, a theory of mixing that is applicable to stellar conditions is still missing. Hitherto stellar evolution theory has resorted to more or less heuristic prescriptions of mixing, the validity of which is often unclear. One of the prominent mechanisms responsible for mixing in stars is the shear or Kelvin-Helmholtz instability. This instability is particularly important since most stars are are known to be at least partly differentially rotating. Therefore, shear mixing is likely to play a significant role in stellar evolution.

Observational evidence seems to suggest that current prescriptions underestimate the efficiency of the mixing processes at work, especially in fast rotating stars. The observational evidence for mixing is increasing rapidly. Herrero et al. (1992) find that all fast rotating O-stars show significant surface He-enrichments. This seems to suggest that mixing is strong enough to transport the nuclearly processed material to the surface in a fraction of the life time of the star. Other observations include the N/C and 13C and 12C enrichments of the Red Giant Branch (e.g. Kraft 1997, Charbonnel 1995), the He and N excesses in OBA supergiants (Fransson et al. 1989), the depletion of boron in most B-type stars (Venn, Lambert & Lemke 1996) and the ratio of the number of blue to red supergiants in galaxies (Langer & Maeder 1995). For a more detailed review of the observational evidence see Maeder (1995), Kraft (1994) and references therein. But the conclusion from these observations is that both, the enrichment in CNO elements and the depletion of fragile elements such as boron can be explained if some form of mixing is introduced (Langer 1997).

However, a convincing theory of mixing under most astrophysical conditions is still missing. One of the pioneering works in this field was performed by Endal & Sophia (1978) who estimated the efficiency of mixing due to several rotationally induced instabilities in massive stars.

The Sun is also known to rotate differentially. Pinsonneault et al. (1989) treated the efficiency of rotational mixing as a free parameter and fitted it to a solar model. In the meantime, helioseismological observations from the spacecraft SOHO and from the GONG network have yielded precise measurements of the rotation rate in the solar interior. In a number of publications the effects of mixing on stellar evolution have been examined, see, e.g. Meynet & Maeder (1997) and Staritsin (1999). Recently, Denissenkov & Tout (2000) investigated the effects of ‘deep mixing’ on the evolution of red giants. However, for want of a valid fundamental theory, the mixing formalisms employed in the aforementioned studies have all been based on more or less heuristic estimates. It is the aim of this work to make steps towards a more fundamental theory of mixing.

In chemically homogeneous, stratified shear flows, the Kelvin-Helmholtz instability occurs when the destabilising effect of the relative motion in the different layers dominates over the stabilising effect of buoyancy (see e.g. Chandrasekhar 1961). The competition between the two effects is described by the Richardson number, Ri. For a non-dissipative, parallel, steady flow a simple linear stability criterion can be derived as follows: Let the horizontal velocity of the fluid be a function of $`z`$ only, i.e. $`U=U(z)`$. Now suppose that two volumes of fluid at heights $`z`$ and $`z+\mathrm{\Delta }z`$ are interchanged where the densities at those two heights are $`\rho `$ and $`\rho +\mathrm{\Delta }\rho `$, and their horizontal velocities $`U`$ and $`U+\mathrm{\Delta }U`$, respectively. The work per unit mass that must be done against gravity is

$$\mathrm{\Delta }W=g\mathrm{\Delta }\rho \mathrm{\Delta }z.$$
(1)
After the exchange both cells have an average velocity of $`\overline{U}=\frac{1}{2}(2U+\mathrm{\Delta }U)`$. Therefore, the kinetic energy available to do this work is

$$\frac{1}{2}\rho [U^2+(U+\mathrm{\Delta }U)^2\frac{1}{2}(2U+\mathrm{\Delta }U)^2]=\frac{1}{4}\rho (\mathrm{\Delta }U)^2.$$
(2)
Hence, a condition for stability is

$$\frac{1}{4}\rho (\mathrm{\Delta }U)^2<g\mathrm{\Delta }\rho \mathrm{\Delta }z,$$
(3)
or

$$\left(\frac{dU}{dz}\right)^2<4\frac{g}{\rho }\frac{d\rho }{dz}.$$
(4)
A more general criterion can be derived (see e.g. Shu 1992):

$$\mathrm{Ri}:=\frac{N^2}{(dU/dz)^2}>\frac{1}{4},$$
(5)
where $`N`$ is the Brunt-Väisälä frequency given by

$$N^2=\frac{g\delta }{H_P}\left[\left(\frac{\mathrm{ln}T}{\mathrm{ln}P}\right)_{\mathrm{ad}}\frac{\mathrm{ln}T}{\mathrm{ln}P}+\frac{\varphi }{\delta }\frac{\mathrm{ln}\mu }{\mathrm{ln}P}\right],$$
(6)
with $`\delta =(\mathrm{ln}\rho /\mathrm{ln}T)_{P,\mu }`$ and $`\varphi =(\mathrm{ln}\rho /\mathrm{ln}\mu )_{P,T}`$, where $`T`$ is temperature, $`P`$ is pressure and $`\mu `$ chemical potential. $`H_P`$ is the pressure scale height, i.e. $`(\mathrm{ln}P/z)^1`$. Ri is the so-called Richardson number. Stability is predicted if the local Richardson number exceeds 1/4 everywhere.

More detailed analytical theories of shear mixing have been developed (Maeder 1995, Canuto 1998) and some will be reviewed briefly in Sec. 2. In Sec. 3, we present the results of three-dimensional simulations of shear instabilities in stratified fluids. These numerical simulations were performed for a simple configuration and a range of initial condition. Our aim is to elucidate the nonlinear dynamics of the shear instability and to study the parametric dependence of the mixing efficiency on the Richardson number. We quantify the mixing of the fluid by introducing ‘tracer’ particles that are advected with the fluid. With conditions in stellar interiors in mind, we only consider subsonic flows. From the motion of these particles a heuristic diffusion constant can be derived which is then compared to the analytic prescriptions described in Sec. 2.

Apart from mixing by shear instabilities, there is another important candidate believed to be reponsible for extra mixing in stars, namely convective overshoot. In models of stellar evolution the fluid is assumed to be mixed whereever some convective instability criterion (such as the Schwarzschild criterion) is fulfilled, and convectively stable where this criterion is not fulfilled. Mathematically the boundary between these two regions is sharp. In reality, however, the convective fluid elements still carry some momentum when they reach this boundary and subsequently penetrate into the convectively “stable” region. This mechanism is called convective overshoot and has been investigated numerically by a number of groups (Freytag, Ludwig & Steffen 1996, Nordlund & Stein 1996, Singh, Roxburgh & Chan 1998) all of which find some degree of overshoot. This conclusion is supported by observations. The most compelling evidence stems from isochrone fitting to stellar clusters and binary systems which suggests that convective overshoot is significant (see Zahn 1991 for a review).

Unlike convective overshoot which provides mixing at the boundaries of convection zones, we investigate mixing by rotational instabilities that operate in convectively stable regions.

## 2 Theories of shear mixing

### 2.1 Linear stability analysis of a shear layer

The stability of a stratified shear layer has been studied extensively using linear theory (see e.g. Chandrasekhar 1961, Howard & Maslowe 1973, Maslowe & Thompson 1971, Hazel 1972, Morris et al. 1990). Recently, Balmforth & Morrison (1998) reviewed the conditions for instability in inviscid shear flows. In astrophysics, Kippenhahn & Thomas (1987) and MacDonald (1983) examined the redistribution of matter and angular momentum in accreting white dwarfs by means of a linear stability analysis. In this section we will present a simple example of a linear stability analysis for a plane shear layer. The results of this analysis will then be compared against our numerical simulations in Sec. 3.

The stability equation for wave-like perturbations in its inviscid and incompressible form can be written as (e.g. Chandrasekhar 1961)

$$\frac{d^2\varphi }{dz^2}\left(k^2+\frac{U^{\prime \prime }}{Uc}\frac{g\beta }{(Uc)^2}\right)\varphi =0,$$
(7)
where $`\beta =\rho ^{}/\rho `$ and $`\varphi (z)\mathrm{exp}[ik(xct)]`$ is the perturbation stream function. Eqn. (7) is sometimes referred to as the Orr-Sommerfeld equation. As an example for the linear analysis of the Kelvin-Helmholtz instability we will cite a case which was first considered by Taylor: It consits of two superposed fluids of different densities separated by a transition layer of intermediate density in which the shear velocity varies linearly from that of the lower layer to that of the upper layer, i.e.

layer 1: $`z>d`$, $`\rho =\rho _0(1ϵ)`$, $`U=U_0`$
layer 2: $`d>z>d`$, $`\rho =\rho _0`$, $`U=U_0z/d`$
layer 3: $`z<d`$, $`\rho =\rho _0(1+ϵ)`$, $`U=U_0`$.

It is easily verified that the Richardson number for this problem is given by

$$\mathrm{Ri}=\frac{gϵd}{U_0^2}.$$
(8)
The solutions of eqn (7) in each of the three regions are

layer 1: $`\varphi _1=A_1e^{kz}`$
layer 2: $`\varphi _2=A_2e^{kz}+B_2e^{kz}`$
layer 3: $`\varphi _3=A_3e^{kz}`$.

Continuity of $`\varphi `$ at the interfaces, and requiring $`\varphi `$ to vanish at $`\pm \mathrm{}`$, gives the characteristic equation:

$$\begin{array}{c}e^{2\stackrel{~}{k}}=\left[1\frac{\stackrel{~}{k}(\stackrel{~}{c}+1)^2}{\mathrm{Ri}+(\stackrel{~}{c}+1)+\frac{1}{2}ϵ\stackrel{~}{k}(\stackrel{~}{c}+1)^2}\right]\hfill \\ \hfill \times \left[1\frac{\stackrel{~}{k}(\stackrel{~}{c}1)^2}{\mathrm{Ri}(\stackrel{~}{c}1)\frac{1}{2}ϵ\stackrel{~}{k}(\stackrel{~}{c}1)^2}\right],\end{array}$$
(9)
where $`\stackrel{~}{k}=2kd`$ and $`\stackrel{~}{c}=c/U_0`$. Eqn (9) can now be solved numerically for $`c`$ for a range of wave and Richardson numbers. When $`c`$ becomes imaginary the perturbation becomes unstable. The parameter region in which this occurs, i.e. the region of instability, is shown in Fig. 1.

The efficiency of mixing is sometimes parametrized by a diffusion coefficient which can be loosely defined as

$$D(\mathrm{Ri})\frac{\mathrm{Im}[c_{\mathrm{max}}(\mathrm{Ri})]}{k_{\mathrm{max}}(\mathrm{Ri})},$$
(10)
where $`\mathrm{Im}[c_{\mathrm{max}}]`$ is the maximal growth rate of the perturbations and $`k_{\mathrm{max}}`$ the corresponding wave number - both evaluated at a given Richardson number. The diffusion coefficient as a function of the Richardson number is shown in Fig. 2. As expected, $`D`$ rapidly decreases with increasing Richardson number.

In this section we presented an example of just one of the many instability analyses found in the literature. The example was chosen for its simplicity. Of course, more realistic velocity and density distributions have been studied since (see e.g. Chandrasekhar 1961, Howard & Maslowe 1973) but qualitatively they all yield the same picture as shown in Fig. 1. Linear stability analyses, however, are only valid for small perturbations. They can be useful for studying the onset of instability, but they are not very reliable in describing the form or the efficiency of mixing itself. This is a problem that requires numerical simulations and a set of such simulations will be presented below.

In the next subsection we will quote a thermodynamic estimate of the diffusion coefficient, before we move on to the numerical simulations in Sec. 3.

### 2.2 Maeder’s prescription

Maeder (1995, 1997) reexamined the Richardson criterion taking into account radiative losses (also see Maeder & Zahn 1998). These authors derive a diffusion coefficient given by

$$D=K\frac{(dU/dz)^2}{N^2},$$
(11)
where $`K`$ is the thermal diffusivity $`K=4acT^3/3\kappa \rho ^2C_P`$, $`a`$ being the radiation density constant, $`c`$ the speed of light, $`C_P`$ the specific heat at constant pressure and $`\kappa `$ the opacity. Maeder’s prescription for mixing due to differential rotation has been used, e.g., by Langer (1992) and Denissenkov (1994) for evolutionary calculations for O and B stars in our Galaxy and the Large Magellanic Cloud. More recently, Denissenkov & Tout (2000) invoked a diffusion constant of the form of Eq. (11) to explain extra mixing in globular cluster red giants.

## 3 Numerical simulations

Here we present the results of numerical simulations of the dynamical shear instability in a stratified fluid. The simulations were obtained using the ZEUS-3D code which was developed especially for problems in astrophysical hydrodynamics (Clarke & Norman 1994). The code uses finite differencing on a Eulerian or pseudo-Lagrangian grid and is fully explicit in time. It is based on an operator-split scheme with piecewise linear functions for the fundamental variables. The fluid is advected through a mesh using the upwind, monotonic interpolation scheme of van Leer. For a detailed description of the algorithms and their numerical implementation see Stone & Norman (1992a, b).

In our simulations we employed an ideal gas equation of state, we ignored the effects of magnetic fields, rotation, nuclear reactions and variations in radiative processes. The simulations were computed on a Cartesian grid and the computational domain was chosen to have the dimension $`210^8`$ cm $`\times `$ $`210^8`$ cm $`\times `$ $`10^8`$ cm (in the $`x`$-, $`y`$\- and $`z`$-direction, where gravity acts in the $`z`$-direction). It was covered by 100 $`\times `$ 100 $`\times `$ 50 grid points. The calculations were performed on a CRAY Jedi parallel-processor and an IBM RS/6000 cluster.

The gravitational acceleration was taken to be that of a point mass of 1 $`M_{}`$ at a distance of 1 $`R_{}`$ from the lower boundary, where the gravitational acceleration acted purely in the $`z`$-direction. An analytic model for an isothermal density distribution under constant gravitational acceleration was relaxed in 1D until hydrostatic equilibrium was attained. The residual velocities after the relaxation were six orders of magnitude smaller than the maximum shear velocities. The initial density distribution is shown in Fig. 3. Then a shear velocity profile was imposed on the fluid. It was chosen to have the form of a hyperbolic tangent in order to minimise the effect of the boundaries onto the shear layer, i.e.

$$U(z)=U_0\mathrm{tanh}[(zz_0)/h],$$
(12)
where $`U_0`$ is the amplitude of the shear velocity, $`z`$ the vertical position of the shear layer, and $`h`$ its extent. In order to keep the shear layer away from the boundaries, $`h`$ was taken to be smaller than the vertical extent of the simulation region. Finally, $`U_0`$ was chosen to yield a range of initial Richardson numbers of 0.05 - 3, where the Richardson number is taken in its original simple definition, i.e. Ri$`=g\rho ^{}/\rho U_{}^{}{}_{}{}^{2}`$, and is measured at $`z=z_0`$. Fig. 3 shows the initial conditions for Ri = 0.1, 0.4, 0.6 and 1.6. The boundary conditions were chosen to be periodic in the $`x`$ and $`y`$ direction and reflecting in the $`z`$ direction.

Figs. 4 and 5 show the $`x`$\- and $`z`$-velocities in a vertical slice through the computational domain at various times, for a flow with an initial Richardson number of 0.4 in the shear layer. The occurrence of turbulence is clearly seen. One can see that early on in the simulation the flow undulates slightly before the instability becomes fully nonlinear. Then one can observe that two vortices form which mix material over a large section of the computational domain. The vortices remain quite stable and increase in size. Their centres coincide with the centre of the shear layer. One can also note that at the centre of the vortex, the flow velocity is very small. For efficient transport the velocity fluctuations in the horizontal and the vertical have to be correlated. This correlation is apparent in the formation of vortices as seen in Fig. 5. More quantitatively it can be expressed through the vorticity. Maps of the vorticity in the direction both perpendicular to the gravitational acceleration and the initial flow are shown in Fig. 6. (for Ri=0.4). The contours in Fig. 6a. which show the vorticity after $`t=1500`$ s span a range of $`310^3`$ s-1 compared to a range of $`10^2`$ s-1 in Fig. 6b. after 9000 s. This example shows the production of vorticity by the shear instability.

The eddies have the form of long “rolls” which are aligned with the $`y`$-axis, i.e. perpendicular to, both, the shear flow and the gravitational acceleration. They are found to stay fairly symmetric in the $`y`$-direction through the entire depth of the computational box. The structures on the big scales, as for example the “rolls” in our simulations, ‘feel’ the symmetry of the background model in the $`y`$-direction and therefore, they retain this symmetry. This has also been observed in laboratory experiments with incompressible fluids. Only on the small scales homogeneous, isotropic 3D turbulence will develop but we do not resolve these scales in our simulations. In the case of the shear instability, mixing is provided primarily by motions on big scales and, consequently, we argue that for our purpose we do not have to resolve these small scales. A further discussion of this issue can be found in Sec. 3.2.

Furthermore, over longer distances the “rolls” will start to feel the curvature of the star which is not accounted for in the slab geometry of our simulations and they may have the effect of breaking up the symmetry in the $`y`$-direction. These issues remain uncertainties in the results presented here and will have to be examined in a later study.

In order to study mixing processes the ZEUS code was modified to follow the motion of 1000 ‘tracer’ particles which are advected with the fluid. The positions of some of the tracer particles at various times (again for the case with Ri = 0.4) are shown in Fig. 7.

The diffusion constant can then be defined as

$$D=\sigma ^2/t,$$
(13)
where $`\sigma ^2=\frac{1}{N}_N[z(N)z_0]^2`$, $`z_0`$ being the original height of the $`N`$th tracer particle and $`z(N)`$ its height after a time $`t`$. As apparent from Fig. 7 the diffusion coefficient varies with time. The diffusion constant as a function of time have been plotted in Fig. 8. In Fig. 8 the Richardson number is 0.4, which is greater than 0.25 and, therefore, the flow should be stable to the shear instability according to the simple Richardson criterion. The simulations reveal that it is not, as has in fact been shown in previous simulations. One can observe that initially $`D`$ rises with time before it eventually approaches a value which remains nearly constant for some time. The remaining scatter is mainly due to the stochastic nature of the turbulent mixing. The transient phase during which $`D`$ rises is longer, the greater the Richardson number. Eventually, dissipation (mainly numerical) becomes noticeable and $`D`$ starts to slowly decrease again. The constant value of $`D`$ to which the curves are converging, is the value that is of greatest interest for the purpose of constructing stellar models. This value for $`D`$ depends on the Richardson number. An understanding of the quantitative dependence of $`D`$ on parameters such as the Richardson number will be useful for constructing stellar models which treat mixing through a diffusion equation.

We have plotted this ‘stationary’ value of $`D`$ as a function of Ri for all the simulations which we have performed (see Fig. 10). The errorbars indicate the residual scatter observed in the simulations. Configurations with Ri $`>`$ 1.6 proved to be stable for times up to 30,000 s.

### 3.1 2D versus 3D simulations

Two-dimensional simulations of convective instabilities are known to introduce numerical artefacts and to yield unrealistic results. In 2D simulations of convection the energy cascades to larger scales where viscous dissipation is less effective, whereas 3D simulations reveal that the energy cascades to smaller scales where it is absorbed. Moreover, 2D simulations often produce ‘rolls’, i.e. vortices in the direction perpendicular to the simulated domain, which, for convection, tend to overestimate the mixing.

In order to investigate the importance of performing these simulations in three dimensions, the runs described above were repeated in two dimensions. The ZEUS code is designed to efficiently compute on 1- and 2-dimensional grids, so the runs could be readily repeated in 2D. The same initial model and the same set of physical parameters as in the 3D simulations were used, and a grid of 100 $`\times `$ 50 points was chosen. Fig. 9 shows a sequence of snapshots of the velocity.

The diffusion coefficient derived from the 2D simulations as a function of time is displayed as a dashed line in Fig. 8. The figures show that in 2D simulations $`D`$ rises more slowly with time than in three dimensions. But eventually $`D`$ converges to about the same value which was obtained in the 3D simulations. This is an encouraging result. Not only does this raise the trust in the credibility and self-consistency of these results, it also suggests that 2D simulations may be sufficient to determine the efficiency of shear mixing. This conclusion may not be too surprising given that the 3D simulations revealed a flow structure consisting of long rolls that are parallel to the $`y`$-axis, i.e. perpendicular to the shear flow. On smaller scales, near the Kolmogoroff scale, the picture would look very different because one would observe the onset of full isotropic 3D turbulence. But here we are only simulating the Kelvin-Helmholtz instability and not full 3D turbulence which is a much more difficult problem and which requires a much finer (and yet unattained) numerical resolution. This will be discussed further in the next section.

### 3.2 Numerical viscosity

One frequently voiced objection to these kinds of direct numerical simulations is that Reynolds numbers as high as those encountered in stars are unattainable on current computers, and that therefore the results are unrealistic. But, as pointed out by Balbus, Hawley & Stone (1996), this criticism is unjustified for simulations of the shear instability. They argue, that in order to simulate the onset of instability in a laminar flow, it is merely necessary that the ‘typical’ wavelength of the instability is resolved by the numerical scheme and that the numerical diffusion at this wavelength is less than the growth rate. This makes the simulation of shear instabilities an easier task than the simulation of viscous instabilities where in theory one would have to resolve everything down to the viscous length scale.

If the flow is unbounded, there is no viscous boundary layer which might interfere with the results. The nonlinear instabilities are fundamentally inviscid in character. Therefore, we only need a resolution capable of resolving a range of wavelengths with numerical diffusion errors less than the growth rates. This view is also supported by experimental observations which suggest that mixing through shear instabilities is dominated by large scale motions.

In Fig. 12 we have plotted the kinetic energy in the vertical ($`z`$-direction) as a function of time. The vertical kinetic energy is expressed in units of the total initial kinetic energy. In the beginning the energy is very small until the instability sets in. Then the kinetic energy jumps to about 10% of its total initial value and slowly decreases amidst big fluctuations. This decrease is mainly due to numerical viscosity and shows yet again that we are not simulating proper 3D turbulence which would show a saturation of the kinetic energy.

Porter & Woodward (1994) have estimated the Reynolds number of hydrodynamical simulations based on a PPM (piecewise parabolic method) code. The Reynolds number depends on the truncation error of the finite-difference algorithm, the Courant number, the background advection and the number of grid points. They found that the effective Reynolds number is proportional to the third power of the number of grid points, with the main dissipation occurring at short wavelengths. Above this critical wavelength the diffusion was found to be small. Since the ZEUS code uses piecewise linear functions instead of piecewise parabolic functions, the truncation errors of ZEUS will be larger than those of a PPM code. But even if they were only proportional to the second power of the number of grid points in one direction, we would still expect a numerical Reynolds number of around $`10^4`$.

The effect of the numerical viscosity is to suppress fluctuations beneath the ‘effective’ viscous length scale. This will become important for high enough Richardson numbers and all numerical simulations will find stability above a certain Richardson number. In our case the flow remained stable for Ri $`>`$1.6. Higher resolutions, and therefore higher Reynolds numbers will push this limiting Richardson number to slightly higher values.

### 3.3 Dependence on grid size

We have investigated the dependence of our results on the resolution of the computational grid by repeating the 2D runs presented above (100 $`\times `$ 50) on grids with (b) 200 $`\times `$ 100 and (c) 400 $`\times `$ 200 grid points. The variation of $`D`$ with time is shown in Fig. 8. by the dotted (b) and dash-dotted (c) lines. It is immediately evident that the behaviour of $`D`$ with time varies with the numerical resolution. For more finely resolved grids, $`D`$ rises more quickly. This was to be expected since higher resolutions produce higher Reynolds numbers and therefore a less viscous flow. For higher effective Reynolds numbers the turbulence grows more rapidly, and this is observed here. But again it merits mention that for all simulations $`D`$ eventually converges to rougly the same value for a given Richardson number, as long as the Richardson number is not too big. For large Richardson numbers, the numerical viscosity becomes important, but here we will restrict the discussion to Richardson numbers less than 2.

Finally, it can be seen in Fig. 5 that the eddies nearly span the entire vertical extent of the simulation domain. Hence, it is interesting to verify whether the vertical size of the computational domain affects the size of the eddies, and therewith the value of the mixing coefficients. We have repeated some 2D simulations in computational domains which had a vertical extent of more than $`50h`$ or ten times the vertical size of the previously shown simulations. We did not find that the eddies became bigger when the size of the computational box was increased.

## 4 Discussion

In this paper results of direct numerical simulations of shear instabilities were presented. The simulations were carried out using the ZEUS hydro code. A chemically homogeneous stratified fluid in hydrostatic equilibrium was set up, onto which a shear flow was superimposed. The mixing was quantified by studying the dispersion of ‘tracer’ particles that are advected by the fluid. Thus, the variation of the diffusion constant with the Richardson number was examined.

In accordance with previous simulations, it was found that efficient mixing occurs for Richardson numbers substantially higher than 0.25, contrary to the simple linear stability criterion. Now the diffusion coefficients that were derived from the motion of the tracer particles can be compared to the analytical estimates discussed in Sec. 2. Applying Maeder’s formalism to our initial model, i.e. Eqn (11), yields substantially lower values for $`D`$ than found in our simulations. For our initial models, Maeder’s diffusion coefficients in the centre of the shear layer lie between values of $`10^6`$ \- $`10^7`$ cm$`{}_{}{}^{2}/`$s. This is more than 3 orders of magnitude smaller than the diffusion coefficients found here. Finally, one might want to compare our results with Fig. 2, even though this figure assumed a different shear velocity distribution than the hyperbolic tangens used in the simulations. In Fig. 11 where the dimensionless diffusion coefficient (in units of $`U_0h`$) is plotted against the Richardson number. This dimensionless presentation is also useful for a wider application of our results. It is interesting to note that now our values for $`D`$ are smaller by almost two orders of magnitude than those predicted by the linear analysis of our simple model of Sec. 2.

It was argued that numerical viscosity has no great effect on our findings as the effective Reynolds number is small on the scales on which the Kelvin-Helmholtz instability operates. This is true as long as the Richardson number is not too big. Furthermore, the differences between 2D and 3D simulations of shear instabilities were studied. We found that the dimensionality of the simulation hardly affects the final ‘steady state’ diffusion coefficient. It only affects the dynamics of the onset of the instability until a quasi-stationary state has been reached. The same qualitative observation was made when the resolution of the numerical grid was varied. It was found that the resolution affects the initial temporal variation of the diffusion coefficient, but hardly affects the final value of $`D`$. However, the similarity between the 2D and 3D simulations is only valid on the big scales which we are resolving. On smaller scales, near the Kolmogoroff scale, one would observe full isotropic 3D turbulence which looks very different in two dimensions.

To make the Richardson number the only controlled parameter in our study, is clearly a simplification of the factors that determine the efficiency of mixing. In reality, the efficiency of mixing will depend on the density stratification and the velocity gradient separately, and not solely on the Richardson number. Moreover, the mixing will depend on the exact shape of the velocity profile and, only to a first approximation, on its first derivative. Nevertheless, if one is interested in the efficiency of mixing as a function of a single parameter, e.g. for stellar evolution studies, the Richardson number is still the most suitable parameter, since in the limit of an inviscid and incompressible shear flow, the stability is governed by the Richardson number.

In this work the effects of the dynamical shear instability in a plane shear layer were investigated. Clearly, this is not the only mechanism responsible for mixing in stars. In rotating fluid bodies other instabilities may occur, such as, e.g., the Solberg-Holand instability, the Eddington-Sweet instability and the Goldreich-Schubert-Fricke instability. However, the shear instability operates on the dynamical timescale whereas the instabilities tied to the rotation occur on the longer Eddington-Vogt timescale. Therefore, these instabilities operate on a much longer timescale and therefore play a lesser role in diffusion processes, even though they can affect the long-term evolution of the star.

The consequences of the newly found diffusion coefficients onto stellar evolution and surface abundances are difficult to foresee. In stars the speed as well as the depth of mixing determine the balance between mixing and nuclear burning. Therefore, the effect of mixing depends sensitively on the conditions prevailing in the star, as for example the positions and extents of its convection and burning shells.

Finally we should mention that a number of factors can inhibit or facilitate mixing such as gradients in the chemical potential of the fluid, diffusion of radiation, magnetic fields and effects pertaining to the spherical geometry. In late-type stars strong chemical composition gradients exist, which will have a stabilising effect on the stratification. Therefore, especially for stellar evolution studies, the chemical composition gradient is an important parameter, which will have to be included in future work. Similarly, poloidal magnetic fields, as they are believed to exist in the Sun, will also suppress mixing processes. These effects are potentially of great importance in stars and will be studied in a forthcoming paper.

## Acknowledgements

A part of the simulations were performed on computers of the Rechenzentrum Garching. We thank Phil Armitage for helpful discussions and the referee for useful comments.