Modulation of fluid temperature fluctuations by particles in turbulence

Abstract Modulation of fluid temperature fluctuations by particles due to thermal interaction in homogeneous isotropic turbulence is studied. For simplicity, only thermal coupling between the fluid and particles is considered, and momentum coupling is neglected. Application of the statistical theory used in cloud turbulence research leads to the prediction that modulation of the intensity of fluid temperature fluctuations by particles is expressed as a function of the Damköhler number, which is defined as the ratio of the turbulence large-eddy turnover time to the fluid thermal relaxation time. Direct numerical simulations are conducted for two-way thermal coupling between the fluid temperature field and point particles in homogeneous isotropic turbulence. The simulation results are shown to agree well with the theoretical predictions.


Introduction
Turbulent flows in a broad range of natural and engineering applications include small particles, such as water droplets in clouds (Chen et al. 2018), dust particles in protoplanetary disks (Ishihara et al. 2018) and pneumatic transport of solid particles (Ebrahimi & Crapper 2017). These small particles are transported in turbulence and interact with and modify the carrier flow by various processes, having significant results. For example, water droplets in clouds modify the surrounding supersaturation (humidity) field by condensation/evaporation (Grabowski & Wang 2013), and small heavy particles in turbulence modify the velocity field of the carrier fluid by the drag force (Elghobashi & Abou-Arab 1983). These turbulence modulations by small particles crucially affect the subsequent development of the system, such as the fluid mean velocity, heat and mass transfer, dispersion and the growth of particles. Moreover, particles in turbulence lead to many non-trivial phenomena, such as inertial clustering of particles (Sundaram & Collins 1997), the sling effect (Falkovich & Pumir 2007) and the clustering of particles in the scalar fronts (Bec, Homann & Krstulovic 2014). These complex two-way interactions between turbulence and particles, as well as associated phenomena, make it difficult to develop a generally applicable theory or model that can predict the behaviour of particle-laden turbulent flows (Poelma, Westerweel & Ooms 2007).
Recently, the development of a model for the particle-turbulence interaction in the research field of cloud turbulence has advanced remarkably (Shaw 2003;Devenish et al. 2012;Grabowski & Wang 2013). In the turbulent environment in natural clouds, cloud droplets experience variations in the surrounding supersaturation and velocity field, resulting in different growth histories for cloud droplets by condensation/evaporation and hence broadening of the drop size distributions and acceleration of rain initiation (Sedunov 1974;Cooper 1989;Korolev & Mazin 2003;Lasher-Trapp, Copper & Blyth 2005). Two-way interactions between turbulence and cloud droplets by condensation/evaporation are modelled by Langevin-type equations. This Langevin model predicts modulation of the supersaturation field by cloud droplets and broadening of drop-size distributions. The validity and usefulness of the Langevin model have been demonstrated in various studies, including numerical simulations (Sardina et al. 2015;Siewert, Bec & Krstulvic 2017;Sardina et al. 2018;Saito, Gotoh & Watanabe 2019a) and laboratory experiments (Chandrakar et al. 2016;Desai et al. 2018;Chandrakar et al. 2020). In addition, a novel type of cloud microphysics parameterization based on a similar Langevin model was developed in order to take into account the effect of supersaturation fluctuations in unresolved scales on the growth of cloud droplets (Grabowski & Abade 2017;Abade, Grabowski & Pawlowska 2018;Thomas, Grabowski & Kumar 2020;Chandrakar et al. 2021). Saito, Watanabe & Gotoh (2019b) considered the case in which the particle-turbulence interaction is due to the drag force. By applying knowledge from cloud turbulence research, they found a new time scale that characterizes turbulence modulation by particles. Based on the results of scaling analysis and direct numerical simulations (DNSs), they demonstrated that turbulence modulation by particles can be expressed as a function of the Damköhler number, which is defined as the ratio of the turbulence large-eddy turnover time to the new time scale.
In order to further extend the study by Saito et al. (2019b), the present study considers the case of thermal coupling, in which particles interact with the fluid temperature field by heat transfer in turbulence. This case has been investigated through numerous DNS studies, a few of which are described below. Bec et al. (2014) conducted DNSs of heavy inertial particles transported by a turbulent flow and demonstrated that particles tend to cluster on the fronts of the temperature field, resulting in anomalous scaling laws. Pouransari & Mani (2018) investigated particle-to-fluid heat transfer in particle-laden turbulence for application to particle-based solar receivers  and demonstrated the importance of the two dimensionless numbers defined based on a turbulence large-eddy turnover time, i.e. the particle Stokes number and the heat-mixing parameter. They also derived the leading-order thermal evolution equation which reveals the ensemble averaged heat transfer between the particle and fluid phases. Carbone, Bragg & Iovieno (2019) investigated the multiscale aspects of particle-turbulence thermal coupling through statistical analysis, using such as probability density functions and structure functions, and demonstrated the dominance of thermal caustics at small scales, at which the particle temperature differences for small separations rapidly increase as the Stokes number and the thermal Stokes number are increased. However, despite these previous studies, a model or theory that can predict the modulation of fluid temperature fluctuations by particles in turbulence is not yet available.
In order to shed new light on this long-standing problem, we apply the Langevin model used in cloud turbulence research and derive an analytical prediction for modulation of fluid temperature fluctuations by particles. We also conduct DNSs to confirm the theoretical prediction. Being aware of the difficulty in understanding the two-way coupling between particles and turbulence, we here make a step forward by considering the simplest case, in which particles are assumed to be monodisperse, spherical point particles that have finite thermal inertia but have no momentum inertia. In other words, particles move as fluid particles and interact with the surrounding fluid only by thermal coupling. Moreover, no compressibility is considered for a fluid due to temperature change. The fluid is assumed to be incompressible with a constant mass density. Through these simplifications, we can straightforwardly apply the Langevin model to obtain the analytical prediction, which provides a fundamental understanding of the mechanism of the two-way thermal coupling between particles and turbulence.
The remainder of the present paper is organized as follows. Section 2 describes the governing equations. Section 3 provides the theory. Important parameters and non-dimensional numbers are described. The Langevin model used in cloud turbulence research is applied and an analytical prediction for modulation of the temperature fluctuation by particles is derived. In § 4, we describe the DNS of the two-way thermal coupling between particles and turbulence and confirm the theoretical prediction. A summary and discussion are provided in § 5.

