1. Introduction
In recent decades, the numerical simulation of turbulent flows has become a prominent tool for the prediction of natural phenomena and the design of efficient industrial applications. In most practical situations, the number of degrees of freedom required to represent turbulence fully is intractable for current computational resources, and turbulence modelling is necessary to yield computationally feasible simulations. One of the most successful modelling approaches is the large-eddy simulations (LES) technique, which resolves explicitly the dynamics of turbulence above a prescribed filter scale, and substitutes the interactions with smaller scales by a subgrid-scale (SGS) model (Meneveau & Katz Reference Meneveau and Katz2000).
This paper focuses on incompressible turbulence away from solid boundaries. It is assumed that the filter scale is small enough compared to the forcing mechanisms of turbulence, and large enough compared to the viscous scale. In these conditions, the flow at the filter scale is statistically isotropic and governed by inertial forces, and the interactions between the resolved and subgrid scales lead to the turbulence energy cascade, a highly non-linear process that results in a net transfer of energy from the resolved to the subgrid scales. To yield a correct kinetic energy distribution in the resolved scales, and to reproduce the correct turbulent structure at the filter scale, SGS models should replicate the dynamics of the inertial cascade, not only in an average sense, but also locally.
The energy cascade goes on average from large to small scales, yet it is believed to be strongly bidirectional in a local sense. This is supported by the observation that the local SGS energy fluxes computed by filtering turbulent velocity fields in the inertial range very often have a direction opposite to the mean energy flux (Leslie & Quarini Reference Leslie and Quarini1979; Piomelli et al. Reference Piomelli, Cabot, Moin and Lee1991; Borue & Orszag Reference Borue and Orszag1998; Cerutti & Meneveau Reference Cerutti and Meneveau1998; Tao, Katz & Meneveau Reference Tao, Katz and Meneveau2002; Aoyama et al. Reference Aoyama, Ishihara, Kaneda, Yokokawa, Itakura and Uno2005; Cardesa et al. Reference Cardesa, Vela-Martín, Dong and Jiménez2015; Carter & Coletti Reference Carter and Coletti2018; Ballouz & Ouellette Reference Ballouz and Ouellette2018). This phenomenon, known as backscatter, is commonly identified with an inverse cascade that transports energy from the small to the large scales, and is widely considered a fundamental feature of inertial dynamics that must be modelled explicitly (Carati, Ghosal & Moin Reference Carati, Ghosal and Moin1995; Ghosal et al. Reference Ghosal, Lund, Moin and Akselvoll1995). However, constructing SGS models that faithfully reproduce backscatter has posed a huge challenge for researchers that has not been accomplished satisfactorily to date. As a consequence, the use of purely dissipative SGS models, which cannot produce backscatter, is very extended. Although some attempts have been made to design models that reproduce backscatter (Bardina, Ferziger & Reynolds Reference Bardina, Ferziger and Reynolds1980; Chasnov Reference Chasnov1991; Germano et al. Reference Germano, Piomelli, Moin and Cabot1991; Ghosal et al. Reference Ghosal, Lund, Moin and Akselvoll1995), these models are either impractical or introduce little advantage with respect to dissipative models at the expense of a higher computational complexity and extra tuning parameters.
Despite the conspicuousness of backscatter in the local SGS fluxes, evidence suggests that its relevance to the energy cascade is limited. Simulations performed with dissipative SGS models agree reasonably well on average with fully resolved direct numerical simulations (DNS) and experiments. For instance, dissipative SGS models reproduce the correct kinetic energy spectrum and decay rate of isotropic turbulence, despite their inability to produce backscatter (see Meneveau & Katz (Reference Meneveau and Katz2000) and references therein). Also, the multifractal scaling of the structure functions in the inertial range is well reproduced by the Smagorinsky (Reference Smagorinsky1963) model, which is purely dissipative (Linkmann, Buzzicotti & Biferale Reference Linkmann, Buzzicotti and Biferale2018). This suggests that backscatter is not necessary to keep the energy balance of the resolved scales or to reproduce fundamental features of the energy cascade, at least in a statistical or average sense. Recent evidence questions the relevance of backscatter also in a local sense. Vela-Martín & Jiménez (Reference Vela-Martín and Jiménez2021) analysed the local SGS fluxes in numerical simulations of isotropic turbulence, and reported that backscatter is not correlated in scale (whereas direct energy transfer is), and that it disappears when the SGS fluxes are averaged over scales of the order of the filter width. These findings are in agreement with Ballouz & Ouellette (Reference Ballouz and Ouellette2018), who studied the efficiency of the interscale energy transfer in the inertial range, and found it to be low. They concluded that a large fraction of the SGS stresses drive spatial transport rather than the net energy cascade. This evidence suggests that backscatter represents conservative spatial fluxes at the filter scale, rather than interscale interactions, and is not part of the effective energetic exchange between the subgrid and resolved scales.
This investigation elaborates further on this idea, providing direct evidence that backscatter is associated not with an inverse cascade of energy, but with conservative spatial fluxes, and that it need not be reproduced to construct a ‘perfect’ model of the energy cascade. Although the local SGS fluxes in the inertial range often have a direction opposite to that of the mean energy cascade, we will see that this is a consequence of the commonly adopted representation of the SGS stress tensor; backscatter can be eliminated almost completely by transforming the SGS stress tensor, but without modifying the effect that the subgrid scales have on the resolved scales. This is possible because the filtered Navier–Stokes equations (NSE) depend not directly on the SGS stress tensor, but on its spatial divergence, and one is free to choose any representation of this tensor provided that its spatial divergence is preserved (Jiménez Reference Jiménez2016; Osawa & Jiménez Reference Osawa and Jiménez2018). This approach is based on the well-known gauge invariance of continuum mechanics (Jackson & Okun Reference Jackson and Okun2001), which is exploited here for modelling purposes. In this investigation, simple numerical optimisation procedures are applied to the SGS stress tensor calculated by filtering DNS of isotropic turbulence, leading to alternative representations of this tensor that satisfy the filtered NSE, but that produce almost no backscatter.
From a physical perspective, these findings are justified by resorting to non-equilibrium statistical mechanics. It is argued that non-equilibrium leads to the statistical irreversibility of inertial-range dynamics, which imposes a strongly unidirectional flow of energy from large to small scales. Although possible, global inverse cascades go against the tendency of non-equilibrium systems to evolve towards absolute equilibrium, and are therefore overwhelmingly less probable than global direct cascades (Vela-Martín & Jiménez Reference Vela-Martín and Jiménez2021). That backscatter can be transformed into a conservative spatial flux suggests that this global limitation also holds locally, restricting the spontaneous generation of local inverse cascades. The results presented in this paper show that the energy cascade is unidirectional enough that it can be modelled locally as an irreversible sink of energy, justifying the extended use of purely dissipative SGS models for LES of turbulence.
This paper is organised as follows. In § 2, we present a review of the LES modelling approach in isotropic turbulence, and its relation to the energy cascade. In § 3, we expose the concept of gauge invariance in the NSE and how it applies to the SGS stress tensor, and in § 4, we propose a methodology to exploit this invariance in order to eliminate backscatter almost completely. This methodology is applied to the SGS fluxes calculated by filtering isotropic turbulent flows. The database and the filtering approach are described § 5, and the results are presented in § 6. Finally, the implications of these results are discussed in § 7, and the conclusions are offered in § 8.
2. Modelling the energy cascade
2.1. The energy cascade in physical space
Let us consider the three-dimensional (3-D) incompressible NSE with periodic boundary conditions,
\begin{equation} \left.\begin{gathered} \partial_t u_i + u_j\,\partial_j u_i ={-}\partial_i p + \nu\,\partial_{kk} u_i + q_i, \\ \partial_i u_i =0, \end{gathered}\right\}\end{equation}
where 
$u_i$ is the 
$i$th component of the velocity vector, 
$p$ is the kinematic pressure, 
$\nu$ is the kinematic viscosity, and 
$q_i$ is a body force per unit mass that acts on the large scales of the flow to sustain turbulence. Repeated indices imply summation.
We consider the classical view of turbulence as a multiscale phenomenon in which energy is injected at an integral scale 
$L$, and then transferred to smaller scales by the turbulence cascade (Richardson Reference Richardson1922; Kolmogorov Reference Kolmogorov1941). At the end of the cascade, the kinetic energy is dissipated at the Kolmogorov scale, 
$\eta =(\nu ^{3}/\varepsilon )^{1/4}$, where 
$\varepsilon$ is the mean kinetic energy dissipation per unit mass.
To study the energy cascade in physical space, we apply an isotropic filter at an inertial scale 
$\ell$, such that 
$L\gg \ell \gg \eta$. The filter operation is denoted by 
$\overline {(\cdot )}$, and serves to separate the dynamics of scales above 
$\ell$ from the dynamics of scales below. The equation for each component of the filtered velocity vector reads
\begin{equation} \left.\begin{gathered} \partial_t \bar{u}_i + \bar{u}_j\,\partial_j \bar{u}_i ={-}\partial_i \bar{p} + \partial_j \tau^{S}_{ij} + \bar{q}_i, \\ \partial_i \bar{u}_i =0, \end{gathered}\right\} \end{equation}where
\begin{equation} \tau^{S}_{ij}=\bar{u}_i \bar{u}_j-\overline{u_i u_j} \end{equation}is the so-called SGS stress tensor, and we have neglected the viscous diffusion of momentum. Note that the SGS stress tensor defined in (2.3) has a sign opposite – but is formally similar – to the usual definition in the literature, and is adopted here to ease the manipulation of the equations.
By taking the dot product of (2.2) with the filtered velocity vector, we obtain the equation of the kinetic energy of the resolved scales, 
$\mathcal {E}=\tfrac {1}{2}\bar {u}_i\bar {u}_i$:
\begin{equation} \partial_t \mathcal{E} + \bar{u}_j\,\partial_j \mathcal{E} ={-}\partial_i (\bar{u}_i\bar{p}) + \bar{u}_i\,\partial_j\tau^{S}_{ij} + \bar{u}_i\bar{q}_i, \end{equation}
where the effect of the scales below 
$\ell$ on the energy of the scales above 
$\ell$ is represented by 
$\bar {u}_i\,\partial _j\tau ^{S}_{ij}$. This term is usually decomposed into
\begin{equation} \bar{u}_i\,\partial_j\tau^{S}_{ij}=D^{S} - T^{S}. \end{equation}The first term on the right-hand side,
\begin{equation} D^{S}=\partial_j (\bar{u}_i \tau^{S}_{ij}), \end{equation}
represents a spatial flux with zero average, 
$\langle D^{S} \rangle =0$, where the angular brackets denote the average over the flow domain. The second term,
\begin{equation} T^{S}=\bar{S}_{ij}\tau^{S}_{ij}, \end{equation}
is a source term usually identified with the local energy fluxes of the turbulence energy cascade (Leonard Reference Leonard1975; Borue & Orszag Reference Borue and Orszag1998; Eyink & Aluie Reference Eyink and Aluie2009; Cardesa et al. Reference Cardesa, Vela-Martín, Dong and Jiménez2015), where 
$S_{ij}$ is the rate-of-strain tensor. When the filter scale is within the inertial range, this term has a positive average, 
$\langle T^{S}\rangle >0$, indicating that energy flows, on average, from large to small scales. Moreover, its space and ensemble average, denoted by a double angular bracket, equates to the kinetic energy dissipation rate and to the power input, 
$\langle \langle T^{S} \rangle \rangle =\langle \langle \bar {u}_i \bar {q}_i \rangle \rangle =\varepsilon$.
The energy cascade in physical space is associated with 
$T^{S}$, rather than with 
$\bar {u}_i\,\partial _j \tau _{ij}^{S}$, because the former quantity is Galilean-invariant whereas the latter is not, i.e. 
$T^{S}$ does not change in a uniformly moving frame of reference (Eyink & Aluie Reference Eyink and Aluie2009). In addition, 
$T^{S}$ represents scale-local energy fluxes (Eyink Reference Eyink1995) consistent with the scale-locality of the energy cascade (Zhou Reference Zhou1993; Domaradzki & Carati Reference Domaradzki and Carati2007), and is free from non-local convective sweeping caused by large scales (Eyink Reference Eyink2005; Aluie & Eyink Reference Aluie and Eyink2009).
The sign of 
$T^{S}$ is commonly associated with the local direction of the energy fluxes: 
$T^{S}>0$ for a direct cascade, and 
$T^{S}<0$ for an inverse cascade. In the inertial range of isotropic turbulence, 
$T^{S}$ fluctuates strongly around its mean, and becomes negative in a large fraction of the domain (Domaradzki, Liu & Brachet Reference Domaradzki, Liu and Brachet1993; Cerutti & Meneveau Reference Cerutti and Meneveau1998; Tao et al. Reference Tao, Katz and Meneveau2002; Cardesa et al. Reference Cardesa, Vela-Martín, Dong and Jiménez2015). These negative values of 
$T^{S}$ constitute the so-called energy backscatter. The amount of backscatter depends strongly on the filter type, and is particularly profuse for sharp-Fourier filters, which produce backscatter in 
${\sim }40\,\%$ of the flow domain (Cardesa et al. Reference Cardesa, Vela-Martín, Dong and Jiménez2015; Buzzicotti et al. Reference Buzzicotti, Linkmann, Aluie, Biferale, Brasseur and Meneveau2018b). For well-behaved filters in physical space, such as the Gaussian filter, backscatter appears in 
${\sim }15\,\%$ of the flow domain (Cardesa et al. Reference Cardesa, Vela-Martín, Dong and Jiménez2015, see also § 5). Also, 
$D^{S}$ and 
$\bar {u}_i\,\partial _j\tau _{ij}$ become negative in a large fraction of the domain (
${\sim }50\,\%$ for Gaussian or sharp filters, see § 6.3), but these negative events are associated with conservative spatial fluxes related to the advection of the SGS stress tensor. Vela-Martín & Jiménez (Reference Vela-Martín and Jiménez2021) showed that when 
$\bar {u}_i\,\partial _j\tau _{ij}^{S}$ is averaged over volumes of the order of the filter width, this term becomes positive almost everywhere and highly correlated to the local average of 
$T^{S}$. This suggests that an important fraction of the negative events in 
$\bar {u}_i\,\partial _j\tau _{ij}^{S}$ comes from 
$D^{S}$, which seems to represent spatial fluxes that cancel over scales of the order of the filter width.
2.2. Modelling the subgrid-scale stress tensor
In the LES technique, the NSE are coarse-grained, usually by spatial filtering, and simulated on a coarse grid, leading to an equation similar to (2.2). This alleviates the computational cost of simulating a fully resolved turbulent flow, but requires that we model the SGS stress tensor, which contains the interactions between the subgrid and resolved scales. The modelled SGS stress tensor should be able to reproduce the energetic exchange between the resolved and subgrid scales.
A very popular class of LES models relies on the idea of the eddy viscosity; the traceless part of the modelled SGS stress tensor is assumed to be proportional to the filtered rate-of-strain tensor,
\begin{equation} \tau^{M}_{ij} - \tfrac{1}{3}\delta_{ij}\tau^{M}_{kk}=2 \nu_s \bar{S}_{ij}, \end{equation}
where 
$\nu _s$ is the eddy viscosity, which depends on the resolved velocity field, and is designed, among other things, to generate the appropriate amount of SGS energy transfer, 
$T^{M}=2\nu _s \bar {S}_{ij}\bar {S}_{ij}$. A very popular eddy-viscosity model is due to Smagorinsky (Reference Smagorinsky1963).
In order to account for backscatter, the eddy viscosity must take negative values, which generally leads to critical numerical instabilities. This is the case, for instance, in the dynamic procedure (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991) in which a local 
$\nu _s$ is computed dynamically following self-similarity arguments. In order to remove numerical instabilities, 
$\nu _s$ is set to zero, where it attains negative values (clipped) or averaged in homogeneous directions, leading to models for which 
$\nu _s>0$ and 
$T^{M}>0$ everywhere (Meneveau & Katz Reference Meneveau and Katz2000). These purely dissipative eddy-viscosity models have been used with considerable success despite their inability to reproduce backscatter, and despite the weak correlation that exists between the filtered rate-of-strain tensor and the SGS stresses calculated in DNS and experiments (Clark, Ferziger & Reynolds Reference Clark, Ferziger and Reynolds1979; Meneveau Reference Meneveau1994; Borue & Orszag Reference Borue and Orszag1998; Tao, Katz & Meneveau Reference Tao, Katz and Meneveau2000).
At variance with the eddy-viscosity models, the so-called similarity models (Bardina et al. Reference Bardina, Ferziger and Reynolds1980) do reproduce to a large extent the structure of the SGS stress tensor, and are able to produce numerically stable backscatter. They are based on the simple argument that the SGS stresses at the filter scale are proportional to the SGS stresses at larger scales. However, these models lead to inaccurate results, partly because they are not dissipative enough. This is why they are usually combined with a purely dissipative eddy-viscosity model (Zang, Street & Koseff Reference Zang, Street and Koseff1993; Meneveau & Katz Reference Meneveau and Katz2000). A similar problem exists for the tensor-diffusivity model, which also produces numerically stable backscatter but is not dissipative enough (Winckelmans et al. Reference Winckelmans, Wray, Vasilyev and Jeanmart2001).
Contrary to what is commonly believed, we will see that backscatter need not be reproduced to faithfully model the SGS interactions, and that the inability of SGS models to reproduce backscatter is not a limitation. Instead, it will be argued in this paper that this is a feature consistent with the strong statistical irreversibility of the energy cascade, which enforces an almost unidirectional flux of energy from large to small scales, even in a local sense.
3. Alternative subgrid-scale stresses
The filtered NSE have the general form of a conservation equation: the temporal variation of the local momentum equates to the divergence of spatial fluxes. Since the early development of continuum mechanics (Jackson & Okun Reference Jackson and Okun2001), it was noted that fluxes are not uniquely defined; there are infinitely many possible representations of these fluxes that produce the same physical effect on the system. As remarked by Jiménez (Reference Jiménez2016), this indeterminacy applies to the momentum fluxes in the NSE, and thus to the SGS stress tensor.
Let us consider (2.2). To 
$\tau ^{S}_{ij}=\bar {u}_i \bar {u}_j-\overline {u_i u_j}$, we can add any tensor field with zero divergence, 
$\partial _j \tau ^{*}_{ij}=0$, so that
\begin{equation} \tau^{A}_{ij}=\tau^{S}_{ij} + \tau^{*}_{ij}, \end{equation}
and 
$\tau ^{S}_{ij}$ produces the same physical effect on the filtered velocity field. This is so because although the SGS stress tensor changes, its divergence is not modified by the transformation in (3.1), i.e.
\begin{equation} \partial_j \tau^{S}_{ij}=\partial_j \tau^{A}_{ij}. \end{equation}
Therefore 
$\tau ^{S}_{ij}$ and 
$\tau ^{A}_{ij}$ are equivalent from the physical point of view. The NSE equations are said to be invariant under the gauge transformation 
$\tau ^{S}_{ij} \rightarrow \tau ^{A}_{ij}$. As is the case of electromagnetism with the Coulomb or Lorenz gauge (Jackson Reference Jackson2002), we are free to choose any gauge for the SGS stresses other than 
$\tau ^{*}_{ij}=0$, at our own convenience. An example applied to the momentum transfer in turbulent channels is Osawa & Jiménez (Reference Osawa and Jiménez2018). In the context of LES modelling, gauge invariance implies that it is necessary to model not the SGS stress tensor, but rather its divergence. In other words, we are free to model any particular representation of the SGS stresses, provided that its divergence produces the required effect on the resolved scales.
From an energetic perspective, any transformation of the SGS stress tensor leads to the same average interscale energy transfer, but, as we will show, produces in general different spatial distributions of the local SGS energy fluxes. From (3.2) follows
\begin{equation} \bar{u}_i\,\partial_j \tau^{S}_{ij}=\bar{u}_i\,\partial_j \tau^{A}_{ij} \end{equation}and
\begin{equation} D^{S} - T^{S} = D^{A} - T^{A}, \end{equation}
where 
$T^{A}=\bar {S}_{ij}\tau ^{A}_{ij}$ and 
$D^{A}=\partial _j(\bar {u}_i\tau ^{A}_{ij})$ are alternative representations of the local SGS fluxes and the local conservative term due to the SGS interactions. Although (3.4) holds for any 
$\tau _{ij}^{A}$, we will see that, in general, 
$T^{A}(\boldsymbol {x})\neq T^{S}(\boldsymbol {x})$ and 
$D^{A}(\boldsymbol {x})\neq D^{S}(\boldsymbol {x})$. Thus by changing the gauge of 
$\tau ^{S}_{ij}$, we are locally modifying both 
$D^{S}$ and 
$T^{S}$, but not 
$D^{S} - T^{S}$. Modifying the SGS stress tensor is equivalent to ‘moving’ the effect of the subgrid scales from a source term into a conservative spatial flux, which produces no net energy transfer to the unresolved scales. Since 
$\langle D^{S}\rangle =\langle D^{A}\rangle =0$, we have
\begin{equation} \langle u_j\,\partial_i \tau^{S}_{ij}\rangle=\langle u_j\,\partial_i \tau^{A}_{ij}\rangle=\langle T^{S}\rangle=\langle T^{A}\rangle, \end{equation}and the average energy flux is not modified by a change of the SGS stress tensor.
In the following, it will be shown that it is possible to find transformations of the SGS stress tensor, obtained by filtering DNS data of isotropic turbulence, that eliminate the backscatter in the SGS fluxes almost completely.
4. Reducing backscatter in the subgrid-scale fluxes
Any physical representation of the SGS stress tensor must fulfil
\begin{equation} \partial_{j}\tau_{ij}=Q_i, \end{equation}
where 
$Q_i=\partial _j \tau _{ij}^{S}=\partial _{j}(\bar {u}_i \bar {u}_j - \overline {u_i u_j})$. This equation, together with the symmetry constraint
\begin{equation} \tau_{ij} - \tau_{ji}=0, \end{equation}defines the linear subspace that contains all the possible representations of the SGS stress tensor. We explore this linear subspace in order to find representations of the SGS fluxes without backscatter.
We search for representations of the SGS stresses in which the fluctuations of 
$T$ around its mean are small enough to eliminate its negative values. This idea is formulated as the minimisation of the square of the standard deviation of 
$T$, 
$\sigma ^{2}=\langle T^{2} \rangle - \langle T \rangle ^{2}$. Since 
$\langle T \rangle$ is invariant under transformations of the SGS stresses, minimising the standard deviation of 
$T$ is equivalent to minimising the average of 
$T^{2}$.
Let us consider the minimisation of the functional
\begin{equation} \mathcal{L} = \int_V [ T^{2} + \gamma_i (\partial_{j}\tau_{ij} - Q_i) + \mu_{ij}(\tau_{ij} - \tau_{ji})]\,\mathrm{d}V, \end{equation}
where 
$\gamma _i$ and 
$\mu _{ij}$ are Lagrange multipliers that impose (4.1) and (4.2). The integral extends over the full computational volume. It is possible to apply calculus of variations to obtain a closed linear system of equations on 
$\tau _{ij}$, 
$\gamma _i$ and 
$\mu _{ij}$, which, together with the constraints, can be solved to obtain the optimal form of 
$\tau _{ij}$. The exact expression of this system of equations is presented in Appendix A. To the knowledge of the author, there is no analytical solution to this problem, and one must resort to numerical procedures.
We apply a projected gradient descent (PGD) method. Let us define a state vector that contains the components of the SGS stress tensor at each point in the domain, 
$\boldsymbol \chi =\{\tau _{11}(\boldsymbol {x}), \tau _{22}(\boldsymbol {x}),\dots,\tau _{ij}(\boldsymbol {x})\}$, and iterate to minimise 
$T^{2}$ by following
\begin{equation} \left.\begin{gathered} \boldsymbol{\chi}^{*}=\boldsymbol{\chi}^{n} - \beta\,\Delta \boldsymbol{\chi},\\ \boldsymbol{\chi}^{n+1}=\mathcal{P}(\boldsymbol{\chi}^{*}), \end{gathered}\right\} \end{equation}
where 
$\beta$ is the step size, 
$\mathcal {P}$ is a linear projector, and 
$\Delta {\boldsymbol {\chi }}$ is the unprojected descent direction, which is given by the gradient of the cost function with respect to the state vector,
\begin{equation} \Delta \boldsymbol{\chi}=\partial{\boldsymbol{\chi}}\left(\int_V T^{2}\,\mathrm{d}V\right)=2T\boldsymbol{\varSigma}. \end{equation}
Here, 
$\boldsymbol {\varSigma }=\{\bar {S}_{11}(\boldsymbol {x}),\bar {S}_{22}(\boldsymbol {x}),\dots,\bar {S}_{ij}(\boldsymbol {x})\}$ contains the components of the filtered rate-of-strain tensor at each point in the domain, and 
$\partial _{\boldsymbol {\chi }}$ represents the partial derivative with respect to 
$\boldsymbol {\chi }$.
The linear operator 
$\mathcal {P}(\boldsymbol {\chi }^{*})$ projects 
$\boldsymbol {\chi }^{*}$ into the closest state vector (under the 
$L^{2}$-norm) that fulfils the linear constraints in (4.1) and (4.2). This projector is derived by minimising a functional using calculus of variations, which leads to a linear set of partial differential equations that can be solved readily with periodic boundary conditions. The details of this procedure are presented in Appendix B. The step size, 
$\beta$, is fixed numerically so that the descent is maximised (Nocedal & Wright Reference Nocedal and Wright2006). The standard representation of the SGS stress tensor, 
$\tau _{ij}^{S}=\bar {u}_i \bar {u}_j - \overline {u_i u_j}$, is chosen as an initial condition of the PGD algorithm.
The PGD algorithm is numerically affordable, and allows us to obtain at each iteration a true physical representation of the SGS stresses, i.e. a representation that fulfils the constraints of the problem to numerical accuracy. This last property is especially advantageous because we are interested not in finding the minimum standard deviation of 
$T$, but in reducing it enough to eliminate backscatter.
5. Local energy fluxes in the inertial range of isotropic turbulence
5.1. Database
The optimisation algorithm presented above is applied to the SGS stresses obtained by filtering turbulent flow fields. We consider an extensive database of isotropic turbulence generated by DNS in a triply periodic cubic box, and forced linearly at the large scales. The database is composed of hundreds of statistically independent fields at three different Reynolds numbers, 
$Re_\lambda =U\lambda /\nu =140$, 
$240$ and 
$324$, which correspond to numerical grids of 
$256^{3}$, 
$512^{3}$ and 
$1024^{3}$ points, respectively. Here, 
$U$ is the root-mean-squared velocity fluctuations, and 
$\lambda =15\sqrt {(\nu /\varepsilon )}U$ is the Taylor microscale. Three datasets have been analysed, one for each 
$Re_\lambda$ value. These datasets contain time-consecutive snapshots from a single history spanning several integral eddy-turnover times, 
$T=L/U$, where 
$L$ is the characteristic large scale of the flow, calculated from the integral of the longitudinal velocity autocorrelation (Pope Reference Pope2001). From each dataset, independent fields separated by at least 
$0.5T$ are selected. In total, 
$290$ independent flow fields are analysed. Table 1 shows a summary of the details of the database. Further details of the simulations and the code can be found in Cardesa et al. (Reference Cardesa, Vela-Martín, Dong and Jiménez2015) and Cardesa, Vela-Martín & Jiménez (Reference Cardesa, Vela-Martín and Jiménez2017).
Table 1. Details of the homogeneous isotropic turbulence database. Here, 
$N$ is the number of grid points of the numerical mesh in each direction, 
$k_{max}$ is the maximum resolved Fourier wavenumber, 
$T_{{eto}}$ is the number of eddy-turnovers of each dataset, and 
$N_{{snap}}$ is the number of snapshots used for the optimisation procedure. The snapshots are separated in time approximately 
$T_{{eto}}/N_{{snap}}$. In the last two columns, we show the maximum and minimum filter scale applied to each simulation normalised with the integral and Kolmogorov scales.
5.2. Filtering operations
For each velocity field, the standard SGS stresses, 
$\tau ^{S}_{ij}$, are calculated by applying a filtering operation,
\begin{equation} \bar{\phi}(\boldsymbol{x})=\int \phi(\boldsymbol \xi)\,G(\boldsymbol{x} - \boldsymbol \xi) \,\mathrm{d}^{3}\boldsymbol{\xi}, \end{equation}where the integral is taken over the full computational domain. The filtering operations are performed as a product in Fourier space,
\begin{equation} \hat{\bar{\phi}}(\boldsymbol{k})=\hat{G}(\boldsymbol{k})\,\hat{\phi}(\boldsymbol{k}), \end{equation}
where the caret denotes the Fourier transform, and 
$\boldsymbol {k}=\{k_i\}$ is the wavenumber vector in Fourier space.
Two different LES filters are applied at each scale 
$\ell$. First, a sharp-Fourier filter (SF) at wavenumber 
${\rm \pi} /\ell$, defined as
\begin{equation} \left.\begin{array}{ll} \hat{G}(\boldsymbol{k})=1, & \text{if } |\boldsymbol{k}|<\dfrac{\rm \pi}{\ell}, \\[6pt] \hat{G}(\boldsymbol{k})=0, & \text{if } |\boldsymbol{k}|\geqslant \dfrac{\rm \pi}{\ell}, \end{array}\right\}\end{equation}and second, a Gaussian filter (GF)
\begin{equation} \hat{G}(\boldsymbol{k})=\exp\left(-\frac{\boldsymbol{k}^{2} \ell^{2}}{2{\rm \pi}^{2}}\right). \end{equation} In addition to the LES filters, a numerical filter is also applied to reproduce the discretisation of the filtered NSE on a numerical grid (Carati, Winckelmans & Jeanmart Reference Carati, Winckelmans and Jeanmart2001). This filtering operation, which is denoted by 
$\tilde {(\cdot )}$, projects a field on a truncated Fourier basis with maximum wavenumber in each direction equal to
\begin{equation} k^{max}_i=\frac{\alpha {\rm \pi}}{\ell}, \end{equation}
where 
$\alpha$ is a parameter. This numerical filter is used for consistency with LES simulations, and to limit the number of degrees of freedom on which the optimisation procedure operates, thus alleviating the computational costs.
The main results of this paper are obtained with 
$\alpha =2$, so that the grid spacing of the numerical mesh is half the filter length, but results for 
$\alpha =3$ are also presented to show that they are not sensitive to the resolution of the numerical grid. For the case of the SF, 
$\alpha =1$ is also used. This is similar to the common case in which the LES filter is imposed by the grid. The standard SGS stresses and the SGS fluxes are effectively calculated as
\begin{equation} \left.\begin{gathered} \tau^{S}_{ij}=\tilde{{\bar{u}}_i{\bar{u}}_j}-\tilde{\overline{u_i u_j}},\\ T=\tau_{ij}\tilde{\overline{S}}_{ij}, \end{gathered}\right\} \end{equation}
where 
$\tau _{ij}$ is any physical representation of the SGS stresses. All operations are performed on the truncated numerical grid with proper aliasing procedures, including the triple products present in the unprojected descent gradient in (4.5). Note that because the numerical filter is a projector, the SGS stresses and the SGS fluxes are Galilean-invariant (Buzzicotti et al. Reference Buzzicotti, Linkmann, Aluie, Biferale, Brasseur and Meneveau2018b). This property also extends to the gradient descent direction in (4.5), and therefore to the alternative SGS stresses produced by the optimisation procedure.
5.3. Statistics of the subgrid-scale energy fluxes
The turbulent flow fields in the database are filtered at scales separated by a factor of 
$2$, from 
$\ell =0.65L$ to the end of the inertial range, 
$\ell \approx 30\eta$, and 
$T^{S}$ is calculated by applying (5.6) with 
$\alpha =2$. The results are very similar for 
$\alpha =3$ and 
$\alpha =1$ (in the case of the SF). The maximum and minimum filter scales used for each Reynolds number are shown in table 1. In figure 1, we show the average energy fluxes at each scale for the three different Reynolds numbers and the two filter types. The average of 
$T^{S}$ is very close to the total energy dissipation for all cases, indicating that the filter scales lie within an inertial range.
Figure 1. Average SGS fluxes normalised by the average energy dissipation for (a) an SF, and (b) a GF, for different filter scales and Reynolds numbers.
In figures 2(a,b), we show the probability density function (p.d.f.) of 
$T^{S}$ at different scales for 
$Re_\lambda =240$ and for the SF and GF. The local SGS fluxes are negative in a large fraction of the domain, particularly for the SF, in agreement with previous investigations (Aoyama et al. Reference Aoyama, Ishihara, Kaneda, Yokokawa, Itakura and Uno2005; Cardesa et al. Reference Cardesa, Vela-Martín, Dong and Jiménez2015; Buzzicotti et al. Reference Buzzicotti, Linkmann, Aluie, Biferale, Brasseur and Meneveau2018b). To quantify the amount of backscatter, we define the volume fraction in which 
$T<0$,
\begin{equation} v^{-}=\frac{1}{V}\int_V \mathcal{H}(T)\,\mathrm{d}V, \end{equation}and the net amount of backscatter normalised with the total energy flux,
\begin{equation} f^{-}={-}\frac{1}{\langle T \rangle V}\int_V T\,\mathcal{H}(T)\,\mathrm{d}V, \end{equation}
where 
$\mathcal {H}(T)$ is a step function that takes values 
$\mathcal {H}=1$ when 
$T<0$, and 
$\mathcal {H}=0$ when 
$T\geqslant 0$, and 
$V$ is the volume of the computational domain.
Figure 2. (a,b) P.d.f.s of 
$T^{S}$ for different filter scales for 
$Re_\lambda =240$, and for (a) an SF, and (b) a GF. The vertical dotted line marks the origin, and the vertical dashed line the average of 
$T^{S}$. (c,d) Volume fractions of backscatter, 
$v^{-}$, and net backscatter, 
$f^{-}$, in 
$T^{S}$ for different filter scales and Reynolds numbers, and for (c) an SF, and (d) a GF. The markers correspond to the average in the database and the error bars to the standard deviation.
These quantities are calculated for each independent flow field in the database, and their average for different filter scales and Reynolds numbers is shown in figures 2(c,d). To characterise how much these quantities change across different flow fields, their standard deviations in the database are shown as error bars. For the SGS fluxes calculated with the SF, backscatter is present in around half of the points in the domain, 
$v^{-}\approx 0.5$, but this fraction reduces to less than one-fifth when a GF is applied, 
$v^{-}\approx 0.15$. The differences between filters is more drastic for the net backscatter, which is close to unity for the SF but decreases to 
$f^{-}\approx 0.05$ for the GF. Both 
$f^{-}$ and 
$v^{-}$ increase with decreasing scale, especially 
$f^{-}$ for filter scales close to the dissipative range, i.e. for the smallest filter at each 
$Re_\lambda$ value. This increment may reflect the reduction of the mean energy flux due to viscous dissipation, and the increased fluctuations of the local SGS fluxes due to intermittency effects (Cerutti & Meneveau Reference Cerutti and Meneveau1998). The standard deviations of these quantities in the database are small compared to their means, indicating that they do not change substantially across different flow fields.
6. Results
6.1. Convergence of the projected gradient descent algorithm
For each independent flow field and filter scale, we apply the PGD algorithm described in § 4 using 
$\tau ^{S}_{ij}$ as an initial condition. At each iteration, the convergence of the algorithm is quantified using the standard deviation of 
$T$ in the domain,
\begin{equation} \sigma(T)=(\langle T^{2} \rangle - \langle T \rangle^{2})^{1/2}, \end{equation}
and the 
$L^{2}$-norm of the projected gradient of the cost function,
\begin{equation} |\Delta \boldsymbol{\chi}_p|=\frac{1}{\beta}|\boldsymbol{\chi}^{n} - \boldsymbol{\chi}^{n-1}|, \end{equation}
which is derived from (4.4). To assess the ability of the algorithm to reduce backscatter, also 
$f^{-}$ and 
$v^{-}$ are calculated at each iteration. In this section, only results for 
$\ell =0.33L$ and 
$Re_\lambda =240$ are shown, but the convergence is qualitatively similar for the rest of the filter scales and Reynolds numbers.
In figures 3(a–d), we show the evolution of 
$\sigma (T)$ and 
$|\Delta \boldsymbol {\chi }_p|$ averaged across the database as functions of the number of steps in the optimisation procedure. In the case of the SF, the PSD algorithm is able to decrease 
$\sigma (T)$ by a factor of 
$10$ before reaching a plateau. For the GF, the initial standard deviation of 
$T$ is approximately half that of the SF, but the algorithm reaches a plateau with a similar value. These results are quantitatively very similar for different values of 
$\alpha$, indicating that the size of the numerical grid on which the optimisation procedure is performed has no effect on the final outcome of the PGD algorithm. For both filter types, the gradient of the cost function decreases approximately by a factor of 
$10$ and reaches a plateau. It is unclear whether the algorithm converges to a local or global minimum.
Figure 3. Evolution of: (a,b) 
$\sigma (T)/\langle T\rangle$; (c,d) 
$|\Delta \boldsymbol {\chi }_p|$ normalised by its value at the first step 
$|\Delta \boldsymbol {\chi }_p|_{n=0}$; (e,f) 
$v^{-}$ and (g,h) 
$f^{-}$, averaged over the database as functions of the iteration step, 
$n$. The figures in (a,c,e,g) correspond to the SF, and the ones in (b,d,f,h) to the GF. The different line styles correspond to different values of the resolution parameter, 
$\alpha$, used for the numerical filter in (5.5). The shaded blue contour marks plus/minus the value of the standard deviation of each quantity in the database for 
$\alpha =2$. The data correspond to 
$\ell =0.33L$ and 
$Re_\lambda =240$.
In any case, the optimisation procedure is effective, and yields alternative SGS fluxes with reduced backscatter. In figures 3(e–h), we show the evolution of 
$f^{-}$ and 
$v^{-}$ averaged across the database as functions of the number of steps in the optimisation procedure. Initially, the standard SGS fluxes calculated with the SF yield 
$f^{-}\approx 0.6$ and 
$v^{-}\approx 0.4$, but after the optimisation, these quantities are reduced approximately by three orders of magnitude to 
$v^{-}\sim 10^{-3}$ and 
$f^{-}\sim 10^{-4}$. The optimisation procedure is also effective for the GF, and reduces 
$f^{-}$ and 
$v^{-}$ to comparable values. Again, we observe only small variations of 
$f^{-}$ and 
$v^{-}$ with the resolution of the numerical grid, indicating that the outcome of the algorithm is not sensitive to the value of 
$\alpha$. Hereafter, we consider 
$\alpha =2$.
In figures 3(a–h), we have also plotted the standard deviation of each quantity with respect to the mean across different fields. In all cases, this deviation is at most of the order of the mean, indicating that the algorithm is effective in reducing backscatter in general. The SGS fluxes obtained with this optimisation procedure are denoted by 
$T^{O}$, and will be referred to as the optimised SGS fluxes (rather than optimal SGS fluxes) to remark that it is not clear whether they are optimal in the sense of (4.3).
6.2. Statistics of the subgrid-scale energy fluxes with reduced backscatter
This section is devoted to characterising the reduction of backscatter achieved by the optimisation algorithm for different Reynolds numbers and filter scales. We will see that the results reported in the previous section for 
$\ell =0.33L$ and 
$Re_\lambda =240$ are quantitatively very similar for all the filter scales and Reynolds numbers tested.
In figures 4(a,b), we show the standard deviation of the optimised SGS fluxes, 
$\sigma (T^{O})/\langle T^{O}\rangle$, for different filter scales and Reynolds numbers. This quantity increases slowly with decreasing filter scale, and is very similar for both filter types. The best least-squares fit to the data, assuming it depends logarithmically on the filter scale, is 
$\sigma = \langle T \rangle (0.25-0.037\log (\ell /L) )$ for the SF, and 
$\sigma = \langle T \rangle (0.25-0.033\log (\ell /L))$ for the GF. This fit shows that the standard deviation of the optimised stresses depends only weakly on the filter scale.
Figure 4. (a,b) The standard deviation, (c,d) the volume fraction, and (e,f) the net amount of backscatter, as functions of the filter scale, 
$\ell /L$, for different Reynolds numbers. The markers denote the average of each quantity in the database, and the bars indicate the standard deviation. In (a,b), the bars are smaller than the markers and are not plotted.
In figures 4(c–f), we show the backscatter volume fraction and the net backscatter as functions of the filter scale for different Reynolds numbers. The backscatter volume fraction is reduced in all cases to 
$v^{-}\approx 10^{-3}$, and the net backscatter to 
$f^{-}\approx 10^{-4}$. These values change only slightly with decreasing filter scale, and are very similar for different Reynolds numbers and filter types.
The deviations of all the quantities across the database are shown as error bars in figures 4(a–f), and are small compared to the means (or at most of the same order as the means), indicating that backscatter can be reduced in general, independently of the filter scale and the Reynolds number.
6.3. The structure of the optimised subgrid-scale fluxes
In view of the results presented above, we now characterise how the optimisation algorithm changes the structure of the SGS fluxes. In figures 5(a,b), the p.d.f.s of the standard and optimised SGS fluxes are compared for 
$Re_\lambda =140$ and 
$\ell =0.33L$. Here, the SGS fluxes are normalised by their ensemble average, denoted by a double bracket. There is a strong reduction of the standard deviation of the SGS fluxes for the two filter types, which leads to a concentration of the p.d.f. around the mean. Since the mean remains invariant under transformations of the SGS stresses, the narrowing of the p.d.f. entails that only a negligible fraction of 
$T^{O}$ reaches negative values. Conversely to 
$T$, 
$D$ is not strongly modified by the optimisation procedure. As shown in figures 5(c,d), the p.d.f.s of 
$D^{O}$ and 
$D^{S}$ normalised by the mean SGS fluxes are very similar. Finally, the p.d.f.s of 
$T^{S}$ and 
$T^{O}$ without their ensemble average (denoted by a prime), and divided by their ensemble standard deviation, 
$\langle \langle T'^{2}\rangle \rangle ^{1/2}$, are shown in figures 5(e,f). For the SF, there is a reasonable collapse of both p.d.f.s for values close to the mean, but differences are noticeable in the tails, which are fatter for 
$T^{S}$ than for 
$T^{O}$. For the GF, the p.d.f.s are substantially different, with 
$T^{S}$ being more positively skewed and flatter that 
$T^{O}$. These results suggest that the optimisation procedure modifies not only the variance of 
$T$, but also its higher-order moments.
Figure 5. Probability density functions of (a,b,e,f) 
$T^{S}$ and 
$T^{O}$, and (c,d) 
$D^{S}$ and 
$D^{O}$. The figures in (a,c,e) correspond to quantities calculated with an SF, and the ones in (b,d,f) with a GF. In (a–d), quantities are normalised by the ensemble average of the SGS fluxes, 
$\langle \langle T \rangle \rangle$, and in (e,f) quantities are normalised by the standard deviation 
$\langle \langle T'^{2} \rangle \rangle ^{1/2}$. The prime denotes quantities without their ensemble average. In (b,f), the circles mark the p.d.f.s of the optimised SGS fluxes calculated with an SF. The data correspond to 
$\ell =0.33L$ and 
$Re_\lambda =140$.
We focus now on the third- and fourth-order moments, the skewness and flatness factors, respectively, of 
$T^{S}$ and 
$T^{O}$, which are defined as
\begin{equation} s(T)=\frac{\langle\langle T'^{3}\rangle\rangle}{\langle\langle T'^{2}\rangle\rangle^{3/2}} \end{equation}and
\begin{equation} \kappa(T)=\frac{\langle\langle T'^{4}\rangle\rangle}{\langle\langle T'^{2}\rangle\rangle^{2}}. \end{equation} In figures 6(a,b), 
$s(T^{O})$ and 
$s(T^{S})$ are shown as functions of the filter scale for the SF and the GF. The skewness of 
$T^{S}$ increases with decreasing scales for the SF and the GF, but is larger for the GF than for the SF, in agreement with the p.d.f.s shown in figures 2(a,b). For the SF, 
$T^{S}$ is rather symmetric, while for the GF the p.d.f. of 
$T^{S}$ displays stretched right tails. After the optimisation, the skewness is not strongly modified for the SF, and 
$s(T^{O})$ is similar to 
$s(T^{S})$. On the other hand, for the GF, the skewness of 
$T^{O}$ is substantially reduced with respect to 
$T^{S}$ and becomes similar to the data of the SF. This reflects the disappearance of the strong positive events in the optimised SGS fluxes.
Figure 6. (a,b) Skewness 
$s$, and (c,d) flatness factor 
$\kappa$, of the optimised and standard SGS fluxes, and of the total SGS interactions, 
$D-T$, as functions of the filter scale, for (a,c) the SF, and (b,d) the GF. For 
$\ell /L=0.64$ and 
$0.32$ the data come from 
$Re_\lambda =140$, for 
$\ell /L=0.16$ from 
$Re_\lambda =240$, and for 
$\ell /L=0.08$ from 
$Re_\lambda =320$. The black diamonds in (d) correspond to the data by Cerutti & Meneveau (Reference Cerutti and Meneveau1998) (CM98) at 
$Re_\lambda \approx 150$, which have been plotted only until 
$\ell /\eta =20$, where 
$\eta$ is the Kolmogorov scale.
In figures 6(c,d), 
$\kappa (T^{O})$ and 
$\kappa (T^{S})$ are shown as functions of the filter scale for the SF and the GF. The flatness factor calculated here is in agreement with the flatness factor of 
$T^{S}$ at different scales calculated by Cerutti & Meneveau (Reference Cerutti and Meneveau1998) for 
$Re_\lambda \approx 150$, which is shown in figure 6(d). For 
$T^{S}$, the flatness grows as a power law with decreasing scale due to intermittency effects. However, for 
$T^{O}$, the flatness remains approximately constant in the range of filter scales used here for both filter types, indicating that the optimised SGS fluxes are free from intermittency effects in their fourth-order moment. Note that this implies not that intermittency can be eliminated from the SGS interactions, only that it need not be modelled as part of the SGS energy fluxes. In figures 6(c,d), we also show the flatness factor of 
$D-T$, which increases with decreasing scale. This quantity represents the total effect of the subgrid scales on the resolved kinetic energy, and is invariant under gauge transformations. This suggests that intermittency should be modelled as a robust feature of the SGS interactions, regardless of the particular definition of the SGS stress tensor. A more elaborate discussion on this idea is offered in the conclusion, § 8.
Despite the strong differences between 
$T^{S}$ when calculated with the SF and the GF, the p.d.f.s of the optimised SGS fluxes are very similar for the two filter types, as shown by their good collapse in figures 5(b,f). We define the correlation coefficient
\begin{equation} \rho(\phi,\psi)=\frac{\langle\langle \phi' \cdot \psi'\rangle\rangle}{\sqrt{\langle\langle \phi'^{2} \rangle\rangle\langle\langle \psi'^{2} \rangle\rangle}}, \end{equation}
and use it to quantify the differences between the optimised SGS fluxes calculated with the SF, 
$T_{{s}}$, and with the GF, 
$T_{{g}}$. For the standard SGS fluxes, 
$\rho (T^{S}_{{s}},T^{O}_{{g}})=0.4$, whereas for the optimised SGS fluxes, the correlation coefficient increases to 
$\rho (T^{O}_{{s}},T^{O}_{{g}})\approx 0.65$. These values are similar for the different Reynolds numbers and filter scales, and suggest that the differences between 
$T^{S}_{{g}}$ and 
$T^{S}_{{s}}$ are partially related to conservative spatial fluxes.
As discussed in § 3, the changes in 
$T^{O}$ with respect to 
$T^{S}$ imply the modification of 
$D$. However, this term changes only slightly due to the optimisation procedure. This is shown by its p.d.f.s in figures 5(c,d), and by the correlation 
$\rho (D^{S},D^{O})\approx 0.9$. On the other hand, 
$\rho (T^{S},T^{O})\approx 0.01$. Again, these values are comparable for different Reynolds numbers, filter scales and filter types. The reduction of backscatter, which implies a strong modification of the structure of 
$T$, corresponds to a slight modification of the conservative spatial fluxes produced by the SGS interactions.
The low correlation between 
$T^{S}$ and 
$T^{O}$ reflects the differences between 
$\tau ^{S}_{ij}$ and 
$\tau ^{O}_{ij}$. We calculate the correlation coefficient between these two tensors (including the diagonal and the trace), and find that it is not high. For the GF and the SF, 
$\rho (\tau ^{S}_{ij},\tau ^{O}_{ij})\approx 0.6$ for 
$\ell /L=0.64$, but it decreases with decreasing scale down to 
$\rho (\tau ^{S}_{ij},\tau ^{O}_{ij})\approx 0.4$ for 
$\ell /L=0.08$. Borue & Orszag (Reference Borue and Orszag1998) reported a similar decay of the correlation between 
$\tau ^{S}_{ij}$ and the Smagorinsky model with scale, although they found values between 
$0.4$ and 
$0.25$. However, they considered only the correlation of the off-diagonal elements of the tensor. Here, when only the traceless part of the tensor – or the off-diagonal elements of the tensor – is considered, the correlations are lower than 
$0.5$ in any case, and particularly low for the SF. Let us remark that 
$\tau ^{S}_{ij}$ and 
$\tau ^{O}_{ij}$ are equivalent physical representations of the SGS stress tensor, yet they are not strongly correlated locally.
6.4. Visualisation of the standard and optimised subgrid-scale fluxes
To gain a qualitative understanding of the changes produced by the optimisation algorithm, we compare two-dimensional (2-D) cuts of 
$T^{S}$ and 
$D^{S}$ with the same cuts of 
$T^{O}$ and 
$D^{O}$. These cuts are shown in figures 7(a–d). For ease of visualisation, only data for 
$Re_\lambda =320$, 
$\ell =0.08L$ and the GF are plotted. We observe that 
$T^{S}$ is particularly active in certain regions of the domain, which also correspond to regions of intense activity in 
$D^{S}$. Intense positive events of 
$T^{S}$ appear frequently surrounded by intense negative events, which resemble the characteristic patterns generated by advection. These clusters of positive and negative events are transformed into a spatial flux by the optimisation procedure. Compared to 
$T^{S}$, 
$T^{O}$ is much more spread across the domain, and has values closer to the mean. On the other hand, the structure of the conservative term remains largely unchanged by the optimisation procedure, and 
$D^{O}$ and 
$D^{S}$ have a very similar structure, in agreement with their high correlation coefficient.
Figure 7. (a–d) Colour plots of the SGS fluxes 
$T$, and of the conservative spatial fluxes 
$D$, for 
$Re_\lambda =320$, 
$\ell =0.08L$, and a GF. The panels correspond to: (a) 
$T^{S}$; (b) 
$T^{O}$; (c) 
$D^{S}$; (d) 
$D^{O}$. (e,f) Isocontours of 
$T=0$ (blue) and 
$T=\langle T \rangle$ (red), for (e) the standard SGS fluxes, and (f) the optimised SGS fluxes. All panels correspond to the same plane in the same flow field.
To conclude this section, we provide a visual demonstration of the substantial reduction of backscatter due to the optimisation procedure. In figures 7(e,f), we show isocontours of 
$T=0$, which enclose backscatter events, for the standard and optimised SGS stresses. For comparison, we also show isocontours of 
$T=\langle T\rangle$. While for 
$T^{S}$ backscatter appears distributed across most of the domain, for 
$T^{O}$ backscatter is reduced to isolated spots of size 
${\sim }\ell ^{3}$ and intensity of the order of 
$-0.1\langle T\rangle$. These backscatter events are too few and too weak to have any effect on the energy balance of the resolved scales.
7. Implications of the backscatterless subgrid-scale stresses
7.1. A purely dissipative ‘perfect’ eddy viscosity
An important consequence of the results presented in the previous sections is that it is not necessary to reproduce backscatter to model the interactions between the subgrid and resolved scales in LES. Here we illustrate this idea with a simple example.
Let us consider an eddy-viscosity model, in which the traceless part of the modelled SGS stress tensor is proportional to the resolved rate-of-strain tensor,
\begin{equation} \tau^{M}_{ij} - \tfrac{1}{3}\tau^{M}_{kk}\delta_{ij}=2\nu_s\bar{S}_{ij}, \end{equation}
and 
$\nu _s$ is the eddy viscosity. Assuming that a real representation of the SGS stress tensor is available, it is possible to choose an eddy viscosity that minimises the difference between the modelled SGS stresses and some representation of the true stresses. This constitutes a sort of ‘optimal’ a priori eddy viscosity. In the simple case in which differences are measured with an 
$L^{2}$-norm, this optimal eddy viscosity is obtained by minimising the functional
\begin{equation} \mathcal{M}=\int_V(2\nu_s\bar{S}_{ij} - \tau_{ij})^{2}\,\mathrm{d}V, \end{equation}
where 
$\tau _{ij}$ is any representation of the SGS stress tensor. Introducing a variation on 
$\nu _s$, 
$\delta \nu _s$, and imposing that the variation of 
$\mathcal {M}$ is zero, leads to
\begin{equation} (2\nu_s\bar{S}_{ij} - \tau_{ij})\,4\bar{S}_{ij}\,\delta \nu_s=0. \end{equation}
The eddy viscosity that satisfies this equation for any 
$\delta \nu _s$ is
\begin{equation} \nu_s^{opt} =\frac{\tau_{ij} \bar{S}_{ij}}{2\bar{S}_{ij}\bar{S}_{ij}}, \end{equation}
for 
$\bar {S}_{ij}\bar {S}_{ij}>0$. With this definition of the eddy viscosity, the SGS model in (7.1) produces exactly the same local SGS fluxes as the real turbulent flow,
\begin{equation} T^{M}=2\nu_s^{opt} \bar{S}_{ij}\bar{S}_{ij}=\tau_{ij} \bar{S}_{ij}. \end{equation}
Note that 
$\nu _s^{opt}$ is proportional to 
$\tau _{ij}$, which is defined up to a gauge transformation. This implies that there are infinitely many possible ‘optimal’ eddy viscosities, one for each different representation of the SGS stresses. For the standard definition 
$\tau ^{S}_{ij}=\bar {u}_i \bar {u}_j-\overline {u_i u_j}$, 
$\nu _s^{opt}$ becomes negative in a large fraction of the domain (see § 5.3), and modelling it leads, most surely, to numerical instabilities (Meneveau & Katz Reference Meneveau and Katz2000). However, as shown in this paper, there are representations of the SGS stress tensor for which 
$\tau _{ij}\bar {S}_{ij}$, and therefore 
$\nu _s^{opt}$, is positive almost everywhere, independently of the filter scale and filter type. This example illustrates the main point of this work: there is at least a ‘perfect’ (from the perspective of the SGS fluxes) eddy-viscosity model that reproduces exactly a representation of the SGS fluxes with almost no backscatter.
The idea of the ‘optimal’ eddy viscosity described here shows that even in the case in which the modelled SGS stress tensor is constrained to be parallel to the filtered rate-of-strain tensor, we still have the freedom to choose the most convenient form of the eddy viscosity. This is done by choosing what representation of the SGS stresses to model. The ‘optimal’ eddy viscosity constructed from the alternative SGS stresses found in § 6 may not be easy to model only with information from the resolved scales or may fail to reproduce other crucial aspects of turbulence dynamics. But it might be possible to find representations of the SGS stresses that yield a 
$\nu _s^{opt}$ that is backscatterless, physical in a general sense and practical from a modelling perspective, i.e. that can be appropriately reproduced using easy-to-compute quantities of the resolved scales, for instance, the filtered velocity gradients (Johnson Reference Johnson2020).
7.2. Non-equilibrium statistical mechanics and the local direction of the energy cascade
The approach followed here to show that backscatter is not an essential feature of the SGS fluxes may seem fundamentally mathematical, and somehow independent of the physics of the cascade, which has not been addressed directly up to now. This section discusses the findings of this paper on physical grounds.
In the context of statistical mechanics, turbulence is a strongly dissipative and non-equilibrium phenomenon. The energy cascade is a fundamental consequence of this non-equilibrium state, and its direction from large to small scales can be justified on the basis of equilibrium statistical mechanics. It can be shown, by analysing the truncated Euler equations, which reproduce inertial-range dynamics under certain conditions (She & Jackson Reference She and Jackson1993; Cichowlas et al. Reference Cichowlas, Bonaïti, Debbasch and Brachet2005), that a direct energy cascade from large to small scales is consistent with the tendency of isolated systems to evolve towards absolute equilibrium (Orszag Reference Orszag1974; Vela-Martín & Jiménez Reference Vela-Martín and Jiménez2021). This equilibrium is achieved when energy is equally distributed across all Fourier modes. Due to the fact that more degrees of freedom (Fourier modes) represent the small than the large scales, the evolution towards equilibrium of turbulence with a Kolmogorov spectrum leads to a direct energy cascade.
This argument is used to justify the direction of the energy cascade, but it is probabilistic and does not exclude the possibility of inverse cascades. In the inertial range, the effect of viscous diffusion is negligible and the dynamics is fully time-reversible. Thus interscale interactions are essentially bidirectional: the same mechanisms that operate to produce a direct energy cascade may operate to produce an inverse cascade. In fact, as shown by Carati et al. (Reference Carati, Winckelmans and Jeanmart2001) and Vela-Martín & Jiménez (Reference Vela-Martín and Jiménez2021), sustained global inverse cascades (in which 
$\langle T\rangle <0$) are possible in the inertial range. Although possible, global inverse cascades can be observed for only a limited time and appear not spontaneously, but due to carefully selected initial conditions (Vela-Martín & Jiménez Reference Vela-Martín and Jiménez2021). This limitation reflects the broken temporal symmetry of inertial-range dynamics and the emergence of statistical irreversibility, a consequence of the chaotic interactions of many degrees of freedom. Inverse energy cascades are highly improbable because they move turbulence away from the path towards absolute equilibrium by arranging the energy in a few large-scale structures. In fact, Vela-Martín & Jiménez (Reference Vela-Martín and Jiménez2021) showed that, unlike direct cascades, global inverse energy cascades require the organisation of the full flow. This is more unlikely the more chaotic the flow and the larger its number of degrees of freedom. As a consequence, only global direct energy cascades are observable empirically.
The above discussion is restricted to isotropic 3-D turbulence, which conserves (in the inviscid case) one positive-definite invariant, the kinetic energy. Other flows with cascades, such as 2-D turbulence (Rutgers Reference Rutgers1998), homochiral turbulence (Biferale, Musacchio & Toschi Reference Biferale, Musacchio and Toschi2012) or magnetohydrodynamic turbulence (Christensson, Hindmarsh & Brandenburg Reference Christensson, Hindmarsh and Brandenburg2001), conserve more positive-definite invariants apart from the energy – enstrophy, sign-definite helicity and magnetic helicity, respectively – and display inverse cascades that can be predicted by similar non-equilibrium arguments (Kraichnan Reference Kraichnan1967; Frisch et al. Reference Frisch, Pouquet, Léorat and Mazure1975; Alexakis & Biferale Reference Alexakis and Biferale2018).
Despite the unlikeliness of global inverse cascades in 3-D turbulence, the time-reversibility of the inertial-range dynamics makes it conceivable that local inverse cascades involving a reduced number of degrees of freedom may be observed for limited times. These local inverse cascades could appear often enough and for sufficiently long times to have practical consequences. This idea was supported in part by the phenomenon of backscatter, which has been often associated with an inverse cascade of energy. However, the results presented here rule out this scenario. That backscatter can be almost completely eliminated suggests that the limitations that operate to prevent global inverse cascades operate also locally on the SGS fluxes. Due to its out-of-equilibrium nature, the energy cascade is unidirectional enough to prevent local inverse energy cascades in time scales of practical relevance. Thus, although possible, local inverse cascade are far too unlikely to require modelling; irreversible dissipative SGS models reproduce appropriately the statistical irreversibility of the energy cascade.
8. Conclusions
In summary, the results presented in this paper show that the so-called backscatter, defined as 
$T<0$ (see § 2.1), does not represent a net inverse flux of energy from the small to the large scales, but conservative spatial fluxes. Thus it need not be modelled to reproduce the local energy exchange between the subgrid and resolved scales in large-eddy simulations (LES) of turbulence.
This has been shown by exploiting the invariance of the filtered Navier–Stokes equations (NSE) to divergence-preserving transformations of the subgrid-scale (SGS) stress tensor. Because the SGS stress tensor enters the filtered NSE as a divergence, we can add to it any divergence-free tensor to obtain an alternative representation the SGS stress tensor that produces the same physical effect on the NSE (Jiménez Reference Jiménez2016). Any alternative representation of the SGS stress tensor is physically valid and produces the same average SGS flux. However, different representations of the SGS stress tensor yield, in general, different local SGS fluxes.
Numerical optimisation procedures have been used to find alternative representations of the SGS stress tensor that produce SGS fluxes with negligible backscatter. In the equation of the resolved kinetic energy, these alternative SGS stresses represent backscatter as conservative spatial fluxes instead of as a source term. These conservative fluxes produce zero net energy transfer, and are not directly related to the energy cascade. By analysing an extensive database at different Reynolds numbers and filter scales, it has been shown that these results are general for the SGS fluxes in the inertial range of isotropic turbulence.
The results presented here justify the extended use of purely dissipative SGS models for LES of turbulence away from solid boundaries, and when the filter scale is small enough not to be strongly affected by large-scale anisotropies. In the case of strongly anisotropic scales in which inverse energy fluxes may be relevant, for instance in shear- or rotation-dominated flows (Campagne et al. Reference Campagne, Gallet, Moisy and Cortet2014; Deusebio et al. Reference Deusebio, Boffetta, Lindborg and Musacchio2014; Dong et al. Reference Dong, Huang, Yuan and Lozano-Durán2020), the approach presented in this paper could provide a better estimate of the exact ratio of direct to inverse energy transfer, as well as a better quantitative local measure of these fluxes and of their relation to the local structure of the flow (Buzzicotti et al. Reference Buzzicotti, Aluie, Biferale and Linkmann2018a; Dong et al. Reference Dong, Huang, Yuan and Lozano-Durán2020). This may help to improve LES models, in particular by reducing as much as possible the amount of unstable backscatter – if any – that needs to be reproduced by the SGS model.
In the same way that it is not necessary to model backscatter, it is not necessary to model the strong positive fluctuations of the SGS fluxes. This idea is already present in the so-called minimum dissipation models (Verstappen Reference Verstappen2011; Rozema et al. Reference Rozema, Bae, Moin and Verstappen2015), which prescribe the minimum eddy viscosity required to keep the local energy balance of the resolved scales. Moreover, the SGS fluxes found in this investigation have a fourth-order moment that is almost constant with the filter scale, suggesting that, as noted by Jiménez (Reference Jiménez2016), the intermittency of the SGS fluxes need not be modelled either.
The approach followed here to reduce backscatter is not based on physical or modelling arguments, but is of a demonstrative nature. Backscatter need not be modelled, but it may as well be. This leeway is advantageous for cases in which backscatter is problematic, such as in eddy-viscosity models. However, eliminating backscatter may not be necessary, or even undesirable, in models that produce numerically stable backscatter, such as the similarity or tensor-diffusivity models.
The same applies to the intense positive fluctuations or the intermittency of the SGS fluxes. The latter case is particularly relevant and illustrative in this discussion. It is widely accepted that intermittency is an essential property of inertial-range dynamics (Sreenivasan & Antonia Reference Sreenivasan and Antonia1997). Yet we have shown that it can be removed from the SGS fluxes, although it remains as part of the total SGS interactions (§ 6.2). This may suggest, erroneously, that intermittency need not be considered in the modelling of the SGS stress tensor. This argument stems from a modelling approach based only on the energetics of the cascade, but reproducing exactly the local energy fluxes may not be enough to produce a ‘good’ SGS model. An example is the eddy-viscosity model in § 7.1, which is ‘perfect’ in the sense that it produces the exact SGS fluxes, but may fail to reproduce other crucial properties of turbulence in the inertial range. We have presented physical reasons why we believe that backscatter may be unnecessary from a modelling perspective (see § 7.2). Conversely, we have no physical reasons to think the same of intermittency, which is a more general property of turbulence dynamics than backscatter.
What other properties a ‘good’ SGS model should reproduce remains an open question beyond the scope of this paper, but this work shows that we have more freedom when looking for ‘good’ models than was commonly believed. This work stresses the practical importance of gauge invariance for LES modelling (Jiménez Reference Jiménez2016). The fact that the SGS stress tensor enters the NSE inside a divergence implies that some of its degrees of freedom are not physical but mathematical, and do not have a directly measurable effect on the filtered-momentum equations. Thus it is ineffective to model a particular representation of the SGS stresses when it is only necessary to model its physical effect on the flow, i.e. its divergence. This also implies that a poor correlation between the standard SGS stress tensor and the modelled SGS stress tensor is not necessarily an indication of poor modelling. In fact, it has been shown here that equivalent physical representations of the SGS stress tensor are not highly correlated (see § 6.3). An SGS model may not reproduce well a particular representation of the SGS stress tensor, yet its physical effect on the flow could be reproduced appropriately. This seems to be the case, to some extent, with eddy-viscosity models (Clark et al. Reference Clark, Ferziger and Reynolds1979; Meneveau Reference Meneveau1994; Borue & Orszag Reference Borue and Orszag1998). The freedom in modelling any physical representation of the SGS stresses can be leveraged to improve and develop SGS models in a variety of turbulent flows, and to conveniently choose well-behaved or practical models, for instance, those that do not produce numerically unstable backscatter.
From a physical perspective, the findings of this paper have been justified by resorting to statistical mechanics. In the context of non-equilibrium dynamics, global inverse cascades contravene the tendency of out-of-equilibrium systems to evolve towards equilibrium, and, although possible, are highly unlikely due to statistical irreversibility (Vela-Martín & Jiménez Reference Vela-Martín and Jiménez2021). That backscatter can be almost completely eliminated from the SGS fluxes suggests that the limitations imposed by statistical irreversibility to global cascades also apply to local inverse cascades, which seem to be unlikely enough to have little practical importance in the modelling of the cascade.
Let us emphasise here the fundamental difference between space-local inverse energy cascades and backscatter. The latter strictly represents negative values of the SGS fluxes, 
$T<0$, as defined in (2.7), while the former imply a space-local inverse energy flux due to interactions across scales. As shown here, backscatter is not, in general, related to a local inverse cascade of energy in fully developed 3-D isotropic turbulence. From a fundamental perspective, backscatter may represent a local inverse cascade in special cases, for instance in the numerical experiments by Carati et al. (Reference Carati, Winckelmans and Jeanmart2001) and Vela-Martín & Jiménez (Reference Vela-Martín and Jiménez2021) discussed in § 7.2. Note that in these numerical experiments, it is not possible to reduce backscatter in 
$T$ following the procedure described here because 
$\langle T\rangle <0$. It seems inconceivable that there exists a representation of 
$T$ that does not fluctuate in space and time. This implies that backscatter cannot always be completely eliminated because fluctuations may always allow for the spontaneous – but unlikely – generation of local or global inverse cascades. This idea is formally described by the fluctuation relations (Evans & Searles Reference Evans and Searles2002), and is consistent with the microscopic reversibility of the energy cascade and with the fact that inertial-range dynamics admits inverse cascades (Vela-Martín & Jiménez Reference Vela-Martín and Jiménez2021). However, from a practical perspective, the fluctuations of 
$T$ leading to space-local inverse cascades are unlikely enough to be of practial relevance.
The limited practical relevance of inverse cascades agrees with the idea of a strongly unidirectional turbulence cascade, in which energy (and probably also ‘information’) flows preferentially from large to small scales (Meneveau & Lund Reference Meneveau and Lund1994; Cardesa et al. Reference Cardesa, Vela-Martín and Jiménez2017; Ballouz, Johnson & Ouellette Reference Ballouz, Johnson and Ouellette2020; Vela-Martín Reference Vela-Martín2021). As a consequence, the interactions between the subgrid and resolved scales can be modelled as an irreversible sink of energy. This idea is akin to the successful modelling of microscopically reversible phenomena, such as momentum diffusion or heat conduction, by irreversible macroscopic models. The unidirectionality of the energy cascade is perhaps the reason why its modelling is possible despite its complex dynamics. It is difficult to conceive that the LES approach would have been so successful if the energy cascade was a strongly bidirectional process.
The results presented here may suggest that the definition of the local SGS fluxes is fraught with ambiguity. However, this is true only when the SGS fluxes are analysed from a strictly local perspective, in which case they are affected strongly by conservative spatial fluxes. As shown by Vela-Martín & Jiménez (Reference Vela-Martín and Jiménez2021), these spatial fluxes can be eliminated by averaging carefully the SGS fluxes, leading to more universal markers of the energy cascade. Modifying the definition of the SGS stresses to eliminate these fluxes has a similar effect. For instance, the strong differences in the SGS fluxes for different filter types (Cardesa et al. Reference Cardesa, Vela-Martín, Dong and Jiménez2015; Buzzicotti et al. Reference Buzzicotti, Linkmann, Aluie, Biferale, Brasseur and Meneveau2018b) disappear when the conservative spatial fluxes are removed from the SGS fluxes by modifying the SGS stress tensor (see figure 5b).
Finally, let us note that we have limited here the analysis of gauge transformations in the SGS stress tensor to a modelling perspective. We suggest that a study of the energy cascade from a purely physical perspective, and beyond modelling purposes, should consider gauge transformations motivated by physical arguments, and in which the properties of the cascade are taken into account (e.g. scale and space locality, or intermittency). In this sense, the standard definition of the SGS fluxes may be an appropriate marker of the local interscale fluxes, particularly after local averaging, although its robustness should be tested against other possible definitions (see Osawa & Jiménez (Reference Osawa and Jiménez2018) for an example).
Here, we have considered gauge transformations only in physical space, not in scale space. The latter should be also considered for descriptions of the cascade that require approximate locality both in scale and in physical space (Meneveau Reference Meneveau1991; Farge Reference Farge1992).
Acknowledgements
The author would like to thank J.I. Cardesa, M.P. Encinar and M. Avila for fruitful discussions that led to the improvement of the manuscript.
Declaration of interest
The author reports no conflict of interest.
Appendix A. Minimising the standard deviation of the subgrid-scale fluxes
Let us consider the functional
\begin{equation} \mathcal{L} = \int_V [ T^{2} + \gamma_i (\partial_{j}\tau_{ij} - Q_i) + \mu_{ij}(\tau_{ij} - \tau_{ji})]\,\mathrm{d}V, \end{equation}
where 
$T=\tau _{ij}\bar {S}_{ij}$, and 
$\gamma _i$ and 
$\mu _{ij}$ are Lagrange multipliers that impose
\begin{equation} \left.\begin{gathered} \partial_j \tau_{ij}=Q_i,\\ \tau_{ij}-\tau_{ji}=0, \end{gathered}\right\} \end{equation}
where 
$Q_i=\partial _j\tau ^{S}_{ij}$. To find the minimum of 
$\mathcal {L}$, we introduce a variation on 
$\tau _{ij}$, 
$\delta \tau _{ij}$,
\begin{equation} \delta\mathcal{L} = \int_V [ 2T\bar{S}_{ij}\,\delta \tau_{ij} + \gamma_i\, \partial_{j}\delta\tau_{ij} + \mu_{ij}(\delta\tau_{ij} - \delta\tau_{ji})]\,\mathrm{d}V, \end{equation}
and impose that the variation of the functional is zero, 
$\delta \mathcal {L}=0$, regardless of the form of 
$\delta \tau _{ij}$. Considering the relations
\begin{equation} \gamma_i\,\partial_{j}\delta\tau_{ij}=\partial_j(\gamma_i \delta\tau_{ij}) - \partial \gamma_i\,\delta \tau_{ij} \end{equation}and
\begin{equation} \mu_{ij}(\delta\tau_{ij} - \delta\tau_{ji})=(\mu_{ij} - \mu_{ji})\delta \tau_{ji}, \end{equation}and imposing periodic boundary conditions, we obtain
\begin{equation} \delta\mathcal{L} = \int_V [ 2T\bar{S}_{ij} - \partial_{j}\gamma_i + (\mu_{ij} - \mu_{ji})]\,\delta\tau_{ij}\,\mathrm{d}V. \end{equation}
For 
$\delta \mathcal {L}=0$ independently of 
$\delta \tau _{ij}$,
\begin{equation} 2T\bar{S}_{ij} - \partial_{j}\gamma_i + (\mu_{ij} - \mu_{ji})=0. \end{equation}
By adding this equation to its transpose, we eliminate 
$\mu _{ij}$, and obtain a set of equations for 
$\tau _{ij}$ and 
$\gamma _{i}$:
\begin{equation} \tfrac{1}{2}(\partial_{j}\gamma_i + \partial_{i} \gamma_j)=2T\bar{S}_{ij}. \end{equation}These six equations, together with the three equations of one of the constraints,
\begin{equation} \partial_j\tau_{ij}-Q_i=0, \end{equation}
determine the three components of 
$\gamma _i$ and the six components of 
$\tau _{ij}$.
Appendix B. The linear projector in the gradient descent algorithm
Here, we derive the linear projector used in the gradient descent algorithm,
\begin{equation} \boldsymbol{\chi}^{P}=\mathcal{P}(\boldsymbol{\chi}^{*}), \end{equation}
where 
$\boldsymbol {\chi }=\{\tau _{11}(\boldsymbol {x}), \tau _{22}(\boldsymbol {x}),\dots,\tau _{ij}(\boldsymbol {x})\}$, the asterisk denotes the test state vector, and the superscript ‘
$P$’ denotes its projection. The projection is accomplished by minimising the functional
\begin{equation} \mathcal{R} = \int_V \left[\tfrac{1}{2}(\tau^{*}_{ij} - \tau^{P}_{ij})^{2} + \gamma_i (\partial_{j}\tau^{P}_{ij} - Q_i) + \mu_{ij}(\tau^{P}_{ij} - \tau^{P}_{ji})\right]\mathrm{d}V, \end{equation}
where the first term in the integral is the 
$L^{2}$ distance between the projected tensor and the test tensor, and 
$\gamma _i$ and 
$\mu _{ij}$ are Lagrange multipliers that impose
\begin{equation} \left.\begin{gathered} \partial_j \tau_{ij}^{P}=Q_i,\\ \tau_{ij}^{P}-\tau_{ji}^{P}=0, \end{gathered}\right\} \end{equation}
where 
$Q_i=\partial _j\tau ^{S}_{ij}$.
Introducing a variation on 
$\tau _{ij}^{P}$, 
$\delta \tau _{ij}^{P}$, and operating, we obtain the variation of the functional
\begin{equation} \delta\mathcal{R} = \int_V [(\tau^{*}_{ij} - \tau^{P}_{ij})\,\delta \tau_{ij}^{P} + \gamma_i\, \partial_{j}\delta\tau^{P}_{ij} + \mu_{ij}(\delta\tau^{P}_{ij} - \delta\tau^{P}_{ji})]\,\mathrm{d}V. \end{equation}Considering the relations
\begin{equation} \gamma_i\,\partial_{j}\delta\tau^{P}_{ij}=\partial_{j}(\gamma_i\delta\tau^{P}_{ij}) - \partial_j \gamma_i\,\delta \tau^{P}_{ij} \end{equation}and
\begin{equation} \mu_{ij}(\delta\tau^{P}_{ij} - \delta\tau^{P}_{ji})=(\mu_{ij} - \mu_{ji})\,\delta \tau^{P}_{ji}, \end{equation}and imposing periodic boundary conditions, we obtain
\begin{equation} \delta\mathcal{R} = \int_V (\tau^{*}_{ij} - \tau^{P}_{ij} -\partial_j \gamma_i + \mu_{ij} - \mu_{ji})\,\delta\tau^{P}_{ij}\,\mathrm{d}V. \end{equation}
For 
$\delta \mathcal {R}$ to be zero independently of 
$\delta \tau _{ij}^{P}$, we have that
\begin{equation} \tau^{*}_{ij} - \tau^{P}_{ij} -\partial_j \gamma_i + \mu_{ij} - \mu_{ji}=0. \end{equation}
We eliminate 
$\mu _{ij}$ by adding this equation to its transpose, and obtain
\begin{equation} \partial_i \gamma_j + \partial_j \gamma_i = 2(\tau^{P}_{ij} - \tau^{*}_{ij}), \end{equation}
where we have assumed that the test SGS stress tensor is symmetric. By taking derivatives of this equation and using the constraint 
$\partial _i \tau ^{P}_{ij}=Q_i$, we arrive at an equation for 
$\gamma _i$:
\begin{equation} \partial_{kk} \gamma_i + \partial_i (\partial_k\gamma_k) = 2(\partial_{i}\tau^{*}_{ij} - Q_i). \end{equation}We now use the identity
\begin{equation} \partial_{kk} \gamma_i = \partial_i (\partial_k \gamma_k) - \epsilon_{ijk}\epsilon_{klm}\,\partial_{j}\partial_{l}\gamma_{m}, \end{equation}
where 
$\epsilon _{ijk}$ is the antisymmetric Levi-Civita symbol. The last term on the right-hand side corresponds to minus the double rotor of 
$\boldsymbol {\gamma }=\{ \gamma _i\}$, which in vectorial notation reads 
$-\boldsymbol {\nabla }\times (\boldsymbol {\nabla }\times {\boldsymbol {\gamma }})$. Substituting this expression in (B11) and taking the divergence leads to an equation for 
$\varGamma =\partial _k \gamma _k$:
\begin{equation} \partial_{kk} \varGamma = \partial_i \varLambda_i, \end{equation}
where 
$\varLambda _i = \partial _j \tau ^{*}_{ij} - Q_i$. We solve this equation imposing periodic boundary conditions, and substitute 
$\varGamma$ in (B10), leading to
\begin{equation} \partial_{kk} \gamma_i = 2\varLambda_i - \partial_i \varGamma. \end{equation}
Solving for 
$\gamma _i$, and substituting in (B9), leads to
\begin{equation} \tau^{P}_{ij} = \tau^{*}_{ij} - \tfrac{1}{2}(\partial_i \gamma_j + \partial_j \gamma_i), \end{equation}which allows us to readily obtain the projected SGS stress tensor.
Open access
