4.1. Mean passive scalar statistics

The mean passive scalar concentration $\overline {\Theta }_c$ is scaled by the mean centreline concentration, which is linearly related to the distance from the orifice

(4.1) \begin{equation} \frac {\overline {\Theta }_c(z)}{\overline {\Theta }_0}=\frac {B_\Theta D}{z-z_0}. \end{equation}

From the present DNS, we deduce a decay constant of $B_\Theta =4.9$ with a virtual origin at $z_0/D=-0.132$ , which means the virtual origin is close to the virtual origin of the velocity $z_0/D=0$ (Nguyen & Oberlack Reference Nguyen and Oberlack2024a ) taken from the classical velocity scaling (Boersma et al. Reference Boersma, Brethouwer and Nieuwstadt1998)

(4.2) \begin{equation} \frac {\overline {U}_{z,c}(z)}{\overline {U}_0}=\frac {B_u D}{z-z_0}. \end{equation}

Figure 1. Mean inverse passive scalar concentration at the centreline over the distance of the orifice: Birch et al. (Reference Birch, Brown, Dodson and Thomas1978) (), Babu & Mahesh (Reference Babu and Mahesh2005) (), Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001) (), present DNS ().

This is compared with Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001), where they found $z_0/D=0.5$ for the concentration and $z_0/D=5.5$ for the axial velocity. Figure 1 showcases a consistent and smoother linear trend of the inverse passive scalar concentration up to $z/D=75$ compared with other experiments and DNS. A collection of the data from the experiments shown in figure 1 can be found in table 1.

Table 1. Various jet parameters of DNS and experiments.

For the velocity the decay constant is $B_u=5.15$ (Nguyen & Oberlack Reference Nguyen and Oberlack2024a ), implying that $\overline {\Theta }_c$ decays faster than the velocity. Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001) have observed the same, but their virtual origins for the velocity and the passive scalar concentration were further apart. The value of $\overline {\Theta }_c$ reaches self-similarity at $z/D=15$ , as depicted in figure 2. This onset occurs earlier than the assumed $z/D=20$ in Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001). Upon comparing the profile of $\overline {U}_{z,c}$ with that of $\overline {\Theta }_c$ , the latter shows a slightly larger spreading rate. The Prandtl number describes the relative thicknesses of the velocity and thermal boundary layers. With $Pr=0.71$ , the thermal boundary layer is thicker than the velocity boundary layer, leading to a larger spreading rate. The half-width of the mean passive scalar concentration, $\eta _{1/2,\Theta }=0.109$ , exceeds that of the velocity, $\eta _{1/2,u}=0.089$ , as reported in Nguyen & Oberlack (Reference Nguyen and Oberlack2024a ), where $\eta =r/(z-z_0)$ is the usual jet-type similarity coordinate, as derived below in (5.20). Remarkably, $\eta _{1/2,\Theta }$ closely aligns with the observation in Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001), where $\eta _{1/2,\Theta }=0.108$ in their DNS.

Figure 2. Mean passive scalar $\overline {\Theta }$ scaled with the centreline passive scalar concentration $\overline {\Theta }_c$ and plotted as a function of the similarity coordinate $\eta$ (see (5.20)) at different distances from the orifice: $z/D=15$ (), $25$  (), $35$ (), $45$ (), $55$ ().

In summary, the characteristics of $\overline {\Theta }_c$ are consistent with those observed by Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001). The profiles exhibit self-similarity and a faster spread compared with the velocity profile. Moreover, the virtual origins of concentration and velocity appear to have shifted closer to each other compared with the DNS by Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001), which is attributed to the turbulent pipe inlet.

4.2. Second-order passive scalar and mixed statistics

The variance of the passive scalar fluctuation $R_{\Theta \Theta }$ is shown in figure 3 for different distances $z$ from the orifice. The value at $\eta =0$ agrees with the results of Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001) and Darisse et al. (Reference Darisse, Lemay and Benaïssa2015), both of which report values around $0.04$ and slightly above. However, the off-axis peak exceeds the value of $0.05$ measured by Darisse et al. (Reference Darisse, Lemay and Benaïssa2015). The peak is at $\eta =0.082$ , which is slightly towards the centre of the inflection point of the passive scalar profile in figure 2 at $\eta =0.089$ . The inflection point in the scalar profile corresponds to the region where the scalar is being transported most effectively from the jet core to the surrounding fluid. In turbulent round jet flows, turbulent eddies begin to form at the edge of the jet inlet. These eddies grow and interact with the scalar field as the jet develops. This process enhances the mixing of the scalar within the shear layer between jet and the surrounding ambient fluid, leading to fluctuations in the scalar field. The mixing begins slightly inside the region where the scalar gradient is steepest which may lead to the peak in fluctuations on the centre-facing side of the inflection point.

Figure 3. Variance of the passive scalar fluctuations $R_{\Theta \Theta }$ scaled with the centreline passive scalar concentration $\overline {\Theta }^{2}_c$ and plotted as a function of the similarity coordinate $\eta$ at different distances from the orifice: $z/D=25$ (), $35$ (), $45$ (), $55$ ().

Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001) notes in figure 4 that the normalised root mean square (r.m.s.) of the passive scalar fluctuation $\sqrt {R_{\Theta \Theta }}$ shows an increasing trend rather than a horizontal line, the former indicating non-self-similarity. This trend can be attributed to the smaller axial length of their computational box, as the present DNS shows a slight increase up to $z/D=35$ . After this point, however, $\sqrt {R_{\Theta \Theta }}$ stabilises around $0.22$ , supporting Dowling & Dimotakis’s (Reference Dowling and Dimotakis1990) conclusion that $\sqrt {R_{\Theta \Theta }}$ lies between $0.23$ and $0.24$ , with the present DNS slightly below this range.

Figure 4. Centreline r.m.s. of the passive scalar fluctuation $\sqrt {R_{\Theta \Theta }}$ scaled with the centreline mean passive scalar concentration over the distance $z$ from the orifice: Darisse et al. (Reference Darisse, Lemay and Benaïssa2015) (), Birch et al. (Reference Birch, Brown, Dodson and Thomas1978) (), Babu & Mahesh (Reference Babu and Mahesh2005) (), Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001) (), present DNS ().

The turbulent diffusive fluxes shown in figures 5 and 6 exhibit a close resemblance to the radial profiles in Darisse et al. (Reference Darisse, Lemay and Benaïssa2015). The measured maximum values for $R_{r\Theta }/(\overline {U}_{z,c}\overline {\Theta }_c)$ and $R_{z\Theta }/(\overline {U}_{z,c}\overline {\Theta }_c)$ in Darisse et al. (Reference Darisse, Lemay and Benaïssa2015) are $2.2$ and $3$ , respectively, consistent with the $2.2$ and $3.2$ values observed in figures 5 and 6. In addition, the value at $\eta =0$ for $R_{z\Theta }/(\overline {U}_{z,c}\overline {\Theta }_c)$ is approximately $0.024$ in Darisse et al. (Reference Darisse, Lemay and Benaïssa2015), closely matching the range of $0.025-0.026$ in the present DNS.