Governing equations
We consider monodisperse, spherical point particles that have finite thermal inertia but no momentum inertia. Particles move as fluid particles and interact with the surrounding fluid only by thermal coupling. The evolution equations for the particle position are given by where x pj is the jth particle position, and u(x pj , t) is the fluid velocity at x pj at time t. The evolution equation for the particle temperature is given by where θ pj is the temperature of the jth particle, and θ(x pj , t) is the fluid temperature at x pj at time t. In addition, τ θ is the particle thermal relaxation time and is given as where r is the particle radius, κ is the thermal diffusivity of the fluid, ρ p and ρ 0 are the particle and fluid densities, respectively, and c p and c 0 are the specific heat capacities of the particle and fluid, respectively, at constant pressure. The thermal equations for particles (2.2) and (2.3) are based on solutions for spherical particles in an infinite domain. The fluid is assumed to be incompressible with a constant mass density. In order to focus on the thermal coupling, the feedback effect from the particles to the fluid is considered only for the fluid temperature and not for the fluid velocity. The velocity field of a fluid is governed by the incompressible Navier-Stokes equations, where p is the pressure and f u represents the external force, which is solenoidal (∇ · f u = 0). The temperature field for the fluid is governed by the advection-diffusion equation, where f θ represents the external force. Here, F p is the feedback term from particles due to two-way thermal coupling, where N p is the total number of particles in the domain, m p = 4ρ p πr 3 /3 is the particle mass and δ() is the Dirac delta function. Note that, since the temperature equations (2.6) and (2.7) are linear for both θ and θ pj , we can regard θ and θ pj as deviations from a constant reference temperature T ref . Namely, actual fluid and particle temperatures are given by (θ + T ref ) and (θ pj + T ref ), respectively, and T ref is set so that θ = 0 in a statistically steady state of turbulence, where the angle brackets indicate an ensemble average over turbulence realizations.

