2.1 Fluid phase

We consider a statistically stationary, homogeneous and isotropic turbulent flow, governed by the incompressible Navier–Stokes equations. The turbulent velocity field advects the fluid temperature field (assumed a passive scalar), together with the inertial particles. In this study, we account for two-way thermal coupling between the fluid and particles, but only one-way momentum coupling. Therefore, the governing equations for the fluid phase are

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\frac{1}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\boldsymbol{f}, & \displaystyle\end{eqnarray}$$

(2.1c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}T+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T=\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FB}^{2}T-C_{T}+f_{T}. & \displaystyle\end{eqnarray}$$

$\boldsymbol{u}(\boldsymbol{x},t)$

$p(\boldsymbol{x},t)$

$\unicode[STIX]{x1D70C}_{0}$

$\unicode[STIX]{x1D708}$

$T(\boldsymbol{x},t)$

$\unicode[STIX]{x1D705}$

$\unicode[STIX]{x1D708}$

$\unicode[STIX]{x1D705}$

$Pr\equiv \unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$

$Pr=1$

$\boldsymbol{f}$

$f_{T}$

$C_{T}$

$\boldsymbol{x}$

When the forcing is confined to sufficiently large scales, it is assumed that the details of the forcing do not influence the small-scale dynamics. Previous experimental evidence seems to confirm this (Sreenivasan Reference Sreenivasan1996), leading to a universal behaviour of the small scales. However, recent studies (Gotoh & Watanabe Reference Gotoh and Watanabe2015) pointed out that this hypothesis of universality is partially violated by the advected scalar fields, whose inertial range statistics exhibit sensitivity to the details of the imposed forcing. Since we aim to characterize temperature and temperature gradient fluctuations in the dissipation range for different inertias of the suspended particles, we employ a forcing that imposes the same total dissipation rate for all the simulations. This produces results which can be meaningfully compared for different parameters of the suspended particles, since the response of the system to the same injected thermal power can be examined. Therefore, we employ the large-scale forcing (Kumar, Schumacher & Shaw Reference Kumar, Schumacher and Shaw2013; Kumar et al. Reference Kumar, Schumacher and Shaw2014),

(2.2a,b ) $$\begin{eqnarray}\hat{\boldsymbol{f}}(\boldsymbol{k},t)=\unicode[STIX]{x1D700}\frac{\hat{\boldsymbol{u}}(\boldsymbol{k},t)}{\displaystyle \mathop{\sum }_{\boldsymbol{k}_{f}\in {\mathcal{K}}_{f}}\Vert \hat{\boldsymbol{u}}(\boldsymbol{k}_{f},t)\Vert 2}\unicode[STIX]{x1D6FF}_{\boldsymbol{k},\boldsymbol{k}_{f}},\quad \hat{f}_{T}(\boldsymbol{k},t)=\unicode[STIX]{x1D712}\frac{\hat{T}(\boldsymbol{k},t)}{\displaystyle \mathop{\sum }_{\boldsymbol{k}_{f}\in {\mathcal{K}}_{f}}|\hat{T}(\boldsymbol{k}_{f},t)|2}\unicode[STIX]{x1D6FF}_{\boldsymbol{k},\boldsymbol{k}_{f}},\end{eqnarray}$$

in the wavenumber space. A hat indicates the Fourier transform and $\boldsymbol{k}_{f}$ is the wavenumber which here belongs to the set of forced wavenumbers, ${\mathcal{K}}_{f}=\{\boldsymbol{k}_{f}:\,\Vert \boldsymbol{k}_{f}\Vert =k_{f}\}$ ; $\unicode[STIX]{x1D700}$ and $\unicode[STIX]{x1D712}$ are the imposed dissipation rates of velocity and temperature variance, respectively. This employed forcing scheme thus allows us to control the overall dissipation rate and, therefore, to control the Stokes number.

The values of the parameters relative to the fluid flow employed in the simulations are given in table 1. Time-averaged energy and temperature spectra, in the absence of particle thermal feedback, are shown in figure 1(a).

Table 1. Flow parameters for the numerical simulations in this study. Dimensional parameters are non-dimensionalized into arbitrary code units. The characteristic parameters of the fluid flow are defined from its energy spectrum $E(k)\equiv \int _{\Vert \boldsymbol{k}\Vert =k}\Vert \hat{\boldsymbol{u}}(\boldsymbol{k})\Vert ^{2}\,\text{d}\boldsymbol{k}/2$ . The dissipation rate of turbulent kinetic energy is $\unicode[STIX]{x1D700}\equiv 2\unicode[STIX]{x1D708}\int k^{2}E(k)\,\text{d}k$ . The Kolmogorov length $\unicode[STIX]{x1D702}\equiv (\unicode[STIX]{x1D708}^{3}/\unicode[STIX]{x1D700})^{1/4}$ , time scale $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\equiv (\unicode[STIX]{x1D708}/\unicode[STIX]{x1D700})^{1/2}$ and velocity scale $u_{\unicode[STIX]{x1D702}}\equiv \unicode[STIX]{x1D702}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ . The Taylor micro-scale is $\unicode[STIX]{x1D706}\equiv u^{\prime }/\sqrt{\langle |\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\rangle }$ . The root mean square velocity is $u^{\prime }\equiv \sqrt{(2/3)\int E(k)\,\text{d}k}$ and the integral length scale $\ell \equiv \unicode[STIX]{x03C0}/(2u^{\prime 2})\int E(k)/k\,\text{d}k$ . Similarly, the spectrum, root mean square value and dissipation rate of the scalar field are $E_{T}(k)\equiv \int _{\Vert \boldsymbol{k}\Vert =k}|\hat{T}(\boldsymbol{k})|^{2}\,\text{d}\boldsymbol{k}/2$ , $T^{\prime }\equiv \sqrt{2\int E_{T}(k)\,\text{d}k}$ , $\unicode[STIX]{x1D712}_{f}\equiv 2\unicode[STIX]{x1D705}\int k^{2}E_{T}(k)\,\text{d}k$ . The small-scale temperature is determined by the viscosity and dissipation rate: $T_{\unicode[STIX]{x1D702}}\equiv \sqrt{\unicode[STIX]{x1D712}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}}$ . Since the Prandtl number is unitary the small scales of the scalar and the velocity field are of the same order.

2.2 Particle phase

We consider rigid, point-like particles which are heavy with respect to the fluid, and small with respect to any scale of the flow. In particular, the particle density $\unicode[STIX]{x1D70C}_{p}$ is much larger than the fluid density $\unicode[STIX]{x1D70C}_{p}\gg \unicode[STIX]{x1D70C}_{0}$ , and the particle radius $r_{p}$ is much smaller than the Kolmogorov length scale $r_{p}\ll \unicode[STIX]{x1D702}$ . With these assumptions (and neglecting gravity) the particle acceleration is described by the Stokes drag law. Analogously, the rate of change of the particle temperature is described by Newton’s law for the heat conduction

(2.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{x}_{p}}{\text{d}t}\equiv \boldsymbol{v}_{p}, & \displaystyle\end{eqnarray}$$

(2.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{v}_{p}}{\text{d}t}=\frac{\boldsymbol{u}(\boldsymbol{x}_{p},t)-\boldsymbol{v}_{p}}{\unicode[STIX]{x1D70F}_{p}}, & \displaystyle\end{eqnarray}$$

(2.3c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D703}_{p}}{\text{d}t}=\frac{T(\boldsymbol{x}_{p},t)-\unicode[STIX]{x1D703}_{p}}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}}. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D70F}_{p}\equiv 2\unicode[STIX]{x1D70C}_{p}r_{p}^{2}/(9\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D708})$

$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}\equiv \unicode[STIX]{x1D70C}_{p}c_{p}r_{p}^{2}/(3\unicode[STIX]{x1D70C}_{0}c_{0}\unicode[STIX]{x1D705})$

$c_{p}$

$c_{0}$

$St\equiv \unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$

$St_{\unicode[STIX]{x1D703}}\equiv \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$

$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$

We consider nine values of $St_{\unicode[STIX]{x1D703}}$ and three values of $St$ in order to explore the behaviour of the system over a range of parameter values. Since we are accounting for thermal coupling, each combination of $St_{\unicode[STIX]{x1D703}}$ and $St$ must be simulated separately, and when combined with the large number of particles in the flow domain, the set of simulations require considerable computational resources. Therefore, in the present study we restrict attention to $Re_{\unicode[STIX]{x1D706}}=88$ , but future explorations should consider larger $Re_{\unicode[STIX]{x1D706}}$ in order to explore the behaviour when there exists a well-defined inertial range in the flow.

In order to obtain deeper insight into the role of the two-way thermal coupling, we perform simulations with (denoted by S1) and without (denoted by S2) the thermal coupling. The particle parameters employed in the simulations are given in table 2. The second-order longitudinal structure functions of the particle velocity are shown in figure 1(b) for all the Stokes numbers investigated.

2.3 Thermal coupling

In the two-way thermal coupling regime, the thermal energy contained in the fluid is finite with respect to the thermal energy of the particles, therefore, when heat flows from the fluid to the particle the fluid loses thermal energy at the particle position. Due to the point-mass approximation, the feedback from the particles on the fluid temperature field is a superposition of Dirac delta functions, centred on the particles. Hence the coupling term in (2.1c ) is given by

