On the interaction between tides and convection

Gordon I. Ogilvie, Geoffroy Lesur

► To cite this version:

Gordon I. Ogilvie, Geoffroy Lesur. On the interaction between tides and convection. Monthly Notices of the Royal Astronomical Society, 2012, 422, pp.1975-1987. 10.1111/j.1365-2966.2012.20630.x. insu-03612428

HAL Id: insu-03612428

https://insu.hal.science/insu-03612428

Submitted on 18 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the interaction between tides and convection

Gordon I. Ogilvie<fnref n="1" />* and Geoffroy Lesur<fnref n="1" />,<fnref n="2" />

Accepted 2012 January 14. Received 2012 January 5; in original form 2011 November 2

## ABSTRACT

We study the interaction between tides and convection in astrophysical bodies by analysing the effect of a homogeneous oscillatory shear on a fluid flow. This model can be taken to represent the interaction between a large-scale periodic tidal deformation and a smaller scale convective motion. We first consider analytically the limit in which the shear is of low amplitude and the oscillation period is short compared to the time-scales of the unperturbed flow. In this limit there is a viscoelastic response and we obtain expressions for the effective elastic modulus and viscosity coefficient. The effective viscosity is inversely proportional to the square of the oscillation frequency, with a coefficient that can be positive, negative or zero depending on the properties of the unperturbed flow. We also carry out direct numerical simulations of Boussinesq convection in an oscillatory shearing box and measure the time-dependent Reynolds stress. The results indicate that the effective viscosity of turbulent convection falls rapidly as the oscillation frequency is increased, attaining small negative values in the cases we have examined, although significant uncertainties remain because of the turbulent noise. We discuss the implications of this analysis for astrophysical tides.

**Key words:** convection – hydrodynamics – turbulence – planets and satellites: general – planet-star interactions – binaries: general.

## 1 INTRODUCTION

Tidal interactions determine the orbital and spin evolution of astro-
physical bodies when they orbit sufficiently close to one another.
Applications include close binary stars, extrasolar planetary systems
and the satellites of Solar system planets.

In many cases of interest the tidally forced bodies are fully or
partly convective. In order to determine the rate of tidal evolution
it is therefore necessary to study the interaction between tides and
convection. The response of a fluid body to tidal forcing generally
consists of two components: one is non-wavelike and of large scale,
while the other involves internal waves of smaller scale. At the
simplest level the convection could be thought of as providing an
effective viscosity that damps the non-wavelike tidal disturbance
and provides a phase shift in the response. As in the theory of Darwin
(1880), the tidal torque and the rate of tidal evolution are then
directly proportional to this effective viscosity. Convection may also
play an important role in dissipating inertial waves, which constitute
the low-frequency wavelike response of rotating convective zones
of stars and giant planets (e.g. Ogilvie & Lin 2004).

As pointed out by Zahn (1966) and Goldreich & Nicholson
(1977), the effective viscosity estimated from mixing-length the-

ory ought to be reduced when the period of the tidal disturbance
is short compared to the characteristic time-scale of the convec-
tive motion. The suppression factor has been the subject of much
debate, informed by observational constraints such as the apparent
rate of circularization of the orbits of close binary stars (Meibom
& Mathieu 2005, and references therein). Zahn (1966) suggested
a less severe suppression, such that the effective viscosity is pro-
portional to the oscillation period for short periods, which is in
better agreement with these observations. Goldreich & Nicholson
(1977) and Goldreich & Keeley (1977) considered the multiscale
nature of the turbulent flow, and argued that the dominant contri-
bution to the effective viscosity at short periods comes from eddies
with a convective time-scale comparable to the oscillation period.
For a Kolmogorov spectrum, this argument gives a more powerful
suppression, such that the effective viscosity is proportional to the
square of the oscillation period for short periods.

While these authors relied on simple physical arguments and
order-of-magnitude estimates, Goodman & Oh (1997), who also
provide a clear review of the controversy, introduced a more formal
procedure for determining the effective viscosity of a convective
flow, by considering the effect of a homogeneous oscillatory strain
on a turbulent velocity field. Analytical progress was limited by the
complexity of the equations. While their method relies on an expan-
sion in powers of the ratio of the oscillation period to the convective
time-scale, it leads to a result in which the dominant contribution

<footnote n="1">¹Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge,
Wilberforce Road, Cambridge CB3 0WA</footnote>

<footnote n="2">²UJF-Grenoble I/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, Grenoble F-38041, France</footnote>

<footnote>*E-mail: gio10@cam.ac.uk</footnote>
comes from eddies for which these time-scales are similar. Their argument does, however, appear to rule out the hypothesis of Zahn (1966).

More recently, numerical simulations of convection have been brought to bear on this question. Penev et al. (2007) and Penev et al. (2009a) applied the procedure of Goodman & Oh (1997) to the velocity fields obtained in numerical simulations of convection in a deep layer to deduce the effective viscosity as a function of the tidal frequency. In Penev, Barranco & Sasselov (2009b) an oscillatory forcing was introduced directly into the convection simulation and the effective viscosity was estimated by measuring the work done by this force. The results of these studies suggest that something closer to the prescription of Zahn (1966) may be appropriate, not for the reasons originally suggested, but possibly because the power spectrum of the convection is less steep than the Kolmogorov spectrum assumed by Goldreich & Nicholson (1977).

While the results of Penev and collaborators are of considerable interest, the uncertainties in these calculations have not been quantified. Owing to the importance of this problem, we are motivated to examine it from an independent viewpoint. We have devised a theoretical framework that allows us to study statistically homogeneous flows; while this local viewpoint omits important global aspects of stellar or planetary convection, it allows us to resolve convection at high Rayleigh numbers and to quantify the uncertainties due to turbulent noise. We also present complementary analytical results which shed new light on previous theoretical discussions.

The plan of this paper is as follows. In Section 2 we introduce the homogeneous sheared coordinate system used in our analytical and numerical calculations, and discuss the equations of fluid dynamics in this system. In Section 3 we derive the response of a fluid flow to oscillatory shear in the limits of low amplitude and high frequency, and calculate the effective elasticity and viscosity in various cases. Direct numerical simulations of convection in an oscillatory shearing box are reported and interpreted in Section 4, followed by a summary and discussion of the results.

# 2 THE OSCILLATORY SHEARING BOX

## 2.1 Motivation

At the present time it is not practical to compute a global model of an astrophysical body with turbulent convection and to measure its tidal response through the direct application of a time-dependent gravitational potential. As noted in Section 1, the response of a fluid body to tidal forcing generally consists of a large-scale non-wavelike motion together with some internal waves. The non-wavelike tide can be computed without difficulty in a linear approximation and consists of a time-dependent spheroidal deformation of the body. From the local perspective of the smaller scale convective motion, this tidal flow appears as a spatially homogeneous, approximately incompressible motion that is oscillatory in time. We are interested, therefore, in determining the response of the convective motion to such a deformation and the additional stresses that result from this interaction.

A general, 3D, spatially homogeneous, incompressible deformation can be represented by a traceless velocity gradient tensor $\nabla u$ that is independent of position. In this paper we consider a time-dependent simple shear of the form $u_y \propto x$. Using a linear combination of simple shears with different orientations, it is possible to compose an arbitrary velocity gradient tensor that is traceless (see Appendix A). Provided that we are working a linear regime, there-

fore, it is sufficient to measure the response to a time-dependent simple shear.

## 2.2 Sheared coordinates

We initially consider an incompressible fluid in two dimensions, satisfying the Navier-Stokes equations:

$$ (\partial_t + u_x \partial_x + u_y \partial_y)u_x = -\partial_x p + v (\partial_x^2 + \partial_y^2) u_x + \tilde{f}_x, \quad (1) $$

$$ (\partial_t + u_x \partial_x + u_y \partial_y)u_y = -\partial_y p + v (\partial_x^2 + \partial_y^2) u_y + \tilde{f}_y, \quad (2) $$

$$ \partial_x u_x + \partial_y u_y = 0, \quad (3) $$

where $u$ is the velocity, $p$ is the pressure divided by the density, $v$ is the kinematic viscosity and $\tilde{f}$ is a body force per unit mass.

We suppose that the system is subject to a homogeneous time-dependent shear. We define the sheared coordinates

$$ x' = x, \quad y' = y - a(t)x, \quad t' = t. \quad (4) $$

The dimensionless quantity $a$ is the shear, and $\dot{a} = da/dt$ is the shear rate. The Jacobian determinant of the transformation is unity. Partial derivatives transform according to

$$ \partial_x = \partial_x' - a \partial_y', \quad \partial_y = \partial_y', \quad \partial_t = \partial_t' - \dot{a} x \partial_y'. \quad (5) $$

We write

$$ u_x = v_x, \quad u_y = v_y + \dot{a}x, \quad (6) $$

so that $v$ represents the velocity relative to the shearing motion. We continue to refer vector components to the original Cartesian basis.

With these substitutions, the Navier-Stokes equations become

$$ [\partial_t' + v_x(\partial_x' - a\partial_y') + v_y\partial_y']v_x = -(\partial_x' - a\partial_y')p + v\Delta v_x + \tilde{f}_x, \quad (7) $$

$$ [\partial_t' + v_x(\partial_x' - a\partial_y') + v_y\partial_y']v_y + \ddot{a}x' + \dot{a}v_x = -\partial_y' p + v\Delta v_y + \tilde{f}_y, \quad (8) $$

$$ (\partial_x' - a\partial_y')v_x + \partial_y' v_y = 0, \quad (8) $$

involving the Laplacian operator in unsheared coordinates,