Characteristic time scales and non-dimensional numbers
The characteristic response time of the particle temperature to the surrounding fluid temperature field is given by the particle thermal relaxation time τ θ in (2.3). On the other hand, the characteristic response time of the fluid temperature field to the temperature of suspended particles is given by the fluid thermal relaxation time, where n p is the particle number density (τ θ f is referred to as the 'gas thermal relaxation time' in Pouransari & Mani (2018)). Here, τ θ f is inversely proportional to the particle number density n p and the particle radius r. This formula is similar to that for the phase relaxation time, an important time scale in cloud physics that characterizes the particle-turbulence interaction by condensation/evaporation (Politovich & Cooper 1988;Korolev & Mazin 2003;Kostinski 2009). Recently, a time scale of similar form was found for the particle-turbulence interaction by the drag force (Saito et al. 2019b).
There are three important non-dimensional numbers that characterize the two-way thermal interaction between particles and fluid. First, the particle-to-fluid heat capacity ratio is expressed as the ratio of τ θ to τ θ f as where C 0 and C p are the heat capacities of fluid and particles in the domain, respectively, φ m is the mass-loading parameter given as and φ v is the volume fraction. Second, the ratio of T to τ θ f , where T is the large-eddy turnover time for the turbulent velocity field, is represented by the Damköhler number as Here, T = L/u rms where L is the integral scale and u rms is the root-mean-square velocity (see Appendix D for definitions of other turbulence parameters). Third, the ratio of τ θ to T is represented by the thermal Stokes number for the large-scale flow as where the subscript L indicates that St L is that for the large-scale eddies of turbulent flow. From (3.4) and (3.5), φ θ is also expressed by Da and St L as Note that the large-eddy turnover time, T, is used in the definitions of (3.4) and (3.5) because the largest eddies are considered to be the most important in the thermal coupling. This conjecture was first proposed for the particle-turbulence interaction by condensation/evaporation and was confirmed by DNSs of cloud turbulence (Lanotte, Seminara & Toschi 2009;Sardina et al. 2015;Götzfried et al. 2017;Kumar et al. 2018;Saito et al. 2019a). Later this assumption was shown to also hold for the particle-turbulence interaction by the drag force (Saito et al. 2019b). Note also that the importance of the Damköhler number was recently demonstrated by Pouransari & Mani (2018) for the particle-to-fluid heat transfer in particle-laden turbulence. They referred to the inverse of the Damköhler number (Da −1 = τ θ f /T) as the heat-mixing parameter, as opposed to the commonly used heat-mixing parameter defined using the particle thermal relaxation time τ θ . In the present study, we use the term 'Damköhler number' for (3.4) in order to emphasize the relation with cloud-turbulence research. Originally, the Damköhler number came from chemical engineering and described the ratio of the turbulent flow time scale and the chemical reaction time scale (Dimotakis 2005). In the cloud physics literature, a similar parameter was introduced to examine the role of cloud entrainment and turbulent mixing in the growth of cloud droplets by condensation/evaporation. Since the growth of cloud droplets has some similarities to a chemical reaction, the Damköhler number term has been used (Shaw 2003;Devenish et al. 2012). Saito et al. (2019b) introduced another similar parameter for momentum coupling and referred to it also as the Damköhler number. In order to distinguish a variety of Damköhler numbers, it is useful to add an additional adjective before it, for example: the supersaturation Damköhler number for cloud physics (as used by Shaw 2003); the inertial Damköhler number for momentum coupling; the thermal Damköhler number for (3.4). In the present study, however, we simply refer to (3.4) as the Damköhler number, for brevity.

Modulation of fluid temperature fluctuations by particles
We apply the Langevin model in order to derive the analytical expression for modulation of fluid temperature fluctuations by particles. We assume that turbulence is statistically homogeneous, is in a statistically steady state and is excited by the external forces f u and f θ in (2.4) and (2.6), respectively, that are Gaussian random, white in time and homogeneous isotropic. We also assume that, both in cases with and without particles, the injection rate of the temperature-fluctuation variance is kept at χ 0 .
For cloud turbulence, the Langevin equation models the time variation of the supersaturation along the Lagrangian trajectory of cloud droplets. The feedback effect from cloud droplets on the supersaturation field is expressed by the term with the phase relaxation time τ c , and the turbulent mixing is expressed by the term with the large-eddy turnover time T. (For more details, readers are referred to, for example, Chandrakar et al. (2016).) In the present system, i.e. (2.1) to (2.7), the fluid temperature θ corresponds to the supersaturation, and the fluid thermal relaxation time τ θ f corresponds to τ c . Therefore, the Langevin model can be straightforwardly applied to the present system.
The situation under consideration is as follows. The number density of particles is n p and particles are uniformly randomly distributed in the domain. As described in Saito et al. (2019b), we divide the domain into small subdomains and assume that each subdomain includes only a single particle (Tagawa et al. 2012). The volume of the subdomain can be estimated approximately as 1/n p . The particle number density n p is assumed to be sufficiently large so that the characteristic length of the subdomain is smaller than the smallest length scale of the fluid temperature fluctuation (the Batchelor scale η B ) and that the fluid temperature in the subdomain can be regarded as uniform. In the subdomain, the fluid temperature is θ and the particle temperature is θ p . As time goes on, the fluid in the subdomain is advected by the turbulent velocity field as a small Lagrangian fluid parcel. Since we neglect the momentum inertia of the particle, the particle moves in the same way as the surrounding fluid parcel. From (2.6), the evolution equation for the temperature of the Lagrangian fluid parcel, θ, is modelled by the Langevin equation as follows: where θ 0 is the standard deviation of the fluid temperature fluctuation without the effect of particles, η(t) is a Gaussian random number with zero mean and unit variance which is statistically independent of θ and θ p , and α and c are constants that are determined depending on the actual situation of the turbulence. The equation for the particle temperature is the same as (2.2), i.e.
Since (3.7) and (3.8) are linear in θ and θ p , we can derive the analytical expressions for the statistically steady state. We first consider the particle temperature θ p . Multiplying (3.8) by θ p and taking an average, we have d dt where p indicates an average over particles. Since the left-hand side is zero in a statistically steady state, we have Multiplying (3.7) by θ p and taking an average p , we have since the feedback terms from particles cancel out due to (3.10) and η(t)θ p p = 0. Multiplying (3.8) by θ and taking an average p , we have Summing (3.11) and (3.12), we obtain d dt In a statistically steady state, we have (3.14) Combining (3.10) and (3.14) and using the definition of (3.5), we finally obtain This equation predicts the variance ratio of the fluctuations of particle and fluid temperatures. We next consider the fluid temperature θ. Here, we consider the dissipation of the temperature fluctuation variance. Since the injection rate is kept at a constant χ 0 with or without particles, for the case without particles, we have and, for the case with particles, we have Substituting (3.10) and (3.15) into the above equation and using the definitions of (3.2)-(3.6), we finally obtain ( 3.19) This equation predicts modulation of the fluid temperature fluctuation variance by particles as a function of the Damköhler number.
We consider two limits of (3.19). First, for the limit Da → 0 with fixed φ θ , we have which indicates that modulation of the variance approaches a constant depending on φ θ . In summary, in a statistically steady state, (3.15) predicts the variance ratio of the particle and fluid temperature fluctuations, and (3.19) predicts modulation of the fluid temperature fluctuation variance by particles.