(2.4) $$\begin{eqnarray}C_{T}(\boldsymbol{x},t)=\frac{4}{3}\unicode[STIX]{x03C0}\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}_{0}}\frac{c_{p}}{c_{0}}r_{p}^{3}\mathop{\sum }_{p=1}^{N_{P}}\frac{\text{d}\unicode[STIX]{x1D703}_{p}}{\text{d}t}\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{p}).\end{eqnarray}$$

Table 2. Particle parameters in dimensionless code units. The Stokes number is $St\equiv \unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ and the thermal Stokes number $St_{\unicode[STIX]{x1D703}}\equiv \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ ; the particle response times are defined in the text. In the simulations, $St_{\unicode[STIX]{x1D703}}$ is varied by varying the particle heat capacity. The different combinations of $St$ and $St_{\unicode[STIX]{x1D703}}$ are simulated including the two-way thermal coupling (simulations S1) and neglecting it (simulations S2).

2.4 Validity and limitations of the model

The physical model in §§ 2.1–2.3, is normally referred to as point-particle model. The two-way thermal coupling regime is considered, that is, the particles can affect the fluid temperature field while the direct particle–particle thermal interaction is neglected. Previous estimations (Elghobashi Reference Elghobashi1991) have shown that particle–particle interactions become relevant at average volume fractions $\unicode[STIX]{x1D719}$ exceeding $10^{-3}$ . In this work, the volume fraction lies in the two-way coupling regime and the average distance between particles exceeds the particle diameter by an order of magnitude. In our simulations $\unicode[STIX]{x1D719}=4\times 10^{-4}$ , which is small enough to neglect particle collisions, but large enough for two-way momentum coupling between the particles and fluid to be important (e.g. Elghobashi Reference Elghobashi1994). Nevertheless, we ignore two-way momentum coupling in the present study. The motivation is that including both two-way momentum and two-way thermal coupling introduces too many competing effects that would compound a thorough understanding of the problem. We therefore employ a reductionistic approach, seeking first to understand the role of two-way thermal coupling in the absence of momentum coupling, and then in a future study will explore their combined effects. Even aside from this methodological point, the results still have physical relevance since the thermal relaxation time is often larger than the momentum relaxation time in many particle-laden flows. For example, $St_{\unicode[STIX]{x1D703}}/St$ ranges from 2 to 6 for many liquid droplets in air ( ${\approx}4$ for water droplets in air). Therefore thermal feedback can be more relevant to the thermal balance than momentum feedback on momentum balance. This is confirmed by the analysis of the effect of momentum coupling and elastic collisions on the temperature statistics, presented in appendix A. An additional set of simulations shows that, in the range of parameters considered in this work, both phenomena have a minimal effect.

The numerical solution of (2.1b,c ) and (2.3) is considered a DNS, insofar as the Kolmogorov scale is resolved (Kuerten Reference Kuerten2016), even though the details of the flow near the surface of each particle are not resolved. This point-particle simplification is formally valid for particles which are smaller than the smallest active scale in the flow, the Kolmogorov microscale (Toschi & Bodenschatz Reference Toschi and Bodenschatz2009). When the particle Reynolds number (based on the slip velocity between the particle and the local fluid) is small, the effect of the stresses on the particle can be described using a Stokes drag force (Maxey & Riley Reference Maxey and Riley1983). Under an analogous set of conditions, the heat transfer between the particle and the fluid is a diffusive process, that has a time scale $r_{p}^{2}/\unicode[STIX]{x1D705}$ . For small point-like particles, this time scale is much smaller than the Kolmogorov time scale, so that the heat transfer is a quasi-steady process which leads to the Newton-like equation for heat transfer in (2.4) (Zonta et al. Reference Zonta, Marchioli and Soldati2008; Bec et al. Reference Bec, Homann and Krstulovic2014). In this work, the ratio between the particle radius $r_{p}$ and the Kolmogorov scale $\unicode[STIX]{x1D702}$ is well below $0.1$ for $St=0.5$ and $St=1$ and approximately $0.1$ at $St=3$ . Therefore finite size effects are negligible up to $St=1$ , while there may be small errors for $St=3$ . The impact of these small errors should be quantified in a future work using a more sophisticated model that resolves the flow around the particle surface. Also the fluid continuity equation should be in principle modified, due to the volume occupied by inertial particles. However, the error introduced in the continuity equation, proportional to the rate of change of the local volume fraction, is of the same order as the error introduced by other approximations in the model which are quite small.

In the point-mass Eulerian–Lagrangian model here employed, the fluid temperature in (2.3c ) should be understood as the temperature of the carrier fluid flow without the local effect of the disturbance due to the presence of the particles (Boivin, Simonin & Squires Reference Boivin, Simonin and Squires1998). Neglecting this disturbance is justified for particles with a diameter much smaller than the Kolmogorov scale and much smaller than the grid spacing. In the case of two-way momentum coupling, the error introduced by neglecting the disturbance can be estimated to be less than $10\,\%$ for the particle parameters we are considering (Horwitz & Mani Reference Horwitz and Mani2016). A similar estimation is expected for the fluid temperature disturbance due to the particle, since the equations for the particle velocity and temperature are analogous. At the largest simulated Stokes number, $St=3$ , the same estimation indicates that an error of approximately $15\,\%$ can be introduced. Therefore, the error is larger for this case, though certainly a sub-leading effect. Some simplified, efficient ways to compute the effect of the particle disturbance are available for the two-way momentum coupling (Horwitz & Mani Reference Horwitz and Mani2016), but equivalent models for the thermal coupling problem are not well developed or tested. Given this, and the fact that for our small particles the corrections due to the disturbance terms would be sub-leading, it is justified to neglect their effect as a first approximation for understanding this complex problem.

We also emphasize that the statistics of the fluid temperature field presented in the paper refer to the carrier, resolved, temperature field, with the disturbance temperature field produced by the particles neglected. The actual (or total) temperature field is the superposition of the carrier temperature field and the disturbance induced by the particles. In the limit of large scale separation between the particle size and Kolmogorov scale and in the dilute regime, the disturbance field produced by the particle can be determined analytically and so the statistics of the actual temperature field can be reconstructed knowing the resolved carrier flow and the analytic solution for the particle temperature disturbance field. This aspect is discussed in appendix B.

Summarizing, because of the marked scale separation between the particle size and the Kolmogorov scale, $r_{p}\ll \unicode[STIX]{x1D702}$ , and low particle volume fraction, $\unicode[STIX]{x1D719}\ll 1$ , which are the conditions explored in our study, the point-particle method provides a good first approximation to the complex problem under consideration. The statistics for the temperature field reported in the paper refer to the carrier temperature field, in which the near-particle temperature changes are excluded, while the statistics of the actual temperature field can be recovered a posteriori on the basis of the modelling hypothesis. Future work should explore the problem using methods where the flow around the particle is resolved, referred to as FRDNS (fully resolved DNS) (Boivin et al. Reference Boivin, Simonin and Squires1998; Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2015). However, even with currently available high performance computational resources, FRDNS is limited to a small number of particles, making it unfeasible at present to explore the problem of interest in this work (Botto & Prosperetti Reference Botto and Prosperetti2012).

2.5 Numerical method

We perform direct numerical simulation of incompressible, statistically steady and isotropic turbulence on a tri-periodic cubic domain. Equations (2.1a ), (2.1b ) and (2.1c ) are solved by means of the pseudo-spectral Fourier method for the spatial discretization. The $3/2$ rule is employed for dealiasing (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988), so that the maximum resolved wavenumber is $k_{max}=N/2$ . The required Fourier transforms are executed in parallel using the P3DFFT library (Pekurovsky Reference Pekurovsky2012). Forcing is applied to a single scale, that is to all wavevectors satisfying $\Vert \boldsymbol{k}\Vert ^{2}=k_{f}$ , with $k_{f}=2$ , and the equations for the fluid velocity and temperature Fourier coefficients are evolved in time by means of a second-order Runge–Kutta exponential integrator (Hochbruck & Ostermann Reference Hochbruck and Ostermann2010). This method has been preferred to the standard integrating factor because of its higher accuracy and, above all, because of its consistency. Indeed, in order to obtain an accurate representation of small-scale temperature fluctuations, it is critical that the numerical solution conserves thermal energy. The same time integration scheme is used to solve particle equations (2.3a ), (2.3b ) and (2.3c ), thus providing overall consistency, since the system formed by fluid and particles is evolved in time as a whole.

The fluid velocity and temperature are interpolated at the particle position by means of fourth-order B-spline interpolation. The interpolation is implemented as a backward non-uniform fast Fourier transform (NUFFT) with B-spline basis: the fluid field is projected onto the B-spline basis in Fourier space through a deconvolution, then transformed into the physical space by means of a inverse fast Fourier transform (FFT). A convolution provides the interpolated field at particle positions (Beylkin Reference Beylkin1995). Since B-splines have a compact support in physical space and deconvolution in Fourier space reduces to a division, this provide an efficient way to obtain high-order interpolation. This guarantees smooth and accurate interpolation and its efficient implementation is suitable for pseudo-spectral methods (van Hinsberget al. Reference van Hinsberg, ten Thije Boonkkamp, Toschi and Clercx2012). The coupling term (2.4) has to be projected on the Cartesian grid used to represent the fields. This is performed by means of the forward NUFFT with B-spline basis (Beylkin Reference Beylkin1995). Briefly, the algorithm works as follows (Carbone & Iovieno Reference Carbone and Iovieno2018). The convolution of the distribution $C_{T}(\boldsymbol{x},t)$ with the B-spline polynomial basis $B(\boldsymbol{x})$ is computed in physical space, so that it can be effectively represented on the Cartesian grid