$$ \begin{aligned} \Delta &= \partial_x^2 + \partial_y^2 \\ &= (\partial_x' - a\partial_y')^2 + \partial_y'^2 \\ &= \Delta' - 2a\partial_x'\partial_y' + a^2\partial_y'^2. \end{aligned} $$

These equations are spatially homogeneous (not involving $x'$ or $y'$ explicitly) except for the term $\ddot{a}x'$, which arises because a spatially inhomogeneous force is required to maintain an accelerating shear. By writing

$$ \tilde{f}_x = f_x, \quad \tilde{f}_y = \ddot{a}x' + f_y, \quad (9) $$

i.e. by subtracting the inhomogeneous force, we restore spatial homogeneity to the equations. In this representation, $f$ is the body force that, if necessary, drives the (possibly turbulent) motion whose response we wish to measure, while $\ddot{a}x'$ is the force that drives the imposed large-scale shear.

The property of spatial homogeneity results from the translational invariance of the homogeneously sheared system, given that a Galilean transformation easily removes the velocity shift associated with a translation in the $x$ direction. It means that statistically homogeneous turbulence, for example, is possible in this system. It also means that, for computational or analytical purposes, we can
impose periodic boundary conditions on **v** and **p** in sheared coordinates. When employing a spectral method, these variables can be expanded in Fourier series in sheared coordinates, using basis functions exp(ik'_x x' + ik'_y y').

An alternative approach, which leads to identical results, is to solve the original equations in unsheared coordinates but to apply modified periodic boundary conditions of the form

$$p(L_x, y, t) = p(0, (y - a(t)L_x) \mod L_y, t), \quad (10)$$

$$p(x, L_y, t) = p(x, 0, t), \quad (11)$$

and similarly for **v** (not **u**), where $L_x$ and $L_y$ are the dimensions of the shearing box. The appropriate basis is then one composed of shearing waves exp(ik'_x(t)x + ik'_y(t)y), having (quantized) time-dependent wavevectors with components $k'_x(t) = k'_x - a(t)k'_y$ and $k'_y = k'_y$. With these definitions, exp(ik'_x(t)x + ik'_y(t)y) = exp(ik'_x x' + ik'_y y').

Note that the oscillatory shear in our model is not imposed by the boundary conditions as such, but is driven by the inhomogeneous body force $\ddot{a}x' e_y$ as described above. Our model therefore provides a self-consistent local representation of any turbulent flow subject to an oscillatory deformation due to an external body force.

Similar considerations apply to the Boussinesq equations in three dimensions, or indeed to the equations of compressible convection. In the Boussinesq case the basic equations in unsheared coordinates are

$$(\partial_t + u_j \partial_j)u_i = -\partial_i p + \nu \Delta u_i + b e_i + \tilde{f}_i, \quad (12)$$

$$(\partial_t + u_i \partial_i)b + N^2 u_i e_i = \kappa \Delta b, \quad (13)$$

$$\partial_i u_i = 0, \quad (14)$$

where b is a buoyancy variable (proportional to the Eulerian entropy perturbation multiplied by the gravitational acceleration), e is a unit vector in the direction opposite to gravity, $N^2$ is the square of the buoyancy frequency (negative in the convectively unstable case) and $\kappa$ is the thermal diffusivity. (In the case of convection, no additional body force f is required to drive the motion.)

The standard shearing box, as employed by, e.g. Rogallo (1981), Wisdom & Tremaine (1988) and Hawley, Gammie & Balbus (1995), corresponds to setting $a(t) \propto t$ and, if appropriate, adding rotation to the box. In this case the shear is inexorable and the shearing coordinates must be remapped periodically for the purposes of numerical simulation.

In our case we consider $a(t)$ to oscillate sinusoidally and no remapping is required. As the frequency of the oscillation is varied, we can choose either to keep the amplitude of $a$ the same, which means that the maximum angle of the shear is fixed, or to scale the amplitude of $a$ inversely with the frequency, which means that the maximum angular velocity of the shear is fixed. For the most part, we are interested in the regime in which the shear is very small (in either sense), but if it is too small, then its effects cannot be measured reliably in a direct numerical simulation.

The quantity to be measured that is of greatest interest is the shear stress associated with the fluid motions, i.e. the Reynolds stress component

$$-R_{xy} = -(v_x v_y), \quad (15)$$

where the angular brackets denote a suitable averaging operation. (The physical Reynolds stress has an additional factor of the density of the fluid.) It is this stress that exchanges energy with the

large-scale shear. The energy equation in sheared coordinates in the absence of buoyancy forces is

$$\partial'_t (\frac{1}{2} v_i v_i) + \partial_i F_i = -\dot{a} v_x v_y - v |\nabla \times \mathbf{v}|^2 + f_i v_i, \quad (16)$$

where **F** is an appropriate energy flux density. If the Reynolds stress $-R_{xy}$ is positively correlated with the shear rate $\dot{a}$, then, like a viscous stress, it extracts energy from the large-scale shear and imparts it to smaller scale motions. However, there is no reason in principle why the correlation should not be negative, in which case work is done on the large-scale shear at the expense of either the body force or the decaying energy of the smaller scale flow.

# 3 HIGH-FREQUENCY RESPONSE

## 3.1 Linearized equations

In this section we present analytical results on the response of fluid motions to high-frequency shear of low amplitude. We omit buoyancy forces, which are probably inessential for this discussion, and consider the Navier-Stokes equations in sheared coordinates,

$$[\partial_t + v_j(\partial'_j - a\delta_{j1}\partial'_y)]v_i + \dot{a}v_x\delta_{i2} = -(\partial'_i - a\delta_{i1}\partial'_y)p + v(\partial'_j - a\delta_{j1}\partial'_y)(\partial'_j - a\delta_{j1}\partial'_y)v_i + f_i, \\ (\partial'_i - a\delta_{i1}\partial'_y)v_i = 0, \quad (17)$$

where $\delta_{ij}$ is the Kronecker delta and the subscripts 1 and 2 refer to the x and y directions involved in the shear. We now carry out a linearization in the shear amplitude. We consider a basic flow that exists in the absence of shear, satisfying the equations

$$(\partial_t + v_j \partial'_j)v_i = -\partial'_i p + \nu \Delta' v_i + f_i, \quad (18)$$

$$\partial'_i v_i = 0. \quad (19)$$

This flow might be freely decaying, or it might for example be a steady (or statistically steady) flow sustained by a body force. (Since we have omitted buoyancy forces, we cannot consider convection as such in this section, but we can study laminar or turbulent flows driven by a body force.) If we assume that $\partial'_i f_i = 0$, since only the solenoidal part of the force drives the flow, then $p$ satisfies the Poisson equation

$$\Delta' p = -(\partial'_i v_j) \partial'_j v_i. \quad (20)$$

The Laplacian operator $\Delta'$ has a unique inverse $\Delta'^{-1}$ if, as we assume, the boundary conditions are periodic and the fields have zero mean. This inverse is also a self-adjoint operator.

The linearized equations are

