Vincent Mons

Claude Cambon

Pierre Sagaut

Antoine Briard

Keywords:

0) ∼ k σ , unlike unstably stratified homogeneous turbulence where γ strongly depends on σ . The MCS model relies on the truncation at the second order of the spectral two-point velocity correlation expansion into spherical harmonics. The expansion is here pursued at the next even order, the fourth one: the noteworthy consequence is that γ is decreased compared to MCS and is thus closer to values obtained in direct numerical simulations and experiments. Finally, some analytical considerations about odd-order contributions in the expansion of polarisation anisotropy are proposed.

In this work, the spectral modelling developed in MCS [Mons, Cambon, Sagaut. A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors. J Fluid Mech. 2016;788:147-182] for shear-driven turbulence is further analysed and then improved. First, using self-similarity arguments, it is shown that the asymptotic kinetic energy exponential growth rate γ is independent of the large scales infrared slope σ , with E(k →

Introduction

Natural turbulent flows represent a great challenge in the turbulence community in terms of modelling. Indeed, in atmospheric turbulent flows, for example, several complex mechanisms strongly interact in an intricate manner, on distinct characteristic time scales, such as rotation, stratification, helicity, and shear. Among these mechanisms, shear flows are probably the most complex, since no symmetry at all survive, unlike rotation and stratification which nevertheless remain statistically axisymmetric.

In the past years, the authors have analysed separately each of these configurations with the help of an adapted eddy-damped quasi-normal Markovian (EDQNM) closure able of handling strongly anisotropic flows, in order to better understand what are the intrinsic properties of each mechanism. A significant amount of results, both theoretical and numerical, were exposed in various publications [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Spectral modelling for passive scalar dynamics in homogeneous anisotropic turbulence[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF][START_REF] Briard | Anisotropic spectral modelling for unstably stratified homogeneous turbulence[END_REF][START_REF] Burlot | Spectral modelling of high Reynolds number unstably stratified homogeneous turbulence[END_REF][START_REF] Briard | Dynamics of helicity in homogeneous skew-isotropic turbulence[END_REF][START_REF] Cambon | Anisotropic triadic closures for shear-driven and buoyancydriven turbulent flows[END_REF][START_REF] Cambon | Third-order statistics and the dynamics of strongly anisotropic turbulent flows[END_REF][START_REF] Gréa | Dynamics and structure of unstably stratified homogeneous turbulence[END_REF], and still a lot of work needs to be done regarding shear flows specifically. Indeed, the anisotropic EDQNM model developed in Mons, Cambon, and Sagaut (MCS) [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF] handles very satisfactorily various straining processes, for instance when the mean-velocity gradient matrix is symmetric, like in axisymmetric contractions or expansions, or in plane distortions. However, for pure plane shear flows where this matrix is not symmetric and has only one non-zero component, the MCS model gives a kinetic energy exponential growth rate γ larger than common values obtained in direct numerical simulations (DNS): see [START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF] for a review. Despite the clear tendency of increasing values of γ in more recent DNS [START_REF] Isaza | On the asymptotic behaviour of large-scale turbulence in homogeneous shear flow[END_REF][START_REF] Sukheswalla | Simulation of homogeneous turbulent shear flows at higher Reynolds numbers: numerical challenges and a remedy[END_REF] which can reach larger Reynolds numbers with a better spatial resolution, a discrepancy still remains. This quantitative issue regarding the value of γ is probably due to a lack of angular information regarding the distribution of anisotropy within our anisotropic EDQNM modelling. Consequently, we focus here on pure plane shear flows, a particular case of shear-driven turbulence. The reader should nevertheless keep in mind that the model is more general and can handle various shear-driven configurations.

In homogeneous incompressible turbulence, the spectral velocity-velocity correlation Rij can be exactly decomposed into directional and polarisation anisotropies [START_REF] Cambon | Anisotropic developments for homogeneous shear flows[END_REF][START_REF] Cambon | Spectral approach to non-isotropic turbulence subjected to rotation[END_REF], which are quantified by the two scalar quantities E(k, t) and Z(k, t), which both depend on time and on the wavevector k. In order to deal with analytical expressions in the modelling of anisotropy, E and Z are expanded into spherical harmonics, and this expansion is further truncated at the first non-trivial order, the second one, the zeroth-order being the isotropic state where Z = 0 and E = E/(4π k 2 ), with E being the kinetic energy spectrum. To further reduce the numerical cost of the simulations, the detailed equations for the directional and polarisation anisotropies are spherically averaged.

The numerical results obtained with such an approximation were thoroughly discussed in the previously mentioned references and compared quantitatively well with both DNS and experiments, in various configurations, from axisymmetric contractions, expansions, and plane distortions [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF], to the transport of a passive scalar field in the presence of a mean gradient with a variable Prandtl number [START_REF] Briard | Spectral modelling for passive scalar dynamics in homogeneous anisotropic turbulence[END_REF][START_REF] Briard | Prandtl number effects in decaying homogeneous isotropic turbulence with a mean scalar gradient[END_REF], along with the case of unstably stratified turbulence [START_REF] Briard | Anisotropic spectral modelling for unstably stratified homogeneous turbulence[END_REF]. However, for shear flows, the model could not recover accurately the anisotropy distribution by investigating, for instance, several components of the global indicator b ij = u i u j /2K -δ ij /3, where K = ∞ 0 E(k) dk is the kinetic energy, nor the asymptotic value of the kinetic exponential growth rate γ as said earlier. Despite these defects, a good agreement was found for the shear parameter S * = 2 KS/ , where S is the shear rate and the kinetic energy dissipation rate, an important quantity for one-point modelling [START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF]. Thus, it was concluded that anisotropy was not well-enough captured, and more precisely polarisation anisotropy [START_REF] Zhu | Rotating shear-driven turbulent flows: towards a spectral model with angle-dependent linear interactions[END_REF], with notably a too large value for b 13 (for S = -∂U 1 /∂x 3 ). So far, it was impossible to quantify the anisotropic information lost by the second-order truncation. Therefore, we propose here to go beyond the second-order truncation and deal analytically with the fourth-order, to know for sure if this allows to decrease the exponential growth rate γ or not. The derivation itself is rather complex and constitutes one of the main results of the present work. The second significant finding is an analytical proof showing an essential difference between the dynamics of shear flows and unstably stratified turbulence, configurations both analysed with the anisotropic EDQNM model, and in which the kinetic energy grows exponentially, but with completely different properties.

It is worth recalling that the present model addresses only second-order moments of the velocity field, so that high-order statistics are not available. It notably implies that we cannot investigate intermittent properties of shear flows [START_REF] Casciola | Scale-by-scale budget and similarity laws for shear turbulence[END_REF]. However, intermittent statistics are often observed for high-order moments of the velocity field [START_REF] Casciola | Scale-by-scale budget and similarity laws for shear turbulence[END_REF][START_REF] Gualtieri | Scaling laws and intermittency in homogeneous shear flow[END_REF][START_REF] Biferale | Isotropy vs anisotropy in small-scale turbulence[END_REF][START_REF] Schumacher | Derivative moments in turbulent shear flows[END_REF] so that it remains compatible with our approach. Indeed, it has been observed with the model in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF] that small scales of the velocity field second-order moments return to isotropy at large Reynolds numbers, a feature consistent with DNS and experiments [START_REF] Pumir | Turbulence in homogeneous shear flows[END_REF][START_REF] Garg | On the small scale structure of simple shear flow[END_REF][START_REF] Shen | The anisotropy of the small scale structure in high Reynolds number (R λ ∼ 1000) turbulent shear flow[END_REF]: nevertheless, some perspectives to study higher-order statistics with similar modelling techniques are discussed in the concluding section.

Finally, let us mention that other tools to analyse anisotropic features of shear flows are possible. For example, our approach is consistent with the irreducible representation of the SO(3) symmetry group that yields expansions in terms of spherical harmonics as well [START_REF] Rubinstein | Scalar and tensor spherical harmonics expansion of the velocity correlation in homogeneous anisotropic turbulence[END_REF]. Applications to shear flows are also developed directly in physical space, with focus on anisotropic structure functions (see [START_REF] Biferale | Anisotropy in turbulent flows and in turbulent transport[END_REF] for a review): note that expansions in terms of the separation vector r in physical space and in terms of k in spectral space can be partly reconciled [START_REF] Cambon | Third-order statistics and the dynamics of strongly anisotropic turbulent flows[END_REF]. In a recent work [START_REF] Clark | Generation of anisotropy in turbulent flows subjected to rapid distortion[END_REF], the basis functions of the SO(3) symmetry group were used to derive an analytical dynamical model in the rapid distortion linear limit, instead of analysing statistics resulting from DNS as usual. This is in the spirit of the present approach where evolution equations, fully non-linear, however, are written for tensors corresponding to these basis functions.

The paper is organised as follows: as a starting point, the MCS model with the secondorder truncation of the spherical harmonics expansion is recalled in Section 2, with some details gathered in Appendix 1. A first new result is presented regarding the independence of the exponential growth rate with the large scales initial conditions. Then, the fourthorder modelling is presented in Section 3 along with the corresponding new production and non-linear terms. Finally, in Section 4, some analytical developments about odd-order terms of the expansion for the polarisation Z are presented, unfortunately inconclusive. Discussions and further perspectives are gathered in Section 5.

The anisotropic EDQNM modelling for homogeneous turbulence

In this part, the main elements of the MCS anisotropic spectral model are recalled for clarity. First, the so-called E-Z decomposition of the spectral two-point second-order velocity-velocity correlation is presented. Afterwards, the modelling of anisotropy through spherically averaged deviatoric descriptors is exposed, along with brief details about the EDQNM procedure used to close the non-linear terms in the evolution equations. Then, the final spherically averaged equations for homogeneous turbulence are derived.

The E-Z decomposition

The spectral two-point second-order velocity-velocity correlation, or spectral Reynolds tensor, is defined as

Rij (k, t)δ(k -p) = û * i (p, t) ûj (k, t) , ( 1 )
where • denotes the Fourier transform, (•) * is the complex conjugate, • is an ensemble average, and k and p are wavevectors. The tensor Rij further verifies the so-called Craya equation, namely

∂ ∂t -A ln k l ∂ ∂k n + 2νk 2 Rij (k) + M in (k) Rnj (k) + M jn (k) Rni (k) = T NL ij (k), ( 2
)
where ν is the kinematic viscosity, T NL ij is the total non-linear transfer, M ij (k) = (δ in -2α i α n )A nj with α i = k i /k, and A ln = ∂U l /∂x n the mean-velocity gradient matrix, whose symmetric and antisymmetric parts are respectively identified with + andsuperscripts in what follows.

Since Rij is a homogeneous, symmetric, and solenoidal second-order tensor, it can be intrinsically and exactly decomposed into three distinct parts [START_REF] Cambon | Anisotropic developments for homogeneous shear flows[END_REF]

Rij (k) = R(iso) ij (k) + R(dir) ij (k) + R(pol) ij (k), (3)
= E(k) 4πk 2 P ij (k) + E(k) - E(k) 4π k 2 P ij (k) + (Z(k)N i (k)N j (k)), ( 4 )
where P ij = δ ij -α i α j is a projector and N i is a helical mode [START_REF] Cambon | Spectral approach to non-isotropic turbulence subjected to rotation[END_REF][START_REF] Waleffe | The nature of triad interactions in homogeneous turbulence[END_REF], perpendicular to the wavevector k. These three parts in the decomposition (4) refer to the isotropic contribution, the directional anisotropy, which is the difference between the energy along the wavevector k and the spherically averaged energy, and the polarisation anisotropy, which reflects the difference of energy between components of Rij . The kinetic energy density E(k, t) is linked to directional anisotropy and is connected to the kinetic energy spectrum E(k, t) through a spherical integration on a sphere S k of radius k

E(k, t) = S k E(k, t) d 2 k = S k Rii (k, t) 2 d 2 k, ( 5
)
whereas the complex-valued scalar Z(k, t) reflects polarisation anisotropy and reads

Z(k, t) = Rij (k, t) 2 N i (-k)N j (-k). ( 6
)
Both E and Z verify exact evolution equations which can be found in several references [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Cambon | Third-order statistics and the dynamics of strongly anisotropic turbulent flows[END_REF], namely

∂ ∂t -A ln k l ∂ ∂k n + 2νk 2 E -A + ij α i α j E + A + ij (ZN i N j ) = T E , ( 7 )
∂ ∂t -A ln k l ∂ ∂k n + 2νk 2 Z + 2iZ CH + A + ij N * i N * j E + A ij N * i N j Z = T Z , (8)
where CH reflects the motion of the time-evolving Craya-Herring frame with respect to a fixed reference [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF], and with

T E = T NL ii /2 and T Z = T NL ij N * i N * j /2.
One could solve these exact evolution equations, but this would be extremely demanding in terms of computational resources. Note that considering the variables (kE) and (kZ) instead of only E and Z simplifies the equations, as done in [START_REF] Cambon | Anisotropic triadic closures for shear-driven and buoyancydriven turbulent flows[END_REF].

Instead, we choose to model anisotropy in two steps: first, the non-linear directional and polarisation transfers T E and T Z are closed by a classical EDQNM procedure. Afterwards, both E and Z are expanded into spherical harmonics to obtain explicit expressions for the directional and polarisation parts of the spectral Reynolds tensor. This expansion was further truncated at the second order for the sake of simplicity in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF]: the immediate drawback is the loss of angular anisotropic information which makes the production terms, linear with the mean-velocity gradient matrix A ln , not exact anymore. The consequences of such an approximation were thoroughly examined and discussed in several references [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Spectral modelling for passive scalar dynamics in homogeneous anisotropic turbulence[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF], so that they are not further argued here. The aim of this study is precisely to improve the modelling of anisotropy through a truncation at the next even order, the fourth-one, in Section 3. Before that, for consistency and clarity, the details of the truncation at the secondorder are recalled in the following part: in addition, a new theoretical feature is presented to better understand the exponential growth rate of kinetic energy obtained numerically with the model.

Truncated expansion of E and Z at the second-order

The expansion into spherical harmonics of E and Z is the key ingredient for the modelling of anisotropy. Such expansions read

E(k) = E 0 (k)(1 + U (dir)2 ij (k)α i α j + U (dir)4 ijpq (k)α i α j α p α q + • • • ), ( 9 )
Z(k) = 1 2 E 0 (k)(U (pol)2 ij (k) + U (pol)4 ijpq (k)α p α q + • • • )N * i (k)N * j (k), (10)
with E 0 = E/(4π k 2 ), and where time dependence has been omitted for clarity. The determination of the dynamics of the symmetric and deviatoric second-order tensors U (dir)2 ij and U

(pol)2 ij
was the topic of [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF], whereas the determination of that of the fourth-order tensors U (dir)4

ijpq and U (pol)4 ijpq constitutes one of the two theoretical contributions of this work, addressed in detail in Section 3. Odd-order terms of the expansion of E are exactly zero, which comes from the symmetry of Rij [START_REF] Cambon | Anisotropic developments for homogeneous shear flows[END_REF][START_REF] Rubinstein | Scalar and tensor spherical harmonics expansion of the velocity correlation in homogeneous anisotropic turbulence[END_REF]; the zeroth-order corresponds to isotropy, configuration in which the polarisation anisotropy Z = 0. The possibility of non-zero odd-order terms for Z is discussed later on in Section 4.

As mentioned before, the spherical average transforms the k-dependence of the twopoint second-order correlations into a k one and thus strongly reduces the computational time of the simulations but also causes a loss of information regarding the distribution of anisotropy of the flow. In [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF], part of this information was nevertheless recovered, thanks to the deviatoric and symmetric tensors H (dir) ij and H (pol) ij , which can be linked to

U (dir)2 ij and U (pol)2 ij through 2E(k, t)H (dir) ij (k, t) = S k R(dir) ij (k, t) d 2 k = - 2 15 U (dir)2 ij (k, t)E(k, t), ( 11
) 2E(k, t)H (pol) ij (k, t) = S k R(pol) ij (k, t) d 2 k = 2 5 U (pol)2 ij (k, t)E(k, t). ( 12
)
Hence, the second-order spectral tensor is given by

φ ij (k, t) = S k Rij (k, t) d 2 k = 2E(k, t) δ ij 3 + H (dir) ij (k, t) + H (pol) ij (k, t) . ( 13
)
With these expressions [START_REF] Sukheswalla | Simulation of homogeneous turbulent shear flows at higher Reynolds numbers: numerical challenges and a remedy[END_REF] and [START_REF] Cambon | Anisotropic developments for homogeneous shear flows[END_REF], it is thus possible to determine the closed expressions of E and Z at the second-order, namely

E(k) = E 0 (k)(1 -15H (dir) ij (k)α i α j ), Z(k) = 5 2 E 0 (k)H (pol) ij (k)N * i (k)N * j (k). ( 14
)
For now, only the truncated expansions at the second-order are considered, so that the modelled spectral Reynolds tensor reads

Rij = E 0 P ij (1 -15H (dir) pq α p α q ) Re2 ij + 5E 0 P in P jm H (pol) mn + 1 2 P ij H (pol) pq α p α q Rz2 ij . ( 15
)
In the following section, the evolution equations of the three main spectra E, EH

(dir) ij
, and

EH (pol) ij
are given within the anisotropic EDQNM framework.

The EDQNM procedure and spherically averaged equations

The final step of the anisotropic EDQNM modelling is to close the non-linear transfers T E and T Z using a classical EDQNM procedure. Details about the fundamentals of EDQNM can be found, for instance, in [START_REF] Lesieur | Turbulence in fluids[END_REF] for HIT, or in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF] with the anisotropic formalism. The modelled equation [START_REF] Zhu | Rotating shear-driven turbulent flows: towards a spectral model with angle-dependent linear interactions[END_REF] of Rij is further injected into the quasi-normal expression of the non-linear transfers. The eddy-damping term is kept isotropic for generality purposes and also to avoid the introduction of new arbitrary constants: this has proven to be relevant (in non-rotating flows) in various comparisons performed against DNS and experiments in different configurations [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Spectral modelling for passive scalar dynamics in homogeneous anisotropic turbulence[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF][START_REF] Briard | Anisotropic spectral modelling for unstably stratified homogeneous turbulence[END_REF].

After some algebra, one obtains the closed expressions of the non-linear transfers T E and T Z in the evolution equations [START_REF] Cambon | Anisotropic triadic closures for shear-driven and buoyancydriven turbulent flows[END_REF] and ( 8) of E and Z

T E = 2 θ kpq kp[(E + X )((xy + z 3 )(E -E) -z(1 -z 2 )( X -X)) + X (1 -z 2 )(x X -y X )] d 3 p, ( 16
)
T Z = 2 θ kpq kp e -2iλ [(E + X )((xy + z 3 )( X -X) -z(1 -z 2 )(E -E) + i(y 2 -z 2 ) X ) + i X (1 -z 2 )(x(E + X) -iy X )] d 3 p, ( 17
) with E(k, t) = E, E(p, t) = E , E(q, t) = E , X = Z(k, t) e 2iλ , X = Z(p, t) e 2iλ
, and X = Z(q, t) e 2iλ . The λ, λ , and λ are rotation angles around k, p, and q which permit to perform the calculations in a frame more adapted than the Craya-Herring one [START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF]. x, y, and z are the cosines of the angles formed by p and q, q and k, and k and p, respectively. And , . . .) were discarded in the non-linear transfers, because they were assumed to be negligible [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF]: this assumption is successfully assessed for the first time in Section 2.5.

θ
Consequently, after spherical averaging, one ends up with three generalised Lin equations that describe the evolution of the isotropic, directional and polarisation parts of the turbulent velocity field, given under the same compact expression

∂ ∂t + 2νk 2 E(k, t) = S L(iso) (k, t) + S NL(iso) (k, t), ( 18
)
∂ ∂t + 2νk 2 E(k, t)H (dir) ij (k, t) = S L(dir) ij (k, t) + S NL(dir) ij (k, t), ( 19
)
∂ ∂t + 2νk 2 E(k, t)H (pol) ij (k, t) = S L(pol) ij (k) + S NL(pol) ij (k, t), ( 20
)
where S NL(iso) , S

NL(dir) ij

, and S NL(pol) ij are the non-linear spherically averaged transfers, and

S L(iso) , S L(dir) ij
, and S L(pol) ij are the production terms resulting from the presence of a meanvelocity gradient. The explicit expressions of these linear and non-linear transfers are given in Appendix 1.

Thanks to the previous spectral anisotropy descriptors

H (dir) ij and H (pol) ij
, one can express the deviatoric normalised part of the Reynolds stress tensor, which refers to the global anisotropy of the flow, according to

b ij (t) = u i u j 2 K - δ ij 3 = 1 K(t) ∞ 0 E(k, t)[H (dir) ij (k, t) + H (pol) ij (k, t)] dk. ( 21
)
The three generalised Lin equations ( 18), [START_REF] Schumacher | Derivative moments in turbulent shear flows[END_REF], and (20) constitute the main result of MCS. This spectral modelling was then applied in [START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF] exclusively for shear flows. Notably, it was shown in the latter reference that:

• The kinetic energy grows exponentially at large St and Reynolds numbers at the rate 0.33 ≤ γ ≤ 0.34, with K ∼ exp(γ St). • This exponential growth rate γ is found numerically to be independent of the large scales initial conditions σ , where E(k < k L ) ∼ k σ , with k L being the integral wavenumber [START_REF] Eyink | Free decay of turbulence and breakdown of self-similarity[END_REF]. Such an independence to large scales was also reported in [START_REF] Clark | A spectral model applied to homogeneous turbulence[END_REF] for shear flows. • This value γ 0.33 is higher than common values obtained in DNS and experiments, as reported in Table 1 of [START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF], even though more recent numerical studies report larger exponential growth rates.

In the following part, we come back on both this value for γ and on its independence with regard to σ with theoretical arguments. In addition, a qualitative comparison is made with the asymptotic regimes of unstably stratified homogeneous turbulence, where the kinetic energy also grows exponentially but with a strong dependence on σ .

The exponential growth rate γ of kinetic energy

In this section, we provide theoretical arguments to explain both the rather large value of the kinetic energy exponential growth rate γ = 0.33 in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF] and its independence with regard to large scales' initial conditions σ : in particular, the latter point means that at sufficiently large St or Reynolds numbers, the growth rate is similar for Saffman (σ = 2) and Batchelor (σ = 4) turbulence.