(2.5) $$\begin{eqnarray}\widetilde{C}_{T}(\boldsymbol{x},t)=C_{T}\ast B=\frac{4}{3}\unicode[STIX]{x03C0}\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}_{0}}\frac{c_{p}}{c_{0}}r_{p}^{3}\mathop{\sum }_{p=1}^{N_{P}}\frac{\text{d}\unicode[STIX]{x1D703}_{p}}{\text{d}t}B(\boldsymbol{x}-\boldsymbol{x}_{p}).\end{eqnarray}$$

The smoothed field $\widetilde{C}_{T}$ is transformed by means of a FFT giving ${\mathcal{F}}[\widetilde{C}_{T}]$ in Fourier space. Finally, the convolution with the B-spline is removed in Fourier space,

(2.6) $$\begin{eqnarray}\widehat{C}(\boldsymbol{k},t)=\frac{{\mathcal{F}}[C_{T}\ast B]}{{\mathcal{F}}[B]}=\frac{{\mathcal{F}}[\widetilde{C}_{T}]}{{\mathcal{F}}[B]}.\end{eqnarray}$$

Since the convolution is removed in Fourier space, increasing the order of B-spline provides higher accuracy, without introducing non-locality, even if the number of grid points influenced by each particle becomes larger in the preliminary convolution. A fourth-order B-spline polynomial is employed in these simulations and a test with the same forcing, same parameters, with a larger resolution ( $256^{3}$ Fourier modes) confirmed the grid independence of the results. Also, the backward and forward transformations are symmetric, that guarantees energy conservation (Sundaram & Collins Reference Sundaram and Collins1996). A detailed description and assessment of the NUFFT in the framework of particles in turbulence can be found in Carbone & Iovieno (Reference Carbone and Iovieno2018).

3.1 Thermal dissipation due to the carrier temperature gradients

Since the flow is isotropic, $\unicode[STIX]{x1D712}_{f}$ is given by

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{f}=3\unicode[STIX]{x1D705}\langle (\unicode[STIX]{x2202}_{x}T)^{2}\rangle .\end{eqnarray}$$

We consider fixed Reynolds number and $Pr=1$ , thus $\unicode[STIX]{x1D705}$ is the same in all the presented simulations, and so $\langle (\unicode[STIX]{x2202}_{x}T)^{2}\rangle$ fully characterizes $\unicode[STIX]{x1D712}_{f}$ . Moreover, given the expected structure of the field $\unicode[STIX]{x2202}_{x}T$ , it is instructive to consider its full probability density function (PDF), in addition to its moments in order to know how different regions of the flow contribute to the average dissipation rate $\unicode[STIX]{x1D712}_{f}$ .

Figure 1. (a) Three-dimensional energy spectrum of the fluid velocity field (open squares) and temperature field (open circles). The temperature field is computed without any feedback from the particles on the fluid flow (simulations S2). (b) Second-order longitudinal structure functions of the particle velocity for various Stokes numbers.

Figure 2. PDF of the carrier flow temperature gradient $\unicode[STIX]{x2202}_{x}T$ from simulations S1, at $St=3$ , for various $St_{\unicode[STIX]{x1D703}}$ (a,b) dissipation rate, $\unicode[STIX]{x1D712}_{f}$ , of the fluid temperature fluctuations, for different $St$ as a function of $St_{\unicode[STIX]{x1D703}}$ . The filled lines indicate the maximum deviations of the dissipation rate occurred in the time interval used to compute averages.

Figure 2(a) shows the normalized PDFs of $\unicode[STIX]{x2202}_{x}T$ for $St=3$ , for various $St_{\unicode[STIX]{x1D703}}$ , where the PDFs are normalized using the standard deviation of the distribution, $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x2202}_{x}T}$ . The distribution is almost symmetric and it displays elongated exponential tails. The largest temperature gradients exceed the standard deviation by an order of magnitude (Overholt & Pope Reference Overholt and Pope1996). Remarkably, the shape of the PDF shows a very weak dependence on $St$ and $St_{\unicode[STIX]{x1D703}}$ , over all the considered range of these parameters, such that the PDF shape scales with $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x2202}_{x}T}$ . Consistently, the kurtosis of the fluid temperature gradients distribution is approximately constant, which confirms that the fluid temperature gradient PDF is approximately self-similar.

The variance of the resolved fluid temperature gradient is proportional to the actual dissipation rate of the temperature fluctuation $\unicode[STIX]{x1D712}_{f}$ (the proportionality factor being $3\unicode[STIX]{x1D705}$ , the same in all the simulations). In contrast to the PDF shape, the suspended particles have a strong impact on $\unicode[STIX]{x1D712}_{f}$ , as shown in figure 2. As $St_{\unicode[STIX]{x1D703}}$ is increased, $\unicode[STIX]{x1D712}_{f}$ decreases. However, this is mainly due to the fact that as $St_{\unicode[STIX]{x1D703}}$ is increased, $\unicode[STIX]{x1D712}_{p}$ increases, and so $\unicode[STIX]{x1D712}_{f}$ must decrease since $\unicode[STIX]{x1D712}=\unicode[STIX]{x1D712}_{f}+\unicode[STIX]{x1D712}_{p}$ is fixed. The influence of the Stokes number on $\unicode[STIX]{x1D712}_{f}$ is very small in the range of parameters considered.

3.2 Thermal dissipation due to the particle dynamics

The dissipation rate due to the particles, $\unicode[STIX]{x1D712}_{p}$ , depends on the difference between the particle temperature and the fluid temperature at the particle position,

(3.4) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{p}=3\unicode[STIX]{x1D705}\unicode[STIX]{x1D719}\left\langle \left(\frac{T(\boldsymbol{x}_{p},t)-\unicode[STIX]{x1D703}_{p}}{r_{p}}\right)^{2}\right\rangle .\end{eqnarray}$$

For notational simplicity, we define $\unicode[STIX]{x1D711}_{p}\equiv \sqrt{3\unicode[STIX]{x1D719}}(T(\boldsymbol{x}_{p},t)-\unicode[STIX]{x1D703}_{p})/r_{p}$ . When $\unicode[STIX]{x1D711}_{p}$ is normalized by its standard deviation, we can relate this to the rate of change of the particle temperature using (2.3c )

(3.5) $$\begin{eqnarray}\frac{\dot{\unicode[STIX]{x1D703}}_{p}}{\unicode[STIX]{x1D70E}_{\dot{\unicode[STIX]{x1D703}}_{p}}}=\frac{\unicode[STIX]{x1D711}_{p}}{\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D711}_{p}}}.\end{eqnarray}$$

Figure 3. PDF of $\dot{\unicode[STIX]{x1D703}}_{p}$ for $St=1$  (a,b) and $St=3$  (c,d), and for various $St_{\unicode[STIX]{x1D703}}$ . Panels (a–c) are from simulations S1, in which the two-way thermal coupling is considered, while panels (b–d) are from simulations S2, in which the two-way coupling is neglected. (e) Dissipation rate $\unicode[STIX]{x1D712}_{p}$ of the temperature fluctuations due to the particles, for different $St$ as a function of $St_{\unicode[STIX]{x1D703}}$ . (f) Kurtosis of the PDF of $\dot{\unicode[STIX]{x1D703}}_{p}$ .

The normalized PDF of $\dot{\unicode[STIX]{x1D703}}_{p}$ for $St=1$ and $St=3$ , and for various $St_{\unicode[STIX]{x1D703}}$ is shown in figure 3. Figure 3(a) shows the normalized PDF of $\dot{\unicode[STIX]{x1D703}}_{p}$ , for $St=1$ for the set of simulations S1, in which the two-way thermal coupling is taken to account. Figure 3(b) shows the corresponding results for simulations S2, in which the two-way thermal coupling is neglected. The normalized PDF of $\dot{\unicode[STIX]{x1D703}}_{p}$ for $St=3$ , with and without the two-way thermal coupling, is shown in figure 3(c,d).

In contrast to the fluid temperature gradient PDFs, the shape of the PDF of $\dot{\unicode[STIX]{x1D703}}_{p}$ is not self-similar with respect to its variance. As $St_{\unicode[STIX]{x1D703}}$ is increased, the normalized PDF becomes narrower. This is due to the fact that as $St_{\unicode[STIX]{x1D703}}$ is increased, the particles respond more slowly to changes in the fluid temperature field, analogous to the ‘filtering’ effect for inertial particle velocities in turbulence (Salazar & Collins Reference Salazar and Collins2012; Ireland et al. Reference Ireland, Bragg and Collins2016a ). The PDF shapes are mildly affected by $St$ , and for larger $St_{\unicode[STIX]{x1D703}}$ , extreme fluid temperature–particle temperature differences are suppressed when the two-way thermal coupling is neglected.