$$\begin{aligned} (\partial'_t + v_j \partial'_j)\delta v_i + (\delta v_j \partial'_j - a v_x \partial'_y)v_i + \dot{a} v_x \delta_{i2} \\ = -\partial'_i \delta p + a \delta_{i1} \partial'_y p + \nu (\Delta' \delta v_i - 2a \partial'_x \delta'_y v_i), \end{aligned} \\ \delta'_i \delta v_i - a \delta'_y v_x = 0, \quad (21)$$

where $\delta v_i$ is the velocity perturbation induced at first order by the shear and $\delta p$ is the accompanying pressure perturbation. In general the linearized equations must be solved numerically.

## 3.2 Asymptotic analysis for high frequencies

We now consider a limit in which the shear is oscillatory with high frequency. By this we mean that the time-scale of the oscillations is small compared to the circulation period of the streamlines of the basic flow and the viscous time-scale. For multiscale flows such as fully developed turbulence, the approximation considered
here applies when the oscillation period is short compared to the convective time-scale of the smallest eddies. Alternatively, it can be viewed as determining the response of those eddies for which the convective time-scale is long compared to the oscillation period.

We use the method of multiple scales and introduce a fast time variable $T' = t'/\epsilon$, where $\epsilon \ll 1$ is a small parameter that characterizes the ratio of time-scales. The rapidity of the shear is expressed by rewriting $a \mapsto a(T')$ and $\dot{a} \mapsto \epsilon^{-1}\dot{a}$, where the new meaning of $\dot{a}$ is $da/dT'$. The basic flow may vary with the ordinary time variable $t'$. We then expand

$$ \delta v_i = \delta v_{i0} + \epsilon \delta v_{i1} + \dots, \quad (22) $$

$$ \delta p = \epsilon^{-1}(\delta p_0 + \epsilon \delta p_1 + \dots), \quad (23) $$

in asymptotic series, where the quantities on the right-hand side depend on $(x', t', T')$. The $\delta p$ expansion is indexed in this way for reasons that will become clear. Since the basic flow does not depend on $\epsilon$, which is a property of the perturbations only, it is not expanded.

At leading order, we find

$$ \partial'_T \delta v_{i0} + \dot{a} v_x \delta_{i2} = -\partial'_i \delta p_0, \quad (24) $$

$$ \partial'_i \delta v_{i0} - a \partial'_y v_x = 0. \quad (25) $$

Roughly speaking, these equations describe an elastic response in which (from the second component of equation 24) $\delta v_{y0} \approx -a v_x$, leading to a shear stress $-\langle v_x \delta v_{y0} \rangle \approx a \langle v_x^2 \rangle$ proportional to the shear and in phase with it. The reality is somewhat more complicated because of the pressure terms and the constraint of incompressibility. The pressure perturbation satisfies a Poisson equation,

$$ \Delta' \delta p_0 = -2 \dot{a} \partial'_y v_x, \quad (26) $$

obtained by eliminating $\delta v$ from the above two equations, and so

$$ \partial'_T \delta v_{i0} = 2 \dot{a} \partial'_y \partial'^{-1} v_x - \dot{a} v_x \delta_{i2}. \quad (27) $$

The linearized shear stress at this order is

$$ -\delta R_{xy0} = -\langle v_x \delta v_{y0} + v_y \delta v_{x0} \rangle \quad (28) $$

and satisfies

$$ \partial'_T (-\delta R_{xy0}) = \dot{a} \langle v_x^2 - 2(v_x \partial'_y + v_y \partial'_x) \partial'_y \Delta'^{-1} v_x \rangle. \quad (29) $$

In terms of the tensor

$$ A_{ijkl} = \langle v_i \partial'_j \partial'_k \Delta'^{-1} v_l \rangle, \quad (30) $$

we have

$$ \partial'_T (-\delta R_{xy0}) = \dot{a} (A_{1jj1} - 2A_{1221} - 2A_{2121}). \quad (31) $$

This tensor has the symmetry property $A_{ijkl} = A_{ikjl}$, which follows immediately from its definition, and the contraction properties $A_{ijkj} = A_{ijkk} = 0$, which follow from $\partial'_i v_i = 0$. The further symmetry property $A_{ijkl} = A_{ljkj}$ follows from integration by parts, if the averaging operation includes a spatial average over a periodic cell. The quantity multiplying $\dot{a}$ on the right-hand side of equation (31) can be regarded as the effective elastic modulus (shear modulus) of the flow (divided by the density of the fluid).

At the next order, we find

$$ \begin{aligned} \partial'_T \delta v_{i1} + (\partial'_t + v_j \partial'_j) \delta v_{i0} &+ (\delta v_{j0} \partial'_j - a v_x \partial'_y) v_i \\ & = -\partial'_i \delta p_1 + a \delta_{i1} \partial'_y p + \nu (\Delta' \delta v_{i0} - 2 a \partial'_x \partial'_y v_i), \end{aligned} \quad (32) $$

$$ \partial'_i \delta v_{i1} = 0. \quad (32) $$

We then obtain the Poisson equation

$$ \begin{aligned} \Delta' \delta p_1 = & -a(\partial'_t + v_j \partial'_j) \partial'_y v_x - 2(\partial'_t v_j) \partial'_j \delta v_{i0} \\ & +a(\partial'_i v_x) \partial'_y v_i + a \partial'_x \partial'_y p + a \nu \Delta' \partial'_y v_x, \end{aligned} $$

which can be inverted in principle. The linearized shear stress at this order is

$$ -\delta R_{xy1} = -\langle v_x \delta v_{y1} + v_y \delta v_{x1} \rangle \quad (33) $$

and satisfies

$$ \begin{aligned} \partial'_T (-\delta R_{xy1}) = & \langle v_x (\partial'_t + v_j \partial'_j) \delta v_{y0} + v_y (\partial'_t + v_j \partial'_j) \delta v_{x0} \\ & + v_x (\delta v_{j0} \partial'_j - a v_x \partial'_y) v_y + v_y (\delta v_{j0} \partial'_j - a v_x \partial'_y) v_x \\ & + (v_x \partial'_y + v_y \partial'_x) \delta p_1 - a v_x \partial'_y p \\ & - v v_x (\Delta' \delta v_{y0} - 2 a \partial'_x \partial'_y v_y) \\ & - v v_y (\Delta' \delta v_{x0} - 2 a \partial'_x \partial'_y v_x) \rangle. \end{aligned} $$

We substitute for $p$ and $\delta p_1$ from the Poisson equations that they satisfy, using the fact that the inverse Laplacian is self-adjoint, and integrating by parts in various places:

$$ \begin{aligned} \partial'_T (-\delta R_{xy1}) = & \langle v_x (\partial'_t + v_j \partial'_j) \delta v_{y0} + v_y (\partial'_t + v_j \partial'_j) \delta v_{x0} \\ & + v_x (\delta v_{j0} \partial'_j - a v_x \partial'_y) v_y + v_y (\delta v_{j0} \partial'_j - a v_x \partial'_y) v_x \\ & - \Delta'^{-1} (\partial'_y v_x - \partial'_x v_y) [-a (\partial'_t + v_j \partial'_j) \partial'_y v_x \\ & - 2 (\partial'_j v_j) \partial'_j \delta v_{i0} + a (\partial'_x v_x) \partial'_y v_i + a \partial'_x \partial'_y p + a \nu \Delta' \partial'_y v_x] \\ & - a (\partial'_j v_j) (\partial'_j v_i) \partial'_y \Delta'^{-1} v_y - v v_x (\Delta' \delta v_{y0} - 2 a \partial'_x \partial'_y v_y) \\ & - v v_y (\Delta' \delta v_{x0} - 2 a \partial'_x \partial'_y v_x) \rangle. \end{aligned} $$

A further replacement of $p$ is required:

$$ \begin{aligned} \partial'_T (-\delta R_{xy1}) = & \langle v_x (\partial'_t + v_j \partial'_j) \delta v_{y0} + v_y (\partial'_t + v_j \partial'_j) \delta v_{x0} \\ & + v_x (\delta v_{j0} \partial'_j - a v_x \partial'_y) v_y + v_y (\delta v_{j0} \partial'_j - a v_x \partial'_y) v_x \\ & - \Delta'^{-1} (\partial'_y v_x + \partial'_x v_y) [-a (\partial'_t + v_j \partial'_j) \partial'_y v_x \\ & - 2 (\partial'_j v_j) \partial'_j \delta v_{i0} + a (\partial'_x v_x) \partial'_y v_i + a \nu \Delta' \partial'_y v_x] \\ & + a \Delta'^{-2} \partial'_x \partial'_y (\partial'_y v_x + \partial'_x v_y) (\partial'_j v_j) (\partial'_j v_i) \\ & - a (\partial'_j v_j) (\partial'_j v_i) \partial'_y \Delta'^{-1} v_y - v v_x (\Delta' \delta v_{y0} - 2 a \partial'_x \partial'_y v_y) \\ & - v v_y (\Delta' \delta v_{x0} - 2 a \partial'_x \partial'_y v_x) \rangle. \end{aligned} $$

Since we have an expression for $\partial'_T \delta v_{i0}$, we consider

$$ \begin{aligned} \partial_T'^2 (-\delta R_{xy1}) = & \langle v_x (\partial'_t + v_j \partial'_j) \partial_T' \delta v_{y0} + v_y (\partial'_t + v_j \partial'_j) \partial_T' \delta v_{x0} \\ & + v_x [(\partial_T' \delta v_{j0}) \partial'_j - \dot{a} v_x \partial'_y] v_y + v_y [(\partial_T' \delta v_{j0}) \partial'_j - \dot{a} v_x \partial'_y] v_x \\ & - \Delta'^{-1} (\partial'_y v_x + \partial'_x v_y) [-\dot{a} (\partial'_t + v_j \partial'_j) \partial'_y v_x \\ & - 2 (\partial'_i v_j) \partial'_j \partial'_T \delta v_{i0} + \dot{a} (\partial'_i v_x) \partial'_y v_i + \dot{a} \nu \Delta' \partial'_y v_x] \\ & + \dot{a} \Delta'^{-2} \partial'_x \partial'_y (\partial'_y v_x + \partial'_x v_y) (\partial'_j v_j) (\partial'_j v_i) \\ & - \dot{a} (\partial'_j v_j) (\partial'_j v_i) \partial'_y \Delta'^{-1} v_y - v v_x (\Delta'_T \delta v_{y0} - 2 \dot{a} \partial'_x \partial'_y v_y) \\ & - v v_y (\Delta'_T \delta v_{x0} - 2 \dot{a} \partial'_x \partial'_y v_x) \rangle. \end{aligned} $$

Substituting for $\partial'_T \delta v_{i0}$, we obtain

$$ \begin{aligned} \partial_T'^2 (-\delta R_{xy1}) = & \dot{a} \langle v_x (\partial'_t + v_j \partial'_j) (2 \partial'_y \partial'_y \Delta'^{-1} v_x - v_x) \\ & + v_y (\partial'_t + v_j \partial'_j) (2 \partial'_x \partial'_y \Delta'^{-1} v_x) \\ & + v_x [(2 \partial'_j \partial'_y \Delta'^{-1} v_x) \partial'_j - 2 v_x \partial'_y] v_y \\ & + v_y [(2 \partial'_j \partial'_y \Delta'^{-1} v_x) \partial'_j - 2 v_x \partial'_y] v_x \\ & - \Delta'^{-1} (\partial'_y v_x + \partial'_x v_y) [- (\partial'_t + v_j \partial'_j) \partial'_y v_x \\ & - 2 (\partial'_j v_j) \partial'_j (2 \partial'_y \partial'_y \Delta'^{-1} v_x - v_x \delta_{i2}) + (\partial'_i v_x) \partial'_y v_i \\ & + \nu \Delta' \partial'_y v_x] + \Delta'^{-2} \partial'_x \partial'_y (\partial'_y v_x + \partial'_x v_y) (\partial'_j v_j) (\partial'_j v_i) \\ & - (\partial'_j v_j) (\partial'_j v_i) \partial'_y \Delta'^{-1} v_y \\ & - v v_x (2 \partial'_y \partial'_y v_x - \Delta' v_x - 2 \partial'_x \partial'_y v_y) \rangle. \end{aligned} $$
This can be rearranged as follows, after integration by parts:

