Radiation transport plays an important role in stellar atmospheres, but the effects of turbulence are being obscured by other effects such as stratification. Using radiative hydrodynamic simulations of forced turbulence, we determine the decay rates of sinusoidal large-scale temperature perturbations of different wavenumbers in the optically thick and thin regimes. Increasing the wavenumber increases the rate of decay in both regimes, but this effect is much weaker than for the usual turbulent diffusion of passive scalars, where the increase is quadratic for small wavenumbers. The turbulent decay is well described by an enhanced Newtonian cooling process in the optically thin limit, which is found to show a weak increase proportional to the square root of the wavenumber. In the optically thick limit, the increase in turbulent decay is somewhat steeper for wavenumbers below the energy-carrying wavenumber of the turbulence, but levels off toward larger wavenumbers. In the presence of turbulence, the typical cooling time is comparable to the turbulent turnover time. We observe that the temperature takes a long time to reach equilibrium in both the optically thin and thick cases, but in the former, the temperature retains smaller scale structures for longer.

An important property of turbulence is the mixing of fields that are advected by the flow. The simplest example is that of a passive scalar, a quantity that does not backreact on the flow. The magnetic field is another popular example, because for weak field strengths, it can be treated as a passive vector field, making the mathematics more straightforward compared to the fully nonlinear case. Even the flow itself is mixed by the turbulence, which is a much harder problem. This leads to turbulent viscosity, which acts as an enhanced molecular viscosity, although there can be additional important effects if the turbulence is anisotropic. Examples of additional effects occur in stratified flows in the presence of rotation. Such flows can become differentially rotating through what is called the Λ effect.1 It is a nondiffusive effect, analogous to the α effect in mean-field dynamo theory.2,3 These nondiffusive effects have led to significant attention in astrophysics. Scalars, active or passive, have received comparatively less attention, because nondiffusive effects are generally less profound, but see Rädler et al.4 for the slow-down of turbulent diffusion in certain compressive flows.

Prandtl suggested that turbulence has a smoothing effect—just like molecular diffusion. The molecular diffusion coefficient is generally proportional to the product of the typical velocity of the molecules, which is essentially the sound speed, and the typical mean-free path between collisions. Prandtl generalized this to turbulence by using the product of the typical velocity of the turbulent eddies and their correlation length, which he referred to as the mixing length. Important applications of turbulent mixing in astrophysics include turbulent convection in the Sun and stars, as well as mixing of chemicals in the galaxy. The latter is a typical case of a passive scalar, while in the former case, the quantity that is being mixed is the specific entropy, which is an active scalar, because it affects the density in the momentum equation and can lead to buoyancy. Furthermore, the resulting turbulent diffusion is an enhancement not of molecular diffusion, but of photon diffusion, which is also referred to as radiative diffusion.

Radiative diffusion comes in two different forms: optically thick and optically thin. Optically thick is the usual case, where the mean-free path of photons is short compared to the typical scales of the flow. Optically thin, by contrast, means that the photons can propagate over large distances before they are absorbed and re-emitted again. Radiative diffusion ceases to exist in this case and we have to deal instead with an essentially nonlocal process. The effect of radiation now decreases with increasing mean-free path of the photons, contrary to the diffusive case where it increases. The relevant process, in this case, is Newtonian cooling,5 where the cooling rate is directly proportional to the radiation energy density rather than the divergence of its gradient. The cooling timescale is the ratio of mean-free path to some relevant photon speed, rather than turbulent diffusion, whose coefficient is proportional to the product of mean-free path and the relevant photon speed.

In astrophysics, one usually thinks of optically thin processes being those that happen above the photosphere of a star, where photons can travel all the way to infinity. However, even below the photosphere, a process can be optically thin if we look at small length scales, because then the photon mean-free path again exceeds the relevant scale of the flow structures. In this paper, we are interested in the effects of turbulence, especially in this limit.

There is actually a curious analogy between the optically thin limit, where cooling becomes less efficient at small length scales,6 and turbulent diffusion, which also becomes less efficient at small length scales,7,8 although this concerns here the length scales of the mean fields. This is because turbulent diffusivity is not just a coefficient, but an integral kernel in a convolution with the mean temperature.9 The Fourier transformation of this kernel falls off with wavenumber approximately like a Lorentzian, which is analogous to the case of radiative transfer in the optically thin case.10 We may therefore ask: how does the combined effect of turbulence and small optical thickness modify turbulent diffusion at small length scales?

We begin by reviewing some basics about the cooling time as a function of mean-free path in Sec. II, following which we describe our numerical simulations in Sec. III, and finally present our results in Sec. IV. Our conclusions and scope for future work are given in Sec. V.

In compressible hydrodynamics, the energy equation can be written in terms of the specific entropy s(x,t) as

ρTDsDt=·Frad+2ρνS2,
(1)

where ρ is the density, T the temperature, D/Dt=/t+u· the advective derivative, Frad the radiative flux, ν the viscosity, and S the traceless rate-of-strain tensor with the components Sij=(iuj+jui)/2δij·u/3. The negative divergence of Frad is calculated as the imbalance of the intensity and the source function integrated over all frequencies ν̃ and all directions,11 i.e.,