Figure 5. Turbulent heat flux $R_{r\Theta }$ normalised with $\overline {U}_{z,c}\overline {\Theta }_c$ at different distances from the orifice: $z/D=25$ (), $35$ (), $45$ (), $55$ ().

Figure 6. Turbulent heat flux $R_{z\Theta }$ normalised with $\overline {U}_{z,c}\overline {\Theta }_c$ at different distances from the orifice: $z/D=25$ (), $35$ (), $45$ (), $55$ ().

4.3. Third-order mixed statistics

There are only quite few experimental data on the radial profiles of the third-order mixed statistics $R_{i\Theta \Theta }$ , due to their sensitivity to disturbances in experimental set-ups. In contrast, the present DNS provides a well-converged data set, allowing a detailed examination of these statistics. Comparable experiments are mainly the helium jet of Panchapakesan & Lumley (Reference Panchapakesan and Lumley1993b ) and the heated jet by Darisse et al. (Reference Darisse, Lemay and Benaïssa2015), which provide a basis for validation. Notably, Antonia, Prabhu & Stephenson (Reference Antonia, Prabhu and Stephenson1975) also contributed data for $R_{i\Theta \Theta }$ , although using a jet with a strong co-flow. Nevertheless, the shape of the data is similar to the experiments of Darisse et al. (Reference Darisse, Lemay and Benaïssa2015), but with different magnitudes. The radial profiles of the normalised third-order mixed statistics of the present DNS are showcased in figures 7 and 8 for $R_{r\Theta \Theta }$ and $R_{z\Theta \Theta }$ , respectively. Qualitatively, the shapes of the aforementioned experiments are similar to those of the present DNS. While the maxima of the present DNS closely match Darisse et al. (Reference Darisse, Lemay and Benaïssa2015), the minima have slightly larger magnitudes, and Panchapakesan & Lumley (Reference Panchapakesan and Lumley1993b ) presents an even larger minimum for $R_{r\Theta \Theta }$ .

Figure 7. Turbulent heat flux $R_{r\Theta \Theta }$ normalised with $\overline {U}_{z,c}\overline {\Theta }_c^{2}$ at different distances from the orifice: $z/D=25$ (), $35$ (), $45$ (), $55$ ().

Figure 8. Turbulent heat flux $R_{z\Theta \Theta }$ normalised with $\overline {U}_{z,c}\overline {\Theta }_c^{2}$ at different distances from the orifice: $z/D=25$ (), $35$ (), $45$ (), $55$ ().

4.4. Probability density functions

In figure 9, the PDFs of the passive scalar concentration on the centreline at different orifice distances exhibit excellent collapse, closely resembling a Gaussian distribution. The PDFs are slightly negatively skewed, which can also be read off in figure 11 at $\eta =0$ . The negative skewness on the centreline for a passive scalar concentration was also observed in the methane jet of Birch et al. (Reference Birch, Brown, Dodson and Thomas1978).

Figure 9. Probability density functions of $\Theta (\eta =0,z)/\overline {\Theta }_{c}(z)$ at $z/D=15$ (), $25$ (), $35$ (), $45$ (), $55$ (), $65$ () compared with a Gaussian (). Here, $\Theta /\overline {\Theta }_{c}(z)$ is also shown in terms of the standard deviation $\sigma =0.215$ .

Figure 10. Probability density functions of $\Theta (\eta ,z)/\overline {\Theta }_{c}(z)$ for $z/D=28$ (), $42$ (), $56$ ().

Figure 10 depicts the radial PDF evolution at $z/D=28, 42, 56$ . As $\eta$ increases, significant deviations from Gaussian behaviour and the emergence of heavy and skewed tails becomes apparent. The PDF approaches a delta distribution for large $\eta$ due to the passive scalar concentration being zero in the ambient region. Additionally, we can extract that the PDFs collapse due to the scaling of $\eta$ . Quantitative insights from figure 11 show the collapse of skewness $S$ and kurtosis $K$ for $z/D=28, 42, 56$ . The Gaussian values, $S=0$ and $K=3$ , serve as references. From figure 11 we can also see that the PDF becomes more sub-Gaussian as $\eta$ increases. This is consistent with the observations in Birch et al. (Reference Birch, Brown, Dodson and Thomas1978), where the distribution becomes broader with increasing $\eta$ . This is expected due to intermittent mixing of the turbulent flow in the inner region and the ambient laminar fluid in the outer region. For larger $\eta$ , the mean value of $\Theta (\eta ,z)/\overline {\Theta }_c(z)$ approaches zero. As a result, the PDF can no longer be approximated by a Gaussian because it would require negative values of $\Theta (\eta ,z)/\overline {\Theta }_c(z)$ to maintain a Gaussian distribution. In Birch et al. (Reference Birch, Brown, Dodson and Thomas1978), the distribution becomes bimodal and they interpret this as intermittency in the flow. This is similar to what is observed in figure 12 at $\eta =0.139$ for $z/D=28$ . As the distance from the centreline increases to $\eta =0.149$ , the probability of zero passive scalar concentration quickly becomes dominant.

Figure 11. Kurtosis $K$ (above) and skewness $S$ (below) of $\Theta (\eta ,z)/\overline {\Theta }_{c}(z)$ for $z/D=28$ (), $42$ (), $56$  (). Here, $K=3$ , $S=0$ (dashed) are the Gaussian values.

Figure 12. Probability density functions of $\Theta (\eta ,z=28)/\overline {\Theta }_{c}(z=28)$ for $\eta =0.139$ (solid) and $\eta =0.149$ (dashed).

5.1. Multi-point correlation equations for passive scalar and velocity moments

In the first step, our goal is to derive the MPSEs and then extend this derivation to the MPVSCEs, which include both the MPMEs and the MPSEs. Similar to the MPMEs, the MPSEs are derived by multiplying (2.4) with $m-1$ passive scalars at $m-1$ different locations $\boldsymbol{x}_{(l)}$ followed by statistical averaging. Both the MPSEs and MPVSCEs are an infinite set of linear equations, analogical to the MPMEs. The resulting $m$ th-order MPSE is expressed as

(5.1) \begin{equation} \frac {\partial H_{\Theta _{\{m\}}}}{\partial t}+\sum _{l=1}^{m}\left [\frac {\partial }{\partial x_{k_{(l)}}} \left [H_{\Theta _{\{m+1\}}[\Theta _{(m+1)}\mapsto k_{(l)}]}\left [\boldsymbol{x}_{(m+1)}\mapsto \boldsymbol{x}_{(l)}\right ]\right ] -\frac {1}{Pe}\frac {\partial ^{2} H_{\Theta _{\{m\}}}}{\partial x_{k_{(l)}}\partial x_{k_{(l)}}}\right ]=0, \end{equation}

for $m=1,\ldots ,\infty ,$ where