The variance of $\dot{\unicode[STIX]{x1D703}}_{p}$ is proportional to the particle dissipation rate $\unicode[STIX]{x1D712}_{p}$ ,

(3.6) $$\begin{eqnarray}\langle \dot{\unicode[STIX]{x1D703}}_{p}^{2}\rangle =\left\langle \frac{(T(\boldsymbol{x}_{p},t)-\unicode[STIX]{x1D703}_{p})^{2}}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}^{2}}\right\rangle =\frac{r_{p}^{2}\unicode[STIX]{x1D712}_{p}}{3\unicode[STIX]{x1D705}\unicode[STIX]{x1D719}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}^{2}},\end{eqnarray}$$

and the results for this are shown in figure 3(e), for various $St$ and $St_{\unicode[STIX]{x1D703}}$ , and for simulations S1 and S2. The results show that as $St_{\unicode[STIX]{x1D703}}$ is increased, $\unicode[STIX]{x1D712}_{p}$ increases. This is mainly because as $St_{\unicode[STIX]{x1D703}}$ is increased, the thermal time correlation of the particle increases, and the particle temperature depends strongly on its encounter with the fluid temperature field along its trajectory history for times up to $O(\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}})$ in the past. As a result, the particle temperature can differ strongly from the local carrier flow temperature. The results also show that $\unicode[STIX]{x1D712}_{p}$ is dramatically suppressed when two-way thermal coupling is accounted for. One reason for this is that as shown earlier, two-way thermal coupling leads to a suppression in the fluid temperature gradients. As these gradients are suppressed, the fluid temperature along the particle trajectory history differs less from the local carrier flow temperature than it would have in the absence of two-way thermal coupling, and as a result $\unicode[STIX]{x1D712}_{p}$ is decreased.

The results for kurtosis of $\dot{\unicode[STIX]{x1D703}}_{p}$ , as a function of $St_{\unicode[STIX]{x1D703}}$ and for various $St$ are shown in figure 3(f). The results show that the kurtosis decreases with increasing $St_{\unicode[STIX]{x1D703}}$ . This is mainly due to the filtering effect mentioned earlier, wherein as $St_{\unicode[STIX]{x1D703}}$ is increased, the particles are less able to respond to rapid fluctuations in the fluid temperature along their trajectory. Further, the kurtosis is typically larger when the two-way thermal coupling is taken into account (simulations S1), and is maximum for $St=1$ . This is due to the particle clustering on the fronts of the fluid temperature field, as will be discussed in § 5.

Our results for the PDF of $\dot{\unicode[STIX]{x1D703}}_{p}$ and its moments differ somewhat from those in Bec et al. (Reference Bec, Homann and Krstulovic2014). This is in part due to the difference in the forcing methods employed by Bec et al. (Reference Bec, Homann and Krstulovic2014) and that in our study. The solution of (2.3c ) may be written as (Bec et al. Reference Bec, Homann and Krstulovic2014)

(3.7) $$\begin{eqnarray}\langle \dot{\unicode[STIX]{x1D703}_{p}}^{2}\rangle =\frac{1}{2\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}^{3}}\int _{0}^{\infty }\langle (\unicode[STIX]{x1D6FF}_{t}T_{p}(t))^{2}\rangle \exp \left(-\frac{t}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}}\right)\,\text{d}t,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{t}T_{p}(t)\equiv T(\boldsymbol{x}_{p}(t),t)-T(\boldsymbol{x}_{p}(0),0)$ . In the regime $St_{\unicode[STIX]{x1D703}}\ll 1$ , the exponential in (3.7) decays very fast in time so that the main contribution to the integral comes from $\unicode[STIX]{x1D6FF}_{t}T_{p}$ for infinitesimal $t$ , with $\unicode[STIX]{x1D6FF}_{t}T_{p}\sim t^{n}$ for $t\rightarrow 0$ . Substituting $\unicode[STIX]{x1D6FF}_{t}T_{p}\sim t^{n}$ into (3.7) we obtain the leading-order behaviour

(3.8) $$\begin{eqnarray}\langle \dot{\unicode[STIX]{x1D703}_{p}}^{2}\rangle \sim \frac{1}{2\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}^{3}}\int _{0}^{\infty }t^{2n}\exp \left(-\frac{t}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}}\right)\,\text{d}t\sim St_{\unicode[STIX]{x1D703}}^{2n-2},\quad St_{\unicode[STIX]{x1D703}}\ll 1.\end{eqnarray}$$

Bec et al. (Reference Bec, Homann and Krstulovic2014) used a white in time forcing for the fluid scalar field, giving $n=1/2$ , and yielding $\langle \dot{\unicode[STIX]{x1D703}_{p}}^{2}\rangle \sim St_{\unicode[STIX]{x1D703}}^{-1}$ for $St_{\unicode[STIX]{x1D703}}\ll 1$ . However, the forcing scheme that we have employed generates a field $T(\boldsymbol{x},t)$ that evolves smoothly in time, so $n=1$ and $\langle \dot{\unicode[STIX]{x1D703}_{p}}^{2}\rangle \sim$ constant for $St_{\unicode[STIX]{x1D703}}\ll 1$ .

For $St_{\unicode[STIX]{x1D703}}\gg 1$ , the integral in (3.7) is dominated by uncorrelated temperature increments, $\unicode[STIX]{x1D6FF}_{t}T\sim t^{0}$ , such that $\langle \dot{\unicode[STIX]{x1D703}_{p}}^{2}\rangle \sim St_{\unicode[STIX]{x1D703}}^{-2}$ . The comparison between figure 4(a) and figure 5 of Bec et al. (Reference Bec, Homann and Krstulovic2014) highlights the different asymptotic behaviour of $\unicode[STIX]{x1D70E}_{\dot{\unicode[STIX]{x1D703}_{p}}}^{2}\equiv \langle \dot{\unicode[STIX]{x1D703}_{p}}^{2}\rangle$ for $St_{\unicode[STIX]{x1D703}}\ll 1$ , but the same behaviour $\langle \dot{\unicode[STIX]{x1D703}_{p}}^{2}\rangle \sim St_{\unicode[STIX]{x1D703}}^{-2}$ for $St_{\unicode[STIX]{x1D703}}\gg 1$ . Further, as expected, our DNS data approach these asymptotic regimes for both the cases with and without two-way thermal coupling.

Figure 4. (a) Variance of the particle temperature rate of change as a function of the thermal Stokes number for different Stokes numbers. The dotted lines represent the expected asymptotic behaviour for $St_{\unicode[STIX]{x1D703}}\ll 1$ and $St_{\unicode[STIX]{x1D703}}\gg 1$ . (b) Normalized PDF of the particle temperature rate of change, $\dot{\unicode[STIX]{x1D703}_{p}}$ at $St_{\unicode[STIX]{x1D703}}=1$ for various Stokes numbers. The dotted line shows a Gaussian PDF for reference. Results obtained neglecting the particle thermal feedback.

Another difference is that in the results of Bec et al. (Reference Bec, Homann and Krstulovic2014), the tails of the PDFs of $\dot{\unicode[STIX]{x1D703}}_{p}$ for $St_{\unicode[STIX]{x1D703}}=1$ become heavier as $St$ is increased, whereas our results in figure 3 show that while the kurtosis of these PDFs increases from $St=0.5$ to $St=1$ , it then decreases from $St=1$ to $St=3$ . In order to examine this further, we performed simulations (without two-way thermal coupling) for $St_{\unicode[STIX]{x1D703}}=1$ and $St\leqslant 0.4$ . The results are shown in figure 4(b), and in this regime we do in fact observe that the tails of the PDFs of $\dot{\unicode[STIX]{x1D703}}_{p}$ become increasingly wider as $St$ is increased. Taken together with the results in figure 3, this implies that in our simulations, the tails of the PDFs of $\dot{\unicode[STIX]{x1D703}}_{p}$ become increasingly wider as $St$ is increased until $St\approx 1$ , where this behaviour then saturates, and upon further increase of $St$ the tails start to narrow. This non-monotonic behaviour is due to the particle clustering in the fronts of the temperature field, which is strongest for $St\approx 1$ (see § 5). While the results in Bec et al. (Reference Bec, Homann and Krstulovic2014) over the range $St\leqslant 3.7$ do not show the tails of the PDFs of $\dot{\unicode[STIX]{x1D703}}_{p}$ becoming narrower, their results clearly show that the widening of the tails saturates (see inset of figure 5 in Bec et al. (Reference Bec, Homann and Krstulovic2014)). It is possible that if they had considered larger $St$ , they would have also begun to observe a narrowing of the tails as $St$ was further increased. Possible reasons why the widening of the tails saturates at a lower value of $St$ in our DNS than it does in theirs include is the effect of Reynolds number ( $Re_{\unicode[STIX]{x1D706}}=315$ in their DNS, whereas in our DNS $Re_{\unicode[STIX]{x1D706}}=88$ ), and differences in the scalar forcing method. This raises the question about Reynolds number dependency. We expect that the strong intermittency of advected passive scalars in high Reynolds number flows may affect the results quantitatively. However, two-way coupled simulations at higher Reynolds numbers are computationally demanding and are left for a future work.