·Frad=0κνρ4π(Iν̃Sν̃)dΩdν̃,
(2)

where κν̃ is the opacity per unit mass, Iν̃(x,t,n̂) is the specific intensity corresponding to the energy that is carried by radiation per unit area, per unit time, in the direction n̂, per unit solid angle dΩ, and Sν̃(x,t) is the source function. Throughout this work, we make the gray approximation and thus work with frequency-integrated quantities, which amounts to dropping the subscript ν̃. In the gray approximation, I(x,t,n̂) obeys the radiative transfer equation,

n̂·I=κρ(IS),
(3)

which is solved along a set of rays in different directions n̂ using the method of long characteristics.12 The source function is here written as S=(σSB/π)T4, where σSB is the Stefan–Boltzmann constant. The photon mean-free path is =(κρ)1 and κ is a suitably averaged “gray” opacity.

By assuming infinitesimally small temperature perturbations, Spiegel13 linearized Eqs. (1) and (3) and found that for perturbations of the form proportional to exp(ik·x+λt), the inverse relaxation (or cooling) time λ of the temperature perturbations, or decay rate, is given by

λ=cγ(1arctan kk),
(4)

where k=|k| is the wavenumber,

cγ=16σSBT3/ρcp
(5)

is the characteristic velocity of photon diffusion,14 and cp the specific heat at constant pressure. The dependence λ(k) is what we call in this paper the cooling curve. A useful form of the above expression can be obtained under the Eddington approximation,15 where one expands Eq. (3) in terms of moments of n̂ under the closure assumption that n̂in̂jIdΩ=13δijIdΩ. This yields a closed equation 13()2J=JS for the mean intensity J=IdΩ/4π. The cooling rate is then

λcγk22/31+k22/3.
(6)

It is convenient to introduce now the radiative diffusivity, χ=cγ/3, which also clarifies why cγ is called the characteristic velocity of photon diffusion.

Equations (4)–(6) apply to the case of three-dimensional (3-D) variations of the temperature. In our 3-D numerical experiments, however, we restrict ourselves to one-dimensional (1-D) variations of the mean temperature profile. In that case, the relevant version of Eq. (6) becomes14 

λcγk22/31+k22(1-Dperturbations).
(7)

The corresponding version of Eq. (4) then takes the form

λ=cγ(1arctan3k3k)(1-Dperturbations).
(8)

In Fig. 1, we compare λ(k) obtained from the exact equation (red curve) with the approximate λ(k) obtained under the Eddington approximation (black curve) for the relevant 1-D Eqs. (7) and (8). Our numerical solution for λ, which is based on only six rays, depends on the choice of weight factors used in the angular integration. The weight factors have been chosen such that our numerical results (plus signs) agree with the Eddington approximated solution.14 The basic question we want to answer is how the cooling curve gets modified in the presence of turbulence. We expect the effective λ to be enhanced, atleast in the optically thick limit, where k1; however, we do not know what to expect in the optically thin case, where k1, and how it depends on the scale of the turbulent eddies. To address these questions, we now perform turbulence simulations. We are particularly interested in the regime of moderate temperatures, where the radiation pressure can be ignored.

FIG. 1.

The analytic cooling curve (black), the cooling curve under the Eddington approximation (red), and the results from numerical simulations (plus signs).

FIG. 1.

The analytic cooling curve (black), the cooling curve under the Eddington approximation (red), and the results from numerical simulations (plus signs).

Close modal

The conditions in the Sun are extremely inhomogeneous owing to tremendous stratification. The density changes by nearly six orders of magnitude across the convection zone and the temperature by a factor of about 300. This fact alone can introduce new phenomena such as the spontaneous formation of magnetic flux concentrations.16,17 At a more elementary level, the addition of gravity leads to convection and thereby to turbulent motions, which have been the subject of numerous simulations for a long time.11,18 Newtonian cooling also plays important role in the atmospheres of planets,19 where turbulence is not always explicitly invoked and therefore the role of turbulence needs to be parameterized.20 

In the Sun, the microphysical viscosity is about twelve orders of magnitude smaller than the estimated turbulent viscosity, and over eight orders of magnitude smaller than the radiative diffusivity near the surface. Numerical simulations have therefore routinely employed numerical tools that allow the simulations to proceed by dissipating sufficient energy in local regions where necessary. This precluded the study of turbulent Newtonian cooling, because the small-scale turbulent motions have already been altered by such numerical modeling.21 An additional complication is partial ionization, which tends to make the transition from the deeper optically thick layers to the surface very sharp in a stratified system. However, one could imagine it to introduce new effects of its own if we arranged the average temperature of the domain such that it lies exactly in the middle between those of a fully ionized and a neutral medium, as has been done in some other recent experiments.22 In those cases, however, Newtonian cooling does not necessarily play any obvious role.

To accomplish our goal of identifying turbulent effects in optically thin and thick turbulent flows, we resort to the study of a minimal system where isotropic homogeneous turbulence is produced by a stirring force instead of convection. The viscosity is kept constant, but we consider different values to assess the dependence of our results on the Reynolds number. Partial ionization effects are ignored and other complications from adopting a realistic equation of state are not included. We also restrict our attention to the study of vortical forcing. The concept of turbulent mixing is likely to be similar also for irrotational forcing, but other poorly understood features of such turbulence such as a pileup of kinetic energy near the dissipative subrange (bottleneck effect) are known to occur in such cases.23 It may play a role in interstellar turbulence,24 which motivates a more extended future study of turbulent Newtonian cooling for irrotational forcing.