(5.2) \begin{equation} H_{\Theta _{\{m\}}}=\overline {\Theta (\boldsymbol{x}_{(1)})\Theta (\boldsymbol{x}_{(2)}) \cdot \ldots \cdot \Theta (\boldsymbol{x}_{(m)})}, \end{equation}

and

(5.3) \begin{align} &H_{\Theta _{\{m+1\}}[\Theta _{(l)}\mapsto k_{(l)}]}\left [\boldsymbol{x}_{(l)}\mapsto \boldsymbol{x}_{(p)}\right ]\nonumber \\ & = \overline {\Theta (\boldsymbol{x}_{(1)})\cdot \ldots \cdot \Theta (\boldsymbol{x}_{(l-1)}) U_{k_{(l)}}(\boldsymbol{x}_{(p)})\Theta (\boldsymbol{x}_{(l+1)})\cdot \ldots \cdot \Theta (\boldsymbol{x}_{(m+1)})}. \end{align}

A detailed derivation is available in Appendix A. In a similar approach, the MPVSCEs can be derived with (2.2) and (2.4) giving

(5.4) \begin{align} &\frac {\partial H_{i_{\{n\}}\Theta _{\{m\}}}}{\partial t}\nonumber \\[10pt] & + \sum _{l=1}^{n} \left (\frac {\partial }{\partial x_{k_{(l)}}}\left [H_{i_{\{n+1\}}\Theta _{\{m\}}\left [i_{(n+m+1)}\mapsto k_{(l)}\right ]}\left [\boldsymbol{x}_{(n+m+1)}\mapsto \boldsymbol{x}_{(l)} \right ]\right] \right. \nonumber\\[10pt]&\left. +\,\frac {\partial I_{i_{\{n-1\}}\Theta _{\{m\}}\left [ l\right ]}}{\partial x_{i_{(l)}}} -\frac {1}{Re}\frac {\partial ^{2} H_{i_{\{n\}}\Theta _{\{m\}}}}{\partial x_{k_{(l)}} \partial x_{k_{(l)}}}\right )\nonumber \\[10pt] & +\, \sum _{l=n+1}^{n+m} \left (\frac {\partial }{\partial x_{k_{(l)}}}\left [H_{i_{\{n+1\}}\Theta _{\{m\}}\left [i_{(n+m+1)}\mapsto k_{(l)}\right ]}\left [\boldsymbol{x}_{(n+m+1)}\mapsto \boldsymbol{x}_{(l)} \right ]\right ]-\frac {1}{Pe}\frac {\partial ^{2} H_{i_{\{n\}}\Theta _{\{m\}}}}{\partial x_{k_{(l)}} \partial x_{k_{(l)}}}\right )=0, \nonumber\\\end{align}

for $n,m=1,\ldots ,\infty$ , where

(5.5) \begin{equation} H_{i_{\{n\}}\Theta _{\{m\}}}=\overline {U_{i_{(1)}}(\boldsymbol{x}_{(1)})\cdot \ldots \cdot U_{i_{(n)}}(\boldsymbol{x}_{(n)})\Theta (\boldsymbol{x}_{(n+1)}) \cdot \ldots \cdot \Theta (\boldsymbol{x}_{(n+m)})} ,\end{equation}

and the pressure correlation term is defined as

(5.6) \begin{equation} I_{i_{\{n-1\}}\Theta _{\{m\}}\left [ l\right ]}= \overline {U_{i_{(1)}}(\boldsymbol{x}_{(1)}) \cdot \ldots \cdot P(\boldsymbol{x}_{(l)}) \cdot \ldots \cdot U_{i_{(n)}}(\boldsymbol{x}_{(n)})\Theta (\boldsymbol{x}_{(n+1)}) \cdot \ldots \cdot \Theta (\boldsymbol{x}_{(n+m)})}. \end{equation}

More details about the derivation can be found in Appendix B. Setting $m=0$ yields the MPMEs from (5.4), while $n=0$ yields the MPSEs. The MPVSCEs (5.4) are a set of linear differential equations, a crucial feature for the occurrence of the two statistical symmetries (Wacławczyk et al. Reference Wacławczyk, Staffolani, Oberlack, Rosteck, Wilczek and Friedrich2014) to be used below.

5.2. Symmetry-based scaling laws for a passive scalar

The methodology of Lie-symmetry analysis uses continuous transformation groups to unify methods for solving differential equations and finds applications in various fields of applied mathematics and theoretical physics. For a comprehensive understanding, see Bluman, Cheviakov & Anco (Reference Bluman, Cheviakov and Anco2010). In the following, the concept of symmetries of differential equations is briefly introduced. Symmetries in differential equations are transformations that leave the form of the equations unchanged. Consider a system of partial differential equations

(5.7) \begin{equation} \boldsymbol{F}\left (\boldsymbol{x}, \boldsymbol{u}, \boldsymbol{u}^{(1)}, \ldots ,\boldsymbol{u}^{(p)}\right )=0, \end{equation}

with a transformation

(5.8) \begin{equation} \boldsymbol{x}^{*} = \boldsymbol{\phi }(\boldsymbol{x},\boldsymbol{u};a), \quad \boldsymbol{u}^{*} = \boldsymbol{\psi }(\boldsymbol{x},\boldsymbol{u};a), \end{equation}

where $\boldsymbol{x}$ , $\boldsymbol{u}$ , $a$ are the independent variables, dependent variables and an arbitrary continuous parameter $a \in \mathbb{R}$ , respectively. The transformation is considered a symmetry transformation if

(5.9) \begin{equation} \boldsymbol{F}\left (\boldsymbol{x}, \boldsymbol{u}, \boldsymbol{u}^{(1)}, \ldots ,\boldsymbol{u}^{(p)}\right )=0 \iff \boldsymbol{F}\left (\boldsymbol{x}^{*}, \boldsymbol{u}^{*}, \boldsymbol{u}^{*(1)}, \ldots ,\boldsymbol{u}^{*(p)}\right )=0, \end{equation}

where $\boldsymbol{u}^{(p)}$ refers to the $p$ th derivative of $\boldsymbol{u}$ . In this case, we consider (2.2) with $Re^{-1}=0$ , which admits the following two-parameter symmetry:

(5.10) \begin{equation} t^\ast = \textrm { e}^{a_{St}} t, \thinspace \boldsymbol{x}^\ast = \textrm { e}^{a_{Sx}} \boldsymbol{x}, \thinspace \boldsymbol{U}^\ast = \textrm { e}^{a_{Sx}-a_{St}} \boldsymbol{U}, \thinspace P^\ast = \textrm { e}^{2(a_{Sx}-a_{St})} P, \thinspace a_{Sx},a_{St}\in \mathbb{R} . \end{equation}

