Amine Boumaza
email: boumaza@loria.fr

Bruno Scherrer
email: scherrer@loria.fr

Convergence and Rate of Convergence of a Simple Ant Model

Keywords: Algorithms, Theory Multiagent Systems::Emergent behavior, Multiagent Systems:: Multiagent planning, Agents::Formal models of agency

We present a simple ant model that solves a discrete foraging problem. We provide simulations and a convergence analysis. We argue that the ant population computes the solutions of some optimal control problems and converges in some well defined sense. We also discuss the rate of convergence with respect to the number of ants: we give experimental and theoretical arguments that suggest that this rate is superlinear with respect to the number of agents.

A SIMPLE ANT MODEL

Consider a set of artifical ants that move on a 2D-grid (the environment) on which they update artificial "pheronomone traces". Each cell s of this grid stores two pheromone traces as two real numbers: V1(s) and V2(s). We consider that there are exactly four types of grid cells: one cell is the nest, one cell is a food source, there are several bad cells (or obstacles) and all the others are free cells.

Each ant can be in two possible states: either it is carrying food (state "has food") or it is carrying nothing (state "has nothing"). Ants move around and their state change according to the following natural rules: when an ant gets to the food source, its state is set to "has food"; When an ant gets to the nest, its state is set to "has nothing" ; if before being set to "has nothing", the ant that got to the nest had food (i.e. it has just brought food to the nest) then we increment a food counter. This food counter will measure how much food is foraged by the whole ant population.

We now describe the dynamics of this ant model. At the AAMAS'07 May 14-18 2007, Honolulu, Hawai'i, USA.

beginning: a) the food counter is set to 0, b) the positions of the ants are initialized arbitrarily, c) all ants are set to the "carry nothing" state and d) the pheromone values are initialized arbitrarily. At each time step, each ant does two things:

• It updates the two local pheromone traces V1(s) and V2(s) of its current cell s using the pheromone values of its four neighbours. Let us write N (s) the set of neighbours. For i ∈ {1, 2}, let us note max i (N (s)) max s ∈N (s) Vi(s ) and avg i (N (s))

1 4
P s ∈N (s) Vi(s ). The pheromones updates go as follows:

V1(s) ← 8 > > < > > : -1 if s is a bad cell 1 if s is the food source cell β `αmax 1 (N (s)) + (1 -α)avg 1 (N (s)) ótherwise V2(s) ← 8 > > < > > : -1 if s is a bad cell 1 if s is the nest cell β `αmax 2 (N (s)) + (1 -α)avg 2 (N (s)) ótherwise
where 0 ≤ α ≤ 1 and 0 ≤ β < 1 with the condition that β < 1 if α = 1.

• It moves to one of its neighboring cell: with probability (0 ≤ ≤ 1) (that we shall call the exploration rate) it moves uniformly at random to one of its neighbors. With probability 1 -, it moves to the neighbors that has the highest "pheromone" value V1 or V2 depending on its state: the ant uses the pheromone values of V1 if it is in state "carry nothing" and the values of V2 if it is in state "carry food".

The β parameter can be seen as some sort of evaporation parameter and is typically set close to 1. The α parameter should typically be set either close to 0 or close to 1. Two simple peculiar instances of our model correspond to the parameter choices (α = 1, 0 < β < 1) and (α = 0, β = 1). In the former choice, the general pheromone equation (the "otherwise" case above) reduces to Vi(s) ← βmax i (N (s)) and in the latter choice to: Vi(s) ← avg i (N (s)) which is a simple linear update. The question whether the ant activities (pheromone values updates and moves) are done synchronously or in some sort of asynchronous way is unimportant.

SIMULATIONS

Let us consider a typical run of the algorithm during which we fixed the parameters as follow: = 0.8, α = 0.7, β = 0.9999 and m the population size is set to 150. Figure 1 shows the evolution of the ants starting at the nest and evolving through the environment searching for food. Once the food source is found by an ant, more ants move toward the source and a trail emerges between the nest and the source. As time goes, a good proportion of the ants follow the trail and the dynamics seems to stabilize. Depending on the value of α, we can observe different forms of paths, or even sometimes no path. Figure 2 illustrates the asymptotic distribution of ants for different values of α: there are paths except for the case α = 0.6 (see [START_REF] Boumaza | Convergence and rate of convergence of a simple ant algorithm[END_REF] for an explanation).

(t = 10000) (t = 11000) (t = 12000) (t = 16000)
Suppose a path emerges between the nest and the food source. It would be interesting to measure something objective that reveals this path. This can be done through the food counter we introduced earlier. We can draw the curve showing the increase of food brought back to the nest over time. A typical such curve is shown in figure 3. After some time, we observe that this food quantity increases linearly: this means that the foraging behavior of the ants has somehow converged.

The food quantity curve experimentally shows that there is convergence. It can also be exploited to measure a convergence rate. To do so, we fit a line to the linear part of the curve (see figure 3) and define the time of convergence as the intersection of this line with the the x-axis. The rate of convergence is then computed as the inverse of the time of convergence. One can compute some statistics on the rate of convergence to evaluate its dependence with respect to the number of ants. Taking α = 0.7, β = 0.9999, = 0.8, we measured the rate of convergence for various population sizes m. imentally that when the number of ants increases linearly, the increase in rate of convergence is superlinear. In other words, if the number of ants is doubled, the rate increases by more than a factor two.

ANALYSIS

The ant algorithm we have just described can be analyzed within the framework of Markov Decision Processes (MDP), that aims to solve optimal control problems [START_REF] Sutton | Reinforcement Learning : An Introduction[END_REF]. Indeed, it can be shown that the pheromone values V1 and V2 each converge to the optimal value function (i.e. the solution) of a control problem: Proposition 1. Consider the ant model described in section 1. If the exploration rate > 0, then the pheromone value V1 (resp. V2) asymptotically converges to the optimal value function of the MDP M1 = S, A, T1, R1 (resp. M2 = S, A, T1, R1 ) with discount factor β where: • The state space S is the set of grid cells plus an extra "end state".

• The action space A is the set of cardinal moves (north, east, south, west).

• The transition T1 (resp. T2) is characterized as follows: 1) when, in a state s corresponding to a free cell or the nest cell (resp. to a free cell or the food cell), one chooses one of the four directions a ∈ A, the probability of actually making the move in the direction a is α

+ 1-α 4 = 3α+1 4
while the probability of getting to each of the 3 other neighbors is 1-α 4 .

2) From all the other states, that is from bad cells, the food and the "end state" (resp. from bad cells, the nest and the "end state"), there is, for every action, a probability 1 to get to the "end state", which is an absorbing state.

• The reward R1 (resp. R2) is defined as follows: For all states s corresponding to a free cell or a nest cell (resp. a free cell or a food cell), the reward is 0. For the state corresponding to the food (resp. the nest), the reward is 1. The reward is -1 for all states corresponding to bad cells and 0 for the "end state".

The proof of the above result simply consists in checking that the pheromone updates of our algorithm are equivalent to some asynchronous version of the value iteration algorithm, a standard algorithm for solving optimal control problems [START_REF] Sutton | Reinforcement Learning : An Introduction[END_REF], for the problems M1 and M2 (see [START_REF] Boumaza | Convergence and rate of convergence of a simple ant algorithm[END_REF] for more details). Asymptotically, the moves of the ants, which are a mixture of random uniform moves and the action that climbs up the pheromone value (cf. section 1), correspond to a mixture of random uniform moves and the optimal corresponding moves. Depending on whether an ant is carrying food or not, it uses one pheromone or the other, i.e. it moves toward the nest or the food source. This makes the ants go back and forth between the food source and the nest, and explains the emergence of paths.

