1. Introduction
1.1. Bypass transition
Transition to turbulence that does not involve linear instability paths, such as Tollmien–Schlichting waves, is called bypass transition (Morkovin Reference Morkovin1969). This type of transition can be triggered by high levels of free-stream turbulence, surface roughness, etc. In the case of free-stream turbulence, which is the triggering mechanism considered in this paper, bypass transition comprises three stages. In the first stage, low-frequency fluctuations from the free stream penetrate inside the boundary layer, forming high- and low-speed streaks, while high-frequency fluctuations remain in the free stream due to shear sheltering (Hunt & Durbin Reference Hunt and Durbin1999; Leib, Wundrow & Goldstein Reference Leib, Wundrow and Goldstein1999; Zaki & Saha Reference Zaki and Saha2009). In the second stage, the streaks break down to intermittent turbulent patches (or spots) due to secondary instability (Andersson et al. Reference Andersson, Brandt, Bottaro and Henningson2001; Vaughan & Zaki Reference Vaughan and Zaki2011), while in the final stage the spots propagate and merge, forming a fully turbulent region. More details can be found in the review papers of Durbin & Wu (Reference Durbin and Wu2007), Zaki (Reference Zaki2013) and Durbin (Reference Durbin2017).
Most previous investigations of the structural details of turbulent spots, such as shape, propagation speed, growth rate, etc., have employed analysis of single-point statistics (see Emmons Reference Emmons1951; Wygnanski, Sokolov & Friedman Reference Wygnanski, Sokolov and Friedman1976; Cantwell, Coles & Dimotakis Reference Cantwell, Coles and Dimotakis1978; Perry, Lim & Teh Reference Perry, Lim and Teh1981; Singer Reference Singer1996; Nolan & Zaki Reference Nolan and Zaki2013; more recently Wang et al. Reference Wang, Choi, Gaster, Atkin, Borodulin and Kachanov2021, Reference Wang, Choi, Gaster, Atkin, Borodulin and Kachanov2022). This type of analysis, however, cannot capture the underlying physical mechanisms that explain the spot growth and the amalgamation process of smaller spots to form larger turbulent patches as transition progresses. In order to study this process in more detail, an analysis of two-point statistics is required. The second-order structure function at point 
$X_i$ is defined as the second moment of the fluctuating velocity difference at points 
$x^\pm _i=X_i\pm \frac {1}{2}r_i$, i.e. 
$\overline {{\rm d} q^2}(X_i,r_i)=\overline {(u'^+_i-u'^-_i)^2}$, where the overbar denotes time averaging. The volume integral of 
$\overline {{\rm d} q^2}$ over a sphere of radius 
$r=|r_i|$ (divided by the volume of the sphere) represents physically the energy of eddies located at 
$X_i$ that have size (or scale) less than 
$r$; this is also known as scale energy (for details, see Davidson (Reference Davidson2015)). This is the appropriate quantity to study in order to better understand the process of spot growth and merging.
The transport equation of 
$\overline {{{\rm d} q^2}}(X_i,r_i)$ is known as the Kármán–Howarth–Monin–Hill (KHMH) equation. It contains all the physical mechanisms that determine the energy contained within eddies of scale less than 
$r$, such as transfer of energy in scale space (i.e. from smaller or larger scales), transfer of energy in physical space, production (due to mean shear), dissipation (due to viscosity), etc. It was first derived by Kármán & Howarth (Reference Kármán and Howarth1938) for homogeneous isotropic turbulence using the two-point velocity correlation tensor, 
$\overline {u'^+_iu'^-_j}$, and later reformulated in terms of structure function 
$\overline {{\rm d} q^2}(r_i)$ by Kolmogorov (Reference Kolmogorov1941). The equation was used to prove the famous 
$-4/5$ law, which links the interscale energy flux, the separation between the two points and the dissipation rate. For homogeneous isotropic turbulence the interscale flux is always negative, i.e. energy is transferred from large to small scales; this is known as forward cascade. The most general form of the KHMH equation, applicable to inhomogeneous and anisotropic flows, was derived directly from Navier–Stokes equations by Hill (Reference Hill2002). The equation was applied recently to transitional boundary layers and demonstrated strong inverse cascade in the transition region, especially in the streamwise direction (see Yao, Mollicone & Papadakis Reference Yao, Mollicone and Papadakis2022). Analysis of instantaneous velocity fields and flux vectors revealed that the inverse cascade was related to the growth of turbulent spots.
The KHMH equation has been applied to several other flow settings, such as homogeneous shear flow (Casciola et al. Reference Casciola, Gualtieri, Benzi and Piva2003), channel flow (Marati, Casciola & Piva Reference Marati, Casciola and Piva2004; Cimarelli, De Angelis & Casciola Reference Cimarelli, De Angelis and Casciola2013; Cimarelli et al. Reference Cimarelli, De Angelis, Schlatter, Brethouwer, Talamelli and Casciola2015b, Reference Cimarelli, De Angelis, Jimenez and Casciola2016), von Kármán flow (Dubrulle Reference Dubrulle2019; Knutsen et al. Reference Knutsen, Baj, Lawson, Bodenschatz, Dawson and Worth2020), temporal planar jets (Cimarelli et al. Reference Cimarelli, Mollicone, Van Reeuwijk and De Angelis2021) and wake behind a square cylinder (Alves Portela, Papadakis & Vassilicos Reference Alves Portela, Papadakis and Vassilicos2017, Reference Alves Portela, Papadakis and Vassilicos2020), and has revealed convoluted paths of interscale energy fluxes with a mixture of forward and inverse cascade. In all these investigations, the points 
$x^\pm _i$ used to define 
$\overline {{\rm d} q^2}$ are immersed within a turbulent flow. This is not the case, however, in transitional flows, where for a fixed streamwise location 
$X_i$ and separation 
$r_i$, the two points will experience different flow conditions as a spot propagates; for example, 
$x^\pm _i$ may be within the laminar region, straddle the laminar/turbulent interface or be within a turbulent patch. It is therefore difficult to explain the origin of the inverse cascade found in Yao et al. (Reference Yao, Mollicone and Papadakis2022). For example, does it arise from the turbulent conditions within the spot? What is the role of the laminar/turbulent interface? Is there competition between the different flow conditions, i.e. do some lead to forward and others to inverse cascade?
To answer these questions in the present paper we perform conditional averaging of the interscale energy fluxes based on the state of the two points, i.e. whether they experience laminar or turbulent flow conditions. This process clearly elucidates the effect of different flow states and the role of the laminar/turbulent interface. As seen later, it even characterises the separate roles of the downstream and upstream interfaces of the spot. We also examine the production term of the KHMH equation conditioned on turbulent events within a turbulent patch, and compare it with the production term when the two points are located within the fully turbulent region. Similar comparisons have been made for single-point quantities, such as turbulent kinetic energy (see Marxen & Zaki Reference Marxen and Zaki2019). We derive also the conditionally averaged form of the KHMH equation, which is analogous to the conditionally averaged turbulent kinetic energy equation.
A few studies in the literature have similarities to, but also differences from, the present work. Cimarelli et al. (Reference Cimarelli, Cocconi, Frohnapfel and De Angelis2015a) employed the spectral enstrophy budget equation to study the enstrophy transfer between scales close to the interface between decaying shear-free turbulence and quiescent fluid. They found anisotropic cascade, where large-scale structures exhibit only a cascade process normal to the interface, thus reducing their thickness while retaining their lengths parallel to the interface. Zhou & Vassilicos (Reference Zhou and Vassilicos2020) studied the energy cascade across the turbulent/non-turbulent interface (TNTI) at one axial position of an axisymmetric turbulent wake. They found that the interscale energy transfer at the vicinity of the interface is from small to large scales (inverse cascade) in directions close to the interface's tangent plane where motions are predominantly stretching, but from large to small scales (forward cascade) in the other directions where motions are predominantly compressive. This reflects the fundamental mechanism that sustains the TNTI, i.e. fluid is entrained from the irrotational region and the wake grows due to turbulent diffusion (see schematic 3(d) in Zhou & Vassilicos (Reference Zhou and Vassilicos2020)). Cimarelli et al. (Reference Cimarelli, Mollicone, Van Reeuwijk and De Angelis2021) studied the interscale energy transfer process in a temporal planar jet using the KHMH equation. They also found reverse energy cascade that is responsible for the generation of long and wide structures in the interface region. All the aforementioned works have focused on the interscale energy transfer mechanism that sustains the TNTI. This mechanism, however, is different compared with the one that determines the growth of spots in a transitional boundary layer. Conditional analysis was not performed either.
The paper is organised as follows. In § 2 the bypass transition case is briefly presented, in § 3 we summarise the derivation of the standard time-averaged form of the KHMH equation, while in § 4 we derive the conditionally averaged form; this is followed by the conditional decomposition of the energy fluxes in § 5. The next two sections present the results: in § 6 maps of the conditionally averaged nonlinear energy fluxes are shown (focusing on the flux across the laminar/turbulent interface), while in § 7 we compare the two-point energy production and total flux (both conditioned on turbulent events) with the corresponding quantities in the fully turbulent region. We conclude in § 8.
2. Details of the test case examined
We consider the transition of a boundary layer developing on a flat plate due to free-stream turbulence. At the inlet of the computational domain, a random velocity field is superimposed on the Blasius velocity profile. In the free stream, the random field follows a von Kármán spectrum with turbulence intensity 
$3.4\,\%$ and integral length scale 
$L_{11}=5L_0$, where 
$L_0 = \sqrt {\nu X_0/U_\infty }$ is the Blasius similarity variable, 
$X_0$ is the distance of the inlet of the domain from the leading edge of the plate, 
$\nu$ is the kinematic viscosity and 
$U_\infty$ is the free-stream velocity. The inlet Reynolds number is 
$Re_{L_0}= 160$ (or 
$Re_\theta = 110$ based on momentum thickness).
The size of the computational domain is 
$(3000 \times 200 \times 150)L_0$, with the number of cells 
$2049\times 192\times 169$ in the streamwise (
$X$), wall-normal (
$Y$) and spanwise (
$Z$) directions. Velocities are denoted as 
$u$, 
$v$, 
$w$ in the 
$X,Y,Z$ directions respectively. This notation is used interchangeably with the indexed notation 
$X_i$ and 
$u_i$ (with 
$i=1,2,3$); for example 
$X_2=Y$ and 
$u_2=v$. The spacing is uniform in the streamwise and spanwise directions, with 
${\rm \Delta} x^+_{max}\approx 11.78$ and 
${\rm \Delta} z^+_{max}\approx 7.14$, where subscript ‘
${max}$’ represents the maximum value (located in the fully turbulent region). In the wall-normal direction, grid spacing increases gradually; 
$y^+$ at the centroid of the first cell close to the wall is around 0.24. The generated direct numerical simulation database contains 350 uncorrelated snapshots. The results have been validated against the T3A experimental data (Roach & Brierley Reference Roach and Brierley1992). More details about the computational method and comparison of velocity profiles (mean and root mean square) against experiments can be found in Yao, Alves-Portela & Papadakis (Reference Yao, Alves-Portela and Papadakis2020).
For future reference, the normalised skin friction coefficient, 
$C_f/\max (C_f)$, and the maximum value of the time- and spanwise-averaged intermittency are plotted in figure 1. For the methodology used to compute the intermittency, the reader is referred to the aforementioned paper and also to § 4.1. Vertical lines indicate the streamwise locations in the laminar (
$LA$), transitional (
$TR1,TR2, TR3$) and fully turbulent (
$TU$) regions, where velocity data are extracted in order to compute the interscale fluxes.
Figure 1. Normalised skin friction coefficient, 
$C_f/\max (C_f)$ (blue solid line), maximum intermittency 
$\gamma _{max}$ (red solid line) and 
$\gamma (\xi )$ from the formula of Narasimha (Reference Narasimha1985), where 
$\xi =(X-X_s)/(X_{\gamma =0.75}-X_{\gamma =0.25})$ and 
$X_s-X_0=1100L_0$ is the location where transition starts (black dashed line), plotted against streamwise distance 
$(X-X_0)/L_0$. The purple vertical lines are located in the laminar (
$LA=540L_0+X_0$), transitional (
$TR1,TR2,TR3=(1215,1515,1815)L_0+X_0$) and fully turbulent (
$TU=2415L_0+X_0$) regions.
3. Standard time-averaged KHMH equation
In this section, the form of the standard time-averaged KHMH equation is presented. The basic steps of the derivation are sketched below; more details can be found in Hill (Reference Hill2002). Similarities to and differences from the conditionally averaged form are presented and discussed in § 4.
We start with the Navier–Stokes equations at two points 
$x_i^+$ and 
$x_i^-$ (see figure 2):
$$\begin{gather} \frac{\partial u^+_i}{\partial t} + u^+_j \frac{\partial u^+_i}{\partial x^+_j} ={-}\frac{\partial p^+}{\partial x^+_i} + \nu \frac{\partial^2 u^+_i}{\partial x_j^+ \partial x_j^+}, \end{gather}$$
$$\begin{gather}\frac{\partial u^-_i}{\partial t} + u^-_j \frac{\partial u^-_i}{\partial x^-_j} ={-}\frac{\partial p^-}{\partial x^-_i} + \nu \frac{\partial^2 u^-_i}{\partial x_j^- \partial x_j^-}, \end{gather}$$and define the time- and spanwise-averaged velocities as usual:
$$\begin{gather} U^+_i(X,Y)=\overline{u^+_i}= \frac{1}{{\rm \Delta} T L_z} \int_0^{{\rm \Delta} T} \int_0^{L_z} u_i^+ \, {\rm d} z \, {\rm d} t, \end{gather}$$
$$\begin{gather}U^-_i(X,Y)=\overline{u^-_i}= \frac{1}{{\rm \Delta} T L_z} \int_0^{{\rm \Delta} T} \int_0^{L_z} u_i^- \, {\rm d} z \, {\rm d} t. \end{gather}$$
In the following, we use an overbar to denote the standard averaging operation in the time and 
$Z$ direction, as defined by (3.2). Velocity fluctuations around 
$U^+_i$ and 
$U^+_i$ are denoted using primes, i.e.
\begin{equation} u'^+_i=u^+_i-U^+_i, \quad u'^-_i=u^-_i-U^-_i,\end{equation}and fluctuating velocity differences are defined as
\begin{equation} {\rm d} u_i' \equiv {\rm d} u_i- {\rm d} U_i,\end{equation}
where 
${\rm d} u_i=u_i^+-u_i^-$ and 
${\rm d} U_i=\overline {{\rm d} u_i}=U^+_i-U^-_i$. It is straightforward to prove that 
${\overline {{\rm d} u_i'}=0}$.
Figure 2. Sketch that shows the position vector 
${\boldsymbol {X}}=(X_1,X_2,X_3)$ of the midpoint and the separation vector 
$\boldsymbol {r}=(r_1,r_2,r_3)$ between the two points 
$\boldsymbol{x}^+=(x^+_1,x^+_2,x^+_3)$ and 
$\boldsymbol{x}^-=(x^-_1,x^-_2,x^-_3)$.
Subtracting equation (3.1b) from (3.1a), multiplying each term by 
$2\, {\rm d} u_i'$ and applying the time- and spanwise-averaging operation defined in (3.2), we obtain
\begin{align}
& \underbrace{ \overline{{2\, {\rm d} u_i'}\frac{\partial \,
{\rm d} u_i}{\partial t}}}_{\text{Transient term}} +
\underbrace{ \overline{2\, {\rm d} u_i'u^+_j \frac{\partial
\, {\rm d} u_i}{\partial x^+_j} + 2\, {\rm d} u_i'u^-_j
\frac{\partial \, {\rm d} u_i}{\partial
x^-_j}}}_{\text{Nonlinear term}}\nonumber\\
&\quad =\underbrace{ \overline{{-2\, {\rm d}
u_i'}\left(\frac{\partial \, {\rm d} p}{\partial x^+_i} -
\frac{\partial \, {\rm d} p}{\partial
x^-_i}\right)}}_{\text{Pressure--velocity
correlation}} + \underbrace{ \overline{{2 \nu \, {\rm d}
u_i'}\frac{\partial^2 \, {\rm d} u_i}{\partial x_j^+
\partial x_j^+} + 2\nu \, {\rm d} u_i'\frac{\partial^2 \,
{\rm d} u_i}{\partial x_j^- \partial x_j^-}}
}_{\text{Viscous term}},
\end{align}
where we have used 
${\partial u^+_i}/{\partial x^-_j}=0$ and 
${\partial u^-_i}/{\partial x^+_j}=0$ (because 
$x^+_j$ and 
$x^-_j$ are independent variables).
We now define the second-order structure function as 
$\overline {{\rm d} q^2} =\overline {({\rm d} u'_i)^2}=\overline {(u'^+_i-u'^-_i)^2}$. This function has six dimensions, three in physical space (
$X_i$) and three in scale space (
$r_i$). In the particular case examined in this paper, due to the homogeneity in the spanwise direction, there are only two dimensions in physical space. The integral of 
$\overline {{\rm d} q^2}$ in a sphere of radius 
$r=|r_i|$ (divided by the volume of the sphere) represents the energy of eddies with size smaller than 
$r=|r_i|$ (see Davidson Reference Davidson2015); thus 
$\overline {{\rm d} q^2}$ is usually referred to as scale energy.
We seek the transport equation of 
$\overline {{\rm d} q^2}$ in the physical and scale spaces. Applying the variable transformation 
$X_i=0.5(x^+_i + x^-_i)$ and 
$r_i = x^+_i - x^-_i$ and the definitions (3.3a,b) and (3.4) into (3.5), after some algebra we get the following standard KHMH equation for 
$\overline {{\rm d} q^2}$:
\begin{align}
& \underbrace{\overline{\frac{\partial \, {\rm d}
q^2}{\partial t} }}_{\text{Transient term}}+
\underbrace{\overline{ U^*_{j}\frac{\partial \, {\rm d}
q^2}{\partial X_j} }}_{\text{Mean flow advection}}+
\underbrace{\overline{ u'^*_{j}\frac{\partial \, {\rm d}
q^2}{\partial X_j} }}_{\text{Turbulent advection}}+
\underbrace{\overline{{\rm d} U_j \frac{\partial \, {\rm d}
q^2}{\partial r_j} }}_{\text{Linear interscale transfer}}+\underbrace{\overline{{\rm d}
u_j'\frac{\partial \, {\rm d} q^2}{\partial
r_j}}}_{\text{Nonlinear interscale transfer}} \nonumber\\
&\quad =\underbrace{\overline{-2\, {\rm d} u_i' \frac{\partial \,
{\rm d} p'}{\partial X_i}
}}_{\text{Pressure--velocity
correlation}} -\underbrace{\overline{2\, {\rm d} u_i'
u'^*_{j}\frac{\partial \, {\rm d} U_i}{\partial X_j}} -
\overline{2\, {\rm d} u_i' \, {\rm d} u_j' \frac{\partial
\, {\rm d} U_i}{\partial r_j}} }_{\text{Production by mean
flow } (=\mathcal{P})}\nonumber\\
&\qquad +\underbrace{\overline{\nu \frac{1}{2} \frac{\partial^2 \,
{\rm d} q^2}{\partial X_j \partial X_j} }}_{\text{Physical
diffusion}} +\underbrace{\overline{2\nu \frac{\partial^2 \,
{\rm d} q^2 }{\partial r_j \partial r_j}}}_{\text{Scale
diffusion}} -\underbrace{ 4\nu \left( \frac{1}{4}
\overline{ \frac{\partial \, {\rm d} u'_i}{\partial X_j}
\frac{\partial \, {\rm d} u'_i}{\partial X_j}} +
\overline{\frac{\partial \, {\rm d} u'_i}{\partial r_j}
\frac{\partial \, {\rm d} u'_i}{\partial r_j} }
\right)}_{\text{Dissipation } (=\epsilon)}.
\end{align} Note that 
$U^*_{j}$ and 
$u'^*_{j}$ denote the midpoint values of the time-average and fluctuating velocities respectively, i.e. 
$U^*_{j}=(U^+_j + U^-_j)/2$ and 
$u'^*_{j}=(u'^+_j + u'^-_j)/2$. The physical meaning of each term is also provided; the production by mean flow and dissipation are denoted by 
$\mathcal {P}$ and 
$\epsilon$ respectively.
Assuming that the transient term is 0, the above equation can be written in divergence form as
\begin{equation} \frac{\partial \boldsymbol{\phi}_{s_i} }{\partial X_i} + \frac{\partial \boldsymbol{\phi}_{r_i}}{\partial r_i} = \mathcal{P}-\epsilon,\end{equation}where
\begin{equation} \boldsymbol{\phi}_{s_i}=\underbrace{\overline{U_i^*d q^2} }_{=\phi_{s_i}^M} + \underbrace{\overline{u'^*_i \, {\rm d} q^2}}_{=\phi_{s_i}^F} + \underbrace{\overline{2 \, {\rm d} u'_i \, {\rm d} p'}}_{=\phi_{s_i}^P} \underbrace{-\frac{1}{2}\nu \frac{\partial \, \overline{{\rm d} q^2}}{\partial X_i}}_{=\phi_{s_i}^V} \end{equation}is the total flux vector in physical space and
\begin{equation} \boldsymbol{\phi}_{r_i}= \underbrace{\overline{{\rm d} u'_i \, {\rm d} q^2}}_{=\phi_{r_i}^F} + \underbrace{\overline{{\rm d} U_i \, {\rm d} q^2}}_{=\phi_{r_i}^M} \underbrace{-2\nu \frac{\partial \, \overline{{\rm d} q^2}}{\partial r_i}}_{=\phi_{r_i}^V} \end{equation}
is the total flux vector in scale space. We use the superscript ‘
$F$’ to denote the nonlinear components, 
$\phi _{s_i}^F$ and 
$\phi _{r_i}^F$, of these vectors respectively. The conditional decomposition of the nonlinear fluxes 
$\phi _{s_i}^F$ and 
$\phi _{r_i}^F$ is examined in § 5. We are now ready to proceed with the derivation of the conditionally averaged KHMH equation.
4. Conditionally averaged KHMH equation
4.1. Definitions
In order to derive the conditionally averaged form of the KHMH equation, we need first to define the conditions under which the averaging is performed. To this end, we employ the local binary function 
$\tau (x_i,t)$ to distinguish between instantaneous laminar and turbulent states at point 
$x_i$ at time 
$t$. More specifically, 
$\tau (x_i,t)$ takes the value of 
$0$ for the former (i.e. laminar) and 
$1$ for the latter (i.e. turbulent) state. This binary function is computed using the standard deviation of 
$D=|v|+|w|$; see Marxen & Zaki (Reference Marxen and Zaki2019) and Yao et al. (Reference Yao, Alves-Portela and Papadakis2020) for more details.
Since the KHMH equation involves two points, the conditions for the averaging operation should be defined using the states at both points. Four combinations are possible:
(i) If both
$x_i^+$ and $x_i^-$ are located within a turbulent patch, this is denoted as a turbulent–turbulent (or $TT$) event, and it is defined by the condition $\tau ^+\tau ^-=1$, where $\tau ^+=\tau (x_i^+,t)$ and $\tau ^-=\tau (x_i^-,t)$.(ii) If both points are within the laminar region, this is a laminar–laminar (or
$LL$) event, and it is defined by the condition $(1-\tau ^+)(1-\tau ^-)=1$.(iii) If
$x_i^+$, $x_i^-$ are within a turbulent and a laminar patch respectively, this is a turbulent–laminar (or $TL$) event, defined by $\tau ^+(1-\tau ^-)=1$.(iv) If
$x_i^+$, $x_i^-$ are within a laminar and turbulent region respectively, this is a laminar–turbulent (or $LT$) event, defined by $(1-\tau ^+)\tau ^-=1$.
In our notation, the first upper-case letter denotes the state of point 
$x_i^+$ and the second the state of 
$x_i^-$. We also use the generic notation 
$AA$ to refer to a general event, i.e. 
$AA=TT$ or 
$LL$ or 
$TL$ or 
$LT$. The four events are shown schematically in figure 3.
Figure 3. Contour plot of wall-normal velocity fluctuations in the transitional region. A turbulent spot is clearly visualised on the right half of the figure. Within the spot 
$\tau =1$, while outside (grey region) 
$\tau =0$. The four different two-point event types, 
$LL$, 
$TT$, 
$TL$ and 
$LT$, are shown. All points are located at the same horizontal plane, i.e. 
$r_2=0$.
We can now define the time- and spanwise-average two-point intermittencies as
$$\begin{gather} \gamma^{(TT)}{(X,Y;r_1,r_2,r_3)} \equiv \frac{1}{{\rm \Delta} T L_z} \int_0^{{\rm \Delta} T} \int_0^{L_z} \tau^+\tau^{-} \, {\rm d} z \, {\rm d} t, \end{gather}$$
$$\begin{gather}\gamma^{(TL)}{(X,Y;r_1,r_2,r_3)} \equiv \frac{1}{{\rm \Delta} T L_z } \int_0^{{\rm \Delta} T} \int_0^{L_z} \tau^+(1-\tau^{-}) \, {\rm d} z \, {\rm d} t, \end{gather}$$
$$\begin{gather}\gamma^{(LT)}{(X,Y;r_1,r_2,r_3)} \equiv \frac{1}{{\rm \Delta} T L_z } \int_0^{{\rm \Delta} T} \int_0^{L_z} (1-\tau^+)\tau^{-} \, {\rm d} z \, {\rm d} t, \end{gather}$$
$$\begin{gather}\gamma^{(LL)} {(X,Y;r_1,r_2,r_3)} \equiv \frac{1}{{\rm \Delta} T L_z } \int_0^{{\rm \Delta} T} \int_0^{L_z} (1-\tau^+)(1-\tau^{-}) \, {\rm d} z \, {\rm d} t. \end{gather}$$
Since 
$\tau ^+\tau ^-+(1-\tau ^+)(1-\tau ^-)+\tau ^+(1-\tau ^-)+(1-\tau ^+)\tau ^-=1$, we have 
$\gamma ^{(TT)}+\gamma ^{(TL)}+\gamma ^{(LT)}+\gamma ^{(LL)}=1$. Two point intermittencies were also defined in Yao et al. (Reference Yao, Alves-Portela and Papadakis2020), where 
$TL$ and 
$TL$ events were amalgamated as a combined 
$TL$ event. Here we consider the two event types separately for reasons that will become clear shortly.
The conditional time averages of the general two-point increment 
${\rm d} Q=Q(x_i^+)-Q(x_i^-)$ are defined as
$$\begin{gather} \, \overline{{\rm d} Q}^{(TT)}(X,Y;r_1,{r_2},r_3)\equiv \frac{1}{{\rm \Delta} T L_z \gamma^{(TT)}} \int_0^{{\rm \Delta} T} \int_0^{L_z} \tau^+\tau^{-} \, {\rm d} Q \, {\rm d} z \, {\rm d} t, \end{gather}$$
$$\begin{gather}\overline{{\rm d} Q}^{(TL)}(X,Y;r_1,{r_2},r_3)\equiv \frac{1}{{\rm \Delta} T L_z \gamma^{(TL)}} \int_0^{{\rm \Delta} T} \int_0^{L_z} \tau^+(1-\tau^{-}) \, {\rm d} Q \, {\rm d} z \, {\rm d} t, \end{gather}$$
$$\begin{gather}\overline{{\rm d} Q}^{(LT)}(X,Y;r_1,{r_2},r_3)\equiv \frac{1}{{\rm \Delta} T L_z \gamma^{(LT)}} \int_0^{{\rm \Delta} T} \int_0^{L_z} (1-\tau^+)\tau^{-} \, {\rm d} Q \, {\rm d} z \, {\rm d} t, \end{gather}$$
$$\begin{gather}\overline{{\rm d} Q}^{(LL)}(X,Y;r_1,{r_2},r_3)\equiv \frac{1}{{\rm \Delta} T L_z \gamma^{(LL)}} \int_0^{{\rm \Delta} T} \int_0^{L_z} (1-\tau^+)(1-\tau^{-}) \, {\rm d} Q \, {\rm d} z \, {\rm d} t. \end{gather}$$This means that the standard time average can be decomposed as
\begin{align} \overline{{\rm d} Q}(X,Y;r_1,r_2,r_3) & = \frac{1}{{\rm \Delta} T L_z} \int_0^{{\rm \Delta} T} \int_0^{L_z} \, {\rm d} Q \, {\rm d} z \, {\rm d} t \nonumber\\ & = \gamma^{(TT)}\, \overline{{\rm d} Q}^{(TT)}+\gamma^{(TL)}\, \overline{{\rm d} Q}^{(TL)}+ \gamma^{(LT)}\, \overline{{\rm d} Q}^{(LT)}+\gamma^{(LL)}\, \overline{{\rm d} Q}^{(LL)}. \end{align} The definitions are similar for the conditionally averaged midpoint variable 
$Q^*=0.5[Q(x_i^+)+Q(x_i^-)]$. Thus,
\begin{align} \overline{Q^*}(X,Y;r_1,r_2,r_3) & = \frac{1}{{\rm \Delta} T L_z} \int_0^{{\rm \Delta} T} \int_0^{L_z} Q^* \, {\rm d} z \, {\rm d} t \nonumber\\ & = \gamma^{(TT)} \overline{Q^*}^{(TT)}+\gamma^{(TL)}\overline{Q^*}^{(TL)}+ \gamma^{(LT)} \overline{Q^*}^{(LT)}+\gamma^{(LL)} \overline{Q^*}^{(LL)}. \end{align} Referring back to figure 3, it is clear that 
$\overline {{\rm d} Q}^{(TL)}$ and 
$\overline {{\rm d} Q}^{(LT)}$ are two-point averages across the laminar/turbulent interface. For 
$r_1=0$ and 
$r_3\ne 0$ (case shown in figure 3), these two conditional averages are taken across the interface in the spanwise direction. For 
$r_1>0$ and 
$r_3=0$ (case shown in figure 4), 
$LT$ averaging is taken across the tail (i.e. the downstream end) of a turbulent spot, while 
$TL$ is taken across the head of the spot (i.e. the upstream end). It is therefore possible to distinguish the different properties of the head or tail of a spot using the appropriate conditionally averaged quantity. This is an important observation and facilitates the physical interpretation of the results presented in § 6. Note that 
$\overline {{\rm d} Q}^{(TL)}(r_1,r_2,r_3)=\overline {{\rm d} Q}^{(LT)}(-r_1,-r_2,-r_3)$, but since we consider only positive separations, we need to retain both 
$\overline {{\rm d} Q}^{(TL)}$ and 
$\overline {{\rm d} Q}^{(LT)}$.
Figure 4. Contour plot of instantaneous streamwise velocity fluctuations. Yellow ovals demarcate two turbulent spots. The purple vertical dotted line represents a fixed streamwise location, 
$X$. Two event types, 
$TL$ and 
$LT$, are shown with 
$r_1 \neq 0$ and 
$r_3 = 0$.
We can now proceed to derive the conditionally averaged KHMH equation.
4.2. Derivation of the conditionally averaged KHMH equation
We start again with the Navier–Stokes equations at two points 
$x_i^+$ and 
$x_i^-$, (3.1), and define the conditional velocity fluctuation difference as
\begin{equation} {\rm d} u_i'^{(AA)} \equiv {\rm d} u_i-{\rm d} U_i^{(AA)},\end{equation}
where 
$AA =TT$ or 
$LL$ or 
$TL$ or 
$LT$ as mentioned earlier, 
${\rm d} u_i=u_i^+-u_i^-$ and 
${\rm d} U_i^{(AA)}=\overline {{\rm d} u_i}^{(AA)}$ (from definition (4.2)). It is straightforward to prove that 
$\overline {{\rm d} u_i'^{(AA)}}^{(AA)}=0$; this is the equivalent of 
$\overline {{\rm d} u_i'}=0$ in standard averaging.
Similarly, the conditional fluctuation velocity at the midpoint is defined as
\begin{equation} {u'^*_i}^{(AA)} \equiv {u^*_i}-{U^*_i}^{(AA)},\end{equation}
where 
$u^*_{i}=0.5(u_i^++u_i^-)$ and 
${U^*_i}^{(AA)}=\overline {u^*_{i}}^{(AA)}$, and again 
$\overline {{u'^*_i}^{(AA)}}^{(AA)}=0$.
Subtracting equation (3.1b) from (3.1a), multiplying each term by 
$2\, {\rm d} u_i'^{(AA)}$ and then applying the 
$(AA)$ averaging operation as defined in (4.2), we obtain
\begin{align} &\underbrace{ \overline{2\, {\rm d} u_i'^{(AA)}\frac{\partial \, {\rm d} u_i}{\partial t}}^{(AA)}}_{\text{Transient term}} +\underbrace{ \overline{2\, {\rm d} u_i'^{(AA)}u^+_j \frac{\partial \, {\rm d} u_i}{\partial x^+_j} + 2\, {\rm d} u_i'^{(AA)}u^-_j \frac{\partial \, {\rm d} u_i}{\partial x^-_j}}^{(AA)}}_{\text{Nonlinear term}} \nonumber\\ &\quad =\underbrace{ \overline{-2\, {\rm d} u_i'^{(AA)}\left(\frac{\partial \, {\rm d} p}{\partial x^+_i} - \frac{\partial \, {\rm d} p}{\partial x^-_i}\right)}^{(AA)} }_{\text{Pressure--velocity correlation}} + \underbrace{ \overline{2 \nu \, {\rm d} u_i'^{(AA)}\frac{\partial^2 \, {\rm d} u_i}{\partial x_j^+ \partial x_j^+} + 2\nu \, {\rm d} u_i'^{(AA)}\frac{\partial^2 \, {\rm d} u_i}{\partial x_j^- \partial x_j^-}}^{(AA)} }_{\text{Viscous term}}, \end{align}
where we have used again 
${\partial u^+_i}/{\partial x^-_j}=0$ and 
${\partial u^-_i}/{\partial x^+_j}=0$. This is the equivalent of (3.5) for standard averaging.
We now define the two-point conditional energy as 
$\overline {{\rm d} q^{2^{(AA)}}}^{(AA)} =\overline {{\rm d} u_i'^{(AA)} \, {\rm d} u_i'^{(AA)}}^{(AA)}$ and seek its transport equation in the physical and scale spaces, similar to (3.6). Applying again the variable transformation 
$X_i=0.5(x^+_i + x^-_i)$ and 
$r_i = x^+_i - x^-_i$ and the definitions (4.5) and (4.6) into (4.7), after some algebra we get the following conditionally averaged KHMH equation for 
$\overline {{\rm d} q^{2^{(AA)}}}^{(AA)}$:
\begin{align}
& \underbrace{\overline{\frac{\partial \, {\rm d}
q^{2^{(AA)}}}{\partial t} }^{(AA)}}_{\text{Transient
term}}+ \underbrace{\overline{ {U^*_j}^{(AA)}\frac{\partial
q^{2^{(AA)}}}{\partial X_j} }^{(AA)}}_{\text{Mean
advection}}+ \underbrace{\overline{
{u'^*_j}^{(AA)}\frac{\partial \, {\rm d}
q^{2^{(AA)}}}{\partial X_j} }^{(AA)}}_{\text{Turbulent
advection}}+ \underbrace{\overline{{\rm d} U_j^{(AA)}
\frac{\partial \, {\rm d} q^{2^{(AA)}}}{\partial r_j}
}^{(AA)}}_{\text{Linear transfer}}+\underbrace{\overline{{\rm d} u_j'^{(AA)} \frac{\partial
\, {\rm d} q^{2^{(AA)}}}{\partial r_j}}^{(AA)}
}_{\text{Nonlinear transfer}}\nonumber\\
&\quad =\underbrace{\overline{-2\,
{\rm d} u_i'^{(AA)} \frac{\partial \, {\rm d}
p'^{(AA)}}{\partial X_i}
}^{(AA)}}_{\text{Pressure--velocity
correlation}} +\underbrace{\nu \frac{1}{2} \overline{
\frac{\partial^2 \, {\rm d} q^{2^{(AA)}}}{\partial X_j
\partial X_j} }^{(AA)}}_{\text{Physical diffusion}}+
\underbrace{2\nu \overline{ \frac{\partial^2 \, {\rm d}
q^{2^{(AA)}} }{\partial r_j \partial
r_j}}^{(AA)}}_{\text{Scale diffusion}}\nonumber\\
&\qquad -\underbrace{\overline{2\, {\rm d} u_i'^{(AA)}
{u'^*_j}^{(AA)}\frac{\partial \, {\rm d}
U_i^{(AA)}}{\partial X_j}}^{(AA)} - \overline{2\, {\rm d}
u_i'^{(AA)} \, {\rm d} u_j'^{(AA)} \frac{\partial \, {\rm
d} U_i^{(AA)}}{\partial r_j}}^{(AA)}
}_{\text{Production}}\nonumber\\
& \qquad - \underbrace{2\nu \overline{ \left( \frac{\partial \, {\rm d} u'^{(AA)}_i
}{\partial x^+_j}\frac{\partial \, {\rm d} u'^{(AA)}_i
}{\partial x^+_j} + \frac{\partial \, {\rm d} u'^{(AA)}_i
}{\partial x^-_j}\frac{\partial \, {\rm d} u'^{(AA)}_i
}{\partial x^-_j} \right ) }^{(AA)}}_{\text{Dissipation}}.
\end{align} The form of (4.8) is similar to the standard form (3.6). This is due to the appropriate definitions of the conditional fluctuating quantities in (4.5) and (4.6) that satisfy 
$\overline {{\rm d} u_i'^{(AA)}}^{(AA)}=0$ and 
$\overline {{u'^*_i}^{(AA)}}^{(AA)}=0$. There is, however, an important difference. Standard averaging commutes with the spatial differentiation operation, for example,
\begin{equation} \frac{\partial \, \overline{{\rm d} Q}}{\partial r_j}= \frac{1}{{\rm \Delta} T L_z} \frac{\partial}{\partial r_j} \left( \int_0^{{\rm \Delta} T} \int_0^{L_z} \, {\rm d} Q \, {\rm d} z \, {\rm d} t \right)= \frac{1}{{\rm \Delta} T L_z} \int_0^{{\rm \Delta} T} \int_0^{L_z} \frac{\partial \, {\rm d} Q}{\partial r_j} \, {\rm d} z \, {\rm d} t =\overline{\frac{\partial \,{{\rm d} Q}}{\partial r_j}}. \end{equation}However, this is not the case for conditionally averaged quantities, for example
\begin{equation} \frac{\partial \, \overline{{\rm d} Q}^{(TT)}}{\partial r_j}= \frac{1}{{\rm \Delta} T L_z} \frac{\partial}{\partial r_j} \left( \frac{1}{\gamma^{(TT)}} \int_0^{{\rm \Delta} T} \int_0^{L_z} \tau^+\tau^{-} \, {\rm d} Q \, {\rm d} z \, {\rm d} t \right), \end{equation}while
\begin{equation} \overline{\frac{\partial \,{{\rm d} Q}}{\partial r_j}}^{(TT)}= \frac{1}{{\rm \Delta} T L_z} \left( \frac{1}{\gamma^{(TT)}} \int_0^{{\rm \Delta} T} \int_0^{L_z} \tau^+\tau^{-} \frac{\partial \,{{\rm d} Q}}{\partial r_j} \, {\rm d} z \, {\rm d} t \right) \ne \frac{\partial \, \overline{{\rm d} Q}^{(TT)}}{\partial r_j}, \end{equation}
because 
$\gamma ^{(TT)}$ depends on 
$r_j$. The same issue appears in the conditionally averaged turbulent kinetic energy equation (Marxen & Zaki Reference Marxen and Zaki2019).
This lack of commutation has two important implications. First, (4.8) cannot be written in conservative form. For example, for the nonlinear interscale energy transfer term, we have
\begin{equation} \overline{{\rm d} u_j'^{(AA)} \frac{\partial \, {\rm d} q^{2^{(AA)}}} {\partial r_j}}^{(AA)} = \overline{ \frac{\partial \, {\rm d} q^{2^{(AA)}} \, {\rm d} u_j'^{(AA)}}{\partial r_j}}^{(AA)} \neq \frac{\partial (\overline{{\rm d} q^{2^{(AA)}} \, {\rm d} u_j'^{(AA)}}^{(AA)} )}{\partial r_j}. \end{equation}
The first equality in the above equation is because 
${\partial {\, {\rm d} u_j'^{(AA)}}}/{\partial r_j}=0$ (this is equivalent to 
${\partial \, \overline {{\rm d} u_j'}}/{\partial r_j}=0$ for standard averaging). It is easy to prove; for example if 
$AA=TT$, we have
\begin{align}
\frac{\partial \, {\rm d} u_j'^{(TT)}}{\partial r_j}&=
\frac{\partial \, {\rm d} u_j}{\partial r_j}-\frac{\partial
\, {\rm d} U_j^{(TT)}}{\partial r_j}=\frac{\partial \, {\rm
d} u_j}{\partial r_j}- \frac{\partial}{\partial r_j} \left(
\frac{1}{{{\rm \Delta} T L_z} \gamma^{(TT)}} \int_0^{{\rm \Delta} T}
\int_0^{L_z} \tau^+\tau^{-} \,{\rm d} U_j \, {\rm d} z \, {\rm d} t \right)\nonumber\\
& =\frac{\partial \, {\rm d} u_j}{\partial r_j}- \frac{\partial}{\partial r_j} \left(
\frac{{\rm d} U_j }{{\rm \Delta} T L_z \gamma^{(TT)}}
\int_0^{{\rm \Delta} T} \int_0^{L_z} \tau^+\tau^{-} \, {\rm d} z
\, {\rm d} t \right)\nonumber\\
&= \frac{\partial \, {\rm d}
u_j}{\partial r_j}-\frac{\partial \, {\rm d} U_j}{\partial
r_j}=0-0=0,
\end{align}
because 
${\rm d} u_j$ and 
${\rm d} U_j$ satisfy the continuity equation in the scale space. We also took into account that 
${\rm d} U_j$ is constant with respect to 
$Z$ and 
$t$ and used (4.1a).
The second implication is computational. For standard averaging, all the terms in the KHMH equation can be evaluated numerically by differentiating locally at points 
$x_i^+$ and 
$x_i^-$, for example
\begin{align} \frac{\partial \, {\rm d} U_i }{\partial r_j}&= \frac{\partial \, \overline{{\rm d} u_i} }{\partial r_j} =\overline{\frac{\partial {\, {\rm d} u_i} }{\partial r_j}} =\frac{1}{2}\left( \overline{\frac{\partial {(u_i^+{-} u_i^-)} }{\partial x^+}}-\overline{\frac{\partial {(u_i^+{-} u_i^-)} }{\partial x_j^-}} \right)=\frac{1}{2} \left( \overline{\frac{\partial {u_i^+} }{\partial x_j^+}}+ \overline{\frac{\partial {u_i^-} }{\partial x_j^-}}\right)\nonumber\\ &=\frac{1}{2} \left( \frac{\partial \overline{u_i^+} }{\partial x_j^+} +\frac{\partial \overline{u_i^-}}{\partial x_j^-} \right)= \frac{1}{2} \left( {\frac{\partial U_i^+ }{\partial x_j^+}} +{\frac{\partial U_i^- }{\partial x_j^-}}\right). \end{align}However, this is not the case for the conditionally averaged velocity difference:
\begin{align} \frac{\partial \, {\rm d} U_i^{(AA)} }{\partial r_j}&=\frac{\partial \, \overline{{\rm d} u_i}^{(AA)} }{\partial r_j} \ne \overline{\frac{\partial {\, {\rm d} u_i} }{\partial r_j}}^{(AA)} \left[ = \frac{1}{2} \left ( \overline{\frac{\partial {(u_i^+{-} u_i^-)}}{\partial x_j^+}}^{(AA)} -\overline{\frac{\partial {(u_i^+{-} u_i^-)}}{\partial x_j^-}}^{(AA)} \right ) \right. \nonumber\\ & =\left. \frac{1}{2} \left (\overline{ \frac{{\partial u_i^+ }}{\partial x_j^+}}^{(AA)} +\overline{ \frac{{\partial u_i^-}}{\partial x_j^-}}^{(AA)} \right ) \right]. \end{align}This means that the derivatives of the conditionally averaged two-point variables must be calculated directly in scale space. The process and validation against standard averaging are presented in Appendix A.
5. Conditional decomposition of nonlinear energy fluxes
In this section we decompose the nonlinear energy fluxes in scale and physical spaces, 
$\phi _{r_j}^F=\overline {{\rm d} q^2\, {\rm d} u'_j}=\overline {({\rm d} u'_i)^2\, {\rm d} u'_j}$ and 
${\phi _{s_j}^F}=\overline {u'^*_j \delta q^2}$ respectively, see (3.9) and (3.8), into conditionally averaged components, i.e. we seek to derive expressions similar to (4.3) and (4.4). The decomposition for the second-order structure function (quadratic quantity) was derived in Yao et al. (Reference Yao, Alves-Portela and Papadakis2020). Here we extend the method to energy fluxes (cubic quantities).
For the nonlinear energy flux in scale space 
$\phi _{r_j}^F=\overline {({\rm d} u'_i)^2\, {\rm d} u'_j}$, we have
\begin{align} \phi_{r_j}^F &=\overline{({\rm d} u'_i)^2\, {\rm d} u'_j}= \overline{({\rm d} u_i-{\rm d} U_i)^2({\rm d} u_j-{\rm d} U_j)}\nonumber\\ &=\overline{ ({\rm d} u_i)^2\, {\rm d} u_j} -\overline{({\rm d} u_i)^2}\, {\rm d} U_j -\overline{2\, {\rm d} u_i\, {\rm d} u_j}\, {\rm d} U_i +2\, {\rm d} U_i\, {\rm d} U_i\, {\rm d} U_j \nonumber\\ &=\gamma^{(TT)}\overline{({\rm d} u_i)^2\, {\rm d} u_j}^{(TT)} +\gamma^{(TL)}\overline{({\rm d} u_i)^2\, {\rm d} u_j}^{(TL)} +\gamma^{(LT)}\overline{({\rm d} u_i)^2\, {\rm d} u_j}^{(LT)}\nonumber\\ &\quad +\gamma^{(LL)}\overline{({\rm d} u_i)^2\, {\rm d} u_j}^{(LL)} -\overline{({\rm d} u_i)^2}\, {\rm d} U_j -\overline{2\, {\rm d} u_i\, {\rm d} u_j}\, {\rm d} U_i +2\, {\rm d} U_i\, {\rm d} U_i\, {\rm d} U_j, \end{align}
where (4.3) was applied to 
$\overline {({\rm d} u_i)^2\, {\rm d} u_j}$.
We now write 
$\overline {({\rm d} u_i)^2\, {\rm d} u_j}^{(AA)}$ (
$AA =TT$ or 
$LL$ or 
$TL$ or 
$LT$) in terms of conditional fluctuations. To do this, we use the definition 
${\rm d} u_i'^{(AA)} \equiv {\rm d} u_i-{\rm d} U_i^{(AA)}$, see (4.5), and express the conditional nonlinear energy flux 
${\phi ^F_{r_j}}^{(AA)}=\overline {({\rm d} u'^{(AA)}_i)^2\, {\rm d} u'^{(AA)}_j}^{(AA)}$ as
\begin{align} {\phi^F_{r_j}}^{(AA)}&=\overline{({\rm d} u'^{(AA)}_i)^2\, {\rm d} u'^{(AA)}_j}^{(AA)} =\overline{ ({\rm d} u_i)^2\, {\rm d} u_j}^{(AA)} -\overline{({\rm d} u_i)^2}^{(AA)}\, {\rm d} U_j^{(AA)} \nonumber\\ &\quad -\overline{2\, {\rm d} u_i\, {\rm d} u_j}^{(AA)}\, {\rm d} U_i^{(AA)} -2\, {\rm d} U_i^{(AA)}\, {\rm d} U_i^{(AA)}\, {\rm d} U_j^{(AA)}. \end{align}
Solving for 
$\overline {({\rm d} u_i)^2\, {\rm d} u_j}^{(AA)}$ and substituting into (5.1), we obtain the desired decomposition:
\begin{equation} \phi^F_{r_j} =\gamma^{(TT)}{\phi^F_{r_j}}^{(TT)}+ \gamma^{(TL)}{\phi^F_{r_j}}^{(TL)}+ \gamma^{(LT)}{\phi^F_{r_j}}^{(LT)}+ \gamma^{(LL)}{\phi^F_{r_j}}^{(LL)}+\phi_{r_j},\end{equation}
where the additional term 
$\phi _{r_j}$ is given by
\begin{align} \phi_{rj} &= \sum_{AA,BB=}^{TT, TL, LT, LL} \gamma^{(AA)}\gamma^{(BB)}\, \overline{{\rm d} u_i\, {\rm d} u_i}^{(AA)}({\rm d} U_j^{(AA)}-{\rm d} U_j^{(BB)})\nonumber\\ &\quad +\sum_{AA,BB=}^{TT, TL, LT, LL} 2\gamma^{(AA)}\gamma^{(BB)}\, \overline{{\rm d} u_i\, {\rm d} u_j}^{(AA)}({\rm d} U_i^{(AA)}-{\rm d} U_i^{(BB)})\nonumber\\ &\quad - \sum_{AA,BB,CC=}^{TT, TL, LT, LL}2 \gamma^{(AA)}\gamma^{(BB)}\gamma^{(CC)} \,{\rm d} U_i^{(AA)}({\rm d} U_i^{(AA)}\, {\rm d} U_j^{(AA)}-{\rm d} U_i^{(BB)}\, {\rm d} U_j^{(CC)}). \end{align} Following the same process, the following decomposition can be obtained for the nonlinear energy flux in physical space 
${\phi _{s_j}^F}=\overline {u'^*_j \delta q^2}=\overline {u'^*_j ({\rm d} u'_i)^2}$:
\begin{equation} \phi^F_{s_j} =\gamma^{(TT)}{\phi^F_{s_j}}^{(TT)}+ \gamma^{(TL)}{\phi^F_{s_j}}^{(TL)}+ \gamma^{(LT)}{\phi^F_{s_j}}^{(LT)}+ \gamma^{(LL)}{\phi^F_{s_j}}^{(LL)}+\phi_{s_j},\end{equation}
where 
${\phi ^F_{s_j}}^{(AA)}=\overline {{u'^*_j}^{(AA)}({\rm d} u'^{(AA)}_i)^2}^{(AA)}$ and the additional term 
$\phi _{s_j}$ is given by
\begin{align} \phi_{sj} &=\sum_{AA,BB=}^{TT, TL, LT, LL} \gamma^{(AA)}\gamma^{(BB)}\, \overline{{\rm d} u_i\, {\rm d} u_i}^{(AA)}({U^*_j}^{(AA)}-{U^*_j}^{(BB)})\nonumber\\ &\quad +\sum_{AA,BB=}^{TT, TL, LT, LL} 2\gamma^{(AA)}\gamma^{(BB)}\, \overline{{\rm d} u_iu^*_j}^{(AA)}({\rm d} U_i^{(AA)}-{\rm d} U_i^{(BB)})\nonumber\\ &\quad - \sum_{AA,BB,CC=}^{TT, TL, LT, LL}2 \gamma^{(AA)}\gamma^{(BB)}\gamma^{(CC)} \,{\rm d} U_i^{(AA)}({\rm d} U_i^{(AA)}{U^*_j}^{(AA)}-{\rm d} U_i^{(BB)}{U^*_j}^{(CC)}). \end{align}It is clear that the additional terms arise due to differences in the conditionally averaged quantities. Similar terms appear in single-point statistics; see for example (2.9) in Marxen & Zaki (Reference Marxen and Zaki2019) or (5.306) in Pope (Reference Pope2000).
Figure 5 presents contours of all the terms appearing in (5.3) for the decomposition of the energy flux vector 
$(\phi ^F_{r_3},\phi ^F_{s_2})$ at location 
$TR2$ on the 
$(r_3, Y)$ plane (the separations in the other two directions are equal to 0, i.e. 
$r_1=0, r_2=0$). In the figure, we combine the results of 
$TL$ and 
$LT$ together. Figure 5(a), denoted as LHS, depicts the results from the direct calculation of the flux vector 
$(\phi ^F_{r_3},\phi ^F_{s_2})$ using standard time average, i.e. 
$(\overline {({\rm d} u'_i)^2\, {\rm d} u'_3}, \overline {({\rm d} u'_i)^2u'^*_{2}})$. Figure 5( f), denoted as RHS, presents the sum of all terms in (5.3). It is clear that the results in the two panels are essentially identical, confirming the validity of the derived decomposition. The flux vector originates from the focal point 
$Y/L_0 \approx 2.5$ and 
$r_3/L_0 \approx 10$ and transfers energy radially to different directions. The strongest flux is located in the region 
$0< Y/L_0<10, r_3/L_0>10$ and is found to be positive, indicating strong inverse cascade. This pattern was observed in Yao et al. (Reference Yao, Mollicone and Papadakis2022).
Figure 5. Nonlinear interscale energy flux vector 
$(\phi ^F_{r_3},\phi ^F_{s_2})$ (denoted LHS) and the constituent components 
$(\gamma ^{(AA)}{\phi ^F_{r_3}}^{(AA)}, \gamma ^{(AA)}{\phi ^F_{s_2}}^{(AA)})$, where 
$AA=TT, TL+LT, LL$, in the 
$(r_3,Y)$ plane at location 
$TR2$. Here, RHS denotes the sum of the terms on the right-hand side of (5.3). Contours represent the magnitude of the flux vectors and the purple horizontal lines indicate the local boundary layer thickness.
Apart from verifying the decomposition (5.3), figure 5 provides important insight into the origin of the aforementioned energy flux pattern. It is interesting to observe that the inverse cascade arises from the 
$TL+LT$, 
$LL$ and 
$\phi$ components; all contribute to the strong positive flux in the region 
$0< Y/L_0<10$, 
$r_3/L_0>10$ with values that are of the same order of magnitude. Outside this region, they have small values. On the other hand, the 
$TT$ component contributes to the forward cascade in the region left of the focal point. The shape therefore of the energy flux vectors on the 
$(r_3, Y)$ plane arises from the superposition of the contribution of 
$TT$ term which is responsible for forward cascade, and all the other terms which are responsible for inverse cascade.
The constituent terms plotted in figure 5 are weighted by the corresponding two-point intermittencies, 
$\gamma ^{(AA)}$, which can have very small values depending on the location examined (see figure 1). In the following section, the evolution of the vectors 
$( {{\phi ^F_{r_3}}^{(AA)}},{{\phi ^F_{s_2}}^{(AA)}} )=(\overline {({\rm d} u'^{(AA)}_i)^2\, {\rm d} u'^{(AA)}_3}^{(AA)}, \overline {({\rm d} u'^{(AA)}_i)^2{u'^*_2}^{(AA)}}^{(AA)})$ at different streamwise locations is examined.
6. Evolution of conditionally averaged nonlinear energy fluxes
In this section we consider the evolution of the nonlinear energy fluxes in the transitional region. To simplify the analysis and interpretation of the results we consider 
$r_2=0$, i.e. points 
$x_i^+$ and 
$x_i^-$ have the same wall-normal distance.
6.1. Fluxes on $(r_3,Y)$ plane
Figure 6 shows the conditionally averaged energy fluxes 
$( {{\phi ^F_{r_3}}^{(AA)}},{{\phi ^F_{s_2}}^{(AA)}} )$, where 
$AA = TT$ or 
$TL+LT$ or 
$LL$ at locations 
$TR1,TR2,TR3$. The last row (marked as ‘Total’) depicts 
$( {{\phi ^F_{r_3}}},{{\phi ^F_{s_2}}} )$. In all plots, 
$r_1=0$ and 
$r_2=0$.
Figure 6. Conditionally averaged flux vectors 
$({{\phi ^F_{r_3}}^{(AA)}},{{\phi ^F_{s_2}}^{(AA)}} )$ (a–i) and standard-averaged flux vector 
$( {{\phi ^F_{r_3}}},{{\phi ^F_{s_2}}} )$ (j–l) at locations 
$TR1$ (a,d,g,j), 
$TR2$ (b,e,h,k) and 
$TR3$ (c, f,i,l). The contours represent the magnitude of the corresponding flux vectors and the purple horizontal lines indicate the local boundary layer thickness.
At the early stages of transition, at location 
$TR1$ (left column), the flux vector 
$( {{\phi ^F_{r_3}}^{(LL)}},{{\phi ^F_{s_2}}^{(LL)}} )$ is dominant and is almost identical to the total flux. We can clearly see strong inverse cascade occurring to the right of the focal point, which is mainly due to 
$LL$ events. The vectors 
$( {{\phi ^F_{r_3}}^{(TT)}}, {{\phi ^F_{s_2}}^{(TT)}} )$ and 
$( {{\phi ^F_{r_3}}^{(TL)}}, {{\phi ^F_{s_2}}^{(TL)}} )+( {{\phi ^F_{r_3}}^{(LT)}},{{\phi ^F_{s_2}}^{(LT)}} )$ are very localised and their behaviour is difficult to interpret. This may be due to the small number of 
$TT$ and 
$TL/LT$ events at this early transition location. Notice, however, the high values of the magnitude of the 
$TT$ vector, 
$({{\phi ^F_{r_3}}^{(TT)}}, {{\phi ^F_{s_2}}^{(TT)}} )$, compared with all the other vectors. However, its overall contribution to the total flux, 
$\gamma ^{(TT)}\boldsymbol {\phi} ^{{F}^{(TT)}}$, is small because the values of intermittency 
$\gamma ^{(TT)}$ are negligible at this location (see figure 1).
At location 
$TR2$ (middle column), 
$TT$ and 
$TL/LT$ events start to play a more significant role as expected, but they act in different areas of the map. Events of the 
$LL$ type maintain their dominant contribution to the inverse cascade, while 
$TL$ and 
$LT$ events also amplify and contribute to the inverse cascade in the area 
$10L_0< r_3<20L_0$ above 
$Y\approx 5L_0$. Notice again that 
$TT$-averaged flux vectors have a much higher magnitude compared with the total flux, and show a mixture of the forward and inverse cascade to the left and right of the focal point respectively (with the forward cascade being slightly stronger). This picture is consistent with figure 5 where the conditionally averaged fluxes are weighted by the two-point intermittency.
At location 
$TR3$ (right column), 
$TT$ events now assume the dominant role because 
$\gamma ^{(TT)}$ approaches 1, 
$LL$ events are localised (in the same way that 
$TT$ events were localised at 
$TR1$ location), while 
$TL$ and 
$LT$ events again give rise to inverse cascade. Interestingly, as intermittency increases, 
$TT$ events show stronger inverse cascade. For comparison, the flux vector 
$( \phi ^F_{r_3}, \phi ^F_{s_2} )$ at the fully turbulent location 
$TU$ is plotted in figure 7. At this location, the inverse cascade is weakened (but it is still visible) and energy flows in the wall-normal direction before looping back to small spanwise length scales.
Figure 7. Standard-averaged flux vector 
$( {{\phi ^F_{r_3}}},{{\phi ^F_{s_2}}} )$ at 
$TU$ location. The contour represents the magnitude of the flux vectors.
The strong inverse cascade found at 
$TR2$ and 
$TR3$ locations is clearly not observed in the fully turbulent region. The origin of the inverse cascade arises mainly due to 
$TL$ and 
$LT$ events. Indeed, in the 
$TR2$ location their magnitude is 5–6 times larger compared with 
$LL$ events (that also contribute to inverse cascade), while in 
$TR3$ it is about 2–3 times larger. The laminar/turbulent interface, therefore, plays a crucial role in the inverse cascade process. In the next section, we focus on the 
$TR2$ location and explore in more depth the cascade process in the three-dimensional 
$(r_1,r_3, Y)$ hyperplane.
6.2. Fluxes on the $(r_1,r_3,Y)$ hyperplane and $(r_1,r_3)$ plane
We project the conditionally averaged fluxes 
${\phi ^F}^{(AA)}$ and 
${\phi ^S}^{(AA)}$ on the 
$(r_1,r_3,Y)$ hyperplane, i.e. plot three-dimensional maps of the vector 
$({\phi ^F_{r_1}}^{(AA)}, {\phi ^F_{r_3}}^{(AA)},{\phi ^F_{s_2}}^{(AA)})$. We select location 
$TR2$, which is in the middle of the transition region, and has single-point intermittency of about 0.4 (see figure 1).
Stream-tubes obtained from the total and the nonlinear flux vectors (standard- or conditionally averaged) are shown in figure 8; see caption for details. Figure 8(a) shows the standard-averaged total energy flux vector 
$(\phi _{r_1},\phi _{r_3},\phi _{s_2})$ that includes the nonlinear, linear, pressure and viscous components (see (3.9)). This plot is very similar to that of figure 8(b) that shows the nonlinear component 
$(\phi ^F_{r_1}, \phi ^F_{r_3}, \phi ^F_{s_2})$, the latter being the dominant component (see Yao et al. (Reference Yao, Mollicone and Papadakis2022) for a detailed discussion on the other components). Both plots depict a dense cluster of stream-tubes with energy flowing to larger 
$r_1$ scales (inverse cascade) with energy originating at 
$Y \approx 5L_0$ and 
$r_3\approx 10L_0$. There is milder inverse cascade in larger 
$r_3$ scales. In another (smaller) cluster, energy flux vectors rotate and bend towards the 
$Y$ axis (
$r_1=0,r_3=0$).
Figure 8. Stream-tubes of (a) the total flux vector 
$({{\phi _{r_1}}},{{\phi _{r_3}}},{{\phi _{s_2}}})$, (b) the standard-averaged nonlinear flux vector 
$({\phi ^F_{r_1}}, {\phi ^F_{r_3}}, {\phi ^F_{s_2}})$, and the conditionally averaged vectors (c) 
$({{\phi ^F_{r_1}}^{(LT)}}, {{\phi ^F_{r_3}}^{(LT)}},{{\phi ^F_{s_2}}^{(LT)}})$, (d) 
$({{\phi ^F_{r_1}}^{(TL)}}, {{\phi ^F_{r_3}}^{(TL)}},{{\phi ^F_{s_2}}^{(TL)}})$, (e) 
$({{\phi ^F_{r_1}}^{(TT)}}, {{\phi ^F_{r_3}}^{(TT)}},{{\phi ^F_{s_2}}^{(TT)}})$ and ( f) 
$({{\phi ^F_{r_1}}^{(LL)}}, {{\phi ^F_{r_3}}^{(LL)}},{{\phi ^F_{s_2}}^{(LL)}})$ at 
$TR2$. The plots are generated by placing a sphere of radius 
$5L_0$ at point 
$(r_1,r_3, Y)=(5L_0, 10L_0, 3L_0)$ and tracing the stream-tubes crossing the sphere. The stream-tubes are coloured according to the sign of the first component, i.e. 
${\phi _{r_1}}$, 
${\phi ^F_{r_1}}$, 
${\phi ^F_{r_1}}^{(LT)}$, 
${\phi ^F_{r_1}}^{(TL)}$, 
${\phi ^F_{r_1}}^{(TT)}$ and 
${\phi ^F_{r_1}}^{(LL)}$ (red for positive, blue for negative, thus indicating inverse or forward cascade in the 
$r_1$ direction respectively). The colour bars also refer to the value of the first component (the minimum/maximum values are the same to facilitate comparison).
The decomposition (5.3) allows us to probe in more detail the origin of the strong inverse cascade in the 
$r_1$ direction and identify the flow events that determine it. In the homogeneous spanwise direction 
$TL$ and 
$LT$ events can be combined together as 
$TL+LT$ (see figure 3). However, the streamwise direction is inhomogeneous and these events have to be considered separately. As mentioned earlier, for a fixed 
$(X,Y)$ spatial location and 
$r_3$ separation, a 
$TL$ event with 
$r_1>0$ takes place across the upstream laminar/turbulent interface of a spot. On the other hand, an 
$LT$ event is taken across the downstream interface (see figure 4).
Figures 8(c) and 8(d) show the stream-tubes for the conditionally averaged flux vectors 
$({\phi ^F_{r_1}}^{(LT)}, {\phi ^F_{r_3}}^{(LT)},{\phi ^F_{s_2}}^{(LT)})$ and 
$({\phi ^F_{r_1}}^{(TL)}, {\phi ^F_{r_3}}^{(TL)},{\phi ^F_{s_2}}^{(TL)})$ respectively. It is very clear that the 
$LT$ flux vector contributes most strongly to the inverse cascade; this corresponds to the downstream laminar/turbulent interface of a spot. Indeed there is a cluster of stream-tubes whose direction indicates the transfer of energy to larger streamwise scales, up to 
$r_1 \approx 50L_0$. The stream-tubes then bend towards smaller scales, 
$r_1,r_3 \to 0$. On the other hand, 
$TL$ events that correspond to energy flux across the upstream interface of a turbulent spot also contribute to inverse cascade, but over a shorter streamwise range, 
$0< r_1<20L_0$. This clearly indicates that interscale energy transfer processes are different at the upstream and downstream interfaces of turbulent spots.
Figures 8(e) and 8( f) show the stream-tubes of the conditionally averaged flux vectors 
$({\phi ^F_{r_1}}^{(TT)}, {\phi ^F_{r_3}}^{(TT)},{\phi ^F_{s_2}}^{(TT)})$ and 
$({\phi ^F_{r_1}}^{(LL)}, {\phi ^F_{r_3}}^{(LL)},{\phi ^F_{s_2}}^{(LL)})$. The former corresponds to interscale transfer within spots; there is weak inverse cascade and the stream-tubes bend towards smaller scales. The latter corresponds to time instants where both points are located within laminar regions. There is strong forward cascade to small scales and then bending and energy transfer to larger scales away from the wall. This is a quite complicated energy flux pattern, which is difficult to interpret physically.
Plots in the three-dimensional hyperplane 
$(r_1,r_3, Y)$ visualise the main features of the energy flux paths, but can hide important detail. To uncover this detail, in figure 9 we plot the flux vectors in the 
$(r_1,r_3)$ plane at the specific height 
$Y=4.5L_0$. The total flux vector and the standard-averaged nonlinear component (figures 9a and 9b respectively) are very similar and show a recirculating pattern with inverse cascade for 
$r_3>10$ over the range of 
$r_1$ examined and forward cascade for 
$r_3<10L_0$. It is very interesting to see that around 
$r_3=10L_0$ the energy flux is negligible; this cannot be easily observed from figure 8. The conditionally averaged fluxes also show detail that cannot be discerned from the three-dimensional plots. For example, the strong inverse cascade of 
${\phi ^F_{r_1}}^{(LT)}$ extends over the whole range of 
$r_3$ and 
$r_1$, while for 
${\phi ^F_{r_1}}^{(TL)}$ it extends only in a specific range of separations, 
$r_3 \approx (5-15)L_0$, depending on 
$r_1$, as can be seen from figures 9(c) and 9(d) respectively. The flux component 
${\phi ^F_{r_1}}^{(LL)}$ (figure 9e) clearly demonstrates forward cascade over the whole 
$r_3$ range examined, but 
${\phi ^F_{r_1}}^{(TT)}$ (figure 9f) is very small around 
$r_3\approx 10$ and increases at the boundaries of the domain. These plots confirm that the strong inverse cascade in the standard-averaged flux is due to 
$LT$ events at the downstream interface of the turbulent spots. Interestingly, the energy fluxes due to 
$TL$ and 
$LL$ events almost cancel out around 
$r_3\approx 10L_0$, and this explains the very small fluxes in this area for 
${\phi ^F_{r_1}}$. Notice also that the 
$TT$ events account for the forward cascade observed in 
${\phi ^F_{r_1}}$ and 
${\phi _{r_1}}$ for small 
$r_3$ separations.
Figure 9. Flux vectors (a) 
$(\phi _{r_1},\phi _{r_3})$, (b) 
$(\phi ^F_{r_1}, \phi ^F_{r_3})$, (c) 
$({\phi ^F_{r_1}}^{(LT)}, {\phi ^F_{r_3}}^{(LT)})$, (d) 
$({\phi ^F_{r_1}}^{(TL)}, {\phi ^F_{r_3}}^{(TL)})$, (e) 
$({\phi ^F_{r_1}}^{(TT)}, {\phi ^F_{r_3}}^{(TT)})$ and ( f) 
$({\phi ^F_{r_1}}^{(LL)}, {\phi ^F_{r_3}}^{(LL)})$ on the 
$(r_1,r_3)$ plane at wall-normal height 
$Y=4.5L_0$ and location 
$TR2$. Contours represent the magnitude of the energy flux vectors in the plane.
The above figures have demonstrated the central role of the downstream laminar/turbulent interface in the interscale transfer and in particular the inverse cascade over a large range of scales. Additionally, 
$TL$ events were localised in a smaller range of spanwise and streamwise separations. We now try to explain physically this behaviour with the aid of the cartoons shown in figure 10. More specifically, we consider a fixed 
$X$ location (for figures 8 and 9, 
$X=X_{TR2}$) and follow a turbulent spot of diamond shape as it propagates to the right and crosses this location. The arrowhead shape at the upstream and downstream ends in figure 10 is a simplified, but rather realistic, approximation. This can be seen from figure 11 where we demarcate the boundaries of two turbulent spots. It is also consistent with experimental spot observations (see figures 12 and 14 in Anthony, Jones & Lagraff (Reference Anthony, Jones and Lagraff2005)). The sharp corners of the spot around the maximum thickness are less realistic; the shape is more rounded in this region as can be seen from figure 18 of Marxen & Zaki (Reference Marxen and Zaki2019). However, analysis of this simplified shape can provide significant physical insight, as is seen next.
Figure 10. A propagating turbulent spot of diamond shape crossing a fixed 
$X$ location; the view shown is in the 
$(X, Z)$ plane. Red, yellow and blue circles represent the 
$x_i^-$, 
$X_i$ and 
$x_i^+$ points respectively. Only points with fixed 
$r_3={\rm \Delta} z$ separation are shown.
Figure 11. Contour plot of instantaneous streamwise velocity fluctuations. The boundaries of two turbulent spots in the transition region (marked by dashed yellow lines) indicate spots of approximately diamond shape.
A compilation of reported leading and trailing edge propagation speeds suggests typical values around 
$0.9 U_\infty$ and 
$0.5 U_\infty$ respectively (see Fransson Reference Fransson2010). The fact that the the leading edge propagates faster results in the spot growth in the streamwise direction and is consistent with the dominant role of the downstream laminar/turbulent interface in inverse cascade. This explains the direction of the nonlinear energy fluxes shown in figure 8, but does not explain the difference in the streamwise and spanwise separations over which TL and LT events are active.
The underlying cause for this difference can be elucidated with the aid of figure 10 that depicts snapshots of the propagating spot at three time instants 
$t_1, t_2$ and 
$t_3$. In all snapshots, we consider a fixed middle point (denoted with a yellow circle) located at the streamwise position, 
$X$. The blue and red circles represent the 
$x_i^+$ and 
$x_i^-$ points respectively. Only points with a fixed spanwise separation, 
$r_3={\rm \Delta} z$, are shown in the figure. At time 
$t=t_1$, the downstream apex of the spot lies exactly at the fixed 
$X$ location. It can be seen that for an 
$LT$ event (the only type of event possible at this time instant), the streamwise separation 
$r_1(t_1)$ is very long, of the order of the spot length. For 
${\rm \Delta} z=0$, 
$r_1(t_1)$ attains a maximum value, equal to twice the spot length.
At time instant 
$t=t_2$, approximately half of the spot has crossed 
$X$. The valid spanwise locations of the middle point are determined by the spreading angle of the front apex. Note that for fixed 
$r_3={\rm \Delta} z$, 
$TL$ and 
$LT$ events coexist, but it is clear that the 
$r_1(t_2)$ separation of an 
$LT$ event is shorter than the one at 
$t=t_1$, i.e. 
$r_1(t_2)< r_1(t_1)$. At 
$t=t_3$, the whole spot has crossed the considered location; thus the rear apex is at 
$X$. At this time instant, only 
$TL$ events are possible. It can be seen that only a narrow range of 
$r_1$ separations is admissible for a given 
$r_3$. The actual range depends on the spreading angle of the rear apex. Thus on average, the valid 
$r_1$ separations corresponding to 
$TL$ events at 
$t_2$ and 
$t_3$ are more narrow compared with 
$LT$ events at 
$t=t_1$ and 
$t_2$. This explains the inverse cascade over a wider range of separations for 
$LT$ events shown in figure 9(c). We also conjecture that the largest admissible 
$r_1$ value mentioned earlier explains why the stream-tubes shown in figure 8(c) reach up to a maximum 
$r_1$ (the exact value depends on point of origin of the stream-tubes in the 
$(r_3, Y)$ plane) and then bend backward towards small scales.
To provide further insight into the observed behaviour of the conditionally averaged fluxes, we examine in more detail the flux vector component in the 
$r_1$ direction 
${\phi ^F_{r_1}}^{(AA)}=\overline {({\rm d} u'^{(AA)}_i)^2\, {\rm d} u'^{(AA)}_1}^{(AA)}$. In figure 12, the probability density functions (PDFs) of 
$({\rm d} u'^{(AA)}_i)^2\, {\rm d} u'^{(AA)}_1$ for 
$AA= LT$, 
$TL$, 
$TT$ and 
$LL$ are plotted for separations 
$r_1=30L_0$ and 
$r_3=10L_0$ at plane 
$Y=4.5L_0$. The PDF of the instantaneous flux 
$({\rm d} u'^{(LT)}_i)^2\, {\rm d} u'^{(LT)}_1$ is asymmetric and skewed to positive values (implying more intense instantaneous inverse cascade events than forward cascade events). This is also the case but it is less evident for 
$({\rm d} u'^{(TL)}_i)^2\, {\rm d} u'^{(TL)}_1$, while 
$({\rm d} u'^{(TT)}_i)^2\, {\rm d} u'^{(TT)}_1$ is almost symmetric and 
$({\rm d} u'^{(LL)}_i)^2\, {\rm d} u'^{(LL)}_1$ is skewed to the left. Note the large positive and negative fluctuations of the instantaneous fluxes compared with the time-average values reported in figure 9. This means that instantaneously energy flows in either direction, and intense fluxes of relatively low probability tip the balance in one direction or another after time-averaging.
Figure 12. The PDFs of 
$({\rm d} u'^{(AA)}_i)^2\, {\rm d} u'^{(AA)}_1$ for 
$AA= LT$, 
$TL$, 
$TT$ and 
$LL$ for 
$r_1=30L_0$ and 
$r_3=10L_0$ at 
$Y=4.5L_0$.
7. Conditionally averaged scale energy production and transfer within a turbulent spot and comparison with fully developed turbulence
The conditional averaging operations defined in § 4.1 allow us to compute the scale energy production and interscale transfer within a turbulent spot and compare them with the corresponding quantities in the fully turbulent region. Similar work has been done for single-point statistics, for example, the turbulent kinetic energy by Park et al. (Reference Park, Wallace, Wu and Moin2012), Nolan & Zaki (Reference Nolan and Zaki2013) and Marxen & Zaki (Reference Marxen and Zaki2019). To the best of our knowledge, this is the first time this type of analysis is extended to two-point statistics.
We start by deriving a decomposition similar to (5.3) for the production term 
$\mathcal {P}$ of the KHMH equation (3.6). This term consists of two components, one due to inhomogeneity of the mean flow in scale space, 
$\mathcal {P}_r=-2\overline {{\rm d} u'_i \, {\rm d} u'_{j}}({\partial \, {\rm d} U_i}/{\partial r_j})$, and the other due to inhomogeneity in physical space, 
$\mathcal {P}_s=-2\overline {{\rm d} u'_i u'^*_{j}} ({\partial \, {\rm d} U_i}/{\partial X_j})$. Here we decompose the former component, 
$\mathcal {P}_r$, which is the dominant one; this is also discussed below. Applying a process similar to that described in § 5, we obtain
\begin{align} &\mathcal{P}_r ={-}2\overline{{\rm d} u'_i \, {\rm d} u'_{j}}\frac{\partial \, {\rm d} U_i}{\partial r_j}\nonumber\\ &\quad ={-}2 (\gamma^{(TT)})^2\overline{{\rm d} u'^{(TT)}_i\, {\rm d} u'^{(TT)}_{j}}^{(TT)}\frac{\partial \, {\rm d} U_i^{(TT)}}{\partial r_j}- 2( \gamma^{(TL)})^2\overline{{\rm d} u'^{(TL)}_i\, {\rm d} u'^{(TL)}_{j}}^{(TL)}\frac{\partial \, {\rm d} U_i^{(TL)}}{\partial r_j}\nonumber\\ &\qquad -2(\gamma^{(LT)})^2\overline{{\rm d} u'^{(LT)}_i\, {\rm d} u'^{(LT)}_{j}}^{(LT)}\frac{\partial \, {\rm d} U_i^{(LT)}}{\partial r_j} -2(\gamma^{(LL)})^2 \overline{{\rm d} u'^{(LL)}_i\, {\rm d} u'^{(LL)}_{j}}^{(LL)}\frac{\partial \, {\rm d} U_i^{(LL)}}{\partial r_j} - \phi, \end{align}
where the additional term 
$\phi$ is given by
\begin{align}
\phi&= 2(\gamma^{(TT)})^2 \,{\rm d} U_i^{(TT)}\, {\rm d}
U_{j}^{(TT)}\frac{\partial \, {\rm d} U_i^{(TT)}}{\partial
r_j} +2(\gamma^{(TL)})^2 \,{\rm d} U_i^{(TL)}\, {\rm d}
U_{j}^{(TL)}\frac{\partial \, {\rm d} U_i^{(TL)}}{\partial r_j}\nonumber\\
&\quad +\,2(\gamma^{(LT)})^2 \,{\rm d}
U_i^{(LT)}\, {\rm d} U_{j}^{(LT)}\frac{\partial \, {\rm d}
U_i^{(LT)}}{\partial r_j} +2(\gamma^{(LL)})^2 \,{\rm d}
U_i^{(LL)} \,{\rm d} U_{j}^{(LL)}\frac{\partial \, {\rm d}
U_i^{(LL)}}{\partial r_j}\nonumber\\
&\quad +\,2\gamma^{(TT)}\, \overline{{\rm d} u_i\, {\rm d}
u_{j}}^{(TT)}\left[ \,{\rm d} U_i^{(TT)}\frac{\partial
\gamma^{(TT)}}{\partial r_j}+\frac{\partial}{\partial r_j}(
\gamma^{(LL)}\, {\rm d} U_i^{(LL)}\right.\nonumber\\
&\quad \left.+\,\gamma^{(TL)}\, {\rm d}
U_i^{(TL)} +\gamma^{(LT)}\, {\rm d} U_i^{(LT)})\vphantom{\frac{\partial
\gamma^{(TT)}}{\partial r_j}}\right]\nonumber\\
&\quad +\,2\gamma^{(TL)}\, \overline{{\rm
d} u_i\, {\rm d} u_{j}}^{(TL)}\left[{\rm d}
U_i^{(TL)}\frac{\partial \gamma^{(TL)}}{\partial
r_j}+\frac{\partial}{\partial r_j}( \gamma^{(LL)}\, {\rm d}
U_i^{(LL)}\right.\nonumber\\
&\quad \left. +\,\gamma^{(TT)}\, {\rm d} U_i^{(TT)}
+\gamma^{(LT)}\, {\rm d} U_i^{(LT)} )\vphantom{\frac{\partial
\gamma^{(TT)}}{\partial r_j}}\right]\nonumber\\
&\quad +\,2\gamma^{(LT)}\, \overline{{\rm d} u_i\, {\rm d}
u_{j}}^{(LT)}\left[{\rm d} U_i^{(LT)}\frac{\partial
\gamma^{(LT)}}{\partial r_j}+\frac{\partial}{\partial r_j}(
\gamma^{(LL)}\, {\rm d} U_i^{(LL)}\right.\nonumber\\
&\quad \left. +\,\gamma^{(TL)}\, {\rm d}
U_i^{(TL)} +\gamma^{(TT)}\, {\rm d} U_i^{(TT)} )\vphantom{\frac{\partial
\gamma^{(TT)}}{\partial r_j}}\right]\nonumber\\
&\quad +\,2\gamma^{(LL)}\, \overline{{\rm
d} u_i\, {\rm d} u_{j}}^{(LL)}\left[{\rm d}
U_i^{(LL)}\frac{\partial \gamma^{(LL)}}{\partial
r_j}+\frac{\partial}{\partial r_j}( \gamma^{(TT)}\, {\rm d}
U_i^{(TT)}\right.\nonumber\\
&\quad \left. \vphantom{\frac{\partial
\gamma^{(TT)}}{\partial r_j}} +\,\gamma^{(TL)}\, {\rm d} U_i^{(TL)}
+\gamma^{(LT)}\, {\rm d} U_i^{(LT)} )\right] -2\, {\rm d} U_i\, {\rm d} U_{j}\frac{\partial \,
{\rm d} U_i}{\partial r_j}.
\end{align} Note that the terms 
$\mathcal {P}_r^{(AA)}=-2\overline {{\rm d} u'^{(AA)}_i\, {\rm d} u'^{(AA)}_{j}}^{(AA)}({\partial \, {\rm d} U_i^{(AA)}}/{\partial r_j})$ that appear in (7.1) are the production terms of the conditionally averaged KHMH equation (4.8). Due to the strong shear in the wall-normal direction, the dominant component of 
$\mathcal {P}_r=-2\overline {{\rm d} u'_i \, {\rm d} u'_{j}}({\partial \, {\rm d} U_i}/{\partial r_j})$ is 
$\mathcal {P}_{r(1,2)}=-2\overline {{\rm d} u'_1 \, {\rm d} u'_{2}} ({\partial \, {\rm d} U_1}/{\partial r_2}) |_{r_2=0}$. It can be easily proved (see Yao et al. Reference Yao, Mollicone and Papadakis2022) that for 
$r_1=0$ this component is equal to 
$-2\overline {{\rm d} u'_1 \, {\rm d} u'_{2}} ({\partial U_1}/{\partial x_2})$. The corresponding production component due to inhomogeneity of the mean flow in physical space, 
$\mathcal {P}_s=-2\overline {{\rm d} u'_1 u'^*_{2}} ({\partial \, {\rm d} U_1}/{\partial X_2})$, is much smaller. Therefore in this section we consider only 
$\mathcal {P}_{r(1,2)}=-2\overline {{\rm d} u'_1 \, {\rm d} u'_{2}} ({\partial \, {\rm d} U_1}/{\partial r_2}) |_{r_2=0}$ and the corresponding conditionally averaged 
$(TT)$ component 
$\mathcal {P}_{r(1,2)}^{(TT)}=-2\overline {{\rm d} u'^{(TT)}_1du'^{(TT)}_{2}}^{(TT)} ({\partial \, {\rm d} U_1^{(TT)}}/{\partial r_2} ) |_{r_2=0}$ within a turbulent spot. We focus at location 
$TR2$ and compare the aforementioned component with 
$\mathcal {P}_{r(1,2)}$ evaluated at the fully turbulent region, 
$TU$.
We also extend the decomposition (5.3) to the total fluxes in physical and scale spaces 
$\boldsymbol {\phi _s}$ and 
$\boldsymbol {\phi _r}$, defined in (3.8) and (3.9) respectively. The resulting expressions are
\begin{equation} \boldsymbol{\phi}_{s_j}= \gamma^{(TT)}\boldsymbol{\phi}_{s_j}^{(TT)} + \gamma^{(TL)}\boldsymbol{\phi}_{s_j}^{(TL)}+ \gamma^{(LT)}\boldsymbol{\phi}_{s_j}^{(LT)}+ \gamma^{(LL)}\boldsymbol{\phi}_{s_j}^{(LL)}+ \phi,\end{equation}where
\begin{align}
\boldsymbol{\phi}_{s_j}^{(AA)}&=\overline{{U^*_j}^{(AA)}\,
{\rm d} q^{2^{(AA)}}}^{(AA)} + \overline{{u'^*_j}^{(AA)}\,
{\rm d} q^{2^{(AA)}}}^{(AA)}\nonumber\\
&\quad + \overline{2 \, {\rm d}
u'^{(AA)}_j \, {\rm d} p'^{(AA)} }^{(AA)} - \frac{1}{2} \nu
\overline{\frac{\partial \, {\rm d} q^{2^{(AA)}} }{\partial
X_j}}^{(AA)}\end{align}and
\begin{equation} \boldsymbol{\phi}_{r_j}= \gamma^{(TT)}\boldsymbol{\phi}_{r_j}^{(TT)}+ \gamma^{(TL)}\boldsymbol{\phi}_{r_j}^{(TL)}+ \gamma^{(LT)}\boldsymbol{\phi}_{r_j}^{(LT)}+ \gamma^{(LL)}\boldsymbol{\phi}_{r_j}^{(LL)}+ \phi,\end{equation}where
\begin{equation} \boldsymbol{\phi}_{r_j}^{(AA)}= \overline{{\rm d} u'^{(AA)}_j\, {\rm d} q^{2^{(AA)}} }^{(AA)} +\overline{{\rm d} U_j^{(AA)} \, {\rm d} q^{2^{(AA)}} }^{(AA)} -2\nu \overline{\frac{\partial \, {\rm d} q^{2^{(AA)}} }{\partial r_j}}^{(AA)}.\end{equation}The full expressions, including the remainder terms, are provided in Appendix B.
7.1. Conditionally and standard-averaged production and fluxes
In figure 13, contours of 
$\mathcal {P}_{r(1,2)}^{(TT)}$ at 
$TR2$ and of 
$\mathcal {P}_{r(1,2)}$ at 
$TR2$ and 
$TU$ are plotted in the 
$(r_3,Y)$ plane for 
$r_1=0$. It can be seen that the production peaks within the turbulent spot and the fully turbulent region are located at approximately the same spanwise separation and wall-normal height, 
$r_3\approx 5L_0$ and 
$Y \approx 1.3L_0$ respectively. On the other hand, the peak of 
$\mathcal {P}_{r(1,2)}$ at 
$TR2$ is found to be at larger 
$r_3$ separation and further away from the wall, 
$r_3 \approx 7.5L_0$ and 
$Y \approx 4.5L_0$. We mark the spanwise scales where the peaks appear, 
$r_3=5L_0$ and 
$7.5L_0$, with vertical dotted lines in figure 13, and plot the variation of the three production terms along these lines in figure 14. Notice the very close matching of the conditionally averaged production 
$\mathcal {P}_{r(1,2)}^{(TT)}$ at 
$TR2$ (dashed line) and the standard-averaged production 
$\mathcal {P}_{r(1,2)}$ at 
$TU$ (solid line) close to the wall (for 
$Y\le 1.3L_0$) while further away the two sets deviate. On the other hand, 
$\mathcal {P}_{r(1,2)}$ at 
$TR2$ shows significantly different behaviour even close to the wall, and of course peaks at a different distance.
Figure 13. Contour plots of conditionally averaged production 
$\mathcal {P}_{r(1,2)}^{(TT)}$ at 
$TR2$ (a), standard-averaged production 
$\mathcal {P}_{r(1,2)}$ at 
$TR2$ (b) and 
$TU$ (c) on the 
$(r_3,Y)$ plane for 
$r_1=0$. The red and black vertical dotted lines are placed at 
$r_3=5L_0$ and 
$7.5L_0$, respectively. The horizontal purple line indicates the local boundary layer thickness.
Figure 14. Variation of production terms along the wall-normal distance at 
$r_3/L_0=5.0$ and 
$7.5$. The variation is along the dotted vertical lines shown in figure 13 that pass through the corresponding production peaks. The plots are normalised by the value at 
$r_3/L_0=5.0$. The solid lines represent 
$\mathcal {P}_{r(1,2)}$ at 
$TU$. The dashed lines denote 
$\mathcal {P}_{r(1,2)}^{(TT)}$ and the circles 
$\mathcal {P}_{r(1,2)}$, both at location 
$TR2$.
Figure 15(a) shows stream-tubes in the 
$(r_1,r_3,Y)$ hyperplane obtained from the conditionally averaged total fluxes 
$({\phi _{r_1}^{(TT)}}, {\phi _{r_3}^{(TT)}},{\phi _{s_2}^{(TT)}})$ at 
$TR2$ together with an isosurface of the conditionally averaged production. In figure 15(b) we plot 
$({{\phi _{r_1}}},{{\phi _{r_3}}},{{\phi _{s_2}}})$ and production in the fully turbulent region. The latter figure reflects the dynamics of near-wall turbulence; energy is extracted from the mean flow at the buffer layer where the production peak is located, then it is transferred away from the wall and towards larger 
$r_1$ scales before bending back to smaller scales (dissipation region). This behaviour is related to the self-sustained turbulence mechanism near the wall (see Cimarelli et al. Reference Cimarelli, De Angelis and Casciola2013). A similar pattern can be discerned in figure 15(a), but the inverse cascade and flow of energy away from the wall is over a smaller range of 
$r_1$ separations (up to 
$r_1\approx 10$); the stream-tubes again bend towards small scales. There are also some deviations between the two plots for larger 
$r_3$ separations. If the centre of the sphere (used for identifying which stream-tubes to trace) is placed at smaller 
$r_3$ and the radius is reduced, the similarity between the two panels is more evident (see figure 16 and caption for details).
Figure 15. Stream-tubes of the conditionally averaged total flux vector 
$({\phi _{r_1}^{(TT)}}, {\phi _{r_3}^{(TT)}},{\phi _{s_2}^{(TT)}})$ at 
$TR2$ (a) and of the standard-averaged vector 
$({{\phi _{r_1}}},{{\phi _{r_3}}},{{\phi _{s_2}}})$ at 
$TU$ (b) in the three-dimensional (
$r_1,r_3,Y$) hyperplane. The stream-tubes are coloured according to the sign of 
$\phi _{r_3}$ (red for positive, blue for negative, thus indicating inverse and forward cascade respectively). The colour bar refers to the value of 
$\phi _{r_3}$. The plots were generated by placing a sphere of radius 
$15L_0$ at 
$(r_1=5L_0, r_3=5L_0, Y=1.3L_0)$ and tracing the stream-tubes crossing the sphere. Isosurfaces of the production term with values 
$0.9\times \max (\mathcal {P}_{r(1,2)}^{(TT)})$ (a) and 
$0.85\times \max (\mathcal {P}_{r(1,2)})$ (b) are shown in yellow.
Figure 16. Same as figure 15, but the sphere is now placed at 
$(r_1=5L_0, r_3=2.5L_0, Y=1.3L_0)$ and has smaller radius, 
$4L_0$. Isosurfaces of the production term with values 
$0.95\times \max (\mathcal {P}_{r(1,2)}^{(TT)})$ (a) and 
$0.85\times \max (\mathcal {P}_{r(1,2)})$ (b) are shown in yellow.
The shorter range of inverse cascade in the 
$TR2$ location compared with 
$TU$ is probably because the spots are still developing, the merging is not yet complete and thus they have a smaller footprint in the streamwise direction. Note also the similarities of figures 15(a) and 16(a) to figure 8(e) that shows only the nonlinear component of the flux vector. This similarity confirms that this is the most important component that determines the overall behaviour.
8. Conclusions
We apply conditional averaging to study the interscale energy transfer process during bypass transition. To this end, we define two-point intermittencies and apply them to decompose the energy fluxes into different components that depend on the local conditions at the two points used to define the flux; the points are both within a laminar region or a turbulent spot or straddle the laminar/turbulent interface. The flux terms are evaluated numerically directly in the scale space because conditional averaging does not commute with the spatial derivative operator.
In the 
$(r_1, r_3, Y)$ hyperplane, strong inverse cascade is found in the 
$r_1$ direction, due to the nonlinear fluxes across the downstream and upstream boundaries of a spot. For the former boundary, the inverse cascade extends over a larger range of 
$r_1$ separations compared with the latter boundary. We explain this finding by considering a propagating spot as it passes across a fixed streamwise location.
We derive also the conditionally averaged KHMH equation, which, although not used directly, inspires the conditional decomposition of the production term and the total energy fluxes when both points are located within a turbulent spot (
$TT$ events). We compare with the corresponding terms in the fully turbulent region and find significant similarities, but also some differences. In both plots, a cluster of stream-tubes originates from the production peak and transfers energy to larger scales before bending back to small scales and the near-wall region. This spiral shape is similar to that found in the fully turbulent region and in channel flow. However, the extent of the spiralling motion is confined to smaller separations, probably because the spots have not yet fully merged. Also, a smaller cluster of stream-tubes transfers energy in the 
$r_1$ and 
$r_3$ directions, which is not found in the fully turbulent region.
The conditional averaging approach for two-point statistics developed in the paper can be applied to other flow configurations that exhibit sharp interfaces, such as wakes and jets, where a TNTI separates the irrotational and vortical regions. Important questions remain to be answered; for example, do the conditionally averaged statistics exhibit self-similarity? How does this develop as the jet/wake expands? Research in this direction is left as future work.
Funding
H.Y. acknowledges financial support from the Department of Aeronautics, Imperial College London and the Imperial College-CSC scholarship. The authors also wish to acknowledge the UK Turbulence Consortium (UKTC) for providing access to the ARCHER high performance computing facility through EPSRC grant EP/R029326/1. G.P. is supported by EPSRC grants EP/X017273/1 and EP/W001748/1.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Calculation of derivatives of two-point quantities directly in scale space
In this appendix, we provide the steps for the numerical evaluation of derivatives of two-point quantities directly in scale space.
A.1. Calculation of ${\partial \, {\rm d} U_i^{(AA)}}/{\partial r_2}|_{r_2=0}$
(i) Consider a cell with centroid at distance
$Y=Y_0$ from the wall. The distances of the centroids of the cells located above and below are ${\rm d} y_1$ and ${\rm d} y_2$ respectively (see figure 17). Recall that the grid is non-uniform in the wall-normal direction, so ${\rm d} y_1 \ne {\rm d} y_2$. Calculate ${\rm d} U_i^{(AA)}(r_2={\rm d} y_1+{\rm d} y_2)=U_i^{(AA)}(Y_0+{\rm d} y_1)-U_i^{(AA)}(Y_0-{\rm d} y_2)$, i.e. the velocity difference between the red $+/-$ markers in figure 17.(ii) Calculate
${\rm d} U_i^{(AA)}(r_2=-({\rm d} y_1+{\rm d} y_2))=-{\rm d} U_i^{(AA)}(r_2={\rm d} y_1+{\rm d} y_2)$, i.e. the velocity difference between the blue $+/-$ markers.(iii) Use the central difference scheme to compute
(A1)and store the value at the midpoint (denoted with a green circle in figure 17).\begin{equation} \left. \frac{\partial \, {\rm d} U_i^{(AA)}}{\partial r_2} \right|_{r_2=0} = \frac{{\rm d} U_i^{(AA)}(r_2={\rm d} y_1+{\rm d} y_2) - {\rm d} U_i^{(AA)}(r_2={-}({\rm d} y_1+{\rm d} y_2))}{2({\rm d} y_1+{\rm d} y_2)}, \end{equation}(iv) Repeat steps (i)–(iii) for cell centroids at different heights
$Y_0$.(v) Because the mesh is non-uniform, the middle point is not located at
$Y_0$, so interpolate values of ${\partial \, {\rm d} U_i^{(AA)}}/{\partial r_2}|_{r_2=0}$ at midpoints to obtain the value at $Y_0$.
Figure 17. Sketch for the calculation of 
${\partial \, {\rm d} U_i^{(AA)}}/{\partial r_2} |_{r_2=0}$.
Comparison with evaluation at points 
$x_i^+$ and 
$x_i^-$ for the standard-averaged streamwise velocity shows that the results are identical (see figure 18.
Figure 18. Contour plot of 
${\partial \, {\rm d} U_1}/{\partial r_2} |_{r_2=0}$ in the 
$(r_3,Y)$ plane. Evaluation directly in scale space (a) and from 
$\frac {1}{2}( {{\partial U_1^+ }/{\partial x_2^+}}+{{\partial U_1^- }/{\partial x_2^-}})$ (b).
A.2. Calculation of ${\partial \, {\rm d} U_i^{(AA)}}/{\partial r_1}|_{r_1=0}$
(i) Calculate
${\rm d} U_i^{(AA)}(r_3,r_1=2{\rm \Delta} x)$, i.e. the velocity difference between the red $+/-$ points shown in the left-hand panel of figure 19, at a fixed height $y=Y_0$.(ii) Calculate
${\rm d} U_i^{(AA)}(r_3,r_1=-2 {\rm \Delta} x)$, i.e. the velocity difference between the blue $+/-$ points shown in the right-hand panel of figure 19, at a fixed height $y=Y_0$.(iii) Use the central difference scheme to compute
(A2)and store the value at the middle point (marked with a green circle in both panels of figure 19).\begin{equation} \left. \frac{\partial \, {\rm d} U_i^{(AA)}}{\partial r_1}\right|_{r_1=0}= \frac{{\rm d} U_i^{(AA)}(r_3,r_1=2 {\rm \Delta} x) -{\rm d} U_i^{(AA)}(r_3,r_1={-}2{\rm \Delta} x)}{2(2 {\rm \Delta} x)}, \end{equation}(iv) Repeat steps (i)–(iii) at different heights
$Y_0$.
Figure 19 shows the four points involved for 
${\partial \, {\rm d} U_i^{(AA)}}/{\partial r_1}|_{r_1=0}$. The process is similar for the evaluation of 
${\partial \, {\rm d} U_i^{(AA)}}/{\partial r_1}$ for different 
$r_1$ values. Comparison with evaluation at points 
$x_i^+$ and 
$x_i^-$ for the standard-averaged streamwise velocity shows identical results (see figure 20).
Figure 19. Sketch for the calculation of 
${\partial \, {\rm d} U_i^{(AA)}}/{\partial r_1}|_{r_1=0}$.
Figure 20. Contour plot of 
${\partial \, {\rm d} U_1}/{\partial r_1} |_{r_1=0}$ in the 
$(r_3,Y)$ plane. Evaluation directly in scale space (a) and from 
$\frac {1}{2}( {{\partial U_1^+ }/{\partial x_1^+}}+{{\partial U_1^- }/{\partial x_1^-}})$ (b).
Appendix B. Conditional decomposition of total fluxes in scale and physical spaces
B.1. Decomposition of the total flux in scale space
\begin{equation} \boldsymbol{\phi}_{r_j}= \gamma^{(TT)}\boldsymbol{\phi}_{r_j}^{(TT)}+ \gamma^{(TL)}\boldsymbol{\phi}_{r_j}^{(TL)}+ \gamma^{(LT)}\boldsymbol{\phi}_{r_j}^{(LT)}+ \gamma^{(LL)}\boldsymbol{\phi}_{r_j}^{(LL)}+ \phi, \end{equation}where
\begin{align} \boldsymbol{\phi}_{r_j}^{(AA)}&= \overline{{\rm d} u'^{(AA)}_j\, {\rm d} q^{2^{(AA)}}}^{(AA)} + \overline{{\rm d} U_j^{(AA)} \, {\rm d} q^{2^{(AA)}} }^{(AA)} - 2\nu \overline{\frac{\partial \, {\rm d} q^{2^{(AA)}} }{\partial r_j}}^{(AA)}, \end{align}
\begin{align}
\phi&=\sum_{AA=}^{TT,TL,LT,LL} \left[ \gamma^{(AA)}(
\overline{({\rm d} u_i)^2}^{(AA)}\, {\rm d} U_j^{(AA)}
+\overline{2\, {\rm d} u_i\, {\rm d} u_j}^{(AA)}\, {\rm d}
U_i^{(AA)}\right.\nonumber\\
&\quad \left.-2\, {\rm d} U_i^{(AA)}\, {\rm d} U_i^{(AA)}\,
{\rm d} U_j^{(AA)} )\vphantom{\frac{\partial \, {\rm d} u_i
}{\partial r_j}} + \gamma^{(AA)}\, {\rm d} U_i^{(AA)}\, {\rm d} U_i^{(AA)}\,
{\rm d} U_j -2\nu\gamma^{(AA)}\right.\nonumber\\
&\quad \times \left.\left(2\, {\rm d}
U_i^{(AA)}\frac{\partial \, {\rm d} U_i^{(AA)} }{\partial
r_j} +2\, {\rm d} U_i^{(AA)}\overline{ \frac{\partial \,
{\rm d} u_i }{\partial r_j}}^{(AA)} -\frac{\partial ({\rm
d} U_i^{(AA)})^2 }{\partial r_j}\right) \right] \nonumber\\
&\quad -\overline{({\rm d} u_i)^2}\, {\rm d} U_j
-\overline{2\, {\rm d} u_i\, {\rm d} u_j}\, {\rm d} U_i
+{\rm d} U_i\, {\rm d} U_i\, {\rm d} U_j\nonumber\\
&\quad -2\nu\left({-}2\,
{\rm d} U_i\frac{\partial \, {\rm d} U_i }{\partial r_j}
-2\, {\rm d} U_i\overline{ \frac{\partial \, {\rm d} u_i
}{\partial r_j}} +\frac{\partial ({\rm d} U_i)^2 }{\partial
r_j} \right).
\end{align}B.2. Decomposition of the total flux in physical space
\begin{equation} \boldsymbol{\phi}_{s_j}= \gamma^{(TT)}\boldsymbol{\phi}_{s_j}^{(TT)} + \gamma^{(TL)}\boldsymbol{\phi}_{s_j}^{(TL)}+ \gamma^{(LT)}\boldsymbol{\phi}_{s_j}^{(LT)}+ \gamma^{(LL)}\boldsymbol{\phi}_{s_j}^{(LL)}+ \phi, \end{equation}where
\begin{align}
\boldsymbol{\phi}_{s_j}^{(AA)}&= \overline{{U^*_j}^{(AA)}\,
{\rm d} q^{2^{(AA)}}}^{(AA)} +\overline{{u'^*_j}^{(AA)}\,
{\rm d} q^{2^{(AA)}}}^{(AA)}\nonumber\\
&\quad +\overline{2 \, {\rm d}
u'^{(AA)}_j \, {\rm d} p'^{(AA)} }^{(AA)} -\frac{1}{2} \nu
\overline{\frac{\partial \, {\rm d} q^{2^{(AA)}} }{\partial
X_j}}^{(AA)},
\end{align}
\begin{align}
\phi&=\sum_{AA=}^{TT,TL,LT,TT} \left[ \gamma^{(AA)}(
\overline{({\rm d} u_i)^2}^{(AA)}{U^*_j}^{(AA)}
+\overline{2\, {\rm d} u_i{U^*_j}}^{(AA)}\, {\rm d}
U_i^{(AA)}\right.\nonumber\\
&\quad -2\, {\rm d} U_i^{(AA)}\, {\rm d}
U_i^{(AA)}{U^*_j}^{(AA)} ) +\gamma^{(AA)}\, {\rm d} U_i^{(AA)}\, {\rm d}
U_i^{(AA)}U^*_j -\frac{1}{2}\nu\gamma^{(AA)}\nonumber\\
&\quad \times \left(2\, {\rm
d} U_i^{(AA)}\frac{\partial \, {\rm d} U_i^{(AA)}
}{\partial X_j} + 2\, {\rm d} U_i^{(AA)}\overline{
\frac{\partial \, {\rm d} u_i }{\partial X_j}}^{(AA)} -
\frac{\partial ({\rm d} U_i^{(AA)})^2 }{\partial
X_j}\right) \nonumber\\
&\quad \left. +\,\vphantom{\frac{\partial \, {\rm d} U_i^{(AA)} }{\partial
X_j}}2\gamma^{(AA)}\, {\rm d} U_i^{(AA)}\, {\rm d} P^{(AA)}
\right] -\frac{1}{2} \nu\left({-}2\, {\rm d}
U_i\frac{\partial \, {\rm d} U_i }{\partial X_j} -2\, {\rm
d} U_i\overline{ \frac{\partial \, {\rm d} u_i }{\partial
X_j}} +\frac{\partial ({\rm d} U_i)^2 }{\partial X_j}\right)\nonumber\\
&\quad -\overline{({\rm d} u_i)^2}U^*_j -\overline{2\, {\rm d} u_iU^*_j}\, {\rm d} U_i
+{\rm d} U_i\, {\rm d} U_iU^*_j-2\, {\rm d} U_i\, {\rm d}
P.
\end{align}
Open access