If we now insert (5.10) into (2.2) under the frictionless assumption, we find the prefactor $\textrm { e}^{2a_{St}-a_{Sx}}$ in front of each of the three terms, which can be cancelled out, and we get $({(\partial \boldsymbol{U}^\ast )}/{(\partial \boldsymbol{U}{t^\ast})}) + (\boldsymbol{U}^\ast \cdot \nabla ^\ast ) \boldsymbol{U}^\ast = - \nabla ^\ast P^\ast$ . This proves the scaling symmetry. Two things are important for the following analyses: (i) the scaling symmetries of (2.2) but also (2.4) transfer to the moment equations; (ii) the scaling symmetry of time characterised by the group parameter $a_{St}$ is preserved for the stationary case of the Euler equation, because the above-mentioned scaling factor does not change for this case, even if the unsteady term is no longer included. This also remains the case if one considers a statistically stationary flow for the present case, because the key variables such as velocity $\boldsymbol{U}$ and pressure $P$ apparently also contain time as a basic variable, and this applies equally to the corresponding statistical moments.

The following set of selected symmetries are relevant for the derivation of the scaling laws of a turbulent round jet, which have their origin in the NSEs and the passive scalar equations. In the limit of $Re \to \infty$ they are transferred to the H-approach, which consist of a translation symmetry in space

(5.11) \begin{align} \overline {T}_{x}: \quad &t^{*}=t, \quad {x}^{*}_i={x}_i+a_{x_i},\quad \overline {\Theta }^{*}= \overline {\Theta }, \nonumber \\[6pt] & H^{*}_{\Theta _{\{m\}}}=H_{\Theta _{\{m\}}},\quad H^{*}_{i_{\{n\}}\Theta _{\{m\}}}=H_{i_{\{n\}}\Theta _{\{m\}}}, \end{align}

and a scaling symmetry in space ( $a_{Sx}$ ), time ( $a_{St}$ ) and of the passive scalar ( $a_{S\theta }$ )

(5.12) \begin{align} & \overline {T}_{S}: \quad t^{*}=\mathrm{e}^{a_{St}}t, \quad {x}^{*}_i=\mathrm{e}^{a_{Sx}}{x}_i,\quad \overline {\Theta }^{*}=\mathrm{e}^{a_{S\theta }} \overline {\Theta }, \nonumber \\[6pt] & H^{*}_{\Theta _{\{m\}}}=\mathrm{e}^{m a_{S\theta }} H_{\Theta _{\{m\}}},\quad H^{*}_{i_{\{n\}}\Theta _{\{m\}}}=\mathrm{e}^{ma_{S\theta }+n(a_{Sx}-a_{St})}H_{i_{\{n\}}\Theta _{\{m\}}}. \end{align}

It is noteworthy that the NSEs exhibit a restricted set of scaling symmetries compared with the Euler equations, in particular $a_{St} = 2a_{Sx}$ in the symmetry above which can be shown by substituting the scaling symmetry (5.10) into the NSEs. For turbulent flows at higher $Re$ , viscosity dominates only on length scales between the Taylor and the Kolmogorov length scales or the corresponding wavenumbers in spectral space. This is supported by Oberlack (Reference Oberlack2000), where a boundary layer type of asymptotic expansion was performed in correlation space. As a result, the multi-point equations for turbulent flows for large scales admit the scaling symmetries of the inviscid form of the equation and, hence, to leading order, no $Re$ dependence is implied. In full analogy, this can also be shown for the $Pe$ number dependence. Additionally, a statistical symmetry on the basis of the linear MPVSCEs in (5.4) is considered

(5.13) \begin{align} \overline {T}_{Ss}: \quad &t^{*}=t, \quad {x}^{*}_i={x}_i,\quad \overline {\Theta }^{*}=\mathrm{e}^{a_{Ss}} \overline {\Theta }, \\[6pt] & H^{*}_{\Theta _{\{m\}}}=\mathrm{e}^{a_{Ss}} H_{\Theta _{\{m\}}},\quad H^{*}_{i_{\{n\}}\Theta _{\{m\}}}=\mathrm{e}^{a_{Ss}}H_{i_{\{n\}}\Theta _{\{m\}}},\nonumber \end{align}

which is a measure of intermittency. In Wacławczyk et al. (Reference Wacławczyk, Staffolani, Oberlack, Rosteck, Wilczek and Friedrich2014) this has been concluded by analysing the PDF formulation of the moment equation known as the Lundgren–Monin–Novikov hierarchy. In this formulation, the group parameter $a_{Ss}$ of an equivalent symmetry scales the shape of a PDF such that intermittency is characterised. This so-called intermittency symmetry does not appear in either the NSEs or the passive scalar equation. For the linear systems (5.1)–(5.6), there exists also the generic symmetry of superposition, which does not play a role in the present subsection, however, and only comes into play in §§ 5.4. The combination of the symmetries (5.11)–(5.13) leads to the characteristic condition for the invariant solution, i.e. turbulent scaling laws of the velocity and passive scalar moments of a spatially evolving turbulent round jet. Although the system (5.1)–(5.6) constitutes an arbitrary multi-point moments hierarchy, in the following, we will only consider one-point statistics and scaling laws, i.e. $\boldsymbol{x}_{(1)}=\boldsymbol{x}_{(2)}=\ldots =\boldsymbol{x}_{(n+m)}$ and therefore $H_{i_{\{n\}}\Theta _{\{m\}}}=\overline {U^n_{[i]}\Theta ^m}$ , where $[\,\cdot \,]$ means that no summation is implied. The related characteristic condition reads

(5.14) \begin{align} \frac {\mathrm{d}r}{a_{Sx}r}&=\frac {\mathrm{d}z}{a_{Sx}z+a{_{z}}}=\frac {\mathrm{d}\overline {U_{[i]}^{n} \Theta ^m}}{[n(a_{Sx}-a_{St})+ma_{S\theta }+a_{Ss}]\overline {U_{[i]}^n \Theta ^m}}\nonumber \\[6pt] &=\frac {\mathrm{d}\overline {\Theta }}{[a_{S\theta }+a_{Ss}]\overline {\Theta }}=\ldots =\frac {\mathrm{d}\overline {\Theta ^m}}{[ma_{S\theta }+a_{Ss}]\overline {\Theta ^m}}. \end{align}

Symmetry breaking is introduced through flow-specific invariants for the turbulent jet, i.e. the thermal energy conservation (Sadeghi, Oberlack & Gauding Reference Sadeghi, Oberlack and Gauding2021)

(5.15) \begin{equation} I_\Theta =\int _{0}^\infty \overline {U_z \Theta }\,r\,\mathrm{d}r ,\end{equation}

and the generalised momentum integral as derived in Nguyen & Oberlack (Reference Nguyen and Oberlack2024b )

(5.16) \begin{equation} I_O=\int _0^\infty \left [ \overline {U_z^{2}}-\frac {1}{2}\left (\overline {U_r^{2}}+\overline {U_\varphi ^{2}}\right )+\frac {1}{2} \frac {\mathrm{d}}{\mathrm{d}z} \left (\overline {U_rU_z}\,r\right )\right ]r\,\mathrm{d}r .\end{equation}