In a sustained shear flow, kinetic energy is continuously produced by the mean-velocity gradient and verifies the evolution equation

dK dt = SR 13 (t) -(t), dU 1 dx 3 = -S, A ij = dU i dx j , ( 22
)
where is the kinetic energy dissipation rate. Kinetic energy eventually grows exponentially [START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF][START_REF] Tavoularis | Asymptotic laws for transversely homogeneous turbulent shear flows[END_REF][START_REF] George | The decay of homogeneous isotropic turbulence[END_REF] at a rate γ = 2b 13 -/KS when the global anisotropy indicators b ij have reached an asymptotic state.

Values of γ obtained in experiments [START_REF] Tavoularis | Experiments in nearly homogenous turbulent shear flow with a uniform mean temperature gradient. Part 1[END_REF][START_REF] Tavoularis | Further experiments on the evolution of turbulent stresses and scales in uniformly sheared turbulence[END_REF][START_REF] Souza | The structure of highly sheared turbulence[END_REF] and in DNS [START_REF] Isaza | On the asymptotic behaviour of large-scale turbulence in homogeneous shear flow[END_REF][START_REF] Sukheswalla | Simulation of homogeneous turbulent shear flows at higher Reynolds numbers: numerical challenges and a remedy[END_REF][START_REF] Gualtieri | Scaling laws and intermittency in homogeneous shear flow[END_REF][START_REF] Pumir | Turbulence in homogeneous shear flows[END_REF][START_REF] Brethouwer | The effect of rotation on rapidly sheared homogeneous turbulence and passive scalar transport. Linear theory and direct numerical simulation[END_REF][START_REF] Pumir | Persistent small scale anisotropy in homogeneous shear flows[END_REF] are smaller than the γ = 0.33 obtained with the present anisotropic EDQNM modelling, and some values are reported in Table 1 of [START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF]. In what follows, we first bring theoretical arguments which were not provided in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF] to justify the value γ = 0.33 obtained with the model. To this end, we use the methodology of Soulard et al. [START_REF] Soulard | Large-scale analysis of self-similar unstably stratified homogeneous turbulence[END_REF] who derived a prediction for the growth rate of the Rayleigh-Taylor mixing zone length as a function of the infrared slope σ . These theoretical arguments were also applied to the framework of unstably stratified homogeneous turbulence (USHT) in [START_REF] Briard | Anisotropic spectral modelling for unstably stratified homogeneous turbulence[END_REF][START_REF] Burlot | Spectral modelling of high Reynolds number unstably stratified homogeneous turbulence[END_REF][START_REF] Soulard | Large-scale analysis of self-similar unstably stratified homogeneous turbulence[END_REF] to determine the exponential growth rate of the kinetic energy: results will be discussed hereafter.