4.1 Fluctuations of the carrier temperature field

Figure 5(a,b) shows the normalized one-point PDF of the carrier flow temperature for $St=1$ and $St=3$ , respectively, and for various $St_{\unicode[STIX]{x1D703}}$ . The PDFs are normalized with the standard deviation of the distribution $\unicode[STIX]{x1D70E}_{T}$ . The PDFs are almost Gaussian for low $St_{\unicode[STIX]{x1D703}}$ , while the tails become wider as $St_{\unicode[STIX]{x1D703}}$ is increased. However, we are unable to explain the cause of this enhanced non-Gaussianity. The temperature PDFs are also not symmetric, and display a bump in the right tail. This behaviour was also reported by Overholt & Pope (Reference Overholt and Pope1996) for the case without particles, and it appears to be a low Reynolds number effect that is also dependent on the forcing method employed.

Figure 5. PDF of the carrier flow temperature for $St=1$  (a) and $St=3$ (b), and for various $St_{\unicode[STIX]{x1D703}}$ . (c) Variance of the carrier flow temperature fluctuations for different $St$ as a function of $St_{\unicode[STIX]{x1D703}}$ . (d) Kurtosis of the carrier flow temperature PDF. These results are from simulations S1 in which the two-way thermal coupling is considered.

The effect of $St$ on $\unicode[STIX]{x1D70E}_{T}$ is striking, whereas we saw earlier in figure 2(c) that $\unicode[STIX]{x1D712}_{f}$ only weakly depends on $St$ . To explain the dependence upon the Stokes number we note that the energy balance (3.2) can be rewritten as

(4.1) $$\begin{eqnarray}\unicode[STIX]{x1D712}=\unicode[STIX]{x1D705}\left[\langle \Vert \unicode[STIX]{x1D735}T\Vert ^{2}\rangle +\frac{2}{3}\frac{\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D702}^{2}}\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}_{0}}\frac{1}{St}\langle (T(\boldsymbol{x}_{p},t)-\unicode[STIX]{x1D703}_{p})^{2}\rangle \right].\end{eqnarray}$$

The factor $\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{p}/(\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D702}^{2})$ is constant in our simulations. Therefore, since our DNS data suggest that $\unicode[STIX]{x1D712}_{f}$ is a function of $St_{\unicode[STIX]{x1D703}}$ only (see figure 2 c), from (4.1) and (2.3c ) we obtain

(4.2) $$\begin{eqnarray}\langle T(\boldsymbol{x}_{p},t)^{2}\rangle -\langle \unicode[STIX]{x1D703}_{p}^{2}\rangle \propto Stf(St_{\unicode[STIX]{x1D703}}).\end{eqnarray}$$

The kurtosis of the fluid temperature fluctuation is shown in figure 5(d), as a function of $St_{\unicode[STIX]{x1D703}}$ and for various $St$ . For small $St_{\unicode[STIX]{x1D703}}$ , the kurtosis of the fluid temperature fluctuation is close to the value for a Gaussian PDF, namely $3$ . However, as $St_{\unicode[STIX]{x1D703}}$ is increased, the kurtosis increases. Furthermore, the kurtosis decreases with increasing $St$ for the range considered in our simulations. The explanation of these trends in the kurtosis is unclear.

4.2 Particle temperature fluctuations

Figure 6(a,b) shows the normalized one-point PDF of the particle temperature with $St=1$ , for various $St_{\unicode[STIX]{x1D703}}$ , and for simulations S1 and S2. Figure 6(c–d) shows the corresponding results for $St=3$ , and the PDFs are normalized by their standard deviations. When the two-way thermal coupling is accounted for, the tails of the particle temperature distribution tend to become wider as $St_{\unicode[STIX]{x1D703}}$ is increased. On the other hand, when the two-way coupling is neglected, the PDF of the particle temperature is very close to Gaussian, and its shape is not sensitive to either $St$ or $St_{\unicode[STIX]{x1D703}}$ .

Figure 6. PDF of the particle temperature for $St=1$  (a,b) and $St=3$  (c,d), for various $St_{\unicode[STIX]{x1D703}}$ . Panels (a–c) are from simulations S1, in which the two-way thermal coupling is considered, while panels (b–d) are from simulations S2, in which the two-way coupling is neglected. (e) Variance of the particle temperature fluctuations for different $St$ numbers as a function of $St_{\unicode[STIX]{x1D703}}$ . (f) Kurtosis of the particle temperature distribution.

The variance of the particle temperature fluctuations monotonically decreases with increasing $St_{\unicode[STIX]{x1D703}}$ , as shown in figure 6(e). The results also show a strong dependence on $St$ , but most interestingly, the dependence on $St$ is the opposite for the cases with and without two-way coupling. To understand this we note that using the formal solution to the equation for $\dot{\unicode[STIX]{x1D703}}_{p}(t)$ (ignoring initial conditions) we may construct the result

(4.3) $$\begin{eqnarray}\langle \unicode[STIX]{x1D703}_{p}^{2}(t)\rangle =\frac{1}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}^{2}}\int _{0}^{t}\int _{0}^{t}\langle T(\boldsymbol{x}_{p}(s),s)T(\boldsymbol{x}_{p}(s^{\prime }),s^{\prime })\rangle \text{e}^{-(2t-s-s^{\prime })/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}}\,\text{d}s\,\text{d}s^{\prime }.\end{eqnarray}$$

If we now substitute into this the exponential approximation

(4.4) $$\begin{eqnarray}\langle T(\boldsymbol{x}_{p}(s),s)T(\boldsymbol{x}_{p}(s^{\prime }),s^{\prime })\rangle \approx \langle T^{2}(\boldsymbol{x}_{p}(t),t)\rangle \exp [-|s-s^{\prime }|/\unicode[STIX]{x1D70F}_{T}],\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{T}$ is the time scale of $T(\boldsymbol{x}_{p}(t),t)$ , then we obtain

(4.5) $$\begin{eqnarray}\langle \unicode[STIX]{x1D703}_{p}^{2}(t)\rangle =\frac{\langle T^{2}(\boldsymbol{x}_{p}(t),t)\rangle }{1+\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}/\unicode[STIX]{x1D70F}_{T}}.\end{eqnarray}$$

This result reveals that the particle temperature variance is influenced by $St$ in two ways. First, $\langle T^{2}(\boldsymbol{x}_{p}(t),t)\rangle$ depends upon the spatial clustering of the inertial particles, and this depends essentially upon $St$ . Second, the time scale $\unicode[STIX]{x1D70F}_{T}$ is the time scale of the fluid temperature field measured along the inertial particle trajectories, and hence depends upon $St$ . For isotropic turbulence, this time scale is expected to decrease as $St$ is increased, which would lead to $\langle \unicode[STIX]{x1D703}_{p}^{2}(t)\rangle$ decreasing as $St$ increases, which is the behaviour observed in figure 6(e). In the presence of two-way coupling, however, $\langle T^{2}(\boldsymbol{x},t)\rangle$ increases with increasing $St$ , as shown earlier. In the two-way coupled regime this increase in $\langle T^{2}(\boldsymbol{x},t)\rangle$ leads to an increase in $\langle T^{2}(\boldsymbol{x}_{p}(t),t)\rangle$ that dominates over the decrease of $\unicode[STIX]{x1D70F}_{T}$ with increasing $St$ , and as a result $\langle \unicode[STIX]{x1D703}_{p}^{2}(t)\rangle$ increases with increasing $St$ .

The kurtosis of the particle temperature increases with increasing $St_{\unicode[STIX]{x1D703}}$ when the two-way thermal coupling is accounted for, as shown in figure 6(f) (simulations S1, filled symbols). Conversely, the kurtosis of the particle temperature remains constant as $St_{\unicode[STIX]{x1D703}}$ is increased when the two-way thermal coupling is ignored (simulations S2, open symbols).

5.1 Particle clustering on the temperature fronts

It is well known that inertial particles in turbulence form clusters (Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007), which may be quantified using the radial distribution function (RDF). As shown in figure 7(a), the particle number density in our simulations at small separations is a order of magnitude larger than the mean density when $St=O(1)$ . Bec et al. (Reference Bec, Homann and Krstulovic2014) showed that inertial particles also exhibit a tendency to preferentially cluster in the fluid temperature fronts where the temperature gradients are large. To demonstrate this, they measured the temperature dissipation rate at the particle positions and showed that this was higher than the Eulerian dissipation rate of the fluid temperature fluctuations. Note that previous works (Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2013, Reference Gualtieri, Picano, Sardina and Casciola2015) have shown that the radial distribution function can be reduced by momentum coupling, which we are neglecting. A future work that includes two-way momentum coupling should consider how this affects the way inertial particles sample high temperature gradients in the flow.

Figure 7. (a) Radial distribution function (RDF) as a function of the separation $r/\unicode[STIX]{x1D702}$ for various $St$ . (b) Particle number density conditioned on the magnitude of the fluid temperature gradient at the particle position, for various $St$ . These results are from simulations S2, in which the two-way thermal coupling is neglected.