The group parameters of the symmetries are constrained through these invariants. The symmetry breaking is induced by implementing the symmetries (5.11)–(5.13) in (5.15) through $\overline {U_z \Theta }=\mathrm{e}^{-(a_{S\theta }+a_{Sx}-a_{St}+a_{Ss})}\overline {U_z \Theta }^{*}$ and $r=\mathrm{e}^{-a_{Sx}}r^{*}$ , which gives

(5.17) \begin{equation} I_\Theta =\mathrm{e}^{-(a_{S\theta }+a_{Sx}-a_{St}+2a_{Sx}+a_{Ss})} \int _{0}^\infty \overline {U_z \Theta }^{*}r^{*}\,\mathrm{d}r^{*}, \end{equation}

and in (5.16) through $\overline {U^{2}_{i}}=\mathrm{e}^{-2(a_{Sx}-a_{St})-a_{Ss}}\overline {U^{2}_{i}}^{*}$ and $\overline {U_{i}}^{*}=\mathrm{e}^{-(a_{Sx}-a_{St})-a_{Ss}}\overline {U_{i}}^{*}$ , we obtain

(5.18) \begin{equation} I_O= \mathrm{e}^{-4a_{Sx}+2a_{St}-a_{Ss}}\int _0^\infty \left [ \overline {U_z^{2}}^{*}-\frac {1}{2}\left (\overline {U_r^{2}}^{*}+\overline {U_\varphi ^{2}}^{*}\right ) +\frac {1}{2} \frac {\mathrm{d}}{\mathrm{d}z^{*}} \left (\overline {U_rU_z}^{*}r^{*}\right )\right ]r^{*}\mathrm{d}r^{*}. \end{equation}

Inferring $I_\Theta$ and $I_O$ are invariants, i.e. constants, results in the symmetry breaking $a_{S\theta }+3a_{Sx}-a_{St}+a_{Ss}=0$ and $4a_{Sx}-2a_{St}+a_{Ss}=0$ to maintain invariance. When solved for $a_{St}$ and $a_{S\theta }$ , we find

(5.19) \begin{equation} a_{St}=2a_{Sx}+\frac {1}{2}a_{Ss}, \quad a_{S\theta }=-a_{Sx}-\frac {1}{2}a_{Ss}. \end{equation}

Under the consideration of the symmetry breaking (5.19), the turbulent scaling laws can be derived by integrating the characteristic condition in (5.14)

(5.20) \begin{align} \eta &=\frac {r}{z-z_0} ,\end{align}

(5.21) \begin{align} \widetilde {\overline {U_i^n\Theta ^m}}(\eta )\mathrm{e}^{c_{i,nm}(n+m)}&= \displaystyle\frac {\overline {U_i^n\Theta ^m}(r,z)}{(z-z_0)^{-(n+m)(1+\frac {1}{2}a_{Ss}^{*})+a_{Ss}^{*}}} ,\end{align}

(5.22) \begin{align} \widetilde {\overline {\Theta ^m}}(\eta )\mathrm{e}^{c_{m}m}&=\displaystyle\frac {\overline {\Theta ^m}(r,z)}{(z-z_0)^{-m(1+\frac {1}{2}a_{Ss}^{*})+a_{Ss}^{*}}}, \end{align}

where $z_0 = -a_z /a_{Sx}$ is the virtual origin and $a_{Ss}^{*}=a_{Ss}/a_{Sx}$ . The variables marked with $\widetilde {(\,\cdot \,)}$ are the invariants i.e. the similarity variables.

In Nguyen & Oberlack (Reference Nguyen and Oberlack2024b ), we found that the statistical symmetry (5.13), giving a measure of intermittency (Wacławczyk et al. Reference Wacławczyk, Staffolani, Oberlack, Rosteck, Wilczek and Friedrich2014), allows for a possible variation of the turbulent decay behaviour of the moments in the $z$ -direction agreeing with the hypothesis of George (Reference George1989). However, this observation stands despite the present DNS data showing that, for the intermittency symmetry (5.13), we have $a_{Ss}=0$ . Nevertheless, the intermittency reappears in the reduced equations of the MPMEs through the scaling laws. Subsequently, it is shown that this also applies to passive scalars as detailed in § 5.4.

The invariant $\eta$ is received by integrating the first two terms of (5.14) and (5.21) results in the integration of the second and third terms and finally (5.22) emerges after integrating the second and the remaining terms in (5.14). The integration constants $c_{i,nm}$ and $c_{m}$ emerge from the integration and move into the exponent after exponentiation. From (5.21), we see not only the turbulent scaling law of the mixed moments but also a generalisation of the velocity and passive scalar scaling laws. For $m=0$ , the velocity scaling laws are received as derived in Nguyen & Oberlack (Reference Nguyen and Oberlack2024b ), while for $n=0$ we obtain the passive scalar scaling law (5.22).

5.3. Validation of the scaling laws

Just as in Nguyen & Oberlack (Reference Nguyen and Oberlack2024b ), $a^\ast _{Ss}$ generates a one-parameter family of scaling laws which, as conjectured by George (Reference George1989), may be induced by a variation of the inflow condition. However, the purely hydrodynamic turbulence DNS in Nguyen & Oberlack (Reference Nguyen and Oberlack2024b ) has already shown that $a^\ast _{Ss}=0$ and we find this result confirmed by the coupling of the NSEs (2.1) and (2.2) with the scalar equation (2.4). Furthermore, it can be inferred that $a^\ast _{Ss}$ does not differ for the passive scalar scaling law compared with the velocity scaling law. As $a_{Ss}^{*}=0$ is dependent on the velocity field and also appears in the mixed velocity–scalar moments, the passive scalar cannot have an effect on $a_{Ss}^{*}$ . Physically, this makes sense, since a passive scalar has no effect on the velocity field.

The validation of the scaling laws (5.20)–(5.22) is first conducted on the centreline, where $r=0$ . To facilitate a comparison with the DNS results, the scaling laws presented in (5.21) and (5.22) are expressed as

(5.23) \begin{equation} \overline {\Theta ^m}(r=0,z)=\frac {\widetilde {\overline {\Theta ^m}}(\eta =0)\alpha _{\Theta ,m}\mathrm{e}^{c_m m}}{(z-z_0)^{m}}, \end{equation}

where $\alpha _{\Theta ,m}$ is extracted so that $\widetilde {\overline {\Theta ^m}}(\eta =0)=1$ and

(5.24) \begin{equation} \overline {U_{i}^n\Theta ^m}(r=0,z)=\frac {\widetilde {\overline {U_{i}^n\Theta ^m}}(\eta =0)\alpha _{i\Theta ,nm}\mathrm{e}^{c_{i,nm} (n+m)}}{(z-z_0)^{n+m}}, \end{equation}

where we again extract $\alpha _{i\Theta ,nm}$ so that the invariant $\widetilde {\overline {U_{i}^n\Theta ^m}}(\eta =0)=1$ . These formulations allow for a direct comparison between the scaling laws and the DNS data at the centreline.