The demonstration essentially relies on two features: (i) the linear dynamics of large scales dominated by the anisotropy production and (ii) the self-similarity of the kinetic energy spectrum. In the self-similar state, the kinetic energy and its dissipation rate evolve as K ∼ ∼ exp(γ St): combining this with dimensional analysis for the integral scale L ∼ K 3/2 / yields L ∼ exp(γ St/2). Let us call γ E the exponential growth rate of the kinetic energy spectrum large scales. Assuming that the main contribution to kinetic energy comes from large scales, one has

K(t) = ∞ 0 E(k, t) dk k L 0 k σ exp(γ E St) dk = k L σ +1 σ + 1 exp(γ E St). ( 23
)
Further using the self-similarity hypothesis for the kinetic energy and the integral scale

k L ∼ 1/L ∼ exp(-γ St/2), one obtains γ E = γ 2 (σ + 3). ( 24
)
This result can also be obtained differently, in a manner analogous to [START_REF] George | The decay of homogeneous isotropic turbulence[END_REF]: assuming that at large scales the kinetic energy is given by K ∼ S 2 L 2 , one can expand the kinetic energy spectrum according to

E(k < k L , t) = K(t)L(t)(kL(t)) σ = S 2 L(t) σ +3 k σ ∼ S 2 exp σ + 3 2 γ St k σ , (25)
which further illustrates that the self-similar evolution of the integral scale is crucial in the demonstration. It is worth noting that these two equations hold for both shear flows and USHT, if one replaces S by the stratification frequency in the latter case.

Then, the exponential growth rate γ E is also given by the largest eigenvalue max of the linear operator of the generalised Lin equations system ( 18), [START_REF] Schumacher | Derivative moments in turbulent shear flows[END_REF], and [START_REF] Pumir | Turbulence in homogeneous shear flows[END_REF]. The linear operator of shear flows verifies, dropping the non-linear and viscous terms and using the expression (25) of self-similarity