We quantify the tendency for inertial particles to cluster in the fluid temperature fronts by computing the single particle position probability density, conditioned upon the norm of the fluid temperature gradient,

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D70C}(\boldsymbol{x}_{p}|\,\Vert \unicode[STIX]{x1D735}T\Vert )=\frac{\unicode[STIX]{x1D70C}(\Vert \unicode[STIX]{x1D735}T\Vert (\boldsymbol{x}_{p}))}{\unicode[STIX]{x1D70C}(\Vert \unicode[STIX]{x1D735}T\Vert )}.\end{eqnarray}$$

This conditioned probability can also be understood as the ratio between the fraction of inertial particles $n_{p}(St;\Vert \unicode[STIX]{x1D735}T\Vert )$ located in a region of a given temperature gradient magnitude $\Vert \unicode[STIX]{x1D735}T\Vert$ and the number of particles $n_{p}(0;\Vert \unicode[STIX]{x1D735}T\Vert )$ which would be located in the same region for $St\rightarrow 0$ , that is, when particles follow fluid trajectories,

(5.2) $$\begin{eqnarray}\unicode[STIX]{x1D70C}(\boldsymbol{x}_{p}\mid \Vert \unicode[STIX]{x1D735}T\Vert )=\frac{n_{p}(St;\Vert \unicode[STIX]{x1D735}T\Vert )}{n_{p}(0;\Vert \unicode[STIX]{x1D735}T\Vert )}.\end{eqnarray}$$

By defining $\Vert \unicode[STIX]{x1D735}T\Vert _{rms}$ as the root mean square (r.m.s.) value of $\Vert \unicode[STIX]{x1D735}T\Vert$ , small values of $\Vert \unicode[STIX]{x1D735}T\Vert /\Vert \unicode[STIX]{x1D735}T\Vert _{rms}$ may be interpreted as corresponding to the large scales, and are associated with the Lagrangian coherent structures in which the temperature field is almost constant. Large values of $\Vert \unicode[STIX]{x1D735}T\Vert /\Vert \unicode[STIX]{x1D735}T\Vert _{rms}$ may be interpreted as corresponding to the small scales, and are associated with fronts in the fluid temperature field. The results for $\unicode[STIX]{x1D70C}(\boldsymbol{x}_{p}|\,\Vert \unicode[STIX]{x1D735}T\Vert )$ are shown in figure 7(b), for the simulations without two-way thermal coupling (the results show only a weak dependence on $St_{\unicode[STIX]{x1D703}}$ when the two-way coupling is included). Due to the clustering on the temperature fronts, $\unicode[STIX]{x1D70C}(\boldsymbol{x}_{p}|\Vert \unicode[STIX]{x1D735}T\Vert )$ is an increasing function of $\Vert \unicode[STIX]{x1D735}T\Vert$ and it is larger than unity in the region of large temperature gradient. The probability of observing a gradient of a certain magnitude (which is proportional to $n_{p}(0;\Vert \unicode[STIX]{x1D735}T\Vert )$ ) decays almost exponentially with increasing $\Vert \unicode[STIX]{x1D735}T\Vert$ , as in figure 2(a). For values of $St$ at which the maximum particle clustering takes place, $n_{p}(St;\Vert \unicode[STIX]{x1D735}T\Vert )$ is up to four times larger than $n_{p}(0;\Vert \unicode[STIX]{x1D735}T\Vert )$ in regions of strong temperature gradients. We expect even higher values at the largest $\Vert \unicode[STIX]{x1D735}T\Vert$ , however it is difficult to obtain statistically relevant results in correspondence with such extreme events. These results therefore support the conclusions of Bec et al. (Reference Bec, Homann and Krstulovic2014) that inertial particles preferentially cluster in the fronts of the fluid temperature field where $\Vert \unicode[STIX]{x1D735}T\Vert /\Vert \unicode[STIX]{x1D735}T\Vert _{rms}$ is large.

5.2 Particle motion across the temperature fronts

To obtain further insight into the thermal coupling between the particles and fluid we consider the properties of the particle heat flux conditioned on $\Vert \unicode[STIX]{x1D735}T\Vert$ . In particular, we consider the following quantity

(5.3) $$\begin{eqnarray}q_{n}(\Vert \unicode[STIX]{x1D735}T\Vert )\equiv (T(\boldsymbol{x}_{p})-\unicode[STIX]{x1D703}_{p})^{n}\boldsymbol{v}_{p}\boldsymbol{\cdot }\boldsymbol{n}_{T}(\boldsymbol{x}_{p})|_{\Vert \unicode[STIX]{x1D735}T\Vert },\end{eqnarray}$$

where $\boldsymbol{n}_{T}$ is the normalized, resolved, temperature gradient

(5.4) $$\begin{eqnarray}\boldsymbol{n}_{T}(\boldsymbol{x}_{p})\equiv \frac{\unicode[STIX]{x1D735}T(\boldsymbol{x}_{p})}{\Vert \unicode[STIX]{x1D735}T(\boldsymbol{x}_{p})\Vert }.\end{eqnarray}$$

The statistics of $q_{n}$ provide a way to quantify the relationship between the particle heat flux and the local carrier temperature gradients in the fluid. Understanding this relationship is key to understanding how the particles modify the properties of the fluid temperature and temperature gradient fields. It is justified to investigate the interaction between the resolved temperature field, in which the near-particle disturbances are not represented, since, in the dilute regime, particles are statistically far enough that a particle rarely finds itself in the disturbance region of another particle. As discussed in appendix B, the norm of the gradient of the perturbation field induced by the particle is proportional to $r_{p}^{-2}$ and, therefore, it is usually large.

Figure 8. (a) Results for $\langle |q_{0}(\Vert \unicode[STIX]{x1D735}T\Vert )|\rangle /u_{\unicode[STIX]{x1D702}}$ , for various $St$ . (b) Results for $\langle |\cos \unicode[STIX]{x1D6FC}_{p}|\rangle$ as a function of $\Vert \unicode[STIX]{x1D735}T\Vert$ , for various $St$ . These results are from simulations S2, in which the two-way thermal coupling is neglected.

The efficiency with which the particles cross the fronts in the carrier flow temperature field is quantified by $\langle |q_{0}|\rangle$ , and our results for this quantity in one-way coupled simulations are shown in figure 8(a). The curves are normalized with the Kolmogorov velocity scale $u_{\unicode[STIX]{x1D702}}$ . At moderate Stokes number, particles tend to accumulate near the front, therefore they cross the front with a small velocity. On the other hand, particles with large inertia slowly respond to a change of the local velocity/temperature and therefore they are less affected by the local value of the temperature gradient, carrying large temperature increments across the temperature field. Accordingly, the velocity magnitude becomes nearly independent of the local value of the temperature gradient as the particle inertia is increased, as shown in figure 8.

It is also important to consider whether the reduction of $\langle |q_{0}|\rangle$ as $\Vert \unicode[STIX]{x1D735}T\Vert$ increases is due to the reduction of the norm of the particle velocity or to the lack of alignment between the particle velocity and the fluid temperature gradient at the particle position. Figure 8(b) displays the average of the absolute value of the cosine of the angle between the particle velocity and temperature gradient

(5.5) $$\begin{eqnarray}\cos \unicode[STIX]{x1D6FC}_{p}\equiv \frac{\boldsymbol{v}_{p}}{\Vert \boldsymbol{v}_{p}\Vert }\boldsymbol{\cdot }\frac{\unicode[STIX]{x1D735}T(\boldsymbol{x}_{p})}{\Vert \unicode[STIX]{x1D735}T(\boldsymbol{x}_{p})\Vert },\end{eqnarray}$$

conditioned on $\Vert \unicode[STIX]{x1D735}T\Vert$ . The results show that as $\Vert \unicode[STIX]{x1D735}T\Vert$ is increased, the particle motion becomes misaligned with the local carrier flow temperature gradient. This then shows that the reduction of $\langle |q_{0}|\rangle$ as $\Vert \unicode[STIX]{x1D735}T\Vert$ increases is due to non-trivial statistical geometry in the system. The results also show that as $St$ is increased, the cosine of the angle between the fluid temperature gradient and the particle velocity becomes almost independent of $\Vert \unicode[STIX]{x1D735}T\Vert$ , and $\langle |\cos \unicode[STIX]{x1D6FC}_{p}|\rangle \approx 1/2$ , the value corresponding to $\cos \unicode[STIX]{x1D6FC}_{p}$ being a uniform random variable. This shows that as $St$ is increased, the correlation between the direction of the particle velocity and the local carrier fluid velocity gradient vanishes.

5.3 Heat flux due to the particle motion across the fronts

We now turn to consider the quantity $\langle q_{1}\rangle$ . When the particle moves from a cold to a warm region of the fluid, the component of the particle velocity along the temperature gradient is positive, $\boldsymbol{v}_{p}\boldsymbol{\cdot }\boldsymbol{n}_{T}(\boldsymbol{x}_{p})>0$ . If the particle is also cooler than the local fluid so that $T(\boldsymbol{x}_{p})-\unicode[STIX]{x1D703}_{p}>0$ , then as it moves into the region where the fluid is warmer, $q_{1}>0$ meaning that the particle will absorb heat from the fluid, and will therefore tend to reduce the local fluid temperature gradient. When the particle moves from a warm to a cold region of the flow, if $T(\boldsymbol{x}_{p})-\unicode[STIX]{x1D703}_{p}<0$ then $q_{1}$ is also positive, so that again the particle will act to reduce the local temperature gradient in the fluid. Therefore, $q_{1}>0$ indicates that the action of the inertial particles is to smooth out the fluid temperature field, reducing the magnitude of its temperature gradients, and $q_{1}<0$ implies the particles enhance the temperature gradients.