Analysing the DNS data shown in figure 13, we conclude that for $n=0$ we obtain $c_m=c=1.755$ and $\alpha _{\Theta ,m}=\alpha _{\Theta }=0.722$ for the purely passive scalar scaling (5.22). Note that both constants are independent of the moment order $m$ . This implies that $c$ is an invariant with respect to $m$ . So far, this cannot be derived directly from symmetry theory, but support is essentially from the DNS data.

Extending this to the scaling law of the mixed moments for $i=z$ in (5.24), it can be observed that the constants in figure 13 span a plane. For $i=r,\varphi$ , a ‘discrete’ plane is formed because the moments are zero for odd $n$ on the centreline in figure 14. This results, as explained in Nguyen & Oberlack (Reference Nguyen and Oberlack2024b ), from the fact that the pure velocity scaling law for $i=r,\varphi$ ( $m=0$ ) is described by a centred Gaussian process. Since the analysis is performed on the centreline, the scaling law (5.24) gives the same values for $i=r,\varphi$ due to symmetry. Therefore, only $i=r$ will be considered in the following. The prefactors of the scaling law in (5.24) can be characterised by an arithmetic weighing of the moments with

(5.25) \begin{equation} \alpha _{i\Theta ,nm}=\frac {n \alpha _{i,n}+m\alpha _\Theta }{n+m}, \end{equation}

and

(5.26) \begin{equation} c_{i,nm}=\frac {n\beta _i+m c}{n+m}, \end{equation}

where $\alpha _{z,n}=\alpha _z=0.673$ , $\beta _z=2.123$ , $\alpha _{r,n}=\sqrt {2\mathrm{e}}\ [1.81(n-1) ]^{{n}/{2}}$ and $\beta _r=-0.5$ are taken from Nguyen & Oberlack (Reference Nguyen and Oberlack2024b ).

Figure 13. The exponential prefactor $\alpha _{z\Theta ,nm}\mathrm{e}^{c_{z,nm} (n+m)}$ () from (5.24) determined with the DNS data for mixed moments up to $n+m=6$ (), for $\Theta$ moments up to $m=10$ () and $U_z$ moments up to $n=10$ () is shown. Additionally, (5.24) is highlighted for $m=0$ () and $n=0$ ().

Figure 14. The exponential prefactor $\alpha _{r\Theta ,nm}\mathrm{e}^{c_{r,nm} (n+m)}$ () from (5.24) determined with the DNS data for mixed moments up to $n+m=6$ (), for $\Theta$ moments up to $m=10$ () and $U_r$ moments up to $n=10$ () is shown. Additionally, (5.24) is highlighted for $m=0$ () and $n=0$ ().

5.4. Gaussian behaviour of the higher passive scalar moments

Interestingly, the passive scalar moments behave similarly to the axial velocity moments. The higher instantaneous $H$ -moments of the passive scalar exhibit a Gaussian-like behaviour in $\eta$ , similar to what is observed for the axial velocity moments. By analysing the scaled higher moments of the passive scalar according to (5.22), we find that a Gaussian function provides a highly accurate representation for these moments. The Gaussian function is expressed as

(5.27) \begin{equation} \frac {\overline {\Theta ^m}}{\overline {\Theta ^m_c}}=\widetilde {\overline {\Theta ^m}}(\eta )=\mathrm{e}^{-\gamma _m \eta ^{2}}, \end{equation}

where $\eta$ represents the dimensionless radial coordinate as defined in (5.20). This Gaussian function captures the characteristic distribution of the passive scalar, showing a rapid decrease in amplitude as the distance from the central region increases. One consequence of representing moments in its $H$ form is the dominance of the mean over the fluctuations  $R$ . While it may be tempting to assume $\overline {\Theta ^m}\approx \overline {\Theta }^m$ , implying that $\gamma _m$ is linear and thus trivial due to the strong dominance of the mean passive scalar concentration $\overline {\Theta }$ over the fluctuations, our analysis demonstrates that $\gamma _m$ does not behave linearly for larger moments, as is already implied by figure 3 for $R_{\Theta \Theta }$ . Employing (5.27) in (1.4), we find that only for $\gamma _2=2 \gamma _1$ could the Gaussian curve be factored out. For this reason, the higher moments of the fluctuations do not have a simple relationship, e.g. $R_{\Theta \Theta }$ has an off-axis peak (see figure 3) instead of behaving like a Gaussian curve. In figure 15, we observe the self-preservation of six selected moments of the passive scalar, ranging up to order ten, at multiple stages along the spatial domain ( $z/D=25{-}55$ ). Except for the central region, there is a collapse of the data when normalised using the scaling law described in (5.23). This collapse suggests a universal scaling behaviour for the moments of the passive scalar, reinforcing the idea of self-preservation. Additionally, the Gaussian-like curves become progressively narrower as the moment order increases. The DNS data up to the tenth axial moment are presented in semi-log scaling in figure 16, demonstrating a parabolic trend, suggesting a strong agreement with (5.27), although the specific value of $\gamma _m$ remains undetermined. The deviation of the DNS data from the parabolic curve at the edge is likely attributed to external intermittency, as seen in figure 10. Further, $\overline {\Theta }^m$ is shown for $m=2$ to $m=10$ , i.e. $\gamma _m=m\gamma _1$ , where the discrepancy discloses that $\gamma _m$ is indeed nonlinear. A comparison between figures 15 and 16 reveals that the turbulent scaling laws hold true within the fully turbulent regime. It is important to note that turbulent scaling laws are usually only valid within the confines of a fully turbulent regime and become less applicable as we move into the ambient environment. The estimation of $\gamma _m$ in $m$ for the passive scalar moments can be obtained from figure 17 through a straightforward nonlinear fitting procedure. The resulting fit yields

(5.28) \begin{equation} \gamma _m=-1.27m^{2} + 37.52m + 27.34, \end{equation}

with a coefficient of determination of $R^{2}=0.99$ , clearly showing that $\overline {\Theta ^m}\not \approx \overline {\Theta }^m$ and therefore strong non-Gaussian fluctuations and hence higher moments in general cannot be ignored. In the present DNS $\gamma _1=73.59$ agrees well with the experiment by Birch et al. (Reference Birch, Brown, Dodson and Thomas1978) at $\gamma _1=73.6$ . Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001) found $\gamma _1=59.1$ , which means a wider passive scalar profile. It is interesting to note that Birch et al. (Reference Birch, Brown, Dodson and Thomas1978) utilised a fully developed pipe profile as an inlet, which is similar to the fully developed turbulent pipe profile that is presently used. On the other hand, Lubbers et al. (Reference Lubbers, Brethouwer and Boersma2001) employed a top-hat velocity profile, implying that inlet conditions affect the width of the passive scalar profile.

Figure 15. The radial profiles of the $m$ th axial moment normalised with the scalings in (5.22) at different distances from the orifice: $z/D=25$ (), $35$ (), $45$ (), $55$ (). The black solid lines indicate the Gaussian from (5.27) using $\lambda _m$ from (5.28).

Figure 16. The radial profiles () of the first (top) up to the tenth (bottom) axial moment and the corresponding Gaussian () from (5.27) shown in a semi-logarithmic plot at $z/D=45$ . Here, $\overline {\Theta }^m$ () is depicted from $m=2$ to $m=10$ .