∂ ∂t ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ E EH (dir) 11 EH (dir) 13 EH (dir) 33 EH (pol) 11 EH (pol) 13 EH (pol) 33 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = S ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 2(σ + 2) 0 0 2 0 0 0 -2(σ +11)
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ E EH (dir) 11 EH (dir) 13 EH (dir) 33 EH (pol) 11 EH (pol) 13 EH (pol) 33 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . (26)
Note that (i) the spanwise component does not appear since H 22 = -H 11 -H 33 because of incompressibility and (ii) each component of EH ij was assumed to behave spatially at large scales like ∼ k σ , which is verified numerically. The maximum real eigenvalue max (σ ) of the linear operator is given by

max (σ ) = 1 21 √ 10 [-1112 + 267σ (1 + σ ) + 49 3664 -168σ -159σ 2 + 18σ 3 + 9σ 4 ] 1/2 . ( 27
)
The variations of max as a function of σ are presented in Figure 1(a) (dashed line): one can note that for common values of the infrared slope, basically from σ = 1 to σ = 4, there is a strong variation of max . Recall that max was obtained by assuming that the dynamics of the large scales is linear: this is further assessed in Figure 1(b) where the budget terms of the evolution equation ( 18) of E(k, t) are displayed: it is clear that at large scales the linear production term is much larger than the non-linear transfer term.

S (tot) ij = 2(δ ij S (iso) /3 + S (dir) ij + S (pol) ij
).

linear to non-linear terms ratios are presented for the components () 11 , () 13 , and () 33 . On the contrary, for k > k S , the non-linear transfers drive the dynamics, as observed by values of the ratios much lower than unity. The final theoretical prediction for the kinetic energy exponential growth rate γ th is obtained by equating the most amplified eigenvalue and the kinetic energy spectrum growth rate

γ E = max (σ ) ⇔ γ th = 2 max (σ ) σ + 3 , ( 28
)
and is presented in Figure 1(a), where it is clear that the true exponential growth rate of the kinetic energy (-curve) accounting for non-linearities is much lower than the linear prediction (--curve), and almost independent of the large scales initial conditions σ . In particular, one has for the integer values of the infrared slope γ th (σ = 1) = 0.358, γ th (σ = 2) = 0.339, γ th (σ = 3) = 0.346, and γ th (σ = 4) = 0.367. The implications of Equation ( 28) and the numerical results in Figure 1(a) are twofold. First, there is a good agreement between the new theoretical prediction (28) -which relies on the linear dynamics at large scales and the self-similarity of the kinetic energy spectrum -and the numerical simulations of [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF], which gave a value around 0.33 ≤ γ ≤ 0.34, as recalled for clarity in Figure 3(a) in grey.

Secondly, one can remark that the predictions for the exponential growth rate hardly depend on σ in Figure 1(a), which is also in agreement with the simulations of [START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF]. This is of great theoretical interest because it strongly differs from USHT, where kinetic energy also grows exponentially, but with a strong dependence on large scales according to the theoretical prediction [START_REF] Burlot | Spectral modelling of high Reynolds number unstably stratified homogeneous turbulence[END_REF]

γ USHT = 4 σ + 3 , ( 29
)
where the kinetic energy exponential growth rate γ USHT varies from 1 to 4/7 for σ = 1 up to σ = 4. It is worth noting that the present model is able to capture the strong dependence of the kinetic energy exponential growth rate in USHT as well [START_REF] Briard | Anisotropic spectral modelling for unstably stratified homogeneous turbulence[END_REF]. The reason why the asymptotic anisotropic states strongly differ between shear flows and USHT is because the linear operator of the production terms does not depend on σ in the latter case. Indeed, the σ -contributions in the linear operator of shear flows appear due to the spatial-derivative terms in ∂/∂k (see Appendix 1), which arise from the mean-flow advection term in the Craya equation ( 2): basically, they balance the σ -dependence of γ E . In conclusion, it has been shown analytically within the anisotropic EDQNM modelling that the kinetic energy exponential growth rate in shear flows is independent of the large scales initial conditions σ , consistently with self-similarity arguments. This independence strongly differs from the asymptotic exponential growth rate in USHT which depends on σ . The model is able to recover both phenomenologies: this is the first theoretical contribution of this paper.

Quadratic contributions of anisotropy

For the modelling of non-linear transfers, it was assumed for simplicity reasons in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF] that quadratic contributions of anisotropy could be discarded, namely contributions like

H (dir) il H (dir) lj , H (dir) il H (pol) lj
, . . . and so on. In this part, we verify this assumption. The four corresponding quadratic non-linear transfers were explicitly computed and are gathered in Appendix 2 because the expressions are rather lengthy: the isotropic one Q NL(iso) , the directional and polarisation ones Q

NL(dir) ij and Q NL(pol) ij
, and the return-to-isotropy one

Q NL(RTI) ij
, which should be respectively added to S NL(iso) , S NL(dir) ij , S NL(pol) ij , and S (RTI) ij , the latter being defined in Appendix 1. The isotropic term Q NL(iso) is conservative, and one has

∞ 0 2 Q NL(iso) δ ij 3 + Q NL(dir) ij + Q NL(pol) ij -Q NL(RTI) ij dk = 0. ( 30
)
In fact, we derived analytically Q

NL(RTI) ij

only to check the latter property, which is a good way to identify calculations errors. The impact of quadratic anisotropic contributions on the non-linear transfers is revealed in Figure 2(a). One can note that the global shape of the transfers is preserved and that there is a rather small change in intensity, mainly at large scales which is expected since this is where anisotropy is dominant. The isotropic and polarisation parts of the transverse component () 33 are more affected by the quadratic contributions than the directional part.

In a recent work by Soulard and Gréa [START_REF] Soulard | Influence of zero-modes on the inertial-range anisotropy of Rayleigh-Taylor and unstably stratified homogeneous turbulence[END_REF], a zero-mode analysis was performed on a distinct anisotropic EDQNM model restricted to axisymmetric configurations. It was notably found that zero modes of higher-order harmonics, which could be compared to the quadratic anisotropic contributions here, for example, become negligible for increasing k. This is fully consistent with what is obtained in Figure 2(a), namely that the quadratic anisotropic contributions are completely negligible at small scales.

Regarding the global anisotropy indicator b ij in Figure 2(b): the asymptotic values are almost the same with and without the quadratic anisotropic contributions. One can observe a slight increase of the streamwise anisotropy and decrease of the transverse one. The kinetic energy exponential growth rate is not modified by the presence of the quadratic anisotropic contributions.

In conclusion, the quite heavy analytical calculations which led to the quadratic anisotropic contributions in the non-linear transfers do not provide significant changes and lengthen the numerical simulations. This fully justifies that they were neglected so far. This analysis also indicates that in order to improve the modelling of anisotropy, it is not efficient to refine the non-linear transfers, which remains a noticeable feature: one should rather improve the production terms, which is the topic of the next part.

Modelling with the fourth-order expansion of E and Z

In this section, we aim at improving the anisotropic EDQNM modelling originally developed in MCS by considering the fourth-order terms in the expansions for E and Z in Equations ( 9) and [START_REF] Isaza | On the asymptotic behaviour of large-scale turbulence in homogeneous shear flow[END_REF], respectively. The objective is rather simple, namely to show that taking into account more spherical harmonics improves the results of the model for shear flows, and more precisely decreases the value of the kinetic energy exponential growth rate γ . To this end, the theoretical developments necessary to the fourth-order expansion are exposed, and then the new spherically averaged equations are derived.

The operators of the fourth-order expansion

The first step is to determine the two fourth-order operators P ijpq and N ijpq which, when contracted with E and Z respectively, extract only the fourth-order tensors U (dir) 4 ijpq and U (pol)4 ijpq of Equations ( 9) and [START_REF] Isaza | On the asymptotic behaviour of large-scale turbulence in homogeneous shear flow[END_REF], and erase the H

(dir) ij and H (pol) ij
contributions. This is similar to what was done for the second-order expansion in [START_REF] Cambon | Anisotropic developments for homogeneous shear flows[END_REF], where the operators are, according to Equations ( 11) and ( 12)

2EH (dir) ij = S k E δ ij 3 -α i α j d 2 k, 2EH (pol) ij = S k (ZN i N j ) d 2 k. ( 31
)
The operator P ijpq is given in the recent study of Rubinstein et al. [START_REF] Rubinstein | Scalar and tensor spherical harmonics expansion of the velocity correlation in homogeneous anisotropic turbulence[END_REF]. On the other hand, N ijpq was not and is more complex to handle since it involves both the normalised wavevector α = k/k and the helical mode N. After some algebra, one gets

P ijpq = α i α j α p α q - 1 7 (δ ij α p α q + 5perm.) + 1 35 (δ ij δ pq + δ ip δ jq + δ iq δ jp ), (32)
N ijpq = (N i N j α p α q + N p N q α i α j + 4perm.) - 1 7 (δ ij N p N q + 5perm.), (33)
where 'perm.' refers to the other permutations of indices i, j, p, and q not written explicitly for the sake of brevity. Thanks to these two operators, one defines the fourth-order spectral anisotropy descriptors H (dir) ijpq and H in Equation [START_REF] Tavoularis | Asymptotic laws for transversely homogeneous turbulent shear flows[END_REF], as