For the purpose of our present study, we restrict ourselves to studying a turbulent flow in a triply periodic domain of size L3 by applying plane wave forcing throughout the domain. We therefore solve the equations,

ρDuDt=p+ρf+·(2ρνS),
(9)
DlnρDt=·u,
(10)

where p is the pressure, u the velocity, and f the forcing function. In Eq. (9), we have ignored the radiation force (ρκ/c)Frad, where c is the speed of light, as mentioned above. This term is unimportant for the temperatures considered in this work. Nevertheless, the coupled set of equations (1), (9), and (10) makes s(x,t) an active scalar, because it is related to p and ρ through

s=cvlnpcplnρ+s0,
(11)

where cv is the specific heat at constant volume and s0 is an irrelevant constant. This equation follows from the first law of thermodynamics, written in the form Tds=de+pd(ρ1), where e=cvT is the internal energy for a perfect gas, and the ideal gas equation relating the temperature to p and ρ through

(cpcv)T=p/ρ.
(12)

We then have ds=cvdlnT(cpcv)dlnρ. Using the differentiated ideal gas equation, dlnT+dlnρ=dlnp, we arrive at Eq. (11) after integration. For the ratio of specific heats, γ=cp/cv, we assume γ=5/3, which is appropriate for a monatomic gas such as fully ionized hydrogen at the temperatures considered here (T40000K).

In our numerical work, we use dimensionful quantities, where length is measured in megameters (Mm), speed in kms1, and temperature in kelvin. We also use the symbol ad=11/γ=0.4, which is the adiabatic value of what is in astrophysics commonly referred to as the double logarithmic temperature gradient, dlnT/dlnp.25 Using this, cp can then be written as cp=R/(μad), where we have used cpcv=R/μ, with R=8.314×107cm2s2K1 being the universal gas constant and μ=0.6 the mean molecular weight. We then find cp=0.035km2s2K1.

To simulate a turbulent flow, we apply nearly monochromatic forcing with an average forcing wavenumber kf. The forcing function changes abruptly from one time step to the next, i.e., f(x,t)f(x,t) is proportional to δ(tt), where δ is the Dirac δ function. The forcing is then said to be δ correlated in time. The smallest wavenumber in the cubic domain of side length L is k1=2π/L. The ratio kf/k1 is the scale separation ratio, for which we consider the values 1.5 and 10.

For the forcing function f, we select randomly at each time step a phase π<φπ and the components of the wavevector k from many possible discrete wavevectors with lengths in a certain range around a given value kf. In this way, the adopted forcing function,

f(x,t)=Re{Nf̃(k,t)exp[ik·x+iφ]},
(13)

is white noise in time and consists of plane waves with average wavenumber kf. Here, x is the position vector and N=f0(cs0kfδt)1/2 is a normalization factor, where δt is the time step and f0 is an amplitude factor. In this formulation, the averaged forcing is independent of δt. To ensure that f̃ is solenoidal, i.e., perpendicular to k, we write is as

f̃(k)=(k×ê)/[k2(k·ê)2]1/2,
(14)

where ê is an arbitrary unit vector that is not aligned with k. Note that |f|2=1km4s4Mm2. The coefficient f0 is chosen such that the velocity is about 10% of the sound speed.

We adopt a sinusoidal temperature perturbation and write T(x,t=0) in the form

T(x,t=0)=T0+T1sinkx,
(15)

where k=k1=1Mm1 is chosen. This implies that L=2πMm6.28Mm. We choose cs0=30kms1, so that T040000K and cγ=3.9kms1. The temperature perturbation is taken to be T1=2000K, and periodic boundary conditions are assumed for all quantities.

We define the Mach and Reynolds numbers as

Ma=urms/cs0,Re=urms/νkf.
(16)

For a forcing amplitude f0=0.01km2s2Mm1, we have urms2.2kms1, so that Ma0.08. Using ν=103Mmkms1 and kf=10Mm1, we have Re230, while for kf=1.5Mm1, we have Re1200. To determine the microphysical Prandtl number in the optically thick regime, we define

Pr0=3νk1/cγ,
(17)

so that the Prandtl number is ν/χ=Pr0/(k1). Here we have used χ=cγ/3 for the radiative diffusivity. Following earlier work,14 we choose for the mean density the value ρ0=4×104gcm3.

We determine the effective λ from the time evolution of the decay of the sinusoidal perturbation, which we monitor by taking the difference between the maximum and minimum temperatures at each time. This turns out to be reasonably accurate and we use a time interval during which the decay is exponential.