Figure 17. Constants $\gamma _{m}$ from (5.27) are shown for each moment up to order $n=10$ determined by fitting to the DNS, yielding the following fit: $\gamma _m=-1.27m^{2} + 37.52m + 27.34$ .

To get a first hint towards the understanding of the emergence of (5.27) as a function of the radial coordinate $\eta$ , the MPSE (5.1) is transformed into cylindrical coordinates. Subsequently, the scaling laws (5.20)–(5.22) are introduced into the cylindrical-transformed MPSE, yielding

(5.29) \begin{align} \frac {\partial \tilde {H}_{\Theta _{\{m\}}}}{\partial \tilde {t}}&+\sum _{l=1}^{m}\left [\frac {1}{\eta _{(l)}}\frac {\partial }{\partial \eta _{{(l)}}}\left [\eta _{(l)} \tilde {H}_{\Theta _{\{m+1\}}[\Theta _{(m+1)}\mapsto r_{(l)}]}\left [\boldsymbol{\tilde {x}}_{(m+1)}\mapsto \boldsymbol{\tilde {x}}_{(l)}\right ]\right ] \right . \nonumber \\[6pt] &\left . -\tilde {x}_{k_{(l)}}\frac {\partial }{\partial \tilde {x}_{k_{(l)}}}\left [\tilde {H}_{\Theta _{\{m+1\}}[\Theta _{(m+1)}\mapsto z_{(1)}]}\left [\boldsymbol{\tilde {x}}_{(m+1)}\mapsto \boldsymbol{\tilde {x}}_{(1)}\right ]\right ] \right ]\nonumber \\[6pt] &+\sum _{l=2}^{m}\left [\frac {1}{\eta _{(l)}}\frac {\partial }{\partial \tilde {\varphi }_{{(l)}}} \left [\tilde {H}_{\Theta _{\{m+1\}}[\Theta _{(m+1)}\mapsto \varphi _{(l)}]}\left [\boldsymbol{\tilde {x}}_{(m+1)}\mapsto \boldsymbol{\tilde {x}}_{(l)}\right ]\right ] \right . \nonumber \\[6pt] &\left .+\frac {\partial }{\partial \tilde {z}_{{(l)}}} \left [\tilde {H}_{\Theta _{\{m+1\}}[\Theta _{(m+1)}\mapsto z_{(l)}]}\left [\boldsymbol{\tilde {x}}_{(m+1)}\mapsto \boldsymbol{\tilde {x}}_{(l)}\right ]\right ] \right ]\nonumber \\[6pt] &-\left [(m+1)\left (1+\frac {1}{2}a_{Ss}^{*}\right )-a_{Ss}^{*}\right ]\tilde {H}_{\Theta _{\{m+1\}}[\Theta _{(m+1)}\mapsto z_{(1)}]}\left [\boldsymbol{\tilde {x}}_{(m+1)}\mapsto \boldsymbol{\tilde {x}}_{(1)}\right ]=0 ,\end{align}

for $m=1,\ldots , \infty$ and $\tilde {x}_{k_{(l)}} \neq \tilde {\varphi }_{(l)}$ with a reference point $\boldsymbol{\tilde {x}}_{(1)}$

(5.30) \begin{equation} \boldsymbol{\tilde {x}}_{(1)}=\left (\eta _{(1)},0,0\right )^\textrm{T} =\left (\frac {r_{(1)}}{z_{(1)}-z_0},0,0\right )^\textrm{T}, \end{equation}

and the points $\boldsymbol{\tilde {x}}_{(l)}$

(5.31) \begin{equation} \boldsymbol{\tilde {x}}_{(l)}=\left (\eta _{(l)}, \tilde {\varphi }_{(l)}, \tilde {z}_{(l)} \right )^\textrm{T}=\left (\frac {r_{(l)}}{z_{(1)}-z_0}, \tilde {\varphi }_{(l)}, \frac {z_{(l)}}{z_{(1)}-z_0}\right )^\textrm{T} ,\end{equation}

for $l=2,\ldots ,\infty$ and where $\tilde {t}$ emerges from (5.12) as

(5.32) \begin{equation} \tilde {t}=\frac {t}{(z-z_0)^{2+\frac {1}{2}a_{Ss}^{*}}}. \end{equation}

Detailed steps for obtaining (5.29) and (5.32) can be taken from Appendix C. It is important to note that $\boldsymbol{\tilde {x}}_{(l)}$ must contain $\tilde {\varphi }_{(l)}$ to ensure their unique definition. Equation (5.29) has been reduced by one independent variable through the scaling laws. After a symmetry reduction, the symmetry properties change in relation to the original equation. In particular, scaling symmetries can disappear, but various symmetries are also inherited. Presently, symmetries admitted by (5.29) include

(5.33) \begin{align} \tilde {\overline {T}}_{S}: \quad &\tilde {t}^{*}=\mathrm{e}^{\tilde {a}_{St}}\tilde {t}, \quad \eta _{(\cdot )} ^{*}=\mathrm{e}^{\tilde {a}_{S}}\eta _{(\cdot )}, \quad \tilde {\varphi }_{(\cdot )}^{*}=\tilde {\varphi }_{(\cdot )}, \quad \tilde {z}_{(\cdot )}^{*}=\tilde {z}_{(\cdot )}, \nonumber \\[6pt] &\tilde {H}_{\Theta _{\{m\}}}^{*}=\mathrm{e}^{m\tilde {a}_{S\theta }-\tilde {a}_{S}} \tilde {H}_{\Theta _{\{m\}}}, \quad \tilde {H}_{r_{\{1\}}\Theta _{\{m\}}}^{*}=\mathrm{e}^{m\tilde {a}_{S\theta }-\tilde {a}_{St}} \tilde {H}_{r_{\{1\}}\Theta _{\{m\}}}, \nonumber \\[6pt] &\tilde {H}_{\varphi _{\{1\}}\Theta _{\{m\}}}^{*}=\mathrm{e}^{m\tilde {a}_{S\theta }-\tilde {a}_{St}} \tilde {H}_{\varphi _{\{1\}}\Theta _{\{m\}}}, \quad \tilde {H}_{z_{\{1\}}\Theta _{\{m\}}}^{*}=\mathrm{e}^{m\tilde {a}_{S\theta }-\tilde {a}_{S}-\tilde {a}_{St}} \tilde {H}_{z_{\{1\}}\Theta _{\{m\}}}, \end{align}

with $\tilde {H}_{\Theta _{\{m+1\}}[\Theta _{(m+1)}\mapsto i_{(l)}]}=\tilde {H}_{i_{\{1\}}\Theta _{\{m\}}}$ and $i_{(l)}=r_{(l)}, \varphi _{(l)}, z_{(l)}$ . However, caution is necessary regarding the scaling symmetries given here, because although inserting (5.33) into (5.29) seems to prove the symmetry properties, the symmetries (5.33) are not necessarily symmetries of the complete system (5.4) in reduced form. Therefore, we are dealing with a kind of partial invariance, which is not to be confused with the mathematical concept of partially invariant solutions (see e.g. Meleshko Reference Meleshko2005). The statistical symmetry ( $\tilde {a}_{Ss}$ ) is inherited from (5.13)