2E(k, t)H (dir) ijpq (k, t) = S k E(k, t) P ijpq (k) d 2 k = 24 945 E(k, t)U (dir)4 ijpq (k, t), ( 34
) 2E(k, t)H (pol) ijpq (k, t) = S k (Z(k, t)N ijpq (k)) d 2 k = 4 21 E(k, t)U (pol)4 ijpq (k, t), ( 35
)
so that the fourth-order expansion of E and Z reads

E(k, t) = E 0 1 -15H (dir) ij (k, t)α i α j + 945 12 H (dir) ijpq (k, t)α i α j α p α q , ( 36
)
Z(k, t) = 1 2 E 0 5H (pol) ij (k, t) + 21 2 H (pol) ijpq (k, t)α p α q N * i N * j . ( 37
)
At this point, some words need to be said about the properties of the new fourth-order tensors and related assumptions. One can remark that at the fourth order, the polarisation part

H (pol)
ijpq is contracted with both normalised wavevectors α p α q and helical modes N * i N * j , unlike the directional part

H (dir)
ijpq which is only contracted with α i α j α p α q . Hence, H

(dir)

ijpq is fully symmetric in its indices, as reported in [START_REF] Rubinstein | Scalar and tensor spherical harmonics expansion of the velocity correlation in homogeneous anisotropic turbulence[END_REF]. The latter property of full symmetry is thus a priori not verified by H (pol) ijpq , but we nonetheless made this assumption. Otherwise, multiple independent variables related to the fourth-order polarisation appear in the calculations which cannot be dealt with in this framework. In addition, we assumed that H (pol) ijpq is trace-free, meaning that any contraction of two indices yields zero. For H (dir) ijpq , this is verified according to [START_REF] Rubinstein | Scalar and tensor spherical harmonics expansion of the velocity correlation in homogeneous anisotropic turbulence[END_REF], but for H (pol) ijpq this is a supplementary assumption. These two assumptions for H (pol) ijpq were used to obtain Equation [START_REF] Souza | The structure of highly sheared turbulence[END_REF] and thus Equation [START_REF] Pumir | Persistent small scale anisotropy in homogeneous shear flows[END_REF].

Consequently, there is room for improvement for the modelling of the fourth-order expansion of Z. Nevertheless, the present assumptions for H (pol) ijpq allow to close the different equations, and the numerical results presented in Section 3.5 show that these hypotheses are acceptable. It is not clear what different rules could be used for H (pol) ijpq and this is an open question: one could think of a restriction with distinct i ↔ j and p ↔ q symmetries; the contractions i = p, i = q, j = p, j = q could also be discussed.

With the fourth-order expansion, the modelled expression (15) of the spectral Reynolds tensor Rij now reads

Rij = E 0 P ij 1 -15H (dir) pq α p α q + 945 12 H (dir) rspq α r α s α p α q Re2 ij + Re4 ij + 1 2 E 0 5H (pol) pq + 21 2 H (pol)
pqrs α r α s (2P ip P jq + P ij α p α q )

Rz2 ij + Rz4 ij . ( 38
)
It is worth noting that the spherical average of Rij remains unchanged whether the fourthorder contributions are considered or not, i.e. Equation ( 13) is still true. The final step of the fourth-order modelling is to determine analytically (i) the linear and non-linear transfers associated with EH (dir) ijpq and EH . For the sake of readability, no intermediate steps are given here since the expressions are very lengthy, but they can be found in [START_REF] Briard | Modelling of transport in homogeneous turbulence[END_REF].

Fourth-order productions terms

In this part, we aim at determining the explicit expressions of the production terms, linear with the mean-velocity gradient matrix A ij , linked to the fourth-order contributions in the expansions [START_REF] Brethouwer | The effect of rotation on rapidly sheared homogeneous turbulence and passive scalar transport. Linear theory and direct numerical simulation[END_REF] and ( 37) of E and Z, respectively. They are of three kinds: the fourthorder contributions in the equations of EH

S L(dir4) ij (k, t) = -3A + ln EH (dir) ijln -A + ln ∂(kEH (dir) ijln ) ∂k + 1 6 A + ln EH (pol) ijln , (39)
S L(pol4) ij (k, t) = -A + ln EH (dir) ijln + 1 2 A + ln EH (pol) ijln + 1 6 A + ln ∂(kEH (pol) ijln ) ∂k , ( 40
)
where the superscripts (dir4) and (pol4) refer to the contributions of the fourth-order expansion. Consequently and for the sake of clarity, the linear transfers of Equations ( 19) and ( 20) are renamed S L(dir2) ij and S

L(pol2) ij

, so that the total linear transfers for EH ijpq , which can also be divided into two terms, resulting respectively from the second-and fourth-order contributions

S L(dir) ijpq = 1 4 S k A ln k l ∂ Rmm ∂k n -2M mn Rnm P ijpq d 2 k, ( 41
) S L(pol) ijpq = 1 4 S k A ln k l ∂ Rrs ∂k n -M rn Rns -M sn Rnr N * r N * s N ijpq d 2 k. ( 42
)
First, we determine the impact of the second-order contributions EH ijpq . For this purpose, it is convenient to introduce two symmetric and trace-free operators H (2,e) ijpq [EH () ] and H (2,z) ijpq [EH () ], which refer to terms involving respectively Re2 ij and Rz2 ij of Equation ( 15), namely 43)