Modulation by particles with fixed temperature
An important limit case of the Langevin model given by (3.7) and (3.8) is the case in which the particle temperature fluctuation θ p is fixed at zero. This physically corresponds to particles with large thermal inertia (St L → ∞). In this case, the Langevin model reduces to which is equivalent to the Langevin model used in cloud turbulence research (except for constants c and α -see (2) in Chandrakar et al. (2016)). In a statistically steady state, we obtain the same result as (3.20), which is the same formula as (3) in Chandrakar et al. (2016).

Particle temperature fluctuation for the one-way coupling case
Another important limit case of the Langevin model (3.7) and (3.8) is the case without the feedback effect of particles on the fluid temperature field, namely, the one-way coupling case. In this case, (3.7) reduces to and the equation for the particle temperature is the same as (3.8), as follows: In a statistically steady state, we obtain the same result as (3.15) for the particle temperature, as follows: whereas θ 2 p = θ 2 0 for the fluid temperature because the coupling is one-way. Note that a relationship similar to (3.25) was recently derived by Carbone et al. (2019) (see their (4.5)).

Simulation set-up
In order to confirm the theoretical prediction obtained in the previous section, we conduct DNS of particles in turbulence using the cloud microphysics simulator, which is a DNS model developed for the large-scale computation of particle-laden turbulent flow (Gotoh, Suehiro & Saito 2016;Saito & Gotoh 2018). The numerical domain is a periodic cubic box with a length L box per side. Particles are regarded as spherical point particles. The fluid velocity and temperature at the particle position, u(x pj ) and θ(x pj ), are calculated from the velocity and temperature fields at the surrounding eight grid points by linear interpolation, respectively. The same weighting is used to project the particle temperature, θ pj , onto the surrounding eight grid points to calculate the term F p in (2.7), i.e. the particle-in-cell method (Sundaram & Collins 1996).
We numerically integrate the evolution equations (2.1), (2.2), (2.4) and (2.6) using the second-order Runge-Kutta scheme and the pseudospectral method for temporal and spatial discretization, respectively. The external forces f u and f θ in (2.4) and (2.6) are homogeneous isotropic, applied to all wavenumber components within a shell of 1 ≤ kL box ≤ 4, and are produced by Gaussian random, white-in-time processes (Gotoh, Fukuyama & Nakano 2002). The details of the external forces are provided in Appendix A. These simple forces are enough for the fundamental study. After a statistically steady state is achieved for each simulation, turbulence and particle statistics are calculated.
Numerical and physical parameters in the present DNS are summarized in table 1. We studied five values of the total grid number (N 3 grid = 32 3 -512 3 ) and 10 values of the particle radius (r = 1.0-1000 μm). The total number of particles N p is set equal to the total grid number N 3 grid , corresponding to the particle number density n p = 125 cm −3 . For the particle-to-fluid heat capacity ratio, we chose three values (φ θ = 1.0, 0.3 and 0.1) and conducted two additional simulations. In one simulation, the particle temperature was fixed at zero (θ pj = 0) and, in the other simulation, the thermal coupling was one-way (the feedback from particles onto the fluid was neglected). In total, we conducted 5 × 10 × (3 + 2) = 250 kinds of simulations. Table 2 summarizes the particle radius r, the fluid thermal relaxation time τ θ f , the volume fraction φ v and the particle-to-fluid density ratio ρ p /ρ 0 for each simulation.
In table 1, the values of c 0 and c p are taken from Saito et al. (2019a) and Sato, Deutsch & Simonin (1998), respectively. From (3.2), the mass-loading parameters corresponding to φ θ = 1.0, 0.3 and 0.1 are φ m = 2.6, 0.78 and 0.26, respectively. Although φ m is an important parameter for turbulence modulation by particles due to the drag force (Squires & Eaton 1993;Saito et al. 2019b), φ m is not important in the present study because momentum coupling is neglected and the particles move as fluid particles.
The important premise in the Langevin model in the previous section is that the largest eddies play the central role in the particle-turbulence interaction by the two-way thermal coupling. We check the validity of this premise using the same strategy as that used by Lanotte et al. (2009) (see also Saito et al. (2019b)). Namely, we fix the turbulence statistics for the smallest eddies and expand the number of grids, the domain size and the sizes of the largest eddies in order to investigate the interaction between the largest eddies and particles. Turbulence statistics without particles are summarized in table 3. Note that fluid temperature statistics are modulated by the effect of particles.