Figure 9. Results for $\langle q_{1}(\Vert \unicode[STIX]{x1D735}T\Vert )\rangle /(u_{\unicode[STIX]{x1D702}}T_{\unicode[STIX]{x1D702}})$ for $St=0.5$  (a,b),  $St=1$  (c,d) and $St=3$  (e,f) and for various $St_{\unicode[STIX]{x1D703}}$ . Panels (a,c,e) are from simulations S1, in which the two-way thermal coupling is considered, while panels (b,d,f) are from simulations S2, in which the two-way coupling is neglected.

The results for $\langle q_{1}\rangle$ are shown in figure 9 for various $St$ and $St_{\unicode[STIX]{x1D703}}$ , including (simulations S1) and neglecting (simulations S2) the two-way thermal coupling. On average we observe $\langle q_{1}\rangle \geqslant 0$ , such that the particles tend to make the fluid temperature field more uniform. The results show that $\langle q_{1}\rangle$ tends to zero as $\Vert \unicode[STIX]{x1D735}T\Vert \rightarrow 0$ . This indicates that the particles spend enough time in the Lagrangian coherent structures to adjust to the temperature of the fluid. However, $\langle q_{1}\rangle$ increases significantly as $\Vert \unicode[STIX]{x1D735}T\Vert$ increases, suggesting that inertial particles can carry large temperature differences across the fronts. In the limit $St_{\unicode[STIX]{x1D703}}\rightarrow 0$ , $\langle q_{1}\rangle \rightarrow 0$ reflecting the thermal equilibrium between the particles and the fluid. As $St_{\unicode[STIX]{x1D703}}$ is increased, the heat flux becomes finite, however, if $St_{\unicode[STIX]{x1D703}}$ is too large, the particle temperature decorrelates from the fluid temperature and the heat exchange is not effective. Hence, $\langle q_{1}\rangle$ can saturate with increasing $St_{\unicode[STIX]{x1D703}}$ . The results show that $\langle q_{1}\rangle$ increases with increasing $St$ , associated with the decoupling of $\boldsymbol{v}_{p}$ and $\boldsymbol{n}_{T}(\boldsymbol{x}_{p})$ discussed earlier. Finally, the results also show that two-way thermal coupling reduces $\langle q_{1}\rangle$ . This is simply a reflection of the fact that since the particles tend to smooth out the fluid temperature gradients, the disequilibrium between the particle and local fluid temperature is reduced, which in turn reduces the heat flux due to the particles.

6.1 Structure functions of the carrier temperature field

The $n$ th-order structure function of the resolved fluid temperature field is defined as

(6.1) $$\begin{eqnarray}S_{T}^{n}(r)\equiv \langle |\unicode[STIX]{x0394}T(r,t)|^{n}\rangle ,\end{eqnarray}$$

where $\unicode[STIX]{x0394}T(r,t)$ it the difference in the carrier temperature field at two points separated by the distance $r$ (the ‘temperature increment’). The results for $S_{T}^{2}$ , with different $St$ and $St_{\unicode[STIX]{x1D703}}$ are shown in figure 10. The structure functions of the actual temperature field differ from the structure functions of the carrier temperature field, due to the near-particle temperature disturbances. As discussed in appendix B, the impact of the local near-particle perturbation is marked at small separation and the carrier flow temperature field can be understood as the actual temperature field filtered at the grid resolution scale. In order to emphasize this fact, the carrier flow temperature structure functions are reported only down to the Kolmogorov scale, which is comparable to the grid spacing.

Figure 10. Results for $S_{T}^{2}$ for different $St_{\unicode[STIX]{x1D703}}$ , for $St=0.5$  (a),  $St=1$  (b) and $St=3$  (c). (d) Scaling exponents of the fluid temperature structure functions at small separation, $r\leqslant 2\unicode[STIX]{x1D702}$ , at $St=1$ . The data are from simulations S1 in which the two-way thermal coupling is considered.

The results show that $S_{T}^{2}$ decreases monotonically with increasing $St_{\unicode[STIX]{x1D703}}$ at all scales when the two-way thermal coupling is taken to account. In the dissipation range, $S_{T}^{2}$ is directly connected to the dissipation rate, and is suppressed in the same way for the three different $St$ considered. Conversely, the suppression of the large-scale fluctuations is stronger as $St$ is reduced, at least for the range of $St$ considered here.

The scaling exponents of the structure functions of the carrier temperature field

(6.2) $$\begin{eqnarray}\unicode[STIX]{x1D701}_{T}^{n}\equiv \frac{\text{d}\text{log}\text{S}_{T}^{n}(r)}{\text{d}\text{log}\,r}\end{eqnarray}$$

are shown in figure 10(d) for $r\leqslant 2\unicode[STIX]{x1D702}$ . The results show that the resolved fluid temperature field remains smooth (to within numerical uncertainty) even when inertial particles are suspended in the flow.

6.2 Particle temperature structure functions

The $n$ th-order structure function of the particle temperature $\unicode[STIX]{x1D703}_{p}(t)$ is defined as

(6.3) $$\begin{eqnarray}S_{\unicode[STIX]{x1D703}}^{n}(r)\equiv \langle |\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}|^{n}\rangle _{r},\end{eqnarray}$$

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)$ is the difference in the temperature of the two particles, and the brackets denote an ensemble average, conditioned on the two particles having separation $r$ . The results for $S_{\unicode[STIX]{x1D703}}^{2}$ for different $St$ and $St_{\unicode[STIX]{x1D703}}$ , with and without two-way thermal coupling, are shown in figure 11.

Figure 11. Results for $S_{\unicode[STIX]{x1D703}}^{2}$ for different $St_{\unicode[STIX]{x1D703}}$ , for $St=0.5$  (a,b),  $St=1$  (c,d) and $St=3$  (e,f). Panels (a,c,e) are from simulations S1, in which the two-way thermal coupling is considered, while panels (b,d,f) are from simulations S2, in which the two-way coupling is neglected.

The results show that $S_{\unicode[STIX]{x1D703}}^{2}$ depends on $St_{\unicode[STIX]{x1D703}}$ in much the same way as the inertial particle relative velocity structure functions depend on $St$ (Ireland et al. Reference Ireland, Bragg and Collins2016a ). This is not surprising since the equation governing $\dot{\unicode[STIX]{x1D703}}_{p}$ is structurally identical to the equation governing the particle acceleration. However, important differences are that $\dot{\unicode[STIX]{x1D703}}_{p}$ depends on both $St$ and $St_{\unicode[STIX]{x1D703}}$ , and also that the fluid temperature field is structurally different from the fluid velocity field, with the temperature field exhibiting the well-known ramp–cliff structure.

To obtain further insight into the behaviour of $S_{\unicode[STIX]{x1D703}}^{2}$ and $S_{\unicode[STIX]{x1D703}}^{n}$ in general, we note that the formal solution for $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)$ is given by (ignoring initial conditions)

(6.4) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)=\frac{1}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}}\int _{0}^{t}\unicode[STIX]{x0394}T(\boldsymbol{x}_{p}(s),\boldsymbol{r}_{p}(s),s)\exp \left(-\frac{t-s}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}}\right)\,\text{d}s,\end{eqnarray}$$

where $\unicode[STIX]{x0394}T(\boldsymbol{x}_{p}(s),\boldsymbol{r}_{p}(s),s)$ is the difference in the fluid temperature at the two particle positions $\boldsymbol{x}_{p}(s)$ and $\boldsymbol{x}_{p}(s)+\boldsymbol{r}_{p}(s)$ . Equation (6.4) shows that $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)$ depends upon $\unicode[STIX]{x0394}T$ along the path history of the particles, and $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)$ is therefore a non-local quantity. The role of the path history increases as $St_{\unicode[STIX]{x1D703}}$ is increased since the exponential kernel in the convolution integral decays more slowly as $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}$ is increased. Since the statistics of $\unicode[STIX]{x0394}T$ increase with increasing separation, particle pairs at small separations are able to be influenced by larger values of $\unicode[STIX]{x0394}T$ along their path history, such that $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)$ can significantly exceed the local fluid temperature increment $\unicode[STIX]{x0394}T(\boldsymbol{x}_{p}(t),\boldsymbol{r}_{p}(t),t)$ . This then causes $S_{\unicode[STIX]{x1D703}}^{2}$ to increase with increasing $St_{\unicode[STIX]{x1D703}}$ , as shown in figure 11. This effect is directly analogous to the phenomena of caustics that occur in the relative velocity distributions of inertial particles at the small scales of turbulence (Wilkinson & Mehlig Reference Wilkinson and Mehlig2005), and which occur because the inertial particle relative velocities depend non-locally on the fluid velocity increments experienced along their trajectory history (Bragg & Collins Reference Bragg and Collins2014b ). In analogy, we may therefore refer to the effect as ‘thermal caustics’, and they may be of particular importance for particle-laden turbulent flows where particles in close proximity thermally interact.