It is also possible to study the rate of convergence with respect to the number of ants. Seeing the environment as a graph, and the ants as parallel random walks on it, it can be shown [START_REF] Boumaza | Convergence and rate of convergence of a simple ant algorithm[END_REF] that the convergence time is related to the time necessary for the walks to visit every node of the graph, which is known as the cover time of the graph by these parallel random walks [START_REF] Aldous | Reversible Markov Chains and Random Walks on Graphs[END_REF]. Though difficult to compute in general, an asymptotic result on cover time [START_REF] Aldous | Reversible Markov Chains and Random Walks on Graphs[END_REF] allows to state a proposition that sheds light on our observation of superlinear rate of convergence: Proposition 2. Consider the ant model described in section 1 on a toric grid environment with n cells, with an exploration rate = 1 and initialize the ants uniformly on the environment. When the size of the environment n → ∞ and while the number m of ants satisfies m ≥ 6 log n, the time for the pheromone values to reduce their distance to their limit by a factor β is bounded by (25+o(1))n 2 log 2 n m 2 that is a function that is inversely quadratic in the number m of ants.

It is known [START_REF] Aldous | Reversible Markov Chains and Random Walks on Graphs[END_REF] that such a bound is sharp when the environment is a n-cycle graph, that is a long (cycled) corridor: in such a case, the rate of convergence of the pheromone values to their limits is quadratic in the number of ants. Though this result is limited to the case = 1 (where the ants never exploit the pheromone values), we believe that the above bound remains close when is close to 1, and this supports our empirical observation of superlinearity.

Our work is close to previous recent works. In [START_REF] Panait | A pheromone-based utility model for collaborative foraging[END_REF] the authors tackle foraging problem with multiple pheromone landscapes by exploiting the optimal control theory (in fact they refer to reinforcement learning [START_REF] Sutton | Reinforcement Learning : An Introduction[END_REF] which is also formulated through MDPs). In [START_REF] Simonin | Construction of numerical potential fields with reactive agents[END_REF], the author presents an algorithm that compute shortest path using an asynchronous implementation of the Bellman-Ford algorithm: the local update is of the form U (x) ← 1 + min x ∈N (x) U (x ). Modulo the variable change U ↔ log(V ) log(β) this is equivalent to our pheromone update with α = 1. The main contribution of our work is to deepen the analytical view on these algorithms: if a convergence analysis of [START_REF] Simonin | Construction of numerical potential fields with reactive agents[END_REF] can be done [START_REF] Simonin | [END_REF], [START_REF] Panait | A pheromone-based utility model for collaborative foraging[END_REF] only provides experimental evidence. Furthermore, to our knowledge, our work is the only one that tries to tackle the rate of convergence analytically. More details and discussion can be found in an extended version of this article [START_REF] Boumaza | Convergence and rate of convergence of a simple ant algorithm[END_REF].

Figure 1 :

 1 Figure 1: A typical dynamics of the ant model.

Figure 4 Figure 2 :Figure 3 :

 423 Figure 2: Limit distribution of the ants for various values of α.

Figure 4 :

 4 Figure 4: Convergence rate with respect to the number of ants.


Eun-Jin Kim
email: e.kim@sheffield.ac.uk

Nicolas Leprovost

Self-consistent theory of turbulent transport in the solar tachocline III. Gravity waves

Keywords: Turbulence -Sun, interior -Sun, rotation -waves

Aims. To understand the fundamental physical processes important for the evolution of solar rotation and distribution of chemical species, we provide theoretical predictions for particle mixing and momentum transport in the stably stratified tachocline.

Methods. By envisioning that turbulence is driven in the tachocline, we compute the amplitude of turbulent flow, turbulent particle diffusivities, and eddy viscosity, by incorporating the effect of a strong radial differential rotation and stable stratification. We identify the different roles that the shear flow and stable stratification play in turbulence regulation and transport.

Results. Particle transport is found to be severely quenched due to stable stratification as well as radial differential rotation, especially in the radial direction with an effectively more efficient horizontal transport.

The eddy viscosity is shown to become negative for parameter values typical of the tachocline, suggesting that turbulence in the stably stratified tachocline leads to a non-uniform radial differential rotation. Similar results also hold in the radiative interiors of stars, in general.

Introduction