Two limiting cases
We first investigate the results for the simulations in which the particle temperature is fixed at zero (θ pj = 0). Figure 1(a) shows the ratio θ 2 p /θ 2 0 , where θ 2 0 and θ 2 p are the variances of fluid temperature fluctuations without and with the effect of particles, respectively. The horizontal axis is Da, i.e. the Damköhler number in (3.4). The orange curve indicates the theoretical prediction (3.20) with c = 1.59 (c is estimated by curve-fitting method). All of the simulation results collapse nicely onto the theoretical curve. The slope of the curve is close to −1 for Da > 1, indicating that the scalar energy injected by the external force is dissipated mainly by particles in this range. We next investigate the results for the simulations with one-way thermal coupling. Figure 1(b) shows the ratio θ 2 p p / θ 2 p , where θ 2 p p is the fluctuation variance of the particle temperature. Note that θ 2 p = θ 2 0 for the one-way coupling case. The horizontal axis is St L , which is the thermal Stokes number in terms of the large-scale flow in (3.5). The orange curve indicates the theoretical prediction (3.25) with α = 1.31 (α is estimated by a curve-fitting method). Again, all of the simulation results collapse nicely onto the theoretical curve.

Two-way coupling cases
Now, we investigate the simulation results for the two-way thermal coupling case. The theoretical prediction (3.15) indicates that the variance ratio of the particle and fluid temperature fluctuations, θ 2 p p / θ 2 p , is the same as that for the one-way coupling case (3.25). This is confirmed by figures 2(a) to 2(c), which summarize θ 2 p p / θ 2 p for the simulations with φ θ = 1, 0.3 and 0.1, respectively. The orange curve indicates the theoretical curve (3.15) with the same value of α as in figure 1(b). Although there is slight dependence on the Reynolds number R λ (or the number of grids N grid ) and the heat capacity ratio φ θ , all of the simulation results collapse onto the theoretical curve fairly well. Figures 3(a) to 3(c) summarize the ratio θ 2 p /θ 2 0 , or the modulation of the fluid temperature variance by particles, based on the simulation results for the two-way coupling case with φ θ = 1, 0.3 and 0.1, respectively. For each panel, all simulation results collapse onto the theoretical curve (dashed grey curve) fairly well. The theoretical curve (dashed grey curve) follows the curve for the fixed-particle-temperature case (orange curve) for Da < φ θ (as shown in (3.20)) and starts to deviate from the fixed-particle-temperature case for Da ∼ φ θ and approaches a value of (1 + cαφ θ ) −1 , which is independent of Da for Da > φ θ (as shown in (3.21)). The physical interpretation of the results in figure 3 is as follows. For Da < φ θ , St L > 1 and τ θ > T from (3.4)-(3.6). The particle thermal relaxation time τ θ is greater than the large-eddy turnover time T, and the particle temperature responds only slowly to the change in the surrounding fluid temperature field due to turbulent mixing. As a result, the particle temperature does not deviate from the reference temperature so much, and the variance θ 2 p p is small in comparison with the fluid temperature variance θ 2 p (as shown in figure 2 for St L > 1). Therefore, modulation of the fluid temperature variance by particles is close to that for the case with a fixed particle temperature (orange curves 10 -1 10 0 10 1 10 2 10 -2 10 -3 10 -1 10 0 10 1 10 2 10 -2 10 -3 10 -1 10 0 10 1 10 2 Figure 2. Same as figure 1(b), except that the particle temperature is not fixed and the particle-to-fluid heat capacity ratios are φ θ = 1, 0.3 and 0.1 for (a-c), respectively. Note that the vertical axis is linear and the horizontal axis is from 10 −3 .
in figure 3, see (3.20)). On the other hand, for Da > φ θ , St L < 1 and τ θ < T. The particle temperature responds increasingly quickly to the change in the surrounding fluid temperature, and the temperature difference between the fluid and particles decreases. As a result, modulation of the fluid temperature variance by particles is limited and smaller than the case with a fixed particle temperature (see (3.21)). We emphasize the importance of using Da and St L , both of which include the large-eddy turnover time T in their definitions, in order to obtain the good collapse of the simulation results in figures 1-3. Such a good collapse cannot be obtained if we use a non-dimensional number defined with the Kolmogorov time τ . This is confirmed in Appendix B, where the thermal Stokes number St θ = τ θ /τ is used as an example. These results support that the largest eddies play the most important role in modulation of the fluid and particle temperature variances. Before closing this section, we briefly mention the time evolution and spatial distribution of fluctuations. Although the main focus of the present study is on the steady-state values of the fluid and particle temperature variances, we also took data of time evolution and spatial distribution of fluctuations. Some typical results are provided in Appendix C for the case with N grid = 512, φ θ = 1 and r = 20 μm.