We perform numerical simulations with the Pencil Code (https://github.com/pencil-code), which is a public MHD code that is particularly well suited for simulating turbulence.26 We solve Eqs. (1), (9), and (10) with sixth-order finite differences.27 Equation (3) is solved with second-order accurate finite differences along the coordinate directions and the diagonals, i.e., altogether 26 rays. The radiation transport has been parallelized in the Pencil Code by splitting the calculation into parts that are local and nonlocal with respect to each processor.28 Two parts are compute-intensive but require no communication, and one part is nonlocal but does not require waiting for any computation to be done and is therefore fast. We use the third-order time-stepping scheme of Williamson.29 

The code's local cooling and heating properties have been verified,14 and its cooling time has been compared with the analytic cooling time obtained by Spiegel.13 It is this cooling time that determines the relevant time step constraint for radiation simulations,6,30 and not some generalization of the usual Courant condition.31 The latter would erroneously imply a limiting time step that is proportional to the mesh spacing, when it is actually quadratic in the mesh spacing in the optically thick limit and independent of mesh spacing in the optically thin limit.

The code has been applied to sunspot simulations,32 and to a range of more idealized problems of atmospheric stability14 and magnetic spot formation.17,22 It should be emphasized that our calculations classify as direct numerical simulations in the sense that the equations are solved as stated, albeit with unrealistically large viscosity and unrealistically small opacities compared with solar conditions. By comparison, most simulations of solar convection are performed using large eddy simulations.18,30,33

The fact that the Pencil Code uses sixth-order finite differences and a third-order time-stepping scheme does not tell us much about the actual accuracy and convergence of our results. For example, the longer a simulation, the more numerical errors should accumulate, but this is not normally seen. This was recently addressed in a study comparing the numerical accuracy of turbulence and waves.34 In that study, it was concluded that the existence of a forward cascade in turbulence prevents the systematic loss of energy at small scales, where discretization errors are the largest. This is not the case for waves, which therefore need to be solved with much more care. Another such example was a recent three-dimensional study of electromagnetic waves.35 In this and the previous case, it proved advantageous to use an exact time integration under the assumption that the turbulent source is unchanged between two time steps. This approximation is justified because at high wavenumbers, the relevant timescale of waves is much shorter than that of turbulence.34 

Furthermore, radiation can introduce short time scales. As already alluded to in the beginning, the cooling time can be very short and severely restrict the relevant time step constraint.6,30 This makes the simulations very costly,31 but we are not aware of any reports on loss of accuracy in such cases. Furthermore, in the present studies, we are probably not affected by this constraint, because turbulent Newtonian cooling only plays a role when the turbulent time scales are short compared with radiative ones. This is here not the case.

In this section, we present the results for λ obtained using various values of the forcing wavenumber kf and the wavenumber k of the initial perturbation. The Reynolds number varies between 20 and 1200, and the number of mesh points, N3, is varied between 643 and 2563; see Table I. The values of Pr0 are then small, as is also expected for the Sun, and they are 8×103 for our runs with 643 mesh points (small Reynolds numbers) and 8×104 for 2563 mesh points (larger Reynolds numbers). For each series of runs, we perform simulations where we vary the opacity κ and thereby ; see Fig. 2 for visualizations of T on the periphery of the computational domain for runs of Series A for k=0.1, 1, and 10 and at different times (in kiloseconds [ks]). We see that the large-scale temperature contrast (wavenumber k = k1) decreases the fastest for k=1, and more slowly for k=0.1 and 10. For k=10, however, which is the optically thin case, the temperature retains smaller scale structures for longer.

TABLE I.

Summary of the parameters for the series of runs presented in this paper. Within each series of runs, κ and are varied.

Seriesk/k1kf/k1urms/cγPr0ReN3
A′ 10 0.50 8×103 20 643 
10 0.59 8×104 230 2563 
1.5 0.62 8×103 160 643 
1.5 0.46 8×104 1200 2563 
1.5 0.46 8×104 1200 2563 
Seriesk/k1kf/k1urms/cγPr0ReN3
A′ 10 0.50 8×103 20 643 
10 0.59 8×104 230 2563 
1.5 0.62 8×103 160 643 
1.5 0.46 8×104 1200 2563 
1.5 0.46 8×104 1200 2563 
FIG. 2.

Temperature on the periphery of the computational domain for k=0.1 (left column), k=1 (middle column), and k=10 (right column), for t=0.3ks,t=1ks,t=2ks, and t=3ks (top to bottom). The (x, y, z) coordinates are indicated in the first panel and the approximate 1Mm scale is shown in the second panel.

FIG. 2.

Temperature on the periphery of the computational domain for k=0.1 (left column), k=1 (middle column), and k=10 (right column), for t=0.3ks,t=1ks,t=2ks, and t=3ks (top to bottom). The (x, y, z) coordinates are indicated in the first panel and the approximate 1Mm scale is shown in the second panel.

Close modal

In spite of the forcing being monochromatic, the resulting turbulence is excited over a broad range of scales. This is demonstrated in Fig. 3, where we plot kinetic energy spectra, EK(k)=|ũ|2k2dΩk, for Series A and C. Here, ũ is the Fourier transformation of u, and dΩk is the solid angle differential in wavenumber space. The spectra are normalized such that EK(k)dk=urms2/2. In the case of Series C, where Re=1200 and kf/k1=1.5, there is a short inertial range k5/3 together with a bottleneck, i.e., a shallower spectrum near the dissipative subrange.36,37 We note that the bottleneck effect is physical, but much weaker in the one-dimensional spectra that are accessible to laboratory and atmospheric measurements.38 It is also seen in the highest resolution turbulence simulations today.39 

FIG. 3.

Kinetic energy spectra for Series A (dashed blue) and Series C (solid red). The k5/3 slope is overplotted for comparison.

FIG. 3.

Kinetic energy spectra for Series A (dashed blue) and Series C (solid red). The k5/3 slope is overplotted for comparison.

Close modal

In the simulations with larger scale separation (Series A), however, the spectrum is more peaked around k=kf. This occurrence of this spike at kf is partially explained by the smaller Reynolds number (Re=230 in this case).

In general, higher scale separation allows us to see more clearly the various mean-field effects. In this connection, we must remember that the standard concept of turbulent diffusion does require sufficient scale separation and that the lack of scale separation requires one to study the full scale dependence, in which case turbulent diffusion corresponds to an integral kernel in real space, or a multiplication with a k-dependent diffusivity in Fourier space.7 For this reason, we also discuss the aspect of scale dependence below.

In Fig. 4, we plot λ vs k and compare with the laminar case shown in Fig. 1. In all the cases, we see that λ is enhanced relative to the laminar curve. Varying the viscosity, and thereby changing Re from 20 (Series A′) to 230 (Series A), has a very minor effect; compare the dotted and solid blue lines for k/kf=0.1 in Fig. 4. Decreasing kf/k1 from 10 to 1.5, that is, increasing k/kf from 0.1 (Series A) to 0.7 (Series B), has a more significant effect, and λ is seen to increase by a factor that is between 4 and 8, depending on the value of k.

FIG. 4.

Effective turbulent decay rate λ vs k for Series A′ (dotted blue), A (solid blue), B (solid orange), C (dotted orange), and D (dashed orange). The short lines on the left y axis, with line types matching those of the curves in the plot, give the theoretical expectations explained in the text. The curve for the laminar case is shown in red.

FIG. 4.

Effective turbulent decay rate λ vs k for Series A′ (dotted blue), A (solid blue), B (solid orange), C (dotted orange), and D (dashed orange). The short lines on the left y axis, with line types matching those of the curves in the plot, give the theoretical expectations explained in the text. The curve for the laminar case is shown in red.

Close modal

Keeping the value of kf unchanged and increasing k/kf (for Series C and D), i.e., making the scale separation poorer, results in a weak decline of λ. Theoretically, we would expect the turbulent decay rate to be λ=χtk2, where χt=χt0urms/3kf is the nominal turbulent diffusivity in the case of perfect scale separation. For poor scale separation, however, we expect χt=χt0/[1+(k/kf)2]. We see from the short lines overplotted on the left y axis of Fig. 4 that the actual decay rates are somewhat larger.

We note that in Fig. 4, the decay rates of lines having different values of k (but the same value of kf) are all separated by factors that are close to k itself. To demonstrate that this is mostly the result of normalizing λ/cγ by k, we show in Fig. 5 the result of normalizing λ/cγ by k1, which is the same for all runs. The lines are now no longer so strongly separated for different values of k. Note that the abscissa is also scaled by k1 instead of k. Consequently, the small peaks in λ values near k1=1 occur at similar positions.

FIG. 5.

Similar to Fig. 4, but the abscissa is scaled with k1 instead of k, and the ordinate is divided by k1 instead of k. Again, Series A′ and A are denoted by dotted and solid blue lines, respectively, and Series B, C, and D are denoted by solid, dotted, and dashed orange lines, respectively. The curve for the laminar case is shown again in red.

FIG. 5.

Similar to Fig. 4, but the abscissa is scaled with k1 instead of k, and the ordinate is divided by k1 instead of k. Again, Series A′ and A are denoted by dotted and solid blue lines, respectively, and Series B, C, and D are denoted by solid, dotted, and dashed orange lines, respectively. The curve for the laminar case is shown again in red.

Close modal

The standard concept of turbulent diffusion with a diffusion operator of the form χt2 requires one to have sufficient scale separation, as is the case for our runs of Series A. If scale separation is poor, the operator χt2 has to be replaced by a convolution in real space.7 This subject continues to attract significant attention, especially in plasma physics40 and astrophysics.41,42

In Fig. 6, we summarize the results for λ/(cγk1/3) as a function of k/kf for k=0.01 (optically thick regime), and k=100 (optically thin regime). Neither of the two regimes exhibits a k2 dependence, as would be expected for a turbulent diffusion process with k/kf1, i.e., when the turbulent diffusivity is approximately scale-independent. In the optically thick case, λ increases approximately linearly with k for small values of k and then reaches a maximum. In the optically thin case, on the other hand, λ increases with k approximately like k1/2.

FIG. 6.

Dependence of λ/(cγk1/3) on k/kf for k=0.01 (blue, optically thick), and k=100 (red, optically thin). The later obeys a k1/2 scaling. For comparison, we also show the linear scaling in k (blue) and the quadratic scaling (green).

FIG. 6.

Dependence of λ/(cγk1/3) on k/kf for k=0.01 (blue, optically thick), and k=100 (red, optically thin). The later obeys a k1/2 scaling. For comparison, we also show the linear scaling in k (blue) and the quadratic scaling (green).

Close modal

In the optically thick regime, k1, where the radiative diffusion approximation should be applicable, we might expect some analogy between turbulent diffusion of active scalars (such as temperature) and passive scalars (such as chemical concentrations). For the latter, the scale dependence has previously been investigated,8 and it was found to be similar to that for magnetic fields7 in that both had a Lorentzian shape. In these papers, the microphysical and turbulent diffusivities were referred to as κ and κt for the passive scalar8 (not to be confused with the opacity in the present paper) and η and ηt for the magnetic diffusivity.7 They have the same meaning as χ and χt in the present paper. In these cases, we plot the timescale on which a large-scale sinusoidal profile of the passive scalar or the magnetic field gets diffused away.

We reproduce in Fig. 7 the result from the test-field method for passive scalars.8 These authors also studied the effects of rotation and magnetic fields, but those results were not used for the present comparison. Their passive scalar diffusivity κt obeyed a Lorentzian fit such that

κt(k)=urms/3kf1+(ak/kf)2,
(18)

where a = 0.62 is an empirical parameter. The corresponding decay rate, κtk2, is normalized by urmsk1/3 and shows a clear quadratic growth for small k and levels off near k=kf, as expected.

FIG. 7.

Dependence of the passive scalar diffusion rate κt(k)k2 on k/kf for scale separations kf/k1=3 and 8.

FIG. 7.

Dependence of the passive scalar diffusion rate κt(k)k2 on k/kf for scale separations kf/k1=3 and 8.

Close modal

Similar Lorentzian fits have been found over a broad range of different applications to turbulent magnetic diffusion: values of a0.5, a0.7, and a0.2 were found for isotropic turbulence,7 anisotropic turbulence with shear,43 and passive scalar diffusivity with shear,44 respectively.

Our work has demonstrated that the concept of turbulent diffusion carries over to radiative turbulent diffusion as well, in both optically thick and thin limits. While this was expected for the optical thick limit, it was not obvious how this would be modified in the optically thin limit, which is not a diffusion process. Instead, the optically thin case is characterized by Newtonian cooling, which then turns into turbulent Newtonian cooling. Both processes are shown to be scale-dependent, i.e., they are really described by integral kernels.

We can now also answer the question regarding the combined effect of decreased turbulence and small optical depth on the cooling at small length scales. As we have seen, turbulence always enhances the microphysical cooling rates. Thus, at small length scales where radiative diffusion is replaced by the much less efficient Newtonian cooling, turbulence speeds up this effect again. Mathematically, this process is still treated like Newtonian cooling, but now with a cooling time that is no longer given by /cγ, but by the turbulent turnover time (urmskf)1. Radiation no longer enters explicitly, except through the condition k>1 for turbulent Newtonian cooling, as opposed to k<1 for turbulent radiative diffusion.

As for the scope of future work, independent verifications of our results would certainly be desirable. In particular, it is conceivable that one can develop a test-field method similar to that employed for passive scalars.8 It would also be useful to study the effects of turbulent radiative diffusion and turbulent Newtonian cooling by comparing direct numerical simulations with mean-field models. This could be particularly insightful in more realistic situations involving stratification, turbulence, and magnetic fields, which could give rise to interesting phenomena such as magnetic spot formation.16,17

We thank Matthias Rheinhardt for useful comments on the manuscript. This work was supported in part by the Swedish Research Council, Grant No. 2019-04234. We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm.

The authors have no conflicts of interest to disclose.

The source code used for the simulations in this study, the Pencil Code,45 is freely available on https://github.com/pencil-code/. The DOI of the code is https://doi.org/10.5281/zenodo.2315093. The simulation setup and the corresponding data46 are freely available on https://doi.org/10.5281/zenodo.4085411.

1.
G.
Rüdiger
, “
Reynolds stresses and differential rotation I. On recent calculations of zonal fluxes in slowly rotating stars
,”
Geophys. Astrophys. Fluid Dyn.
16
,
239
261
(
1980
).
2.
H. K.
Moffatt
,
Magnetic Field Generation in Electrically Conducting Fluids
(
Cambridge University Press
,
Cambridge
,
1978
).
3.
F.
Krause
and
K.-H.
Rädler
,
Mean-Field Magnetohydrodynamics and Dynamo Theory
(
Pergamon Press
,
Oxford
,
1980
).
4.
K.-H.
Rädler
,
A.
Brandenburg
,
F. D.
Sordo
, and
M.
Rheinhardt
, “
Mean-field diffusivities in passive scalar and magnetic transport in irrotational flows
,”
Phys. Rev. E
84
(
4
),
046321
(
2011
).
5.
E.
Schatzman
,
Cosmic Gas Dynamics: Part 1. Cosmic Gas Dynamics
(
Wiley-Interscience
,
New York
,
1974
), p.
50
.
6.
A.
Brandenburg
and
U.
Das
, “
The time step constraint in radiation hydrodynamics
,”
Geophys. Astrophys. Fluid Dyn.
114
,
162
195
(
2020
).
7.
A.
Brandenburg
,
K.-H.
Rädler
, and
M.
Schrinner
, “
Scale dependence of alpha effect and turbulent diffusivity
,”
Astron. Astrophys.
482
,
739
746
(
2008
).
8.
A.
Brandenburg
,
A.
Svedin
, and
G. M.
Vasil
, “
Turbulent diffusion with rotation or magnetic fields
,”
Mon. Not. R. Astron. Soc.
395
,
1599
1606
(
2009
).
9.
K.-H.
Rädler
, “
Mean-Field Magnetohydrodynamics as a Basis of Solar Dynamo Theory
,” in
Basic Mechanisms of Solar Activity, Proceedings from IAU Symposium No. 71 Held in Prague, Czechoslovakia
, edited by
V.
Bumba
and
J.
Kleczek
(
D. Reidel Publishing Company Dordrecht
,
1976
), pp.
323
344
.
10.
W.
Unno
and
E. A.
Spiegel
, “
The Eddington approximation in the radiative heat equation
,”
Publ. Astron. Soc. Jpn.
18
,
85
95
(
1966
).
11.
Å.
Nordlund
, “
Numerical simulations of the solar granulation I. Basic equations and methods
,”
Astron. Astrophys.
107
,
1
10
(
1982
).
12.
E.-J.
Rijkhorst
,
T.
Plewa
,
A.
Dubey
, and
G.
Mellema
, “
Hybrid characteristics: 3D radiative transfer for parallel adaptive mesh refinement hydrodynamics
,”
Astron. Astrophys.
452
,
907
920
(
2006
).
13.
E. A.
Spiegel
, “
The smoothing of temperature fluctuations by radiative transfer
,”
Astrophys. J.
126
,
202
207
(
1957
).
14.
A.
Barekat
and
A.
Brandenburg
, “
Near-polytropic stellar simulations with a radiative surface
,”
Astron. Astrophys.
571
,
A68
(
2014
).
15.
J. M.
Edwards
, “
Two-dimensional radiative convection in the Eddington approximation
,”
Mon. Not. R. Astron. Soc.
242
,
224
234
(
1990
).
16.
A.
Brandenburg
,
I.
Rogachevskii
, and
N.
Kleeorin
, “
Magnetic concentrations in stratified turbulence: The negative effective magnetic pressure instability
,”
New J. Phys.
18
,
125011
(
2016
).
17.
B.
Perri
and
A.
Brandenburg
, “
Spontaneous flux concentrations from the negative effective magnetic pressure instability beneath a radiative stellar surface
,”
Astron. Astrophys.
609
,
A99
(
2018
).
18.
R. F.
Stein
and
Å.
Nordlund
, “
Simulations of solar granulation: I. General properties
,”
Astrophys. J.
499
,
914
933
(
1998
).
19.
Z.
Zhu
,
Y.-F.
Jiang
,
H.
Baehr
,
A. N.
Youdin
,
P. J.
Armitage
, and
R. G.
Martin
, “
Global 3D radiation hydrodynamic simulations of proto-Jupiter's convective envelope
,”
Mon. Not. R. Astron. Soc.
(published online,
2021
).
20.
L.
Soucasse
,
B.
Podvin
,
P.
Rivière
, and
A.
Soufiani
, “
Low-order models for predicting radiative transfer effects on Rayleigh-Bénard convection in a cubic cell at different Rayleigh numbers
,”
J. Fluid Mech.
917
,
A5
(
2021
).
21.
J.
Leenaarts
, “
Radiation hydrodynamics in simulations of the solar atmosphere
,”
Living Rev. Sol. Phys.
17
,
3
(
2020
).
22.
P.
Bhat
and
A.
Brandenburg
, “
Hydraulic effects in a radiative atmosphere with ionization
,”
Astron. Astrophys.
587
,
A90
(
2016
).
23.
A. J.
Mee
and
A.
Brandenburg
, “
Turbulence from localized random expansion waves
,”
Mon. Not. R. Astron. Soc.
370
,
415
419
(
2006
).
24.
C.
Federrath
,
J.
Roman-Duval
,
R. S.
Klessen
,
W.
Schmidt
, and
M.-M.
Mac Low
, “
Comparing the statistics of interstellar turbulence in simulations and observations. Solenoidal versus compressive turbulence forcing
,”
Astron. Astrophys.
512
,
A81
(
2010
).
25.
R.
Kippenhahn
and
A.
Weigert
,
Stellar Structure and Evolution
(
Springer
,
Berlin
,
1990
).
26.
Pencil Code Collaboration
,
A.
Brandenburg
,
A.
Johansen
,
P. A.
Bourdin
,
W.
Dobler
,
W.
Lyra
,
M.
Rheinhardt
,
S.
Bingert
,
N. E. L.
Haugen
,
A.
Mee
,
F.
Gent
,
N.
Babkovskaia
,
C.-C.
Yang
,
T.
Heinemann
,
B.
Dintrans
,
D.
Mitra
,
S.
Candelaresi
,
J.
Warnecke
,
P. J.
Käpylä
,
A.
Schreiber
,
P.
Chatterjee
,
M. J.
Käpylä
,
X.-Y.
Li
,
J.
Krüger
,
J. R.
Aarnes
,
G. R.
Sarson
,
J. S.
Oishi
,
J.
Schober
,
R.
Plasson
,
C.
Sandin
,
E.
Karchniwy
,
L. F. S.
Rodrigues
,
A.
Hubbard
,
G.
Guerrero
,
A.
Snodin
,
I. R.
Losada
,
J.
Pekkilä
, and
C.
Qian
, “
The Pencil Code, a modular MPI code for partial differential equations and particles: Multipurpose and multiuser-maintained
,”
J. Open Source Software
6
,
2807
(
2021
).
27.
A.
Brandenburg
, “
Computational aspects of astrophysical MHD and turbulence
,” in
Advances in Nonlinear Dynamos (the Fluid Mechanics of Astrophysics and Geophysics, Vol. 9)
, edited by
A.
Ferriz-Mas
and
M.
Núñez
(
Taylor & Francis
,
London/New York
,
2003
), pp.
269
344
.
28.
T.
Heinemann
,
W.
Dobler
,
Å.
Nordlund
, and
A.
Brandenburg
, “
Radiative transfer in decomposed domains
,”
Astron. Astrophys.
448
,
731
737
(
2006
).
29.
J. H.
Williamson
, “
Low-storage Runge-Kutta schemes
,”
J. Comput. Phys.
35
,
48
56
(
1980
).
30.
B.
Freytag
,
M.
Steffen
,
H.-G.
Ludwig
,
S.
Wedemeyer-Böhm
,
W.
Schaffenberger
, and
O.
Steiner
, “
Simulations of stellar convection with CO5BOLD
,”
J. Comput. Phys.
231
,
919
959
(
2012
).
31.
S. W.
Davis
,
J. M.
Stone
, and
Y. F.
Jiang
, “
A radiation transfer solver for Athena using short characteristics
,”
Astrophys. J.
199
,
9
(
2012
).
32.
T.
Heinemann
,
Å.
Nordlund
,
G. B.
Scharmer
, and
H. C.
Spruit
, “
MHD simulations of penumbra fine structure
,”
Astrophys. J.
669
,
1390
1394
(
2007
).
33.
M.
Rempel
,
M.
Schüssler
, and
M.
Knölker
, “
Radiative magnetohydrodynamic simulation of sunspot structure
,”
Astrophys. J.
691
,
640
649
(
2009
).
34.
A.
Roper Pol
,
A.
Brandenburg
,
T.
Kahniashvili
,
A.
Kosowsky
, and
S.
Mandal
, “
The timestep constraint in solving the gravitational wave equations sourced by hydromagnetic turbulence
,”
Geophys. Astrophys. Fluid Dyn.
114
,
130
161
(
2020
).
35.
A.
Brandenburg
and
R.
Sharma
, “
Simulating relic gravitational waves from inflationary magnetogenesis
,”
Astrophys. J.
(in press); arXiv:2106.03857 (
2021
).
36.
G.
Falkovich
, “
Bottleneck phenomenon in developed turbulence
,”
Phys. Fluids
6
,
1411
1414
(
1994
).
37.
Q.
Zheng
,
J.
Wang
,
M. M.
Alam
,
B. R.
Noack
,
H.
Li
, and
S.
Chen
, “
Transfer of internal energy fluctuation in compressible isotropic turbulence with vibrational non-equilibrium
,”
J. Fluid Mech.
919
,
A26
(
2021
).
38.
W.
Dobler
,
N. E. L.
Haugen
,
T. A.
Yousef
, and
A.
Brandenburg
, “
Bottleneck effect in three-dimensional turbulence simulations
,”
Phys. Rev. E
68
,
026304
(
2003
).
39.
C.
Federrath
,
R. S.
Klessen
,
L.
Iapichino
, and
J. R.
Beattie
, “
The sonic scale of interstellar turbulence
,”
Nat. Astron.
5
,
365
371
(
2021
).
40.
A.
Brandenburg
and
L.
Chen
, “
The nature of mean-field generation in three classes of optimal dynamos
,”
J. Plasma Phys.
86
,
905860110
(
2020
).
41.
O.
Gressel
and
D.
Elstner
, “
On the spatial and temporal non-locality of dynamo mean-field effects in supersonic interstellar turbulence
,”
Mon. Not. R. Astron. Soc.
494
,
1180
1188
(
2020
).
42.
A. B.
Bendre
and
K.
Subramanian
, “
Non-locality of the turbulent electromotive force
,”
Phys. Rev. Lett.
(submitted); arXiv:2107.10625 (
2021
).
43.
D.
Mitra
,
P. J.
Käpylä
,
R.
Tavakol
, and
A.
Brandenburg
, “
Alpha effect and diffusivity in helical turbulence with shear
,”
Astron. Astrophys.
495
,
1
8
(
2009
).
44.
E. J. M.
Madarassy
and
A.
Brandenburg
, “
Calibrating passive scalar transport in shear-flow turbulence
,”
Phys. Rev. E
82
,
016304
(
2010
).
45.
A.
Brandenburg
and
W.
Dobler
, “
Pencil Code
,” Astrophysics Source Code Library, ascl:1010.060, Dataset http://ui.adsabs.harvard.edu/abs/2010ascl.soft10060B; .
46.
A.
Brandenburg
and
U.
Das
(
2020
), “Turbulent radiative diffusion and turbulent newtonian cooling,”
Zenodo
, V.2020.10.13, Dataset