Since the formation of the radiative core, which marked the beginning of its journey on the main sequence, the sun has slowed down significantly due to the loss of angular momentum from its surface (e.g. see [START_REF] Stix | The sun: an introduction[END_REF][START_REF] Schatzman | The stars[END_REF]. The angular momentum transport must have been very efficient during its spin-down in order for the sun to have rotational profile as

Send offprint requests to: E. Kim observed today (see, e.g. Charbonneau et al 1998). Vigorous turbulence in the convection zone and possibly thermal wind (Miesch, Brun & Toomre 2006) can readily provide a mechanism for efficient radial momentum transport, thereby eradicating radial differential rotation therein. Such turbulent transport is however considered to be absent in the stably stratified radiative interior, which has also spun-down during the solar evolution, presently rotating uniformly at the rate roughly the same as the mean average rotation rate on the solar surface. Whichever mechanism is responsible for momentum transport in the interior (which itself is an important problem), it should be closely related to the transport in the tachocline through which the surface spindown is communicated to the interior. Transport in the tachocline also plays a crucial role in the overall mixing of light elements (lithium, beryllium, etc) (see e.g. [START_REF] Schatzman | The stars[END_REF]Brun, Turck-Chiéze, & Zahn 1999), thereby determining the level of their surface abundances on the sun. Therefore, it is essential to understand physical mechanisms for transport in the tachocline and then to formulate a consistent theory starting from first principles based on those processes. This is particularly true since virtually all the previous theoretical modelling heavily relies on a simple parameterization of transport process, which is then adjusted to obtain the agreement with observations.

We have initiated the development of a consistent theory of turbulent transport in the tachocline in the previous papers, by taking into account the crucial effect of the large-scale shear flows, provided by a strong radial differential rotation (Kim 2005) as well as latitudinal differential rotation (Leprovost & Kim 2006). By envisioning that the tachocline is perturbed externally [e.g. by plumes penetrating from the convection zone above (e.g. see Gilman 2000;Brummell, Clune & Toomre 2002;Rogers & Glazmaier 2005)], or by instabilities (e.g. see Watson 1981;Charbonneau, Dikpati & Gilman 1999), we demonstrated how turbulence level and transport are reduced via shearing in a non-trivial manner (see also Burrell 1997;Kim & Diamond 2003;Kim 2004;2006) in a simplified three dimensional (3D) hydrodynamic turbulence. In particular, turbulent transport of chemical species and angular momentum are shown to become strongly anisotropic with effectively much more efficient transport in the horizontal (latitudinal) direction than the vertical (radial) direction due to shear stabilization by strong radial shear. The resulting anisotropic momentum transport was shown to reinforce a strong radial shear (Leprovost & Kim 2006), with a positive feedback on the confinement of the tachocline (Spiegel & Zahn 1992) while chemical species are predicted to have latitudinal dependent mixing due to the variation of radial shear (Kim 2005). Furthermore, the results indicate that the turbulence regulation by a shear flow (i.e. differential rotation) leads to weak turbulence and mixing in turbulent tachocline.

The purpose of this paper is to investigate how a stable stratification in the tachocline modifies the predictions obtained in these studies. We again envision that the tachocline is turbulent, driven by a forcing as in Kim (2005). In the presence of a stable stratification, the turbulence in the tachocline is no longer completely random as the stable stratification provides a restoring force against radial displacement of fluid elements, supporting the propagation of internal gravity waves [START_REF] Lighthill | Waves in Fluids[END_REF]. These waves tend to increase the memory of otherwise random turbulent fluid motion and can reduce the overall transport due to turbulence. We shall show that the turbulent transport due to shear stabilization found in the previous studies is further enhanced in the presence of stable stratification (gravity waves). Thus, gravity waves acting together with shear stabilization can lead to a weak mixing in the tachocline, as required for the surface depletion of lithiums (e.g. Pinsonneault et al. 1989).

It is important to contrast our approach to those adopted in most previous works, which focus on the momentum transport by gravity waves themselves through dissipative processes (e.g. Plumb 1978;Kim and MacGregor 2001;2003;[START_REF] Talon | [END_REF]. For instance, gravity waves are considered to be generated in the convection zone (Press 1981) and deposit their momentum into background shear flow via radiative damping as they propagate through the tachocline and the interior (Kim and MacGregor 2001;2003). Relying crucially on radiative damping, the momentum transport by gravity waves would occur on a long time scale. [Recall that the transport by waves requires molecular dissipation, and is thus a slow process.] Similarly, weak mixing due to damped gravity waves was also suggested as a mechanism for a modestly enhanced mixing of light elements (lithium) (Press & Rybicki 1981;Garc ía López & Spruit 1991). In these works, a certain level of background turbulence could not be ruled out and was invoked to provide an enhanced viscosity for the evolution of a mean flow (Kim and MacGregor 2001;2003) and the gravity waves (Charbonnel & Talon 2005;[START_REF] Talon | [END_REF]. The value of the effective viscosity is often arbitrarily taken to be a positive constant, being much larger than molecular viscosity, not affected by shear flow nor by gravity waves, or it is simply parameterized.

In this paper, we shall treat turbulence and gravity waves on an equal footing and consistently compute the values of turbulent eddy viscosity and particle diffusivity, by incorporating the effect of a shear flow (provided by a radial differential rotation), and identify different roles that gravity waves and shear flow play in turbulent transport. Specifically, we shall show that unlike shear flows, which reduce both turbulence level and transport, a stable stratification can suppresses turbulent transport without much effect on turbulence level. We shall further demonstrate that the stratification favors a negative eddy viscosity and thus tends to sharpen the radial gradient of large-scale shear flows rather than smoothing it out. This tendency was also found in Kim & MacGregor (2001;2003). We note that the elucidation of the effects of gravity waves is essential for understanding the momentum transport in the radiative interior as they have often been advocated as a mechanism to explain a uniform rotation in that region. Stratified turbulence with a shear flow is also important for the transport in radiative interiors and/or envelops of stars, in general, as well as in geophysical systems. In particular, it has actively been studied in geophysical systems (e.g. see Jacobitz, Sarkar & van Atta 1997;Stacey, Monismith, & Burau 1999), where stable stratification was shown to inhibit turbulent transport in the direction of a background density gradient, leading to the two dimensional (2D) turbulence.

The remainder of the paper is structured as follows. We first investigate the effect of gravity waves on turbulence in Sect. 2. In Sect. 3, we incorporate the effect of strong radial differen-tial rotation and study how the gravity waves modify the overall turbulent transport. Section 4 contains the conclusion and discussions. To simplify the problem, we shall consider incompressible fluid with local cartesian coordinates x, y, and z for radial, azimuthal, and latitudinal directions, respectively (see Fig. 1), and use Boussinesq approximation to capture the effect of stratification. We assume that the fluid is stirred by a forcing on small scales, giving rise to fluctuations. In the absence of stratification, this fluid forcing will drive turbulence on small scales and maintain it at the level at which the injected energy is balanced by dissipation in the system. As the stratification increases, the forcing will generate not only random turbulent motion but also coherent (gravity) waves. Alternatively, some of turbulent motion will be turned into packets of gravity waves. Therefore, we consider both turbulence and gravity waves as fluctuations on scales much smaller than those associated with mean density or mean background shear flows. Specifically, we express the total mass density ρ = ρ 0 (x) + ρ ′ where ρ 0 (x) and ρ ′ are the background and fluctuating mass densities, respectively; the total velocity u = U 0 + v where U 0 and v are a large-scale shear flow due to differential rotation and small-scale fluctuations; the total particle density of chemical elements n = n 0 (x) + n ′ where n 0 (x) and n ′ are mean and fluctuating components. Then, the main governing equations for fluctuations v, ρ ′ , and n ′ , involving both turbulence and gravity waves, are as follows (see, e.g., Kim & MacGregor 2003;[START_REF] Moffatt | Magnetic field generation in electrically conducting fluids[END_REF]:

Internal Gravity Waves

(∂ t + U 0 • ∇)v = -∇p -gρ ′ x + ν∇ 2 v + f , (1)
∇ • v = 0 , (2)
(∂ t + U 0 • ∇)ρ ′ = ρN 2 g v x + µ∇ 2 ρ ′ , (3)
(∂ t + U 0 • ∇)n ′ = -∂ x n 0 v x + D∇ 2 n ′ . (4)
Here, ν, µ and D are molecular viscosity, thermal diffusivity, and particle diffusivity, respectively; f in Eq. ( 1) is the small-scale forcing driving turbulence; ρ = ρ 0 (x = 0) is the constant background density [measured at the bottom of the convection zone (e.g. see Kim and MacGregor 2003)],

and N = -g(∂ x ρ 0 + ρg/c 2 s )/ρ is the Brunt-Väisälä frequency, where c s is the sound speed. Note that the typical values of ν, D, µ, and N in the tachocline are 10 2 cm 2 s -1 ,10 2 cm 2 s -1 , 10 7 cm 2 s -1 , and 3 × 10 -3 s -1 , respectively.

In this section, we first ignore a background shear flow and study how the turbulence and transport are modified due to stable stratification in the tachocline. To obtain the overall turbulent transport, we solve coupled equations ( 1)-( 4) in terms of Fourier transform for fluctuating quantities φ ′ :

φ ′ (x, t) = 1 (2π) 3 d 3 k φ(k, t) exp {i(k x x + k y y + k z z)} . (5)
Equations ( 1)-( 3) then give us an equation for ρ = e νk 2 t ρ in the form

∂ tt ρ + (µ -ν)k 2 ρ + N 2 ρ = N 2 g e νk 2 t h1 . (6)
Here, N 2 = γN 2 /(γ+a 2 ). Since N ≫ µk 2 is valid on a broad range of reasonable scales l ≫ 10 4 ∼ 10 5 cm in the tachocline, we can find the solutions to leading order in (νk 2 /N ≪) µk 2 /N ≪ 1 as:

ṽx ∼ 1 γ + a 2 dt ′ G µ (t, t ′ ) cos N(t -t ′ ) h1 (t ′ ) , (7) ṽy
∼ - a γ ṽx - β γ dt ′ G(t, t ′ ) h2 (t ′ ) , (8)
ṽz ∼ -

aβ γ ṽx + 1 γ dt ′ G(t, t ′ ) h2 (t ′ ) , (9) ρ
∼ N g γ(γ + a 2 ) dt ′ G µ (t, t ′ ) sin N(t -t ′ ) h1 (t ′ ) . (10)
Here, G(t, t ′ ) = e -νk 2 (t-t ′ ) and G µ (t, t ′ ) = e -(µ+ν)k 2 (t-t ′ )/2 are the Green's functions; a = k x /k y , β = k z /k y , and γ = 1 + β 2 ; h1 and h2 in Eqs. ( 6) and ( 7)-( 10) are the forcing terms which are related to fi as

h1 = (1 + β 2 ) fx -a fy -aβ fz , (11) h2
= -β fy + fz . (12)
Equations ( 7)-( 9) show that the vertical motion ṽx , involving G µ only, is always subject to radiative damping µ while the horizontal motions ṽy and ṽx , containing both G µ and G, are much less affected by µ.

For simplicity, we assume the forcing to be homogeneous in space with a short correlation

time τ f : fi (k 1 , t 1 ) f j (k 2 , t 2 ) = τ f (2π) 3 δ(t 1 -t 2 )δ(k 1 + k 2 )ψ i j (k 2 ) . ( 13
)
Note that in the case where the forcing f due to plumes induces a stronger turbulence towards the bottom of the convection zone, f becomes inhomogeneous with the power spectrum ψ i j = ψ i j (x, k) depending on the radial coordinate x.

A long but straightforward algebra by using Eqs. ( 7)-( 13) then gives us the following results on turbulence level and particle diffusivities defined by

n ′ v i = -D i j T ∂ j n 0 : v 2 x ∼ τ f d 3 k (2π) 3 φ 11 (k) 2(γ + a 2 ) 2 1 (µ + ν)k 2 , (14) v 2 y ∼ τ f d 3 k (2π) 3 a 2 φ 11 (k) 2γ 2 (γ + a 2 ) 2 1 (µ + ν)k 2 + β 2 φ 22 (k) γ 2 1 2νk 2 , ( 15
)
v 2 z ∼ τ f d 3 k (2π) 3 a 2 β 2 φ 11 (k) 2γ 2 (γ + a 2 ) 2 1 (µ + ν)k 2 + φ 22 (k) γ 2 1 2νk 2 , ( 16
)
D xx T ∼ τ f N 2 d 3 k (2π) 3 φ 11 (k) 4γ(γ + a 2 ) , ( 17
)
D yy T ∼ τ f d 3 k (2π) 3 γ 2         a 2 φ 11 (k) 4N 2 (γ + a 2 ) 2 + β 2 φ 22 (k) 2ν(ν + D) 2 k 4         , (18) D zz T ∼ τ f d 3 k (2π) 3 γ 2         a 2 β 2 φ 11 (k) 4N 2 (γ + a 2 ) 2 + φ 22 (k) 2ν(ν + D) 2 k 4         . (19)
Here, we assumed spatial symmetries φ i j (k y ) = φ i j (-k y ) and φ i j (k y ) = φ i j (-k y ); k inside the integrals in Eqs. ( 14)-( 19) is the wavenumber of the forcing; φ i j (k) (i = 1, 2) is the power spectrum of the forcing defined by: hi

(k 1 , t 1 ) h j (k 2 , t 2 ) = τ f (2π) 3 δ(t 1 -t 2 )δ(k 1 + k 2 )φ i j (k 2 ) . (20)
In this paper, we shall consider an incompressible forcing (∇ • f = 0) for simplicity, in which case φ i j in Eq. ( 20) is related to ψ i j in Eq. ( 13) as:

φ 11 = k 4 k 4 y ψ 11 , φ 12 = 1 k 4 y (k 2 k x k z ψ 11 + k 2 k 2 H ψ 13 ) , φ 22 = 1 k 4 y (k 2 x k 2 z ψ 11 + k 4 H ψ 33 + 2k x k z k 2 H ψ 13 ) . (21)
Note that φ 11 and φ 22 can signify strong radial forcing (e.g. by plumes) and horizontal forcing [e.g. due to the instability of the latitudinal shear in the tachocline (see also Kim 2005)]. For instance, in the case of strong radial forcing by plumes (ψ 11 ≫ ψ 33 ) with k H ≫ k x (see Fig. 2 in Leprovost & Kim 2006), φ 11 will dominate over φ 22 . Note further that in the 2D limit with k z = 0 and ψ 33 = 0, φ 22 vanishes.

We discuss some of the important implications of Eqs. ( 14)-( 19). First, all three components of fluctuating velocity amplitude in Eqs. ( 14)-( 16) are independent of N, clearly showing that the amplitude of the turbulent flow is not influenced by stable stratification in the case of the forcing with a short correlation time. The exact level of fluctuations is determined by the characteristics of the external forcing (φ i j ). We examine this in the simple 2D limit where k z = 0 (γ = 1) and ψ 33 = 0. In this limit, the substitution of φ 22 = 0 in Eqs. ( 14)-( 16)

gives v 2 x / v 2 y ∼ 1/a 2 = k 2 y /k 2
x and v 2 z = 0. This is an expected result for the 2D incompressible fluid in xy domain. If the forcing f contains only gravity modes, the wave number of the forcing would satisfy the local dispersion relation k x = ±k y N 2 /(ω -U 0 k y ) 2 -1 (here, ω is the frequency of the gravity waves). Furthermore, if these gravity modes are generated by the overturning of fluids in the convection zone (e.g. Press 1981) with U 0 ∼ 0, they would have a strong power at low frequencies ω ≪ N, thereby giving k x /k y ∼ N/ω ≫ 1, and thus v 2

x / v 2 y ≪ 1. Note that this is the situation normally considered in the previous works (e.g. Kim & MacGregor 2001;2003;[START_REF] Talon | [END_REF], where the main focus was on the momentum deposition by such gravity waves due to radiative damping after entering the tachocline from the bottom of the convection zone. In contrast, in this paper, the forcing f is not restricted to gravity modes, but is taken to be general, including any form of perturbation with arbitrary values of k x /k y . For instance, in the case of a strong radial forcing due to plumes with φ 11 ≫ φ 22 and a

= k x /k y ≪ 1 (see Fig 2 in Leprovost & Kim 2006), v 2
x / v 2 y ≫ 1 with a stronger radial fluctuation than horizontal one. Furthermore, Eqs. ( 14)-( 16) show that the level of anisotropy measured by the ratio of the turbulence amplitude depends on ν and µ. For instance, for an isotropic forcing with φ 11 ∼ φ 22 , v 2

x / v 2 y ∝ ν/µ ≪ 1 becomes small as the radiative damping µ increases. This is because a large µ tends to decrease (vertical) turbulence level by introducing large (thermal) dissipation. Thus, without a shear flow, a stronger horizontal turbulence (level) than vertical one can be caused by large thermal diffusivity µ (but not by stratification) in the case of a temporally short correlation forcing. Note that the anisotropy in turbulence level could be related to the Peclet number Pe = vl/µ as v 2

x / v 2 y ∝ Pe. Here, v and l are the characteristic velocity and length scale of the forcing (which is fixed).

Second, turbulent transport in Eqs. ( 17)-( 19) are strongly affected by stratification as N 2 increases, in contrast to fluctuation levels in Eqs. ( 14)-( 16). Since N ≫ µk 2 on reasonable scales l ≫ 10 4 ∼ 10 5 cm, the vertical (radial) mixing is severely quenched as the stratification increases (i.e. as N -2 ), in qualitative agreement with Brun, Turck-chiéze & Zahn (1999). In comparison, the part of the horizontal (latitudinal) mixing due to the radial forcing φ 11 is reduced proportional to N -2 while the one due to the horizontal forcing φ 22 is independent of N. The comparison of these two contributions (or alternatively D xx and D yy ) in the case of φ 11 ∼ φ 22 gives us a cutoff scale l c ∼ √ ν/N above which vertical (radial) mixing is strongly reduced compared to the horizontal (latitudinal) mixing. That is, in the presence of both radial and horizontal forcings of comparable strength, stable stratification mainly reduces the vertical transport without much effect on horizontal transport on scales l > l c . For parameter values typical of the tachocline ν ∼ 10 2 cm 2 s -1 and N ∼ 3 × 10 -3 s -1 , l c ∼ 10 2 cm. Thus, the stratification is very likely to play an important role over a broad range of physically reasonable scales. In the limit of strong radial forcing with φ 11 ≫ φ 22 and a ∼ 0, the stratification influences both the radial and horizontal transports to the same degree while the incompressibility renders

D xx T /D yy T ∼ D xx T /D zz T ∼ 1/a 2 ∼ v 2
x / v 2 y ≫ 1, with an effectively more efficient radial transport. It is interesting to compare our result D xx ∝ N -2 with the mixing due to radiatively damped waves (e.g. Garc ía López and Spruit 1991; [START_REF] Talon | [END_REF]. For instance, Garc ía López and Spruit (1991) estimated the vertical mixing due to gravity waves to be proportional to µ/N 2 . While the reduction in vertical mixing for large N is in qualitative agreement with our results (∝ N -2 ), the increase in vertical mixing in Garc ía López and Spruit (1991) is due to the fact that damped waves are necessary for wave transport. Finally, we note that without a shear flow, momentum transport vanishes (i.e. v x v y = v x v z = 0) for an isotropic forcing. Scaling of turbulence amplitude, turbulent viscosity (ν T ), and turbulent diffusivity D T are summarised in Table 1 for an isotropic forcing.

A large thermal diffusivity (µ = 10 5 ν) in the tachocline is often considered to reduce the stabilizing effect of stable stratification via the weakening of the buoyancy restoring force. To highlight this effect, it is illuminating to consider the extreme limit of strong thermal diffusion where density fluctuation becomes stationary with ∂ t ρ ′ = N 2 v x /gµk 2 ρ ′ = 0 in Eq. (3). In this limit, by following a similar analysis as previously (with D = ν), we can easily obtain the vertical and horizontal particle diffusivities as follows:

D xx ∼ τ f 2α 2 d 3 k (2π) 3 φ 11 (k) (γ + a 2 ) 2 , ( 22
)
D yy ∼ τ f d 3 k (2π) 3 γ 2 a 2 φ 11 (k) γ 2 (γ + a 2 ) 2 1 2α 2 + β 2 φ 22 (k) (2νk 2 ) 2 , ( 23
)
where α = γN 2 /µk 2 (γ + a 2 ). Equation ( 22) shows that the vertical mixing

D xx ∝ µ 2 /N 4 ∝ Pe -2
decreases for large N while increasing for large µ. This is because the reduction in the vertical mixing due to buoyancy force is weaken by a strong radiative damping µ. The comparison of Eqs.

( 22)-( 23) further shows that the reduction in vertical mixing relative to horizontal mixing is given by a factor of (νk 2 ) 2 /α 2 for an isotropic forcing and is thus weaker than that in the case of weak radiative damping µk 2 ≪ N [see Eqs. ( 17)-( 19)]. Furthermore, this reduction appears on scales larger than the critical scale l cµ = (µν) 1/4 /N 1/2 ∼ 10l c . Here, l c = √ ν/N is the critical scale in the case of µk 2 ≪ N. These results thus show that a strong thermal diffusion weakens the buoyancy effect and makes the effect of stratification become important on larger scales, compared to the case of a weak thermal diffusion. This result is thus consistent with the expectation employed in previous works.

To summarize, a stable stratification can dramatically quench turbulent transport with a more effective mixing in the horizontal directions orthogonal to the background density gradient. It does not however affect the amplitude of the turbulent flow, the ratio of which is found to depend only on ν/µ (for a temporally short-correlated forcing).

Consistent theory

The results in Sec. 2 showed that for a temporally short correlated forcing, a stable stratification reduces turbulent transport only, leading to anisotropic turbulent transport, without much effect on turbulence level. In this section, we study how these results are modified by a stable background shear flow (differential rotation in the tachocline). In particular, we will show that a shear flow not only inhibits the vertical mixing further, enhancing the anisotropic transport, but also reduces turbulence levels anisotropically, thereby leading to effectively stronger horizonal turbulence. For simplicity, we ignore the latitudinal differential rotation compared with the radial differential rotation since it is weaker in the tachocline due to thin tachocline (h < 0.03 ∼ 0.05 of the solar radius R). The inclusion of the latitudinal differential rotation would introduce a small correction term in our results. For instance, for turbulence amplitude, this correction term is of order (h/R) 2 ≪ 1, as shown in Leprovost and Kim (2006). Note that the latitudinal shear is crucial for non-vanishing horizontal momentum transport in Leprovost & Kim (2006). Again,

A = 0, N 0 A 0, N = 0 A 0, N 0 v 2 x µ -1 A -1 A -1 v 2 y ∼ v 2 z ν -1 A -1 ξ -1/3 ν A -1 ξ -1/3 ν ν T 0 A -2 -A -2 D xx T N -2 A -2 N -2 ξ D ξ -2/3 µ D yy T ∼ D zz T ν -1 (ν + D) -1 A -2 ξ -2/3 ν A -2 ξ -2/3 ν v 2 y / v 2 x µ/ν ξ -1/3 ν ξ -1/3 ν D yy T /D xx T N 2 ξ -2/3 ν (N/A) 2 (µ/ν) 2/3 ξ -1 ν Table 1.
Scaling of turbulence amplitude, turbulent viscosity (ν T ), and turbulent diffusivity D T for an isotropic forcing with φ 11 ∼ φ 22 in the case N ≫ µk 2 and D ∼ ν ≪ µ. ξ µ = µk 2 y /A ≪ 1, ξ ν = νk 2 y /A ≪ 1, and ξ D = Dk 2 y /A ≪ 1 are small parameters representing strong shear limit. The second (A = 0 and N 0), the third (A 0 and N = 0), and the fourth (A 0 and N 0) columns contain the results for stratified unsheared case in Sec. 2, for unstratified sheared case (Kim 2005), and for stratified sheared case in Sec. 3, respectively.

we envision that turbulence is maintained in the tachocline by an external forcing while gravity waves are excited due to this external forcing in the stably stratified tachocline. We shall then compute the overall turbulent transport consistently by taking into account the interaction among turbulence, shear flow and gravity waves, instead of simply assuming a (large) constant value of turbulent viscosity for mean shear flow (and gravity waves). Note that this treatment is essential when there is no clear scale separation between gravity waves and turbulence, in which case turbulence cannot be considered to give an enhanced value of viscosity for gravity waves (c.f.

Charbonnel & Talon 2005).

For the evolution of fluctuations, we approximate the radial differential rotation by a linear shear flow with U 0 = -xAŷ to keep the analysis tractable. Here, A is the shearing rate which we assume to be positive without loss of generality. As done in previous papers (Kim 2005;Leprovost & Kim 2006), to capture the effect of shearing due to radial differential rotation (A ∼ 3 × 10 -6 s -1 for the tachocline) non-perturbatively, we use the special Fourier transform for fluctuating quantities φ ′ :

φ ′ (x, t) = 1 (2π) 3 d 3 k φ(k, t) exp {i(k x (t)x + k y y + k z z)} . ( 24
)
Here, k x = k x (t) is the time dependent [unlike constant k x in Eq. ( 5)], satisfying an eikonal equation

∂ t k x (t) = k y A . ( 25
)
Equation ( 25) implies that k x linearly increases in time as k x (t) = k x (0) + k y At, manifesting the main effect of shearing by a shear flow U 0 (x)ŷ, i.e., generation of fine scales in the x direction due to tilting and distortion of fluid eddies (e.g. see Burrell 1997;Kim 2004;Kim 2005). The efficient generation of fine scales by shearing leads to the breakup of eddies and enhancement of the overall dissipation, thereby reducing turbulence amplitude and transport (e.g. see Kim 2005).

A similar effect by a shear flow is expected to persist for a more realistic radial shear [e.g. for an error function used in helioseismic inversions (Kosovichev 1996;Corbard et al 1999)] since the basic mechanism of shearing (e.g. see Fig. 2 in Kim 2005) is the same regardless of the details of the profile of radial shear. The efficiency of the shearing could depend on the details of the profile, possibly leading to a slightly different scaling.

It is interesting to note that a gravity wave with an initially positive value of k x (0) can change its sign after the time interval k x (0)/k y A (for k y > 0) due to shearing. Since the local (radial) group velocity of gravity waves is given by v gx = -k 2 y k x N 2 /k 4 (ω -U 0 k y ) (e.g., see Kim & MacGregor 2003), the gravity wave thus alters its propagation direction as k x flips its sign. Therefore, a gravity wave which initially propagates downward to the interior from the convection zone with a negative vertical group velocity (i.e. v gx < 0) can propagate upwards when the vertical group velocity becomes positive (v gx > 0) due to shearing.

For parameter values typical of the tachocline N ∼ 3 × 10 -3 s -1 and A ∼ 3 × 10 -6 s -1 , Ri = N 2 /A 2 ≫ Ri c ≃ 1/4, satisfying the stability criterion [START_REF] Lighthill | Waves in Fluids[END_REF]. We thus assume that the radial shear flow is stable with large value of Ri = N 2 /A 2 ≫ 1 in the remainder of the paper (see also [START_REF] Schatzman | [END_REF]. Note that the buildup of chemical composition gradient (the so-called µ gradient in the solar interior [e.g. see Michaud and Zahn (1998) and references therein] would further increase the values of N and Ri, making the radial shear flow more stable, although this effect could be counteracted by a radiative damping, as shown in Sec 2. The shearing rate A ∼ 3 × 10 -6 s -1 due to this radial differential rotation is larger than the dissipation rate due to radiative damping µk 2 y on a broad range of scales l(= 1/k y ) > 10 6 ∼ 10 7 cm ∼ 10 -3 H 0 , where H 0 is the pressure scale height ∼ 6 × 10 9 cm. Thus, we focus on the strong shear limit in the following by using ξ µ = µk 2 y /A ≪ 1 as a small parameter.

For Ri ≫ 1 and ν ∼ D ≪ µ, a long but straightforward algebra can give us the solutions to Eqs. ( 1)-( 4), as shown in Appendix A. By using these solutions and the correlation functions of forcing defined in Eq. ( 20), we obtain the following results for turbulence level and transport coefficient defined by n ′ v i = -D i j T ∂ j n 0 for i = 1, 2 and 3, and momentum flux v x v y = -ν T ∂ x U 0 in the strong shear limit ξ µ = µk 2 y /A ≪ 1 (see Appendix A for details):

v 2 x ≃ τ f A d 3 k (2π) 3 φ 11 (k) 2γ 3/2 , ( 26
)
v 2 y ≃ τ f A d 3 k (2π) 3 φ 11 (k) 2γ 5/2 | ln ξ µ | + β 2 φ 22 (k) γ 2 G 0 , (27)
v 2 z ≃ τ f A d 3 k (2π) 3 β 2 φ 11 (k) 2γ 5/2 | ln ξ µ | + φ 22 (k) γ 2 G 0 , (28)
D xx T ≃ τ f 2N 2 d 3 k (2π) 3 φ 11 (k) γ γ + a 2 ξ D G 1 , (29)
D yy T ≃ τ f 2 d 3 k (2π) 3 γ 2 φ 11 (k)Γ( 4 3 ) 2N 2 γ γ + a 2 D µ 3 ξ µ 1 3 + β 2 φ 22 (k)Γ( 5 3 ) A 2 3 2ξ ν 2 3 , (30)
ν T ≃ - τ f 2 d 3 k (2π) 3 φ 11 (k) γ(γ + a 2 ) 1 A 2 + 1 12N 2 (γ + a 2 ) . (31)
Here,

ξ ν = νk 2 y /A ≪ 1, ξ D = Dk 2 y /A ≪ 1, ξ µ = µk 2 y /A ≪ 1, G 0 = 1 3 Γ (1/3) (3/2ξ ν ) 1/3 and G 1 = 1 3 Γ (2/3) 3/ξ µ 2/3
. Note again that ξ µ ≪ 1 is valid on a broad range of scales l(= 1/k y ) > 10 6 ∼ 10 7 cm in the tachocline and also that ξ µ ≪ 1 guarantees that ξ D ∼ ξ ν ≪ 1 since

D ∼ ν ≪ µ.
The spectrum φ i j in Eqs. ( 26)-( 31) are given by Eqs. ( 20), which are related to power spectrum of forcing ψ i j in Eq. ( 13) in the incompressible case.

Equations ( 26)-( 31) reveal the following interesting features. Turbulence levels given in Eqs.

(26)-( 28) are again independent of stratification, similarly to the case without a shear flow [Eqs.

( 14)-( 16)] while they are reduced for strong shear A. This indicates that waves and shear flows play different roles in turbulence regulation -waves do not necessarily quench fluctuation levels while shear flows can reduce them through enhanced dissipation via shearing. The turbulence regulation by shearing gives us horizontal velocity fluctuations in Eqs. ( 27) and (28) which are effectively higher than vertical one in Eq. ( 26). While a similar tendency was also found in the absence of gravity waves (Kim 2005), the exact value of the ratio of vertical to horizontal turbulence levels is not the same. For example, v 2 x / v 2 y ∝ (| ln ξ µ |) -1 for a strong radial forcing φ 11 with φ 22 = 0 while v 2

x / v 2 y ∝ ξ 1/3 ν for an isotropic forcing with φ 11 ∼ φ 22 . Note that ξ ν ≪ ξ µ ≪ 1 are small parameters in our problem, representing the strong shear limit. These results are to be compared with v 2 x / v 2 y ∝ ξ 1/3 ν in the unstratified medium (Kim 2005). The results for the isotropic forcing in various cases are summarised and compared in Table 1.

Transport properties in Eqs. ( 29)-( 31) however exhibit a very different behaviour, with both vertical and horizontal mixing being inhibited in a non-trivial manner by strong stratification as well as by shearing. First, the vertical transport is reduced as

D xx T ∝ ξ D ξ -2/3 µ N -2 ∝ Dµ -2/3 A -1/3 N -2
, becoming small as either stratification or shearing increases. Note that the decrease in D xx T for large N agrees with Miesch (2003). Interestingly, D xx T ∝ µ -2/3 ∝ Pe 2/3 decreases as the radiative damping µ increases. This is because large radiative damping increases thermal dissipation, thereby mainly inhibiting the vertical mixing, as noted in Sec. 2. Thus, compared to the case with N = 0 where D xx T ∝ A -2 (Kim 2005), the vertical mixing is much more quenched by a factor of

ξ D ξ -2/3 µ (A/N) 2 ≪ ξ D ξ -2/3 µ = (D/µ)ξ 1/3
µ ≪ ξ µ ≪ 1 since A/N < 1 and D/µ ∼ 10 -5 in the tachocline (see also Table 1). This clearly shows that shear flow (orthogonal to radial density gradient), stable stratification, and radiative damping can all inhibit the radial transport.

Second, D yy T is less affected by stratification since it involves the two parts -the one from φ 11 is proportional to ξ D ξ -4/3 µ N -2 while the other from φ 22 is proportional to A -2 ξ -2/3 ν , independent of N. In the simplest case of a strong radial forcing with φ 22 = 0, D xx T /D yy T ∼ ξ 2/3 µ ≪ 1, independent of N. This is similar to the case without stratification where D xx T /D yy T ∼ ξ 2/3 ν (Kim 2005). In the general case where φ 11 0 and φ 22 0, D yy T exhibits a non-trivial, complex interplay among stratification and shear flow in determining the overall transport. To appreciate this, we compare Eq. ( 30) with the result obtained without shear flow (18) to find that shear flow enhances the contribution from φ 11 by a factor ξ D ξ -4/3 µ = (D/µ)ξ -1/3 µ while it reduces the part from φ 22 by a factor of (ν/A) 2 ξ -2/3 ν ∝ ξ 4/3 ν ≪ 1. Since D/µ ∼ 10 -5 ≪ 1 and ξ µ ≪ 1, the horizontal transport driven by a strong radial forcing φ 11 can be either inhibited or enhanced by a shear flow (due to φ 11 ) depending on the parameter values. On the other hand, the horizontal transport due to φ 22 is reduced by the shear flow by a factor proportional to ξ 4/3 ν ≪ 1. The end result can easily be shown to be the enhancement of the effect of stratification. To see this, we compare the two contributions from φ 11 and φ 22 to D yy T and obtain a characteristic scale l * = (D/µ) 2 (ν/N) 1/2 (A/N) 5/2 above which the stratification is important, being mainly responsible for quenching transport.

Since l * is smaller than l c = (ν/N) 1/2 obtained without a shear flow, the effect of stratification becomes important compared to the case without shear flows. It is also worth comparing the contributions to D yy T from φ 11 and φ 22 separately to those in the case of N = 0 in Kim (2005), where D yy T ∝ A -2 ξ -2/3 ν from both φ 11 and φ 22 . That is, the contribution from φ 11 is further reduced by a factor of (A/N) 2 ξ 1/3 ν (ν/µ) 4/3 by stratification, while the one from φ 22 is not affected. We now examine how the stratification affects the anisotropy in transport. Equations ( 29) and ( 30) for a strong radial forcing with φ 11 ≫ φ 22 and a ∼ 0 gives D xx T /D yy T ∼ ξ 2/3 µ ≪ 1, independent of N. This is analogous to the result D xx T /D yy T ∼ ξ 2/3 ν ≪ 1 obtained in the unstratified case (Kim 2005). In this case, the anisotropy in the transport is solely caused by shear stabilization. For an isotropic forcing with φ 11 ∼ φ 22 and a ∼ 1, D xx T /D yy T ∼ (A/N) 2 (ν/µ) 2/3 ξ ν ≪ 1. In other words, the anisotropy in turbulent transport depends on stratification, shear, and ν/µ, as clearly shown in Table 1. Since (A/N) 2 ≪ 1 and ν/µ ≪ 1 in the tachocline, D xx T /D yy T is much smaller than ξ 2/3 ν obtained in the unstratified case (Kim 2005). That is, the anisotropy in transport is further enhanced due to stratification. We emphasize that the anisotropy in turbulent transport is much stronger than that in turbulence amplitude, discussed previously (see also Table 1). This result also shows the reduction in D xx T for large µ and again highlights the importance of the radiative damping in reducing the vertical transport, thereby increasing the anisotropy in the transport. Finally, Eq. ( 31) demonstrates one of the most important effects of a stable stratification, which is to drive a system away from a uniform rotation with a negative eddy viscosity. Recall that in the absence of stratification the eddy viscosity is positive in 3D while negative in 2D limit (e.g. see Kim 2005). In contrast, the eddy viscosity in Eq. ( 31) is negative in both 2D (β = 0 and γ = 1) and 3D cases. This behavior was also found in recent numerical simulation by Miesch (2003) of a stably stratified turbulence with an imposed shear driven by penetrative convection, who found anti-diffusive radial and diffusive latitudinal momentum transport, thereby offering a mechanism for a proper transition from latitudinal differential rotation in the convection zone to solid body rotation in the radiative interior. Anti-diffusive momentum transport is a generic feature of a strongly stratified medium (i.e. a geophysical system). It is interesting to note that a negative viscosity was also found in the previous work on momentum transport due to radiatively damped gravity waves (e.g. Kim and MacGregor 2001;2003). In addition, the result (31) shows that the (anti-diffusive) momentum transport becomes less efficient for strong stratification (large N), in agreement with Miesch (2003).

To summarise, our results show that the shearing effect by radial differential rotation together with gravity waves is an important mechanism for turbulence regulation, leading to a weak turbulent transport and anisotropic turbulence and transport in the tachocline. Furthermore, we have consistently derived the values of turbulence level, particle mixing and momentum transport starting from the first principle, clearly identifying the different roles of gravity waves and shearing in transport. For instance, in comparison with Chaboyer and Zahn (1992) [or Spiegel and Zahn (1992)] which start with the assumption of a strong horizontal mixing, we identified the source of such an anisotropic turbulence. In particular, we have made a clear distinction between the anisotropy in the turbulence level and turbulent transport.

Discussion and conclusions

We have studied turbulent transport in the stably stratified tachocline with a strong radial differential rotation, when turbulence is driven and maintained by a forcing (e.g. due to plumes penetrating from the convection zone or due to instability). We have assumed that both turbulence and gravity waves are on small scales (with no clear scale separation between the two) and treated the interaction among gravity waves, turbulence and shear flow consistently. Unlike a shear flow which regulates both turbulence level and transport, a stable stratification is shown to mainly inhibit turbulent transport, leading to a further reduction in transport compared to the unstratified case (Kim 2005). Specifically, for parameter values typical of the tachocline (N/A ≫ 1), particle transport due to a strong radial forcing φ 11 (with φ 22 = 0) is reduced as

ξ D ξ -2/3 µ N -2 , becoming much smaller than horizontal transport (∝ ξ D ξ -4/3 µ N -2 ) by a factor of ξ 2/3 µ ≪ 1. Here, ξ ν = νk 2 y /A ≪ 1, ξ D = Dk 2 y /A ≪ 1
, and ξ µ = µk 2 y /A ≪ ξ ν are the small parameters characterizing a strong shear limit (see the main text for more details). Note that in this case D xx T /D yy T ∝ ξ 2/3 µ ≪ 1 depends only on shearing but not on stratification, indicating that the anisotropy in particle transport is mainly governed by radial differential rotation. A similar scaling (∝ ξ 2/3 ν ≪ 1) was also found in the unstratified case (Kim 2005). However, in the case of an isotropic forcing with φ 11 ∼ φ 22 and k x /k y ∼ 1, the horizontal mixing is much less reduced with D yy T ∼ A -2 ξ -2/3 ν (with no effect of stratification), leading to a stronger anisotropy in transport with D xx T /D yy T ∼ (A/N) 2 (ν/µ) 2/3 ξ ν ≪ ξ 2/3 ν ≪ 1. Note that the anisotropy becomes stronger for larger µ.

Furthermore, the vertical momentum transport was shown to be anti-diffusive with a negative eddy viscosity. That is, small-scale turbulence influenced by gravity waves accentuates the gradient in a radial differential rotation rather than makes it uniform. This is similar to the tendency obtained in the case of momentum deposition by gravity waves due to radiative damping (Kim & MacGregor 2001;2003). The sharpening of the gradient of the radial shear due to the

̟(τ) = α - 1 8N 2 γ γ + τ 2 , α = 1 - 1 8N 2 .
(A.10)

By following a similar algebra used in Kim (2005) and by using the correlation function of the forcing given in Eq. ( 20), we can obtain the following correlation fictions v i v j to leading orders in 1/N 2 ≪ 1 and ξ µ = µk 2 y /A ≪ 1: ,a) (γ + τ 2 ) 3/2 1 + cos [2(ϕ(τ) -ϕ(a))] , ,a) φ 11 (k) 2 γ + a 2 τ 2 γ 2 (γ + τ 2 ) 3/2 + e -2ξ ν Q(τ,a) β 2 φ 22 (k) ,a) φ 11 (k) 2 γ + a 2 τ 2 β 2 γ 2 (γ + τ 2 ) 3/2 + e -2ξ ν Q (τ,a) Here, a = k x /k y , Q(τ, a) = (γτ + τ 3 /3) -(γa + a 3 /3); special symmetries ψ i j (k y ) = ψ i j (-k y ) and ψ i j (k z ) = ψ i j (-k z ) were used. We then compute τ integrals in Eq. (A.11) in the strong shear limit ξ ν ≪ ξ µ ≪ 1 and use v x v y = -ν xx T ∂ x U 0 = ν xx T A to obtain Eqs. ( 26)-( 28) and (31). Next, to compute particle transport, we integrate Eq. (A.7) in time to obtain ñ(k(t), t) = -∂ i n 0 dt 1 d 3 k 1 ĝ(k, t; k 1 , t 1 )e -DQ(t,t 1 ) ṽi (k 1 , x, t 1 ) .

v 2 x ∼ τ f 2A d 3 kdτ (2π) 3 φ 11 (k) γ + a 2 e -ξ µ Q(τ
v 2 y ∼ τ f A d 3 kdτ (2π) 3 e -ξ µ Q(τ
γ 2 , v 2 z ∼ τ f A d 3 kdτ (2π) 3 e -ξ µ Q(τ
(A.12)

Here, Q(t, t 1 )

= t t 1 dt ′ [k 2 x (t ′ ) + k 2 H ] = [k 3 x -k 3 1x ]/3k y A + k 2 H (t -t 1 ); k 2 H = k 2 y + k 2
z is the amplitude of wave number in the horizontal plane; k 2 = k 2 H + k 2 x ; ĝ is the Green's function given by ĝ(k, t; k 1 , t 1 ) = δ(k y -k 1y )δ(k z -k 1z )δ k x -k 1x -k 1y (tt 1 )A . (A.13) A similar analysis using Eqs. (A.9), (A.13), ( 11)-( 13), and (20) and nv i = -D i j ∂ j n 0 gives (τ,a) γ + a 2 (γ + τ 2 ) 1/2 , ,a) φ 11 (k)τ 2 γ + τ 2 2N 2 γ + a 2 + e -2ξ ν Q(τ,a) β 2 φ 22 (τa) 2 (γ + τ 2 )

D xx T = τ f 2(2π) 3 A 2 N 2 d 3 kdτ ξ D φ 11 (k)e -ξ µ Q
D yy T = τ f (2π) 3 A 2 d 3 kdτ ξ D γ 2         e -ξ µ Q(τ
        . (A.14)
Here, ξ D = Dk 2 y /A ∼ ξ ν ≪ 1. Finally, the evaluation of τ integrals in Eq. (A.14) gives us Eqs. ( 29)-(30) in the main text.

Fig. 1 .

 1 Fig. 1. The configuration of our model.

  γ + a 2 ̟(a) 2 e -ξ µ Q(τ,a) τ (γ + τ 2 ) 5/2 × -τ 2 sin [ϕ(τ) -ϕ(a)] + N γ + τ 2 ̟(τ) cos [ϕ(τ) -ϕ(a)]

Acknowledgements. This work was supported by the UK PPARC grant PP/B501512/1.

negative viscosity could eventually lead to time variation in the tachocline (similarly to Kim & MacGregor 2001) or instability (e.g., see [START_REF] Petrovay | [END_REF], causing a rapid radial mixing and thus reducing the anisotropy.

Even in the stably stratified radiative interiors of stars as well as in the tachocline, background turbulence has often been assumed to be present to enhance the value of effective viscosity. While the clarification of the source of turbulence responsible for such an enhanced eddy viscosity is an interesting problem, our results demonstrate how this residual turbulence interacts with gravity waves and shear flow, providing the prediction for the values of eddy viscosity as well as the particle diffusivity which depend on physical quantities like A and N. In particular, the

indicates that turbulent transport of particles can be inhibited due to stable stratification, shear flow and large radiative damping. This finding can have interesting implications for the surface depletion of lithium in the Sun and other stars (i.e., Pinsonneault 1997), and will be studied in a future publication. Furthermore, if a similar physical process operates in the bulk of the radiative interior of the sun, a negative viscosity that we obtain implies that a radial differential rotation which is created during its spin-down would not be eliminated. Note however that Charbonnel & Talon (2005) have shown an efficient momentum transport in the solar radiative interior due to the cumulative effect of large scale meridional flow, shear instability and gravity waves.

The tachocline is believed to possess a strong toroidal magnetic fields of strength 10 4 ∼ 10 5 G, which could have an important influence on turbulent transport in that region. This is an important problem since the presence of a weak poloidal magnetic field in the radiative interior together with a strong toroidal magnetic field in the tachocline (acting as a boundary layer between the radiative interior and convection zone) could offer a mechanism for a uniform rotation in the interior as well as for the tachocline confinement (e.g. Rüdiger & Kichatinov 1996;MacGregor & Charbonneauu 1997;Gough & McIntyre 1998) [see however Brun & Zahn (2006) for a negative result on this scenario]. In this case, the values of effective diffusivity of magnetic fields as well as eddy viscosity play a crucial role in determining the thickness of the tachocline. Therefore, a consistent computation of magnetic diffusivity and eddy viscosity with their dependences on physical quantities such as the strength of magnetic fields, the Brunt-Väisälä frequency, and shearing rate, would be of primary interest (Kim and Leprovost 2007). It is an interesting question, in general, whether a negative eddy viscosity, favored in a stably stratified medium, remains as a robust feature in the presence of magnetic field in view of the forward energy cascade (i.e.

positive eddy viscosity) in MHD turbulence. Furthermore, composition gradients (discussed in Sec. 3) and meridional flows would contribute to the transport in the tachocline. In particular, meridional flows are expected to enhance the radial transport of chemical species by advection as discussed in Kim (2005) although this enhancement would be reduced for stronger horizontal turbulent [e.g., see Kim (2005) and Chaboyer and Zahh (1992)]. These issues are however outside the scope of this paper and will be addressed in future publication.

Appendix A:

In this Appendix, we show how to derive Eqs. ( 26)-(31). To this end, we use U 0 (x)ŷ = -xAŷ in Eqs. ( 1)-( 4), and introduce the transforms ŵ and ˆŵ for w = v i , n, p, h i and f as follows:

where

We then use a new time variable τ = k x /k y + At to rewrite Eqs. ( 1)-( 4) as:

A∂ τ vy -Av x = -ik y p + fy , (A.3)

For a strong shear limit where µ = µk 2 y /A ≪ 1, Eqs. (A.2)-(A.6) can be combined to yield the equation for ˆρ as A.8) to leading order in ξ µ = µk 2 y /A ≪ 1. Here, ĥ1 (τ) = γ ĥxτ ĥyτβ ĥz ; γ = 1 + β 2 and β = k z /k y . In the limit of N 2 /A 2 ≫ 1, Eq. (A.8) can be solved with the following solutions valid up to O(N -2 ):