Summary and discussion
In the present study, we investigated modulation of the fluid temperature fluctuations by particles due to thermal interaction in homogeneous isotropic turbulence. For the benefit of simplicity, only thermal coupling between the fluid and particle was considered, and momentum coupling was neglected. Particles were assumed to move as fluid particles and interact with the surrounding fluid only by heat transfer. Moreover, the fluid is assumed to be incompressible and to have a constant mass density. By virtue of these simplifications, we straightforwardly applied the Langevin model used in cloud turbulence research and obtained an analytical prediction for modulation of the intensity of the fluid temperature fluctuations by particles as a function of the Damköhler number, which is defined as the ratio of the turbulence large-eddy turnover time to the fluid thermal relaxation time. We also conducted DNS of the two-way thermal coupling between particles and turbulence and confirmed that the simulation results agree well with the theoretical predictions.
Modulation of fluid temperature fluctuations by particles due to heat transfer obtained in the present study is qualitatively very similar to turbulence modulation by particles due to the drag force obtained by DNS in Saito et al. (2019b) (compare their figure 3 with our figure 3). It should also be noted that, while Saito et al. (2019b) considered momentum coupling and hence included the effects of inertial clustering and associated phenomena, the present study does not include those effects since momentum coupling is neglected and particles are almost uniformly randomly distributed in the domain. Despite this difference, we obtained very similar results. This similarity strongly suggests that turbulence modulation due to momentum and thermal couplings share a similar underlying mechanism and also that the Langevin model used in cloud turbulence research should be a useful tool for clarification of turbulence modulation by particles. Since a theoretical prediction of modulation (such as (3.21)) has not yet been achieved for the case of the drag force, the application of the Langevin model might be a promising method. This possibility will be explored in a future study.
Since the present study assumed several simplifications of the particle-turbulence interaction, the interpretation of the present results needs some care. Especially, since particle momentum and thermal relaxation times are usually similar in natural and engineering applications, it is important to include momentum coupling. In such a case, the fluid velocity field is also modulated by particles due to the drag force, and the large-eddy turnover time T can no longer be taken as a constant parameter in the Langevin model. In addition, due to the momentum inertia of the particle, the trajectory of the particle deviates from that of the Lagrangian fluid particle. As a preliminary test, we conducted several additional simulations in Appendix E in order to investigate the effect of momentum coupling on modulation of fluid temperature fluctuations. The set-up is the same as that for DNSs with N grid = 128 and φ θ = 0.1 in § 4, except that either oneor two-way momentum coupling is included. For the one-way momentum coupling case, modulation of the fluid temperature variance is almost unchanged as compared with the case without momentum inertia. On the other hand, for the two-way momentum coupling case, the fluid temperature variance is significantly increased by particles as compared with the case without momentum inertia. However, this is only due to attenuation of turbulent kinetic energy by particles which in turn leads to the increase in the large-eddy turn over time, and not due to the effects of inertial clustering and associated phenomena. Appendix E also confirms that the theoretical prediction (3.19) is still useful for the study of cases with momentum coupling. Clearly, more comprehensive simulations are required to be conclusive, but this is left as future work.
Another important task is to clarify the effect of the forcing mechanism used for the injection of kinetic energy and scalar variance in DNS. Although the present DNS used 931 A6-15 a Gaussian white-in-time random force both for velocity and scalar fields, many other forcing mechanisms have been proposed and one of them is the forcing based on a uniform mean scalar gradient: f θ = Γ u 3 , where u 3 is the vertical component of u and Γ is a constant. This form is commonly used in DNSs of passive scalars (Gotoh & Watanabe (2015), Yasuda et al. (2020) and references therein) and is also relevant to cloud physics where the temperature and supersaturation fluctuations are excited by the vertical air motion in the stratified environment (Thomas et al. 2020;Grabowski & Thomas 2021). From our recent analysis of the Langevin model in cloud physics (Saito, Watanabe & Gotoh 2021), we expect that a theoretical prediction of modulation similar to (3.21) can be obtained for the case with this type of scalar forcing. We are aware that the difference in the forcing mechanism substantially affects turbulence modulation by particles. This will be reported in a subsequent study.
Declaration of interests. The authors report no conflict of interest.