(5.34) \begin{equation} \tilde {\overline {T}}_{Ss}: \boldsymbol{\tilde {x}}_{(\cdot )} ^{*}=\boldsymbol{\tilde {x}}_{(\cdot )}, \tilde {H}_{\Theta _{\{m\}}}^{*}=\mathrm{e}^{\tilde {a}_{Ss}} \tilde {H}_{\Theta _{\{m\}}}, \tilde {H}_{i_{\{1\}}\Theta _{\{m\}}}^{*}= \mathrm{e}^{\tilde {a}_{Ss}} \tilde {H}_{i_{\{1\}}\Theta _{\{m\}}}, \end{equation}

and the superposition symmetry of linear systems is presently given by

(5.35) \begin{align} \tilde {T}_{S}:\quad & t^{*}=t, \quad \boldsymbol{\tilde {x}}_{(\cdot )} ^{*}=\boldsymbol{\tilde {x}}_{(\cdot )}, \quad \tilde {H}_{\Theta _{\{m\}}}^{*}=\tilde {H}_{\Theta _{\{m\}}}+\tilde {H}^{\prime}_{\Theta _{\{m\}}},\nonumber \\[6pt] & \tilde {H}_{i_{\{1\}}\Theta _{\{m\}}}^{*}= \tilde {H}_{i_{\{1\}}\Theta _{\{m\}}}+\tilde {H}^{\prime}_{i_{\{1\}}\Theta _{\{m\}}}, \end{align}

where the prime quantities are any additional independent solution of (5.29). Here, the superposition symmetry (5.35) plays a crucial role, as will be explained later. With the symmetries (5.33)–(5.34), a characteristic condition is, once again, derived

(5.36) \begin{equation} \frac {\mathrm{d}\eta }{\tilde {a}_{S} \eta }=\frac {\mathrm{d}\tilde {H}_{\Theta _{\{m\}}}}{(m\tilde {a}_{S\theta } -\tilde {a}_{S}+\tilde {a}_{Ss})\tilde {H}_{\Theta _{\{m\}}}}=\ldots , \end{equation}

where non-relevant terms have been omitted. Integrating the characteristic condition (5.36) leads to

(5.37) \begin{equation} \tilde {H}_{\Theta _{\{m\}}}=C_{m}\eta ^q, \end{equation}

where $C_{m}$ is the integration constant and

(5.38) \begin{equation} q=m\frac {\tilde {a}_{S\theta }}{\tilde {a}_{S}}+\frac {\tilde {a}_{Ss}}{\tilde {a}_{S}}-1. \end{equation}

For the connection of (5.27) to (5.37), the superposition symmetry (5.35) reveals its importance, allowing for a superposition of (5.37). The Taylor series of (5.27)

(5.39) \begin{equation} \mathrm{e}^{-\gamma _m \eta ^{2}}=1-\gamma _m\eta ^{2}+\frac {1}{2}\left (\gamma _m\eta ^{2}\right )^{2}-\frac {1}{6}\left (\gamma _m\eta ^{2}\right )^3+\ldots , \end{equation}

shows that (5.37) are the basis functions of (5.27). Therefore, the sum of the basis functions (5.37) can be written as

(5.40) \begin{equation} \tilde {H}_{\Theta _{\{m\}}}=\sum _{i=0}^{\infty } C_{mi}\eta ^{q_i}, \end{equation}

where

(5.41) \begin{equation} C_{mi}=\frac {(-\gamma _m)^i}{i!}, \quad q_i=2i, \end{equation}

showing the importance of the intermittency symmetry and the superposition symmetry, which are both only admitted by the equations in the $H$ -formulation and not by the Euler or NSEs.

5.5. Gaussian behaviour of the higher mixed moments

We have gathered data for higher mixed moments up to $n+m=6$ , aiming to demonstrate their composition as a combination of axial velocity and passive scalar moments, similar to § 5.3.

Figure 18. The radial profiles (blue, dashed) of $m=1$ (top) up to the $n+m=6$ (bottom) axial mixed moment $\overline {U_z^n\Theta ^m}$ and the corresponding Gaussian (black) from (5.27) shown in a semi-logarithmic plot at $z/D=45$ .

As illustrated in figure 18, we observe that the mixed moments also exhibit Gaussian behaviour at high precision. This characteristic is attributed to both the higher axial velocity and passive scalar moments displaying similar trends. In our recent study (Nguyen & Oberlack Reference Nguyen and Oberlack2024b ) for the axial velocity moments, we established

(5.42) \begin{equation} \frac {\overline {U_z^n}}{\overline {U_{z,c}^n}}=\widetilde {\overline {U_z^n}}(\eta )=\mathrm{e}^{-\gamma _n \eta ^{2}}, \end{equation}

where the Gaussian coefficient $\gamma _n$ was approximated by

(5.43) \begin{equation} \gamma _n=-1.54n^{2}+61.95n+27.10. \end{equation}

From (5.27) and (5.42) we may invoke

(5.44) \begin{equation} \widetilde {\overline {U_z^n \Theta ^m}}(\eta )=\mathrm{e}^{-\gamma _{nm}\eta ^{2}}. \end{equation}

The subsequent approximation is a nonlinear $n$ - $m$ surface fit to the DNS data, as shown in figure 18, and is

(5.45) \begin{equation} \gamma _{nm}=-1.51n^{2}-1.29m^{2}-0.03nm+61.2n+37.34m+30.75, \end{equation}

presented in figure 19 with a coefficient of determination of $R^{2}=0.99$ . The Gaussian coefficients $\gamma _n$ and $\gamma _m$ in (5.43) and (5.28) respectively may be compared with $\gamma _{nm}$ in (5.45) by tentatively assuming $\widetilde {\overline {U_z^n}}\,\widetilde {\overline {\Theta ^m}}\approx \widetilde {\overline {U_z^n \Theta ^m}}$ , which implies $\gamma _{n}+\gamma _{m}\approx \gamma _{nm}$ . Comparing the Gaussian coefficients of (5.28), (5.43) and (5.45) reveals an approximate halving of the constant terms and the emergence of a small cross-product term, while the remaining terms only differ by $2\, \%$ . This observation suggests that the moments are coupled and demonstrate some degree of cross-correlation. Nevertheless, the approximation (5.45) reduces approximately to the special forms (5.28) and (5.43) in the limit cases $n=0$ and $m=0$ .

Figure 19. Constants $\gamma _{nm}$ from (5.44) are shown for pure moments up to order $n,m=10$ and mixed moments up to order $n+m=6$ determined by fitting to the DNS yielding the following fit: $\gamma _{nm}=-1.51n^{2}-1.29m^{2}-0.03nm+61.2n+37.34m+30.75$ .