H (2,e) ijpq [EH () ] = E 8 5 A + ln H () ln (δ ij δ pq + δ ip δ jq + δ iq δ jp ) -4[δ ij (A + lp H () lq + A + lq H () lp ) + δ pq (A + li H () lj + A + lj H () li ) + A + lp (H () lj δ iq + H () li δ jq ) + A + lq (H () lj δ ip + H () li δ jp ) + H () lp (A + jl δ iq + A + il δ jq ) + H () lq (A + jl δ ip + A + il δ jp )] +14(A + ij H () pq + A + pq H () ij + A + ip H () jq + A + iq H () jp + A + jp H () iq + A + jq H () ip ) , (
H (2,z) ijpq [EH () ] = E[4A + ln H () ln (δ ij δ pq + δ ip δ jq + δ iq δ jp ) -10[δ ij (A + lp H () lq + A + lq H () lp ) + δ pq (A + li H () lj + A + lj H () li ) + A + lp (H () lj δ iq + H () li δ jq ) + A + lq (H () lj δ ip + H () li δ jp ) + H () lp (A + jl δ iq + A + il δ jq ) + H () lq (A + jl δ ip + A + il δ jp )] + 35(A + ij H () pq + A + pq H () ij + A + ip H () jq + A + iq H () jp + A + jp H () iq + A + jq H () ip )]. ( 44
)
After some algebra, one eventually obtains

S L(dir2) ijpq = 1 441 - 1 6
H (2,e) ijpq [EH (pol) ] + 4H (2,e) ijpq [EH (dir) ] -H (2,e) ijpq [∂ k (kEH (dir) )] , (45)

S L(pol2) ijpq = 1 441
(6H (2,z) ijpq [EH (dir) ] -4H (2,z) ijpq [EH (pol) ] + H (2,z) ijpq [∂ k (kEH (pol) )]).

(46)

Consistently with the previous calculations, the superscripts (dir2) and (pol2) refer to the contributions of the second-order expansion. Finally, we proceed similarly to determine the impact of the fourth-order contributions on the linear terms of EH (dir) ijpq and EH (pol) ijpq . It is also convenient to define a final operator H (4) ijpq [EH () ] which reflects the terms coming from Re4 ij and Rz4 ij of Equation ( 38):

H (4) ijpq [EH () ] = E A + li H () jlpq + A + lj H () ilpq + A + lp H () ijlq + A + lq H () ijlp - 2 7 A + ln (δ ij H () lnpq + δ pq H () lnij + δ ip H () lnjq + δ iq H () lnjp + δ jp H () lniq + δ jq H () lnip ) . ( 47
)
After intricate algebra involving the spherical integration of terms containing 10 normalised wavevectors α, one finally gets

S L(dir4) ijpq = 1 11 2H (4) ijpq [∂ k (kEH (dir) )] -H (4) ijpq [EH (dir) ] + 3 5
H (4) ijpq [EH (pol) ]

+ E(A - il H (dir) jpql + A - jl H (dir) ipql + A - pl H (dir) ijlq + A - ql H (dir) ijlp ), ( 48
) S L(pol4) ijpq = 1 11 4 5 H (4) ijpq [∂ k (kEH (pol) )] - 2 5
H (4) ijpq [EH (pol) ] + 54H (4) ijpq [EH (dir) ]

+ 3 5 E(A - il H (pol) jpql + A - jl H (pol) ipql + A - pl H (pol) ijlq + A - ql H (pol) ijlp ). ( 49
)
The total fourth-order directional and polarisation linear production terms are consequently S L(dir)

ijpq = S L(dir2) ijpq + S L(dir4) ijpq
, and similarly for S L(pol) ijpq .

Fourth-order non-linear transfer terms

The fourth-order non-linear directional and polarisation transfers S NL(dir) ijpq and S NL(pol) ijpq are defined as

S NL(dir) ijpq (k, t) = 1 2 S k T E (k, t)P ijpq (k) d 2 k, ( 50
) S NL(pol) ijpq (k, t) = 1 2 S k (T Z (k, t)N ijpq (k)) d 2 k. ( 51
)
Discarding the quadratic contributions of anisotropy as for the second-order expansions, one obtains the important fact that S NL(dir)

S NL(dir) ijpq = 2 k θ kpq π 2 k 2 p 2 q(xy + z 3 )E 0 [E 0 (H (dir) ijpq (35z 4 -30z 2 + 3) + H (dir)
ijpq (35y 4 -30y 2 + 3)) -E 0 (8H (dir) ijpq + H (dir) ijpq (35y 4 -30y 2 + 3))] dp dq

+ 2 3 k θ kpq π 2 k 2 p 2 qE 0 [(xy + z 3 )(1 -y 2 )(7y 2 -1)(E 0 -E 0 )H (pol) ijpq + z(1 -z 2 ) 2 (1 -7z 2 )E 0 H (pol) ijpq ] dp dq, ( 52
)
S NL(pol) ijpq = 4 k θ kpq π 2 k 2 p 2 qE 0 [(xy + z 3 )(E 0 H (pol) ijpq (1 -6z 2 + 7z 4 ) -4E 0 H (pol) ijpq ) + z(5 -7z 2 )(y 2 -z 2 )E 0 H (pol) ijpq + (E 0 -E 0 )H (pol) ijpq z(z 2 -1)(1 -6y 2 + 7y 4 ) + xy(5 -7y 2 )(1 -z 2 )E 0 H (pol) ijpq ] dp dq + 60 k θ kpq π 2 k 2 p 2 qE 0 z(1 -z 2 ) × [(1 -8y 2 + 7y 4 )(E 0 -E 0 )H (dir) ijpq + (1 -8z 2 + 7z 4 )E 0 H (dir) ijpq ] dp dq, ( 53
)
where k is the domain where k, p, and q are the lengths of the sides of the triangle formed by the triad, and compact notations are used, E 0 = E(k)/(4πk 2 ), E 0 = E(p)/(4πp 2 ), and E 0 = E(q)/(4π q 2 ).

Final spherically averaged equations

As a consequence of the previous calculations, the evolution equations of the fourth-order anisotropic descriptors EH (dir) ijpq and EH (pol)

ijpq read ∂ ∂t + 2νk 2 E(k)H (dir) ijpq (k) = S L(dir2) ijpq (k) + S L(dir4) ijpq (k) + S NL(dir) ijpq (k), ( 54
)
∂ ∂t + 2νk 2 E(k)H (pol) ijpq (k) = S L(pol2) ijpq (k) + S L(pol4) ijpq (k) + S NL(pol) ijpq (k). ( 55
)
Moreover, the original equations ( 19) and [START_REF] Pumir | Turbulence in homogeneous shear flows[END_REF] of EH (dir) ij and EH

(pol) ij derived in MCS are modified accordingly into At this point, it is crucial to understand that the retro-action of the fourth-order contributions on the second-order ones is uniquely done through the linear production terms S

∂ ∂t + 2νk 2 E(k, t)H (dir) ij (k) = S L(dir2) ij (k) + S L(dir4) ij (k) + S NL(dir) ij (k), ( 56
)
∂ ∂t + 2νk 2 E(k)H (pol) ij (k) = S L(pol2) ij (k) + S L(pol4) ij (k) + S NL(pol) ij (k). ( 57
)
L(dir4) ij and S L(pol4) ij
, and that the impact of the second-order contributions on the fourth-order ones is uniquely done through the linear terms S For future reference, we call this approach the advanced anisotropic EDQNM modelling. In the next section, some numerical results using the fourth-order expansions are presented.

Numerical results using the fourth-order expansion

In this part, some numerical results are proposed to illustrate the interest of considering the fourth-order expansion of E and Z in spherical harmonics. The numerical set-up is the same as in [START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF]: k min = 10 -10 k L (t = 0) and k max = 10 6 k η , where k L and k η = ( /ν 3 ) 1/4 are, respectively, the integral and Kolmogorov wavenumbers. For numerical integration of the five Lin equations of the advanced anisotropic EDQNM modelling, a logarithmic discretisation in wavenumbers is used, such that k i+1 = rk i with r = 10 1/f , f = 17 points per decade, with a third-order implicit Runge-Kutta scheme for the viscous term. The initial Reynolds number based on the Taylor scale is approximately Re λ (0) 10, and an initial isotropic fully developed kinetic energy spectrum, given by [START_REF] Pope | Turbulent flows[END_REF], is chosen.

The main consequence of the fourth-order contributions is that the kinetic energy exponential growth rate is decreased in Figure 3(a), from γ = 0.33 with the MCS model, to γ = 0.28 with the present advanced anisotropic modelling. Even though the value of γ remains rather large, the significant decrease by 15% with the fourth-order expansion of E and Z perfectly proves that taking into account more spherical harmonics goes into the good direction, i.e. diminishes γ towards smaller values, consistently with DNS.

The joint result is, in Figure 3(b), the decrease of b 13 from 0.215 to 0.18, which is a noteworthy feature as well. One can further remark that on the contrary, the fourth-order contributions increase |b 11 | and |b 33 |, which is expected. Indeed, taking into account more harmonics reduces the loss of information due to the spherical integration by restoring part of the anisotropic angular information. As a consequence, the strong anisotropy of the shear flow between the streamwise and transverse directions is better captured.

To better understand the impact of the fourth-order contributions, we investigate the production terms in Figure 4, with the directional and polarisation parts of the streamwise, transverse and cross components in Saffman turbulence. In this figure, the black curves represent simulations with the fourth-order contributions, at Re λ = 9×10 3 . Whereas the grey curves indicate simulations with only the second-order expansion, as in MCS. Since the Reynolds number increases faster in the latter case, the results are presented at St = 43 where Re λ 9×10 3 as well.

One can remark that the effects are different for the directional and polarisation parts: indeed, the fourth-order contributions tend to decrease the intensity of the directional linear terms for the streamwise () 11 and transverse () 33 components, while increasing it for the () 13 component. The opposite happens for the linear polarisation terms. In particular, the strongest difference is observed for the transverse directional transfer S L(dir) [START_REF] Tavoularis | Experiments in nearly homogenous turbulent shear flow with a uniform mean temperature gradient. Part 1[END_REF] , which is positive without the fourth-order contributions, and becomes mostly negative with them.

As a conclusion for the fourth-order contributions, the main consequence is a 15% decrease of the kinetic energy exponential growth rate, namely from γ = 0.33 to γ = 0.28, which is quite significant, without rendering the numerical simulations more costly. It further justifies a posteriori the assumptions made to establish the fourth-order expansion of E and Z.

Odd-order terms in the expansion of Z

In this section, we discuss the possibility of odd-order terms in the spherical harmonics expansion of the polarisation anisotropy Z. Recall that because of the symmetry of Rij , the kinetic energy density E has only even-order contributions [START_REF] Cambon | Anisotropic developments for homogeneous shear flows[END_REF][START_REF] Rubinstein | Scalar and tensor spherical harmonics expansion of the velocity correlation in homogeneous anisotropic turbulence[END_REF]. There are two reasons why we wish to further consider odd-order terms in the expansion of Z: (i) Odd-order expansions could improve the modelling of the 2iZ CH term in the evolution equation [START_REF] Cambon | Third-order statistics and the dynamics of strongly anisotropic turbulent flows[END_REF] of Z; (ii) Recent results show that the main difference between the MCS model (and thus the truncated expansions of E and Z) with an exact treatment of linear terms lies in the polarisation anisotropy [START_REF] Zhu | Rotating shear-driven turbulent flows: towards a spectral model with angle-dependent linear interactions[END_REF]. Unfortunately, these developments are not conclusive yet because they cannot be sustained numerically. Nevertheless, they are presented here for information purposes, and to underline that perhaps different methods than simple tensorial tools are required to model higher-order anisotropy features.

Determination of the third-order expansion of Z

Up to the fourth order, the expansion of Z can be written as

Z(k) = E 0 2 5H (pol) ij (k) + iU (pol)3 ijk (k)α k + 21 2 H (pol) ijpq (k)α p α q N * i (k)N * j (k), (58)
where

U (pol)3 ijk
is a tensor which is assumed to verify, for simplicity reasons and consistency with previous developments, full symmetry under any change of indices, and to be zero when two indices are equal. Note that the expression of the third-order contribution differs from the one in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF] because of the imaginary number i. The latter is crucial, otherwise the third order of Z never contributes in the equations. Two features are needed to prove this statement. First, one requires the following result when computing the polarisation part R(pol

) ij = [ZN i N j ]: N i N j N * p N * q = [P ip P jq + P iq P jp -P ij P pq ] - 1 2 iα a [P jq ipa + P ip jqa + P iq jpa + P jp iqa ], ( 59
)
where ijk is the Levi-Civita permutation tensor. The second one is that the spherical average of an odd number of normalised wavevectors α i is zero [START_REF] Pope | Turbulent flows[END_REF]. Consequently, without the i, [ZN i N j ] has an odd number of α i so that it vanishes with the spherical average. Therefore, the present expansion (58) of Z corrects the Equation (3.15) of [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF]. Then, as for the second-and fourth-order terms, one needs an operator which gives only the third-order contribution U (pol)3 ijk and erases the others. This operator reads N ijk = α k N i N j + α j N i N k + α i N j N k , and we further define