Appendix A. Details of the forces
The velocity forcing f u satisfies where k is the wavenumber vector (k = |k|), P is the projection operator (P ij = δ ij − k i k j /k 2 ), and the spectrum of the force F(k) is constant (= c f ) for 1 ≤ kL box ≤ 4 and zero otherwise. Similarly, the temperature forcing f θ satisfies where F θ (k) is constant (= c θ ) for 1 ≤ kL box ≤ 4, and zero otherwise. For each set of simulations with the same total grid number (N 3 grid = 32 3 -512 3 ), we first conducted simulations without particles and tuned two parameters c f and c θ so that the dissipation rates of kinetic energy and temperature variance at the statistically steady state are and χ in table 1, respectively. We then used the same values of c f and c θ for other simulations with particles. The above kind of random force has been commonly used in simulation studies of turbulence with the assumption that, if the force is confined to sufficiently large scales, the details of the force do not affect the statistical nature of the small-scale dynamics.

Appendix B. On the importance of small-scale eddies
In the present study, we mainly focused on the role of the large-scale eddies and used the large-eddy turnover time T for the non-dimensional numbers as in (3.4) and (3.5). On the other hand, it is also interesting to examine the role of the small-scale eddies, since their importance in particle-turbulence interaction has been demonstrated by the previous studies (Yoshimoto & Goto 2007).
In order to examine the importance of the small-scale eddies, we use the thermal Stokes Here, we show several typical examples of the time evolution of variances. We also show an example of a snapshot of the fluid temperature field. We consider the case with N grid = 512, φ θ = 1 and r = 20 μm in tables 1 and 3. Figure 5(a) shows the time evolution of the turbulent kinetic energy. The statistically steady state is reached around t ∼ 15 s. The large-eddy turnover time at the statistically steady state is T = 2.2 s for this case. Figure 5(b) compares the time evolution of the fluid temperature variances, θ 2 p , for three cases: two-way coupling (black dots); one-way coupling (red triangles); and θ p fixed at zero (blue squares). For all three cases, the variance first decays almost exponentially and reaches the statistically steady state. The time scale for this exponential decay is close to α −1 T for the one-way coupling case (compare with the dashed line), while it is close to [α(c/τ θ f + 1/T)] −1 for the case with θ p fixed at zero (compare with the solid line). These decay rates are consistent with the predictions from the Langevin model ((3.22) and (3.23)). On the other hand, for the two-way coupling case, θ 2 p decays in a similar manner to the one-way coupling case but with smaller amplitude. Figure 5(c) compares the time evolution of the particle temperature variances, θ 2 p p , for the two-way and one-way coupling cases. For both cases, the time scale for the initial exponential decay is close to τ θ (compare with the solid line). Figure 5(d) shows a two-dimensional snapshot of the fluid temperature field at t = 50 s for the two-way coupling case. In this case, the fluid temperature variance is damped by approximately 50 % due to particles.