The results in figure 11 also reveal a strong effect of $St$ , and one way that $St$ affects these results is through the spatial clustering and preferential sampling of the fluid temperature field by the inertial particles. There is, however, another mechanism through which $St$ can affect $S_{\unicode[STIX]{x1D703}}^{2}$ . In particular, since, due to caustics, the relative velocity of the particles increases with increasing $St$ at the small scales, then the values of $\unicode[STIX]{x0394}T(\boldsymbol{x}_{p}(s),\boldsymbol{r}_{p}(s),s)$ that may contribute to $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)$ become larger. This follows since if their relative velocities are larger, then over the time span $t-s\leqslant O(\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}})$ the particle pair can come from even larger scales where (statistically) $\unicode[STIX]{x0394}T(\boldsymbol{x}_{p}(s),\boldsymbol{r}_{p}(s),s)$ is bigger. This effect would cause $S_{\unicode[STIX]{x1D703}}^{2}$ to increase with $St$ for a given $St_{\unicode[STIX]{x1D703}}$ , further enhancing the thermal caustics, which is exactly what is observed in figure 11. The results also show that the thermal caustics are stronger for $St_{\unicode[STIX]{x1D703}}\geqslant O(1)$ when the two-way thermal coupling is ignored. This is mainly due to the reduction in the fluid temperature gradients due to the two-way thermal coupling described earlier, noting that in the limit of vanishing fluid temperature gradients, the thermal caustics necessarily disappear.

At larger scales where the statistics of $\unicode[STIX]{x0394}T$ vary more weakly with $r$ , the non-local effect weakens, the thermal caustics disappear, and a filtering mechanism takes over which causes $S_{\unicode[STIX]{x1D703}}^{2}$ to decrease with increasing $St_{\unicode[STIX]{x1D703}}$ . This filtering effect is directly analogous to that dominating the large-scale velocities of inertial particles in isotropic turbulence, and is associated with the sluggish response of the particles to the large-scale flow fluctuations due to their inertia (Ireland et al. Reference Ireland, Bragg and Collins2016a ).

The particle temperature structure functions $S_{\unicode[STIX]{x1D703}}^{n}$ behave as power laws at small separation, $\log S_{\unicode[STIX]{x1D703}}^{n}(r)\approx \unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}}^{n}\log r+a^{n}$ , and the associated scaling exponents $\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}}^{n}$ are shown in figure 12. The exponents are obtained by fitting the logarithm of the structure function in the dissipation range according to ordinary least squares. To reduce statistical noise, we estimate $\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}}^{n}$ by fitting the data for $S_{\unicode[STIX]{x1D703}}^{n}$ over the range $0.2\unicode[STIX]{x1D702}\leqslant r\leqslant 2\unicode[STIX]{x1D702}$ . Over this range, $S_{\unicode[STIX]{x1D703}}^{n}$ do not strictly behave as power laws, and hence the exponents measured are understood as average exponents. The error bars indicate the largest deviations from the mean exponent observed in the considered range. The results in figure 12 reveal that particle temperature increments exhibit a strong multifractal behaviour. This multifractality is due to the non-local thermal dynamics of the particles and the formation of thermal caustics, described earlier. In particular, there exists a finite probability to find inertial particle pairs that are very close but have large temperature differences because they experienced very different fluid temperatures along their trajectory histories. As with the thermal caustics, the multifractality is enhanced as $St$ is increased. Most interestingly, the results for $\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}}^{n}$ are only weakly affected by the two-way thermal coupling, despite the fact that we observed a significant effect of the coupling on $S_{\unicode[STIX]{x1D703}}^{2}$ . This suggests that the two-way coupling affects the strength of the thermal caustics, but only weakly affects the scaling of the structure functions in the dissipation range.

Figure 12. Scaling exponent of the structure functions of the particle temperature at small separation, $0.2\unicode[STIX]{x1D702}\leqslant r\leqslant 2\unicode[STIX]{x1D702}$ , for various thermal Stokes numbers $St_{\unicode[STIX]{x1D703}}$ , at fixed Stokes number $St=0.5$  (a,b) and $St=1$  (c,d). Panels (a,c) are from simulations S1, in which the two-way thermal coupling is considered, while panels (b,d) are from simulations S2, in which the two-way coupling is neglected.

Figure 13. Second-order mixed structure functions of the carrier flow temperature field, for different thermal Stokes numbers of the suspended particles, at $St=0.5$  (a) and $St=1$  (b). The data refer to the set of simulations S1, with thermal particle back reaction included.

6.3 Mixed structure functions

We turn to considering the behaviour of the flux of the temperature increments across the scales of the flow, which is associated with the mixed structure functions

(6.5) $$\begin{eqnarray}S_{Q}(r)\equiv \langle (\unicode[STIX]{x0394}T(r,t))^{2}\unicode[STIX]{x0394}u_{\Vert }(r,t)\rangle ,\end{eqnarray}$$

where $\unicode[STIX]{x0394}u_{\Vert }$ is the longitudinal relative velocity difference. The results for $S_{Q}$ , for different $St$ and $St_{\unicode[STIX]{x1D703}}$ are shown in figure 13. Just as we observed for the fluid temperature structure functions, $-S_{Q}$ decreases monotonically with increasing $St_{\unicode[STIX]{x1D703}}$ , as was also observed for the fluid temperature dissipation rate $\unicode[STIX]{x1D712}_{f}$ . The mixed structure functions of the carrier temperature field are reported to separation down to the Kolmogorov scale, that is the scale at which the carrier temperature field is resolved, as discussed in appendix B.

Figure 14. Second-order mixed structure functions of the particle temperature, for different thermal Stokes numbers, at $St=0.5$  (a,b) and $St=1$  (c,d). Panels (a,c) refer to the set of simulations S1, in which the thermal particle back reaction is included. Panels (b,d) refer to the set of simulations S2, in which the thermal particle back reaction is neglected.

To consider the flux of the particle temperature increments, we begin by considering the exact equation that can be constructed for $S_{\unicode[STIX]{x1D703}}^{n}$ using PDF transport equations. In particular, if we introduce the PDF ${\mathcal{P}}(\boldsymbol{r},\unicode[STIX]{x0394}\unicode[STIX]{x1D703},t)\equiv \langle \unicode[STIX]{x1D6FF}(\boldsymbol{r}_{p}(t)-\boldsymbol{r})\unicode[STIX]{x1D6FF}(\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)-\unicode[STIX]{x0394}\unicode[STIX]{x1D703})\rangle$ and the associated marginal PDF $\unicode[STIX]{x1D71A}(\boldsymbol{r},t)\equiv \int {\mathcal{P}}\,\text{d}\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$ , where $\boldsymbol{r}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$ are time-independent phase-space coordinates, then we may derive for a statistically stationary system the result (see Bragg & Collins (Reference Bragg and Collins2014a ), Bragg et al. (Reference Bragg, Ireland and Collins2015b ) for details on how to derive such results)

(6.6) $$\begin{eqnarray}\langle [\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)]^{2}\rangle _{\boldsymbol{ r}}=\langle \unicode[STIX]{x0394}T(\boldsymbol{x}_{p}(t),\boldsymbol{r}_{p}(t),t)\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)\rangle _{\boldsymbol{r}}-\frac{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}}{2\unicode[STIX]{x1D71A}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{r}}\boldsymbol{\cdot }\unicode[STIX]{x1D71A}\langle [\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)]^{2}\boldsymbol{w}_{p}(t)\rangle _{\boldsymbol{r}},\end{eqnarray}$$

where $\boldsymbol{w}_{p}(t)\equiv \unicode[STIX]{x2202}_{t}\boldsymbol{r}_{p}(t)$ . The first term on the right-hand side is the local contribution that remains when there exist no fluxes across the scales, and this term determines the behaviour of $\langle [\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)]^{2}\rangle _{\boldsymbol{r}}$ at the large scales of homogeneous turbulence where the statistics are independent of $\boldsymbol{r}$ . The second term on the right-hand side is the non-local contribution that arises for $St_{\unicode[STIX]{x1D703}}>0$ , and it is this term that is responsible for the thermal caustics discussed earlier. It depends on the spatial clustering of the particles through $\unicode[STIX]{x1D71A}$ (which is proportional to the RDF), and the flux $\langle [\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)]^{2}\boldsymbol{w}_{p}(t)\rangle _{\boldsymbol{r}}$ which, for an isotropic system, is determined by the longitudinal component

(6.7) $$\begin{eqnarray}S_{Q_{p}}(r)\equiv \frac{\boldsymbol{r}}{r}\boldsymbol{\cdot }\langle [\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{p}(t)]^{2}\boldsymbol{w}_{p}(t)\rangle _{r}.\end{eqnarray}$$

The results for $S_{Q_{p}}$ from our simulations are shown in figure 14, and they show that without two-way coupling, $-S_{Q_{p}}$ monotonically increases with increasing $St_{\unicode[STIX]{x1D703}}$ at the smallest scales. However, with two-way coupling, $-S_{Q_{p}}$ is maximum for intermediate values of $St_{\unicode[STIX]{x1D703}}$ , and this occurs because as shown earlier, as $St_{\unicode[STIX]{x1D703}}$ is increased, the fluid temperature fluctuations are suppressed across the scales.