H (pol) ijk as 2E(k, t)H (pol) ijk (k, t) = S k [Z(k, t)N ijk (k)] d 2 k = 2 7 E(k, t)U (pol)3 ijk (k, t). ( 60
)
The third-order expansion of Z does not modify the spectral tensor φ ij given in Equation [START_REF] Cambon | Spectral approach to non-isotropic turbulence subjected to rotation[END_REF] Note that using some properties of the helical modes [START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF], the third-order expansion 2Z = 7iH

(pol) pql α l N * p N * q can also be written 2Z = 7H

(pol) pql α l α r prs N * q N * s . Now that the third-order expansion of Z has been properly defined, one needs to compute the linear production terms and non-linear transfers associated to EH . For simplicity, we do not consider in this work the linear coupling between the third and fourth orders.

Remark: It is worth noting that the third-order expansion of polarisation can be related to the stropholysis tensor [START_REF] Kassinos | One-point turbulence structure tensors[END_REF] defined as

Q ijk (t) = ipq α p α k Rjq (k, t) d 3 k = - ijp α p α k E(k, t) d 3 k + α k (Z(k, t)N i (k)N j (k)) d 3 k. ( 62
)
The expression of Q ijk as a function of the second-order anisotropic tensors can be found in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF]. Interestingly, the symmetric stropholysis Q * ijk = (Q ijk + Q ikj + 4 perm.)/6 does not contain directional anisotropy (because ijp =jip ) and depends only on polarisation. More specifically, only odd-order terms of the Z expansion can contribute in Q * ijk , which further justifies the interest of odd-order contributions.

Third-order productions terms

In this part, and following the procedure detailed for the fourth-order contributions, we aim at determining the explicit expressions of the production terms linked to the thirdorder contributions in the expansion (58) of Z. As before, there are three different kinds of terms: the third-order contributions in the equations of EH (dir) (k,t

+ 2EA - ln H () np E[δ ij klp + δ ik jlp + δ jk ilp ]. (67)
After some algebra, one gets

S L(pol2) ijk = - 1 7 H (2,+) ijk [EH (dir) ] + 1 21 [H (2,-)
ijk [EH (pol) ] -H (2,+) ijk [EH (pol) ]]

+ 1 42
H (2,+) ijk [∂ k (kEH (pol) )]. ( 68
)
Now, we proceed similarly to determine the linear term of EH (pol) ijk as a function of the third-order contributions. Eventually, one gets

S L(pol3) ijk = 1 3 [A - il H (pol) jlk + A - jl H (pol) ilk + A - kl H (pol) ijl ], (69)
where notably the symmetric part of the mean-velocity gradient matrix does not intervene. Finally, the total third-order polarisation transfer is S

L(pol) ijk = S L(pol2) ijk + S L(pol3) ijk .

Third-order non-linear transfer

The third-order polarisation non-linear transfer is defined as

S NL(pol) ijk (k, t) = 1 2 S k (T Z (k, t)N ijk (k)) d 2 k. ( 70
)
For the sake of simplicity, quadratic anisotropic contributions are still discarded, as for the second and fourth orders: therefore, only the third-order terms contribute in S NL(pol) ijk .

Furthermore, because H (pol)

ijk is symmetric and trace-free, it follows that third-order expansions vanish in S NL(dir) ij and S NL(pol) ij . Then, injecting the third-order expansion into the expression [START_REF] Gualtieri | Scaling laws and intermittency in homogeneous shear flow[END_REF] of T Z gives the spherically averaged non-linear polarisation transfer

S NL(pol) ijk = 4π 2 k θ kpq k 2 p 2 qE 0 [-4(xy + z 3 )E 0 H (pol) ijk + H (pol) ijk (1 -z 2 )(2z(1 -2y 2 )(E 0 -E 0 ) + xy(1 -3y 2 )E 0 ) + H (pol) ijk E 0 (2(xy + z 3 )(2z 2 -1) -z(3z 2 -1)(y 2 -z 2 ))] dp dq.
(71)

Final spherically averaged equations

The evolution equation of the third-order anisotropic descriptor EH (pol)

ijk reads ∂ ∂t + 2νk 2 E(k)H (pol) ijk (k) = S L(pol2) ijk (k) + S L(pol3) ijk (k) + S NL(pol) ijk (k). (72)
Moreover, the original equations ( 19) and [START_REF] Pumir | Turbulence in homogeneous shear flows[END_REF]

∂ ∂t + 2νk 2 E(k, t)H (dir) ij (k) = S L(dir2) ij (k) + S L(dir3) ij (k) + S NL(dir) ij (k), (73)
∂ ∂t + 2νk 2 E(k)H (pol) ij (k) = S L(pol2) ij (k) + S L(pol3) ij (k) + S NL(pol) ij (k). ( 74
)
Note that the retro-action of the third-order contributions on the second-order ones is uniquely done through the linear transfers S

L(dir3) ij and S L(pol3) ij
, and that the impact of the second-order contributions on the third-order ones is uniquely done through the linear transfers S L(dir2) ijk and S

L(pol2) ijk

. This is completely similar to what was done for the fourth order.

These developments for the third order of the polarisation anisotropy Z are not conclusive yet because the numerical code cannot handle satisfactorily the corresponding equations (72), (73), and (74). This is very likely because of the assumptions we made to define H (pol) ijl , namely the full symmetry between indices. Similar rules were nevertheless successfully used previously for H (dir) ijpq and H (pol) ijpq in Section 3. This suggests that more complex features should be taken into account for odd-order terms with more sophisticated tools.

Conclusions, perspectives, and summary

This work focuses on the spectral modelling of shear flows using an adapted eddy-damped quasi-normal Markovian (EDQNM) closure for homogeneous anisotropic turbulence developed by Mons, Cambon, and Sagaut (MCS) [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF] and further applied in [START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF]. The present paper should be considered as the continuity of the two previous publications. The main objectives were twofold: (i) to explain why in MCS the exponential growth rate of kinetic energy γ is numerically independent of the large scales initial conditions, namely the infrared slope σ of the kinetic energy spectrum E(k < k L ) ∼ k σ . And (ii) to show that γ can actually be reduced by improving the modelling of anisotropy. Two secondary aspects were analysed as well: (iii) to quantify the effects of quadratic anisotropic contributions in the non-linear transfers, an assumption which was not verified in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF]. And (iv) to address the modelling of odd-order terms in the polarisation anisotropy expansion, a work for now inconclusive but still ongoing.

For the first point, it was erroneously thought by the authors in [START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF] that such an independence could be a consequence of the spherical average. It has been shown here analytically in Section 2 using self-similarity arguments that γ is almost independent of σ because the most amplified eigenvalue of the linear operator also depends on σ . This balances the σ -dependence of the self-similar expression of E(k, t) at large scales. Interestingly, the anisotropic EDQNM model recovers well for unstably stratified homogeneous turbulence (USHT) the strong dependence of the kinetic energy exponential growth rate with the infrared slope σ [START_REF] Briard | Anisotropic spectral modelling for unstably stratified homogeneous turbulence[END_REF]: this is because the most amplified eigenvalue of the USHT linear operator does not depend on σ . This analysis consequently shows nicely that the mechanisms at the origin of the exponential growth in shear and USHT are intrinsically different. This is the first important finding of this work.

For the second point, we aimed at improving the MCS modelling for pure shear flows, whose main issue is that the kinetic energy exponential growth rate γ is larger than common value reported in DNS and experiments [START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF]. Consequently, the objective here was to show that γ could be reduced by better capturing the complex distribution of anisotropy. This was done by computing analytically the fourth-order terms in the spherical harmonics expansion of the kinetic energy density E and polarisation Z (in MCS the truncation is at the second order). These developments led to the definitions of two new spherically averaged anisotropic descriptors H (dir) ijpq and H (pol) ijpq whose evolution equations are similar to the ones of the second-order descriptors H . For simplicity reasons, we chose some symmetry rules for H (pol) ijpq whose consequences are that the coupling between the second-order and fourth-order equations is purely linear, i.e. through the production terms.

Numerically, this advanced modelling for the anisotropy is satisfactory because it decreases the exponential growth rate of kinetic energy from γ = 0.33 (MCS) to γ = 0.28, without increasing the computational time. Conjointly, the cross-component of the global anisotropy indicator goes down to |b 13 | = 0.18. Hence, we indeed showed that taking into account more spherical harmonics in the expansion reduces the kinetic energy exponential growth rate by better capturing the anisotropy distribution. This is the second theoretical contribution of this work.

Regarding the third point, quadratic anisotropic contributions were explicitly presented in Appendix 2 and were shown to be numerically negligible in Section 2.5. This notably shows that to improve the description of anisotropy, it is more efficient to refine the linear production terms than the non-linear ones: a useful trail for future works.

Finally, the fourth and last aspect of this work was to address the modelling of oddorder terms in the expansion of polarisation anisotropy Z. All the analytical calculations were done but are inconclusive in terms of numerical simulations. This is probably due to assumptions we made regarding the symmetry properties of the third-order polarisation tensor H (pol) ijk . It was nevertheless shown that odd-order contributions of Z can be related to the stropholysis [START_REF] Kassinos | One-point turbulence structure tensors[END_REF] and as such deserve further investigations.