Appendix D. Definitions of turbulence parameters
The kinetic energy and temperature variance are defined by where u rms and θ rms are the root-mean-square velocity and temperature fluctuations, respectively, and E(k) and E θ (k) are the spectra of the kinetic energy and temperature variance, respectively. The mean dissipation rates for E(k) and E θ (k) are defined by The Kolmogorov and Batchelor length scales are defined by The integral scale and Taylor microscale are defined by The large-eddy turnover time and Kolmogorov time are defined by The Taylor microscale Reynolds number is defined by (a) ( b) ( c) Figure 6. Particle distribution obtained by a simulation with two-way momentum coupling. Particle radius is 1000, 60 and 5 μm from (a) to (c), respectively. Only particles within a slice − x/2 < z < − x/2 are shown. spatial distribution of particles deviates from a uniform-random distribution, especially when the Stokes number St = τ p /τ is close to unity. For simulation with two-way momentum coupling, both fluid velocity and temperature fields are modulated by particles. However, we here focus on modulation of fluid temperature fluctuations. Various parameters at the statistically steady state of turbulence for additional simulations are summarized in table 4. Both one-and two-way momentum coupling cases include runs with the Stokes number ranging from St 1 to St 1. Figure 6 shows the snapshots of particle distribution for the two-way momentum coupling case. Due to momentum coupling, inertial clustering of particles is clearly observed for St = 1.1 (figure 6b, r = 60 μm). Similar particle distributions were obtained for the one-way momentum coupling case. Figure 7(a) shows the ratio θ 2 V /θ 2 0 at the statistically steady state. Note that, for cases with momentum coupling, the average over particles ( p ) is generally different from that over fluid volume ( V ). Since our focus is on modulation of fluid temperature fluctuations, we use the average over fluid volume. The results for the one-way momentum coupling case (red triangles) are very similar to those for the case without momentum inertia (blue squares). On the other hand, the results for the two-way momentum coupling case (black dots) show significantly greater values than the other two cases. Several black dots are even (c) but for one-and two-way momentum coupling cases (red triangles and black dots, respectively) with N grid = 128 and φ θ = 0.1. Note also that the vertical axis is θ 2 V /θ 2 0 instead of θ 2 p /θ 2 0 . Blue squares are the same as in figure 3(c) where the momentum inertia of the particle is neglected. (b) Same as (a) for the two-way momentum coupling case except that the vertical axis is θ 2 V /θ 2 c , where θ 2 c is defined by (E7). See discussion in the main text.
greater than unity, indicating that the fluid temperature variance is increased by particles for these runs. The increase in the fluid temperature variance for the two-way momentum coupling case is caused by attenuation of turbulent kinetic energy by particles. In scalar turbulence, the cascade of scalar variance to smaller scales is mainly driven by the large-scale eddies, and its time scale is roughly estimated by the large-eddy turnover time T. Since the injection rate is kept at a constant χ 0 in the present study, the fluid temperature variance without particles, θ 2 0 , is determined by the scaling relationship (3.16). On the other hand, when the turbulent kinetic energy is attenuated by particles due to two-way momentum coupling, the time scale of the cascade becomes longer, which leads to an increase in the fluid temperature variance. This increase can be estimated by the following scaling argument. If modulation of the fluid temperature variance is solely caused by the change in the large-eddy turnover time (and not by two-way thermal coupling), the variance is expected to become θ 2 c which is determined by where T 1 and T are the large-eddy turnover time of turbulence with and without the effect of two-way momentum coupling, respectively (T 1 > T). By equating (E6) with (3.16), we obtain the relationship between θ 2 c and θ 2 0 as follows: which indicates θ 2 c > θ 2 0 since T 1 > T. This nominal value θ 2 c can be used to estimate the relative effect of turbulence attenuation on modulation of the fluid temperature variance. Figure 7(b) shows the same as figure 7(a) for the two-way momentum coupling case, except that the vertical axis is replaced by the ratio θ 2 V /θ 2 c . By using θ 2 c instead of θ 2 0 , we removed the contribution of the increase in the large-eddy turnover time to modulation of the fluid temperature variance. All results in figure 7(b) are now smaller than unity and agree well with the theoretical curve (grey dashed curve), as for the one-way momentum coupling case and the case without momentum inertia in figure 7(a) (red triangles and blue squares). These results indicate that the increase in the variance for the two-way momentum coupling case shown in figure 7(a) (black dots) is mostly caused by the attenuation of turbulent kinetic energy by particles which in turn leads to the increase in the large-eddy turnover time, and also that inertial clustering and associated phenomena are not so important in the modulation of the fluid temperature variance. In addition, these results also confirm that the theoretical prediction (3.19) (grey dashed curves in figure 7) is still useful for the study of cases with momentum coupling.