$$
\begin{align*}
\partial_t^2 (-\delta R_{xy1}) &= \dot{a}((\partial'_t v_x)(-v_x + \partial'_y \partial'_y \Delta'^{-1} v_x + \partial'_x \partial'_y \Delta'^{-1} v_y) \\
&\quad + v_x(v_j \partial'_j)(2\partial'_y \partial'_y \Delta'^{-1} v_x - v_x) \\
&\quad + v_y(v_j \partial'_j)(2\partial'_x \partial'_y \Delta'^{-1} v_x) \\
&\quad + v_x[(2\partial'_j \partial'_y \Delta'^{-1} v_x)\partial'_j - 2v_x \partial'_y]v_y \\
&\quad + v_y[(2\partial'_y \partial'_y \Delta'^{-1} v_x)\partial'_j - 2v_x \partial'_y]v_x \\
&\quad - \Delta'^{-1}(\partial'_y v_x + \partial'_x v_y)[-(v_j \partial'_j)\partial'_y v_x \\
&\quad - 2(\partial'_i v_j)\partial'_j(2\partial'_y \partial'_y \Delta'^{-1} v_x) + 3(\partial'_y v_j)\partial'_j v_x] \\
&\quad + \Delta'^{-2}\partial'_x \partial'_y(\partial'_y v_x + \partial'_x v_y)(\partial'_i v_j)(\partial'_j v_i) \\
&\quad - (\partial'_i v_j)(\partial'_j v_i)\partial'_y \Delta'^{-1} v_y \\
&\quad - v v_x(\partial'_y \partial'_y v_x - \Delta'^{-1} v_x - 3\partial'_x \partial'_y v_y).
\end{align*}
$$

We define the tensors

$$
B_{ijkl} = \langle (\partial'_i v_l) \partial'_j \partial'_k \Delta'^{-1} v_l \rangle, \quad (34)
$$

$$
C_{ijkl} = -\nu \langle v_i \partial'_j \partial'_k v_l \rangle, \quad (35)
$$

$$
D_{ijkl} = \langle v_i v_j \partial'_k v_l \rangle, \quad (36)
$$

$$
D_{ijklmn} = \langle v_i v_j \partial'_k \partial'_l \partial'_m \Delta'^{-1} v_n \rangle, \quad (37)
$$

$$
D_{ijklmnpq} = \langle v_i v_j \partial'_k \partial'_l \partial'_m \partial'_n \partial'_p \Delta'^{-2} v_q \rangle, \quad (38)
$$

$$
E_{ijkl} = \langle v_m (\partial'_m \partial'_n \Delta'^{-1} \partial'_i v_j) \partial'_n \Delta'^{-1} \partial'_k v_l \rangle, \quad (39)
$$

which satisfy $B_{ijkl} = B_{ikjl}$, $B_{ijkj} = 0$, $C_{ijkl} = C_{ikjl} = C_{lji}$, $C_{ijkj} = 0$ and various other identities for D and E. We then find (after integration by parts)

$$
\begin{align*}
\partial_t^2 (-\delta R_{xy1}) &= -\dot{a}(B_{1jj1} - B_{1221} - B_{1122} \\
&\quad -C_{1221} + C_{1jj1} + 3C_{1122} \\
&\quad -2D_{1jj221} - 2D_{2jj121} - 3D_{1jj221} - 3D_{1jj212} \\
&\quad -D_{ijj1221} - D_{ijj1212} + D_{ijj22} \\
&\quad +4E_{2121} + 4E_{2112}).
\end{align*}
$$

3.3 Interpretation

The two results obtained so far can be written in the forms

$$
\partial_t' (-\delta R_{xy0}) = \dot{a} G_0, \quad (40)
$$

$$
\partial_t^2 (-\delta R_{xy1}) = -\dot{a} G_1, \quad (41)
$$

where $G_0$ and $G_1$ are two coefficients, each of which could be
positive, negative or zero. Given that this is a linear analysis, we may
consider a complex shear $a \propto \exp(-i\omega t)$ with angular frequency
$\omega = O(\epsilon^{-1})$ and deduce that the shear stress in the limit of high
frequency is

$$
-\delta R_{xy} = a \left[ \mathcal{G}_0 - \frac{i\mathcal{G}_1}{\omega} + O(\epsilon^2) \right]. \quad (42)
$$

The term $\mathcal{G}_0$ represents an ideal elastic response, while the term $\mathcal{G}_1$
represents an imperfection of the elasticity associated with dis-
sipation. For comparison, the shear stress of a viscous fluid is
$v\dot{a} = -i\omega va$. Therefore the effective kinematic viscosity of the
flow at high frequencies is $\mathcal{G}_1/\omega^2$.

The calculation of Goodman & Oh (1997) corresponds only to
the leading order of the above expansion. They do not mention the
elastic stress but focus on the dissipation rate at leading order in
a periodic strain, given by their equation (25). As they point out,
this expression evaluates to zero; this is consistent with our analysis
because there is no dissipation associated with a perfect elastic
stress. They obtain a non-zero result at this order by manipulating
the apparent singularity in their expression at zero frequency. While
this procedure may produce a meaningful result, it is not really
justified because the preceding steps have assumed a separation of
time-scales between the tide and the convection; the zero-frequency
pole signals the breakdown of their approximation scheme and its
resolution is not straightforward. In contrast, our calculation of G₁
and the effective viscosity in the high-frequency limit is based on a
systematic asymptotic expansion.

**3.4 Evaluation for statistically isotropic flows**

Although convective flows are naturally anisotropic, analytical
progress is easiest when the flow is assumed to have isotropic prop-
erties when spatially averaged. In this case A<sub>ijkl</sub> must be an isotropic
tensor, i.e.

$$
A_{ijkl} = A_1 \delta_{ij} \delta_{kl} + A_2 \delta_{ik} \delta_{jl} + A_3 \delta_{il} \delta_{jk}. \quad (43)
$$

The symmetry and contraction properties then imply $A_2 = A_1$ and $A_3 = -(d+1)A_1$, where $d$ is the number of spatial dimensions (of course $d=3$, but the case $d=2$ is also of at least theoretical interest), and so

$$
A_{ijkl} = \frac{2K}{d(d-1)(d+2)} [(d+1)\delta_{il}\delta_{jk} - \delta_{ij}\delta_{kl} - \delta_{ik}\delta_{jl}], \quad (44)
$$

where the mean kinetic energy density K is given by

$$
2K = \langle v_i v_i \rangle = A_{ijji} = -d(d-1)(d+2)A_1. \quad (45)
$$

In this isotropic case we then obtain

$$
G_0 = \frac{2(d-2)(d+1)}{d(d-1)(d+2)} K. \quad (46)
$$

The effective elastic modulus is positive for a 3D flow but vanishes
in two dimensions under the assumption of isotropy.

In the isotropic case it can be shown (see Appendix B) that
the triple-correlation tensor *D* vanishes identically. *E* also vanishes
when the identity *E*<sub>ijkl</sub> = −*E*<sub>klij</sub> (which follows from integration by
parts) is combined with the general form of a fourth-rank isotropic
tensor. The tensors *B* and *C* have the form

$$
B_{ijkl} = B[(d+1)\delta_{il}\delta_{jk} - \delta_{ij}\delta_{kl} - \delta_{ik}\delta_{jl}], \quad (47)
$$

$$
C_{ijkl} = C[(d+1)\delta_{il}\delta_{jk} - \delta_{ij}\delta_{kl} - \delta_{ik}\delta_{jl}]. \quad (48)
$$

Furthermore, the energy equation of the basic flow,

$$
\langle (\partial'_i v_i) v_i \rangle = v \langle v_i \Delta' v_i \rangle + \langle f_i v_i \rangle, \quad (49)
$$

implies

$$
B_{ijji} + C_{ijji} = \langle f_i v_i \rangle = P, \quad (50)
$$

the power input per unit volume (or area), and so

$$
d(d-1)(d+2)(B+C) = P \quad (51)
$$

for isotropic statistics. In this case,

$$
G_1 = B(d^2 - 2) + C(d^2 - 6). \quad (52)
$$

This corresponds to a viscous response, with the effective
viscosity being inversely proportional to frequency-squared, for
high-frequency oscillatory shear. The effective viscosity coefficient is $G_1/\omega^2$. In the case $B=0$, when the flow is maintained steadily against viscous dissipation by the body force, the effective viscosity is positive in three dimensions and negative in two dimensions. In the case $P=0$, when the flow is freely decaying ($-B = C > 0$), the effective viscosity is negative in either three or two dimensions.

## 3.5 Evaluation for ABC flows

A widely studied class of incompressible fluid flows in a periodic domain is provided by the ABC flow (named after Arnol'd, Beltrami and Childress; Arnol'd 1965; Galloway & Frisch 1987)

$$v = \begin{pmatrix} A \sin kz' + C \cos ky' \\ B \sin kx' + A \cos kz' \\ C \sin ky' + B \cos kx' \end{pmatrix} \quad (53)$$

in a cube of length $2\pi/k$. This velocity field has the Beltrami property $\nabla' \times v = kv$, so non-linearity is absent; $v \cdot \nabla'v$ is balanced by a pressure gradient. If the flow is unforced, $A, B$ and $C$ decay proportionally to $\exp(-\nu k^2 t)$. Alternatively the flow can be maintained against dissipation by supplying a body force $f = \nu k^2 v$. The most widely studied example has $A = B = C$.

The response coefficients are easily evaluated as

$$G_{0\frac{1}{2}}(A^2 - C^2), \quad (54)$$

$$G_1 = \frac{1}{2} A (\dot{A} + \nu k^2 A). \quad (55)$$

Therefore the elasticity can be positive, negative or zero depending on the anisotropy of the flow and its orientation relative to the shear. The effective viscosity $G_1/\omega^2$ at leading order vanishes for a freely decaying flow but is positive (and inversely proportional to frequency-squared) for a forced flow, assuming that $A \neq 0$.

In fact, probably all analytical examples of 3D fluid flows lack genuine non-linearity, in the sense that, if they are expanded in a Fourier basis with wavenumbers $k$, there are no non-empty triads of interacting components satisfying $k_1 + k_2 + k_3 = 0$. Such triads tend to produce a cascade of energy to larger wavenumbers in the manner of hydrodynamic turbulence. If there are no non-empty triads, the triple-correlation tensors $D$ and $E$ again vanish, and the only contributions to $G_1$ come from the time dependence of the flow (the tensor $B$) and the viscous terms (the tensor $C$). Both of these effects may be small if the viscosity is small.

The ABC flow is stable only for sufficiently small Reynolds number. More realistically, in a typical flow at high Reynolds number there will be a turbulent cascade involving strong triad interactions, meaning that the tensors $D$ and $E$ may be significant (although apparently not in isotropic turbulence). The tensors $B$ and $C$ will also be enhanced by the turbulent cascade. Numerical simulations are required to access this regime.

We note that the 'eddy viscosity' of the $A = B = C$ flow has been calculated, as a function of Reynolds number, by Wirth, Gama & Frisch (1995); see also references therein for related studies. However, the nature of their calculation is different; using multiscale techniques, they determine the behaviour of very slow large-scale deformations of the cellular flow, and deduce that a large-scale instability occurs through the appearance of a negative eddy viscosity. In contrast, our analysis determines the response of a flow to an imposed large-scale deformation of high frequency.

# 4 DIRECT NUMERICAL SIMULATIONS OF CONVECTION IN AN OSCILLATORY SHEARING BOX

## 4.1 Numerical set-up

We have implemented the oscillatory shearing box in the snoopy code, a 3D spectral code solving the equations of incompressible or Boussinesq (magneto)hydrodynamics using a Fourier representation (Lesur & Longaretti 2005, 2007; Lesur & Ogilvie 2010). We assume that an external force creates an oscillatory shear with $\dot{a} = -S \cos(\Omega t)$, where $S$ is the maximum shear rate and $\Omega$ is the tidal frequency. As in previous sections, the shearing motion is in the $y$ direction, with a linear dependence on $x$, i.e. $-S \cos(\Omega t)x e_y$, and the force that drives it is $\Omega S \sin(\Omega t)x e_y$.

In addition to this imposed shearing motion, convection is driven by applying uniform gravity and an unstable entropy gradient, within the Boussinesq approximation. (No additional body force is required to drive the motion.) By setting $N^2 = -1$, we adopt a unit of time related to the unstable stratification. We test two different configurations: convection in the shearwise direction ($e = e_x$) and convection in the spanwise direction ($e = e_z$, where $e$ is the direction of stratification defined in Section 2).

The aspect ratio of the periodic box is adjusted according to the direction of stratification, with $L_x \times L_y \times L_z = 2 \times 1 \times 1$ when $e = e_x$ and $1 \times 1 \times 2$ when $e = e_z$. Having a box elongated in the direction of stratification allows 'elevator modes' to break up more easily through secondary instabilities; otherwise, these tend to dominate the convection at moderate Rayleigh numbers (Lesur & Ogilvie 2010). The resolution is 128 collocation points per unit of length.

To control dissipation on small scales, we introduce an explicit kinematic viscosity $\nu$ and thermal diffusivity $\kappa$, with Prandtl number $\nu/\kappa = 1$ for simplicity. The value of the diffusion coefficients is set according to the Rayleigh number $Ra = |N^2|L_y^4/\nu\kappa$. In our set-up, convection starts when $Ra > (2\pi)^4 \approx 1559$. In the following we will consider simulations exhibiting fully turbulent convection, with $Ra = 4 \times 10^6$.

In order to avoid any artefact of the initial conditions, we initiate our simulations with noise at the largest scales, and we let turbulence evolve without any shear for 100 turnover times (i.e. from $t = -100|N|^{-1}$ to 0). The spectrum of the turbulence we obtain and a typical snapshot are shown in Fig. 1. Once a quasi-stationary turbulent state is reached, we switch on a weak oscillatory shear and start measuring the Reynolds stress $-R_{xy}(t) = -(v_x v_y)$, where $(\cdot)$ denotes a volume average over the box. Such a simulation has to be continued for an integration time $\Delta T$ of hundreds of turnover times in order to reduce the impact of turbulent noise on the measurements.

## 4.2 Turbulent viscosity: definitions and simple models

Traditionally, the turbulent viscosity is associated with a simple closure formula relating the Reynolds stress to the rate of strain; in our system,

$$-R_{xy} = \nu_t \dot{a}. \quad (56)$$

This expression is based on the assumptions that the relationship between stress and rate of strain is linear and instantaneous. In the case of rapid oscillatory shear, the assumption of linearity is reasonable provided that the shear is of sufficiently small amplitude. However, the assumption of instantaneity is generally not valid and
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Figure 1. Simulation snapshot of thermal fluctuations (top) and energy spectra (bottom) at Ra = 4 × 10⁶ in the case without any shear. We show the kinetic energy spectrum (E_k) and the temperature fluctuation spectrum (b_k).</figcaption></figure>

we should allow the Reynolds stress to be linearly related to $\dot{a}(t)$ through an integral over its previous history. In the Fourier domain, this relationship (a convolution) reduces to a multiplication:

$$-\tilde{R}_{xy}(\omega) = v_t(\omega)\tilde{a}(\omega), \quad (57)$$

where $\tilde{\cdot}$ denotes the Fourier transform in time. In this expression, $v_t(\omega)$ is a complex function of frequency, rather than a real constant, indicating that the Reynolds stress and the rate of strain are generally out of phase, and that the relationship is frequency-dependent. Since the relationship involves real-valued functions in the temporal domain, $v_t$ has the Hermitian symmetry $v_t(-\omega) = [v_t(\omega)]^*$ for real $\omega$. With our choice of $\dot{a} = -S \cos(\Omega t)$, we have

$$\tilde{R}_{xy}(\omega) = v_t(\omega)S\pi[\delta(\omega - \Omega) + \delta(\omega + \Omega)]. \quad (58)$$

This makes it possible to measure $v_t(\Omega)$ from the time series $R_{xy}(t)$. In reality, the delta function is replaced by a peak of finite height and non-zero width because of the finite integration time; alternatively, it becomes a Kronecker delta in a discrete Fourier transform. Furthermore, noise is present because there are turbulent fluctuations

in $R_{xy}$ even in the absence of shear. The integration time must be sufficiently long to allow the signal to be detected above the noise.

A simple closure model for the Reynolds stress that can be compared with the results of numerical simulations is a viscoelastic model, which reflects the fact that the turbulent stress cannot instantaneously produce a viscous response to a time-dependent shear rate. [Viscoelastic models for magnetohydrodynamic turbulence in astrophysical discs have been discussed by Ogilvie (2001).] We suppose that the convective turbulence contains a number of viscoelastic components, having effective elastic moduli $c_j$ and relaxation times $1/\gamma_j$, giving rise to the complex turbulent viscosity

$$v_t(\omega) = \sum_j \frac{c_j}{\gamma_j - i\omega}. \quad (59)$$

In this expression, $c_j$ and $\gamma_j$ are real, but we do not require $c_j > 0$, allowing for the possibility that negative Re[$v_t(\omega)$] may appear for some ranges of frequency. For example, the case of constant shear rate ($\omega = 0$) is known to produce negative viscosity for low values of Ra in protoplanetary disc convection, which involves rotation as well as shear (Lesur & Ogilvie 2010). However, in order to have a physical relaxation for all possible input signals, we require $\gamma_j > 0$, indicating that the system we describe is stable for any time-dependent shear. In other words, since $v_t(\omega)$ in equation (57) derives from a causal integral relationship, it should be analytic in the upper half-plane.

This model shares several properties with the analytical results derived in Section 3. In the high-frequency limit $\omega \to \infty$, the model gives

$$v_t = i\omega^{-1} \sum_j c_j + \omega^{-2} \sum_j c_j \gamma_j + O(\omega^{-3}). \quad (60)$$

This behaviour is equivalent to equation (42) obtained for the high-frequency response of an arbitrary flow, if $G_0 = \sum_j c_j$ (the high-frequency elastic modulus) and $G_1 = \sum_j c_j \gamma_j$. (On the other hand, the viscosity obtained in the low-frequency limit $\omega \to 0$ is $\sum_j c_j / \gamma_j$.) As we will see later, the same type of dependency is found in numerical simulations.

In this paper, we will assume that only two viscoelastic components are present. This can be seen as a computational limitation, since numerical simulations cannot probe very high frequencies which are usually associated with scales below the grid size. In the following, we will therefore consider the simplified model

$$v_t(\omega) = \frac{c_1}{\gamma_1 - i\omega} + \frac{c_2}{\gamma_2 - i\omega}. \quad (61)$$

## 4.3 Measuring the turbulent viscosity

A typical example of the time series of the Reynolds stress and its Fourier transform is shown in Fig. 2. In order to reduce the aliasing of low frequencies into the high-frequency domain due to the finite integration time, we have applied a Hanning window to the time series before computing the Fourier transform. As is evident from the time series, the turbulent convection produces large fluctuations in the Reynolds stress which dominate the response to the oscillatory shear. Nevertheless, it is possible to extract useful information by looking at the temporal spectrum of the stress. In particular, a spike is clearly visible at $\omega = \Omega$, which indicates that the turbulent flow is producing a detectable response to our forcing. With the complete Fourier transform of the Reynolds stress, it is possible to measure
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time history of the Reynolds stress (left) and Fourier transform of the Reynolds stress (right) with an oscillatory shear of frequency $\Omega = 16$ and amplitude $S = 0.1$. The periodic signal cannot be seen in the time series. The Fourier transform exhibits a peak at the excitation frequency $\omega = \Omega$.</figcaption></figure>

the turbulent viscosity at the forcing frequency using equation (58). Moreover, evaluating the noise level in the vicinity of the spike allows us to estimate the uncertainty in the measurement of the turbulent viscosity.

It should be noted that the signal-to-noise ratio (defined as the ratio of the amplitude of the spectral spike to the amplitude of the surrounding noise; see Fig. 2, right-hand panel) should decay as $\Delta T^{-1/2}$. Therefore, if the oscillatory signal (the spectral spike) is too weak to be detected in the Fourier transform, one can in principle increase the integration time to reduce the signal-to-noise ratio, but this is a computationally expensive procedure.

## 4.4 Numerical results

The results of the simulations performed in this work are shown in Table 1. They are also plotted in Fig. 3 for shearwise convection and in Fig. 4 for spanwise convection. Despite the presence of significant noise in some of the results, several deductions can be made from these simulations.

(i) At sufficiently high frequency, we find $Re(v_t) < 0$ in every case. Note, however, that noise dominates many high-frequency measurements of this quantity. To substantiate the high-frequency behaviour, we have performed a simulation at lower Rayleigh

<table><caption>**Table 1.** Parameters and results of the simulations discussed in this paper. The error in $v_t$ is estimated based on the rms turbulent noise in the Fourier transform at the oscillation frequency. $K$ is the mean kinetic energy density of the fluctuations.</caption><thead><tr><th>Ra</th><th>e direction</th><th>ΔT</th><th>Ω</th><th>S</th><th>K</th><th>Re(v<sub>t</sub>)</th><th>Im(v<sub>t</sub>)</th><th>Error</th></tr></thead><tbody><tr><td>4 × 10<sup>6</sup></td><td>x</td><td>500</td><td>0.25</td><td>0.1</td><td>2.1 × 10<sup>-1</sup></td><td>7.9 × 10<sup>-2</sup></td><td>2.7 × 10<sup>-2</sup></td><td>1.2 × 10<sup>-2</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>x</td><td>300</td><td>0.5</td><td>0.1</td><td>2.2 × 10<sup>-1</sup></td><td>6.8 × 10<sup>-2</sup></td><td>3.4 × 10<sup>-2</sup></td><td>3.0 × 10<sup>-2</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>x</td><td>300</td><td>1.0</td><td>0.1</td><td>2.1 × 10<sup>-1</sup></td><td>4.5 × 10<sup>-2</sup></td><td>4.1 × 10<sup>-2</sup></td><td>2.2 × 10<sup>-2</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>x</td><td>300</td><td>2.0</td><td>0.1</td><td>2.2 × 10<sup>-1</sup></td><td>2.1 × 10<sup>-2</sup></td><td>3.2 × 10<sup>-2</sup></td><td>1.3 × 10<sup>-2</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>x</td><td>300</td><td>4.0</td><td>0.1</td><td>2.1 × 10<sup>-1</sup></td><td>3.1 × 10<sup>-3</sup></td><td>1.6 × 10<sup>-2</sup></td><td>5.2 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>x</td><td>300</td><td>8.0</td><td>0.1</td><td>2.2 × 10<sup>-1</sup></td><td>-1.2 × 10<sup>-4</sup></td><td>8.6 × 10<sup>-3</sup></td><td>2.1 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>x</td><td>300</td><td>16.0</td><td>0.1</td><td>2.2 × 10<sup>-1</sup></td><td>-5.6 × 10<sup>-5</sup></td><td>3.4 × 10<sup>-3</sup></td><td>3.6 × 10<sup>-4</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>x</td><td>300</td><td>32.0</td><td>0.1</td><td>2.2 × 10<sup>-1</sup></td><td>-7.2 × 10<sup>-5</sup></td><td>1.9 × 10<sup>-3</sup></td><td>4.5 × 10<sup>-5</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>x</td><td>300</td><td>64.0</td><td>0.1</td><td>2.2 × 10<sup>-1</sup></td><td>-1.7 × 10<sup>-5</sup></td><td>9.7 × 10<sup>-4</sup></td><td>4.0 × 10<sup>-6</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>300</td><td>1.0</td><td>0.1</td><td>2.3 × 10<sup>-1</sup></td><td>1.7 × 10<sup>-2</sup></td><td>1.6 × 10<sup>-2</sup></td><td>5.6 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>300</td><td>2.0</td><td>0.1</td><td>2.1 × 10<sup>-1</sup></td><td>1.4 × 10<sup>-2</sup></td><td>1.3 × 10<sup>-2</sup></td><td>4.6 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>300</td><td>4.0</td><td>0.1</td><td>2.2 × 10<sup>-1</sup></td><td>6.3 × 10<sup>-3</sup></td><td>1.3 × 10<sup>-2</sup></td><td>3.4 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>300</td><td>8.0</td><td>0.1</td><td>2.2 × 10<sup>-1</sup></td><td>-6.5 × 10<sup>-4</sup></td><td>5.7 × 10<sup>-3</sup></td><td>1.4 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>300</td><td>16.0</td><td>0.1</td><td>2.2 × 10<sup>-1</sup></td><td>-1.6 × 10<sup>-4</sup></td><td>2.8 × 10<sup>-3</sup></td><td>4.0 × 10<sup>-4</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>700</td><td>0.25</td><td>0.05</td><td>2.2 × 10<sup>-1</sup></td><td>2.3 × 10<sup>-2</sup></td><td>8.8 × 10<sup>-3</sup></td><td>8.9 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>700</td><td>0.5</td><td>0.05</td><td>2.2 × 10<sup>-1</sup></td><td>2.1 × 10<sup>-2</sup></td><td>6.8 × 10<sup>-3</sup></td><td>9.4 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>600</td><td>1.0</td><td>0.05</td><td>2.2 × 10<sup>-1</sup></td><td>1.0 × 10<sup>-2</sup></td><td>1.9 × 10<sup>-2</sup></td><td>8.4 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>600</td><td>2.0</td><td>0.05</td><td>2.2 × 10<sup>-1</sup></td><td>8.9 × 10<sup>-3</sup></td><td>1.3 × 10<sup>-2</sup></td><td>7.2 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>600</td><td>4.0</td><td>0.05</td><td>2.2 × 10<sup>-1</sup></td><td>4.2 × 10<sup>-3</sup></td><td>9.6 × 10<sup>-3</sup></td><td>4.0 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>300</td><td>8.0</td><td>0.05</td><td>2.2 × 10<sup>-1</sup></td><td>1.3 × 10<sup>-3</sup></td><td>5.8 × 10<sup>-3</sup></td><td>2.9 × 10<sup>-3</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>300</td><td>16.0</td><td>0.05</td><td>2.1 × 10<sup>-1</sup></td><td>-1.7 × 10<sup>-4</sup></td><td>2.4 × 10<sup>-3</sup></td><td>6.7 × 10<sup>-4</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>300</td><td>32.0</td><td>0.05</td><td>2.2 × 10<sup>-1</sup></td><td>-8.0 × 10<sup>-5</sup></td><td>1.3 × 10<sup>-3</sup></td><td>9.0 × 10<sup>-5</sup></td></tr><tr><td>4 × 10<sup>6</sup></td><td>z</td><td>300</td><td>64.0</td><td>0.05</td><td>2.2 × 10<sup>-1</sup></td><td>-4.0 × 10<sup>-6</sup></td><td>6.2 × 10<sup>-4</sup></td><td>6.5 × 10<sup>-6</sup></td></tr><tr><td>1 × 10<sup>5</sup></td><td>z</td><td>3600</td><td>64.0</td><td>0.1</td><td>2.6 × 10<sup>-1</sup></td><td>-5.7 × 10<sup>-6</sup></td><td>7.0 × 10<sup>-4</sup></td><td>3.9 × 10<sup>-7</sup></td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Turbulent viscosity of shearwise convection versus the angular oscillation frequency in units of $|N|$. The numerical measurements are shown in blue (open circles) and the best fit of the closure model (61) is shown in green (triangles). Negative values are connected by a dashed line. Uncertainties due to turbulent noise are shown as a shaded region.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Turbulent viscosity of spanwise convection, versus the angular oscillation frequency in units of $|N|$. The numerical measurements are shown in blue (open circles; $S = 0.1$) and red (crosses; $S = 0.05$) and the best fit of the closure model (61) is shown in green (triangles). Negative values are connected by a dashed line. Uncertainties due to turbulent noise are shown as a shaded region.</figcaption></figure>

number ($Ra = 10^5$) for more than 3000 turnover times in order to reduce the effect of turbulent noise. This simulation exhibits clearly $Re(v_t) < 0$ at $\Omega = 64$, indicating that the trend observed in the other configurations is real.

(ii) $Im(v_t) > 0$. This conclusion is rather strong since the noise level is smaller than the measured value of $Im(v_t)$ in most cases. It means that the effective elasticity is positive.

(iii) Larger values of $|v_t|$ are found when the stratification is in the $x$ direction (shearwise convection). This suggests that the full turbulent viscosity tensor is anisotropic.

(iv) Changing $S$ does not change significantly the measured turbulent viscosity, suggesting that the numerical simulations are in the linear regime assumed in writing down equation (57). However, the results obtained with the smaller value of $S$ are subject to significant uncertainty.

(v) The asymptotic behaviour in the high-frequency limit is $Re(v_t) \propto \omega^{-2}$ and $Im(v_t) \propto \omega^{-1}$.

We have also fitted the simple closure model (61) to the numerical results for shearwise and spanwise convection. The best fits are shown as green curves in Figs 3 and 4, and the coefficients we obtained are shown in Table 2. The major viscoelastic component has an effective elastic modulus $c_1$ that is comparable to, but somewhat less than, the mean kinetic energy density $K$ (cf. equation 46, which

<table><caption>**Table 2.** Best-fitting coefficients for the closure model (61).</caption><thead><tr><td>Coefficient</td><td>Shearwise convection</td><td>Spanwise convection</td></tr></thead><tbody><tr><td>c<sub>1</sub></td><td>0.0978</td><td>0.06</td></tr><tr><td>γ<sub>1</sub></td><td>1.07</td><td>1.62</td></tr><tr><td>c<sub>2</sub></td><td>-0.0286</td><td>-0.0178</td></tr><tr><td>γ<sub>2</sub></td><td>6.15</td><td>24.25</td></tr></tbody></table>

suggests a ratio of 4/15 in the isotropic case) and a relaxation rate comparable to the nominal convective turnover rate $|N|$, which is 1 in our system of units. Note that the contribution to $Re(v_t)$ from the first component is $(c_1/\gamma_1)[1 + (\omega/\gamma_1)^2]^{-1}$, which is compatible with the form of the frequency-dependent viscosity used by Ogilvie & Lin (2007). The suggested presence of a second viscoelastic component with a negative $c_2$ and a faster relaxation indicates that, at high frequency, the real part of the viscosity changes sign, as is observed in the numerical results. With more data we might be led to introduce a broader spectrum of viscoelastic components.

In order to identify which spatial scales contribute to the turbulent viscosity, we have computed the spatial 'spectrum' of the turbulent viscosity. This is derived from the instantaneous spatial spectrum of the Reynolds stress,

$$ \hat{R}_{xy}(k, t) = 2 \operatorname{Re} [\hat{v}_x \hat{v}_y^*], \qquad (62) $$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Spatial spectrum of the turbulent viscosity for the simulation of spanwise convection with $\Omega = 16$ and $S = 0.1$. The real part has two components: a negative contribution from large scales and a positive contribution from smaller scale ($k/2\pi \sim 10$).</figcaption></figure>

where $\hat{v}_i$ is the 3D instantaneous spatial Fourier transform of $v_i$ and the overbar denotes a shell-integration procedure in spectral space, such that $R_{xy}(t) = \int_0^\infty \hat{R}_{xy}(k, t) dk$. $\hat{R}_{xy}(k, t)$ is then used in a relation similar to (58) which defines the turbulent viscosity spectrum $\hat{v}_i(k, \omega)$:

$$ \tilde{R}_{xy}(k, \omega) = \hat{v}_i(k, \omega) S \pi [\delta(\omega - \Omega) + \delta(\omega + \Omega)]. \quad (63) $$

We have applied this procedure to the simulation with $\Omega = 16$, $S = 0.1$ and with stratification in the $z$ direction. The corresponding spectra are shown in Fig. 5. As can be seen, much of the turbulent viscosity (both real and imaginary) is due to the largest scales of the system. However, we also observe a component with Re($v_t$) > 0 at 'small' scales ($k/2\pi \sim 10$), indicating that there may be an important contribution from scales whose turnover time is of the order of the tidal frequency. It should be stressed that these spectra are strongly polluted by turbulent noise, especially at low wavenumbers. Therefore, this result should be seen as a plausible trend. Longer simulations having a larger resolution should be used to confirm this finding.

## 5 SUMMARY AND DISCUSSION

In this paper we have studied the interaction between tides and convection in astrophysical bodies by analysing the effect of a homogeneous oscillatory shear on a fluid flow. This model can be taken to represent the interaction between a large-scale periodic tidal deformation and a smaller scale convective motion. We first considered analytically the limit in which the shear is of low amplitude and the oscillation period is short compared to the time-scales of the unperturbed flow. In this limit there is a viscoelastic response, and we obtained expressions for the effective elastic modulus and viscosity coefficient. The effective viscosity is inversely proportional to the square of the oscillation frequency, with a coefficient that can be positive, negative or zero depending on the properties of the unperturbed flow. We also carried out direct numerical simulations of Boussinesq convection in an oscillatory shearing box and measured the time-dependent Reynolds stress. The results indicate that the effective viscosity falls rapidly as the oscillation frequency is increased, attaining small negative values in the cases we have examined.

Our methods and findings differ significantly from those of other authors. The hypothesis of Zahn (1966) that the effective viscosity

of large eddies is reduced by only a linear factor at high frequencies is incompatible with our analytical and numerical results. We find a quadratic reduction factor at high frequencies, similarly to Goldreich & Nicholson (1977) and Goldreich & Keeley (1977), but not necessarily for the same reason. Our analytical study, which goes beyond that of Goodman & Oh (1997), implies that large eddies generically provide a quadratically reduced viscosity [which Goldreich & Nicholson (1977) considered as an upper bound], but with a coefficient that can be positive, negative or zero. Our numerical results also indicate a quadratically reduced viscosity, with a tendency towards negative values at high frequencies in the cases we have examined, an effect due to the largest scales in the turbulent flow.

The appearance of a negative effective viscosity may be alarming. However, there is no reason in principle why the effective viscosity of a convective fluid need be positive. The consequence of a negative effective viscosity would be that tidal evolution proceeds in the reverse direction to that usually assumed. Angular momentum would be transferred in an 'antifrictional' direction, from the less rapidly rotating component to the more rapidly rotating one. For example, if a small satellite orbits a body with a negative effective viscosity, and if the larger body spins more slowly than the orbit, then the orbit would expand and the larger body would spin down. Energy would be added to the spin-orbit system. (Similarly, orbital eccentricity could be caused to increase in situations where circularization would usually be expected.) While this effect might be called tidal 'antidissipation', it does not contradict the laws of thermodynamics. Within limits, work can indeed be done by the convection on the tidal flow; the energy comes from the buoyancy forces that drive the convection, and ultimately from nuclear or gravitational energy. Since we are considering a regime in which the tidal strain is small and the tidal frequency is high (implying that the negative effective viscosity is much smaller in magnitude than the values estimated from mixing-length theory), only a very small fraction of the energy budget of the star would be diverted into the spin-orbit system. The consequences could still be important, because the nuclear energy content of the Sun, for example, is about five orders of magnitude larger than the orbital energy of (say) a solar-type binary star with an orbital period of 10 days.

A well-known example of a related phenomenon is the negative effective viscosity of convection at low or moderate Rayleigh number in an accretion disc (Lesur & Ogilvie 2010, and references
therein), which would lead to antifrictional angular momentum transport and antidiffusion of the surface density of the disc. This example also serves as a reminder that the Coriolis force present in a rotating system can affect the transport of (angular) momentum; it would therefore be of interest to include rotation in the calculations presented in this paper. Antifrictional processes are also well known in the context of the Earth's atmosphere (e.g. McIntyre 2000).

There are several reasons why this reversed tidal evolution may not, in fact, be found in nature. First, it may be that when a more realistic tidal velocity field is considered and the actual anisotropies of the convection and the Coriolis force are taken into account, the net effective viscosity turns out to be positive. Secondly, in a solar-type star where the convective time-scale becomes short near the surface, the net effective viscosity may be dominated by these more rapid convective elements and work out to be positive. Thirdly, it may be that tidal dissipation is dominated by internal waves and that convection makes a relatively small contribution. Fourthly, the negative effective viscosities that we find might be due to explicit (molecular) viscous effects which are totally negligible in astrophysical objects. Indeed, we have not found a negative turbulent viscosity that exceeds the explicit (molecular) viscosity in magnitude. Nevertheless, it is of interest to note the theoretical possibility of reversed tidal evolution and to be alert to any observational evidence of such an effect.

Penev and collaborators have attempted to address the issue of the reduction of the effective viscosity at high frequencies through simulations of convection in a deep layer. While their results suggest only a linear reduction factor (and a positive numerical value), there are many differences between their work and ours. Our numerical study involves a local model of convection, which is in some ways more limited but is also more controlled because a uniform unstable stratification is imposed. We measure the effective viscosity directly rather than relying on a perturbative calculation or an indirect measurement. Penev et al. (2009a) do not find good agreement between different methods, and their results depend on the shear amplitude. Our work covers a wider range of oscillation frequencies and quantifies the uncertainty due to turbulent noise. It should also be pointed out that our work uses a very small forcing compared to Penev and collaborators. In particular, the rms velocities with and without forcing are identical in all our runs, whereas Penev et al. (2009a) have forced velocities of the order of the rms

velocity in their case of weak forcing, and even five times larger in the case of strong forcing. We believe that such strong forcing could significantly affect the measurement of turbulent viscosity. Several tests that we have performed indicate that doubling the forcing amplitude (to $S = 0.2$) is already enough to modify our results.

In our opinion more work is required to extend both analytical and numerical calculations of this type before the frequency-dependent effective viscosity of a realistic stellar or planetary convective zone can be reliably estimated and the efficiency of tidal dissipation in astrophysical bodies can be better understood.

## ACKNOWLEDGMENT

This research was supported by STFC.

## REFERENCES

* Arnol'd V. I., 1965, Comptes Rendus Acad. Sci. Paris, 261, 17

* Darwin G. H., 1880, Philos. Trans. R. Soc. London, 171, 713

* Galloway D., Frisch U., 1987, J. Fluid Mech., 180, 557

* Goldreich P., Keeley D. A., 1977, ApJ, 211, 934

* Goldreich P., Nicholson P. D., 1977, Icarus, 30, 301

* Goodman J., Oh S. P., 1997, ApJ, 486, 403

* Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742

* Lesur G., Longaretti P.-Y., 2005, A&A, 444, 25

* Lesur G., Longaretti P.-Y., 2007, MNRAS, 378, 1471

* Lesur G., Ogilvie G. I., 2010, MNRAS, 404, L64

* McIntyre M. E., 2000, in Batchelor G. K., Moffatt H. K., Worster M. G., eds., Perspectives in Fluid Dynamics. Cambridge Univ. Press, Cambridge, p. 557

* Meibom S., Mathieu R. D., 2005, ApJ, 620, 970

* Ogilvie G. I., 2001, MNRAS, 325, 231

* Ogilvie G. I., Lin D. N. C., 2004, ApJ, 610, 477

* Ogilvie G. I., Lin D. N. C., 2007, ApJ, 661, 1180

* Penev K., Sasselov D., Robinson F., Demarque P., 2007, ApJ, 655, 1166

* Penev K., Sasselov D., Robinson F., Demarque P., 2009a, ApJ, 704, 930

* Penev K., Barranco J., Sasselov D., 2009b, ApJ, 705, 285

* Rogallo R. S., 1981, NASA Technical Memorandum 81315

* Wirth A., Gama S., Frisch U., 1995, J. Fluid Mech., 288, 249

* Wisdom J., Tremaine S., 1988, AJ, 95, 925

* Zahn J. P., 1966, Ann. Astrophys., 29, 489

# APPENDIX A: COMPOSITION OF AN ARBITRARY TRACELESS VELOCITY GRADIENT TENSOR FROM SIMPLE SHEARS

We consider an arbitrary incompressible velocity gradient tensor, being a 3 × 3 matrix with vanishing trace. Each of the off-diagonal elements represents a simple shear, and therefore the off-diagonal part of the matrix is a linear combination of simple shears. To compose the diagonal elements, we first consider the matrix

$$ \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad (A1) $$

which is a linear combination of simple shears. By a rotation through $\pi/4$ about the z axis, we bring this matrix into the diagonal form

$$ \begin{pmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (A2) $$

In a similar way, the matrix

$$ \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix} \qquad (A3) $$
can be constructed from simple shears. Using a linear combination of these two diagonal matrices, the three diagonal elements of the arbitrary velocity gradient tensor can be composed, subject to the constraint that the sum of the diagonal elements vanishes.

APPENDIX B: VANISHING OF CERTAIN TRIPLE VELOCITY CORRELATIONS FOR STATISTICALLY ISOTROPIC FLOWS

Let $v(\mathbf{x})$ be a solenoidal vector field in $d \ge 2$ dimensions satisfying periodic boundary conditions, and let $\langle \cdot \rangle$ be a spatial average over a periodic cell. (We omit primes on $\mathbf{x}$ and $\nabla$ in this appendix.) We assume that $v(\mathbf{x})$ has zero mean and is statistically isotropic so that any averages of tensor products of $v$ and its derivatives are isotropic tensors.<fnref n="1" /> The inverse Laplacian operator $\Delta^{-1}$ is well defined and self-adjoint for periodic fields of zero mean.

We first consider the tensor

$$D_{abcd} = \langle v_a v_b \partial_c v_d \rangle. \qquad (B1)$$

As an isotropic tensor, this must be a linear combination of products of Kronecker deltas:

$$D_{abcd} = t_1 \delta_{ab} \delta_{cd} + t_2 \delta_{ac} \delta_{bd} + t_3 \delta_{ad} \delta_{bc}. \qquad (B2)$$

The symmetry and contraction conditions $D_{abcd} = D_{bacd}$ and $D_{abcc} = 0$ imply $t_2 = t_3$ and $dt_1 + t_2 + t_3 = 0$. Thus

$$D_{abcd} = t_4 [d(\delta_{ac}\delta_{bd} + \delta_{ad}\delta_{bc}) - 2\delta_{ab}\delta_{cd}]. \qquad (B3)$$

However, we also have

$$D_{abba} = \langle v_a v_b \partial_b v_a \rangle = \langle \partial_b (\frac{1}{2} v_a v_b) \rangle = 0 = t_4 d(d-1)(d+2), \qquad (B4)$$

and so

$$D_{abcd} = 0. \qquad (B5)$$

We next consider

$$D_{abcdef} = \langle v_a v_b \partial_c \partial_d \partial_e \Delta^{-1} v_f \rangle. \qquad (B6)$$

Isotropy and symmetry imply

$$D_{abcdef} = t_5(\delta_{ab}\delta_{cd}\delta_{ef} + \delta_{ab}\delta_{ce}\delta_{df} + \delta_{ab}\delta_{cf}\delta_{de}) + t_6(\delta_{ac}\delta_{bd}\delta_{ef} + \delta_{ac}\delta_{be}\delta_{df} + \delta_{ad}\delta_{bc}\delta_{ef} + \delta_{ad}\delta_{be}\delta_{cf} + \delta_{ae}\delta_{bc}\delta_{df} + \delta_{ae}\delta_{bd}\delta_{ef}) \\ + t_7(\delta_{ac}\delta_{bf}\delta_{de} + \delta_{ad}\delta_{bf}\delta_{ce} + \delta_{ae}\delta_{bf}\delta_{de} + \delta_{af}\delta_{bc}\delta_{de} + \delta_{af}\delta_{bd}\delta_{ce} + \delta_{af}\delta_{be}\delta_{cd}). \qquad (B7)$$

(There are three types of term here: those in which the two undifferentiated $v$ are paired, those in which each such $v$ is paired with a $\partial$, and those in which one of them is paired with the differentiated $v$. Symmetry demands that the coefficients of terms of the same type are equal.)

We require the contraction $D_{abcced}$ (and those related by symmetry) to vanish. Thus

$$(d+2)t_5 \delta_{ab} \delta_{cd} + t_6(d+2)(\delta_{ac} \delta_{bd} + \delta_{ad} \delta_{bc}) + 2t_7(\delta_{ab} \delta_{cd} + \delta_{ac} \delta_{bd} + \delta_{ad} \delta_{bc}) = 0, \qquad (B8)$$

which implies $t_5 = t_6$ and $(d+2)t_5 = -2t_7$, and so

$$D_{abcdef} = t_8[(d+2)(\delta_{ac}\delta_{bf}\delta_{de} + \delta_{ad}\delta_{bf}\delta_{ce} + \delta_{ae}\delta_{bf}\delta_{cd} + \delta_{af}\delta_{bc}\delta_{de} + \delta_{af}\delta_{bd}\delta_{ce} + \delta_{af}\delta_{be}\delta_{cd}) \\ -2(\delta_{ab}\delta_{cd}\delta_{ef} + \delta_{ab}\delta_{ce}\delta_{df} + \delta_{ab}\delta_{cf}\delta_{de} + \delta_{ac}\delta_{bd}\delta_{ef} + \delta_{ac}\delta_{be}\delta_{df} + \delta_{ad}\delta_{bc}\delta_{ef} + \delta_{ad}\delta_{be}\delta_{cf} + \delta_{ae}\delta_{bc}\delta_{df} + \delta_{ae}\delta_{bd}\delta_{cf})]. \qquad (B9)$$

The contraction $D_{abcdff}$ (and those related by symmetry) produces a Laplacian, so $D_{abcdff} = D_{abcf}$, which we have already shown to vanish. Thus

$$(d+4)t_8[d(\delta_{ac}\delta_{bf} + \delta_{af}\delta_{bc}) - 2\delta_{ab}\delta_{cf}] = 0, \qquad (B10)$$

which implies $t_8 = 0$ and so

$$D_{abcdef} = 0. \qquad (B11)$$

<footnote n="1">¹ It is debatable whether this assumption is compatible with the periodic boundary conditions. While the cubic symmetry imposed by the boundary conditions of a periodic cube is compatible with isotropy for tensors of second rank, this is not true for higher ranks. One could argue that the flow can be nearly statistically isotropic if it is dominated by scales smaller than the box.</footnote>
The last tensor to consider is

$$D_{ijklm} = \langle v_i v_j \partial_k \partial_l \partial_m \partial_n \Delta^{-2} v_n \rangle. \qquad (B12)$$

We give an abbreviated argument. Isotropy and symmetry imply

$$D_{ijklm} = t_9(\delta_{ab}\delta_{cd}\delta_{ef}\delta_{gh} + \delta_{ab}\delta_{cd}\delta_{eg}\delta_{fh} + \dots) + t_{10}(\delta_{ac}\delta_{bd}\delta_{ef}\delta_{gh} + \delta_{ac}\delta_{bd}\delta_{eh}\delta_{fg} + \dots) + t_{11}(\delta_{ah}\delta_{bc}\delta_{df}\delta_{eg} + \delta_{ah}\delta_{bc}\delta_{dg}\delta_{ef} + \dots). \qquad (B13)$$

Requiring the contraction $D_{ijklm}$ to vanish implies $t_9 = t_{10}$ and $(d+4)t_9 = -t_{11}$. The contraction $D_{ijklm}$ should produce $D_{ijklm}$, which we have already shown to vanish. This is of the form given above, with $t_5 = (d+4)t_9 + 2t_{10}$, $t_6 = (d+6)t_{10}$ and $t_7 = 2t_{10} + (d+4)t_{11}$. Therefore $t_9 = t_{10} = t_{11} = 0$ and so $D_{ijklm} = 0$.

This paper has been typeset from a TeX/LA TeX file prepared by the author.