As perspectives, it could be of great theoretical interest to determine, or at least estimate, at which order one should perform the expansion into spherical harmonics to properly capture all the anisotropic features of a given homogeneous anisotropic turbulent flow. This task could be rather complex, and even if it is out of the scope of the present work, the authors would like to emphasis some aspects.

First, one would have to solve numerically the linear production terms of Equations ( 7) and ( 8) and to project E and Z on an appropriate basis of deviatoric tensors such as H () ij , H () ijk , H () ijpq , . . .: in particular, this would require to determine explicitly the higher-order projectors which may be extremely tedious given the complexity of the algebra involved already at the fourth-order.

Second, it is very likely that this required number of spherical harmonics would depend on both the mean-velocity gradient matrix A ij and the accumulated anisotropy St. Indeed, for the former, it appeared in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF][START_REF] Briard | Decay and growth laws in homogeneous shear turbulence[END_REF] that the MCS model is accurate enough when A ij is symmetric, such as in an axisymmetric expansion or contraction, or in a plane distortion, so that the second-order expansion is sufficient. However, in shear flows, the fact that A ij is not symmetric anymore seems to imply that more harmonics are needed. Furthermore, in shear flows, the MCS model works fine at small St (see the comparison with rapid distortion theory in [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically-averaged descriptors[END_REF]), which means that the number of spherical harmonics needed would increase with St. Consequently, these two aspects make the task of predicting the order at which one should truncate the expansion of E and Z into spherical harmonics quite intricate. An alternative solution could be to combine the anisotropic EDQNM modelling for the nonlinear transfers with an exact treatment of the linear production terms, which is currently the task of Zhu et al. in [START_REF] Zhu | Rotating shear-driven turbulent flows: towards a spectral model with angle-dependent linear interactions[END_REF].

Finally, we propose in Figure 5 a synthetic view of the work accomplished so far, the various achievements, the methodology, the assumptions made, and the remaining objectives as well. Notably, Figure 5 shows that the short-term objective of our work is to improve the modelling of linear production terms with in particular the third-order expansion of Z. This will provide the turbulence community a model which can handle accurately secondorder moments statistics in shear flows, with very few numerical resources and at large Reynolds numbers.

and A = 0.36. The production terms are, with A ij being the mean-velocity gradient matrix where () + and () -denote the symmetric and antisymmetric parts, ). (A9)

S L(iso

Appendix 2. Explicit quadratic contributions of anisotropy for non-linear transfers

Here are the complete expressions of the quadratic anisotropic contributions in the non-linear transfers for the anisotropic EDQNM modelling with truncation at the second order (MCS). The following convenient notation is used

H (dir

Figure 1 .

 1 Figure 1. (a) Theoretical prediction for the kinetic energy exponential growth rate γ th = 2 max (σ )/(σ + 3), where max is provided in Equation (27). The numerical values in grey represent the expected γ th for the classical integer values of the infrared slope σ , whereas the dash-dot line represents the value obtained by the anisotropic EDQNM modelling γ = 0.33. (b) Budget terms of the equation of E(k, t) given in Equation (18), for σ = 2 and Re λ (St = 50) = 2.4×10 4 , normalised by , along with the integral, shear and Kolmogorov wavenumbers k L , k S , and k η . (c) Ratio of the linear to non-linear transfers, for the same configuration as (b). The total transfer is given by S

Figure 2 .

 2 Figure 2. Quadratic anisotropic contributions in the non-linear transfers at St = 50, with σ = 2. Grey curves represent simulations without the quadratic anisotropic contributions. (a) Budget terms along with the integral and Kolmogorov wavenumbers k L and k η . (b) Global anisotropy indicator b ij .

  (pol) ijpq , and (ii) the retro-action of the fourth-order contributions on the equations of E, EH (dir) ij , and EH (pol) ij

  ij, in the equations of EH

  ijpq , and the second-order contributions in the equations of EH(dir) ijpq and EH (pol) ijpq . From the properties of H (dir) ijpq and H (pol) ijpq , it follows that there are no fourth-order contributions in the equation for E. First, the contributions of the fourth-order expansion in the equations of EH (dir) ij and EH (pol) ij are, after some algebra,

  define the linear directional and polarisation production terms S

  ijpq and S NL(pol) ijpq only depend on the fourth-order contributions. In addition, since H (dir) ijpq and H (pol) ijpq are symmetric and trace-free, it follows that fourth-order contributions vanish in S NL(dir) ij and S NL(pol) ij . The calculations yielding to the final expressions of S NL(dir) ijpq and S NL(pol) ijpq are similar to the ones for S NL(dir) ij and S NL(pol) ij , but somehow lengthier:

Figure 3 .

 3 Figure 3. Effects of the fourth-order expansion on the growth of the kinetic energy K(t) and the global anisotropy indicator b ij . (a) K(t) for σ = 2 and σ = 4. (b) b ij for σ = 2. In both figures, the grey curves indicate a simulation without the fourth-order contributions for comparison.

  -the two new for EH (dir) ijpq and EH (pol) ijpq and the two adapted for EH (dir) ij and EH (pol) ij -represent the second important theoretical contribution of this work.

Figure 4 . 33 .

 433 Figure 4. Effects of the fourth-order expansion on the linear terms of EH (dir) ij and EH (pol) ij , for σ = 2,

  (pol) ijk , along with its retro-action on the evolution equations of EH (dir) ij and EH (pol) ij

  ijk , and finally the second-order contributions in the equations of EH (pol) ijk . There are no third-order contributions in the equation for E, and for simplicity reasons, we further discard the third-order contributions in the evolution equations of EH (dir) ijpq and EH (pol) ijpq , and the fourth-order contributions in the equation of EH (pol) ijk . First, the retro-action of the third-order expansion in the equations of EH (dir) ij and EH (pol) ij reads S L(dir3) ij

Figure 5 .

 5 Figure 5. Diagram presenting the methodology of the modelling and summarising the different levels of assumptions, the various achievements, and the ongoing studies as well. The 'path' toward the numerical simulations of the full E-Z equations with non-linear terms closed by EDQNM is broken since this would be almost as costly as DNS in terms of numerical resources.

  kpq is the characteristic time of the third-order correlations, defined more precisely in Appendix 1. Furthermore, quadratic contributions of anisotropy (terms like H

				(dir) il	H (dir) jl	,
	H	(dir) il	H	(pol) jl

  , which is still expressed as a function of H(dir)

			ij	and H ij (pol)	only, and gives a new
	contribution in the modelled spectral Reynolds tensor Rij :
	Rz3 ij (k, t) =	7 2	E 0 (k, t)H pql (k, t)α l α n ( ipn P jq (k) + jqn P ip (k)). (pol)	(61)

  of EH

	(dir) ij	and EH	(pol) ij	derived in MCS are
	modified accordingly into

  ij = {H (dir) , H (dir) } ij .The quadratic anisotropic isotropic transfer term is conservative and readsQ NL(iso) (k, t) = 20 k θ kpq π 2 k 2 p 2 qE 0 [2xy(1z 2 )H

	li	)	H lj (dir)	+ H	(dir) lj	H li (dir)	-	2 3	H (dir) ln	(dir) ln H
											(pol) ln	(E 0 H ln (pol)	-E 0 H ln ) (pol)
			+ 6(xy + z 3 )(2H	(dir) ln	((3x 2 -1)E 0 H ln (dir)	-(3y 2 -1)E 0 H ln ) (dir)
			-H	(pol) ln
											(pol) ln )
			-6H (dir)

δ ((1x 2 )E 0 H (dir) ln -(1y 2 )E 0 H (dir) ln )) + z(z 2 -1)(H (pol) ln ((1 + x 2 )E 0 H (pol) ln -(1 + y 2 )E 0 H ln ((1x 2 )E 0 H (pol) ln -(1y 2 )E 0 H (pol) ln ))] dp dq.

(A10)

Appendices

Appendix 1. Linear and non-linear transfers for the second-order expansion

In this appendix, the explicit expressions of the spherically averaged transfer terms are given. The non-linear ones are

From the total non-linear transfer, a return to isotropy (RTI) term can be extracted

All the non-linear transfers involve the characteristic time θ kpq of the triple correlations, which appears within the EDQNM approximation and is defined as

where μ(k, t) is the (isotropic) eddy-damping factor [START_REF] Lesieur | Turbulence in fluids[END_REF]

The quadratic anisotropic directional transfer term is

The quadratic anisotropic polarisation transfer term is

The quadratic anisotropic RTI term is (dir) , H (dir) } ij -(6y(x + yz) + 2z(1z 2 )){H (dir) , H (pol) } ij + (2y(x + yz) + 6z(1z 2 )){H (pol) , H (dir) } ij + (z(1z 2 )y(x + yz)){H (pol) , H (pol) } ij ) + k(xy + z)E 0 × ((6y(y + xz) -4(1z 2 )){H (dir) , H (pol) } ij + y(y + xz){H (pol) , H (pol) } ij )] dp dq. (A13)