DISSIPATION EFFICIENCY IN TURBULENT CONVECTIVE ZONES IN LOW-MASS STARS

We follow the procedure outlined in GO and assume an external perturbing velocity field given by GO's Equation (8):

$$\begin{equation} \mathbf {V}_t=\mathbf {A}(t)\cdot \mathbf {x}, \end{equation} \tag{ 3 }$$

where A(t) is some matrix that depends only on time, and not space. This is appropriate when the spatial dependence of the perturbation is on scales much larger than our simulation box (e.g., tides). The matrix A(t) is assumed symmetric, since the antisymmetric part corresponds to rotation, and is not expected to contribute to the energy dissipation.

Introducing this velocity field will modify the convective flow (v₀). The time evolution of the change in the turbulent velocity (δv) due to the presence of the above external field can be written as in GO, Equation (19):

$$\begin{eqnarray} \partial _t \delta \mathbf {v}(\mathbf {x},t) &=& {-}2\mathbf {A}(t) \cdot \mathbf {v_0}(\mathbf {x}, t) - \mathbf {v}_0\cdot \mathbf {\nabla }\delta \mathbf {v} - \delta \mathbf {v}\cdot \mathbf {\nabla }\mathbf {v_0}\nonumber\\ && - (\hbox{pressure term}). \end{eqnarray} \tag{ 4 }$$

GO assumed incompressibility, and hence the pressure term simply maintains that ∇ · δv = 0. We use the output of fully compressible simulations, so for us the pressure term should be much more complex. However, the only type of compressibility that we can reasonably incorporate in the analysis is that due to the stratification of the convective layer, and even that we approximate by assuming a constant density scale height. This is discussed in more detail in Section 2.3.

We then follow GO in writing Equation (4) in Fourier space. However, since we are dealing with discretely sampled data, we use discrete Fourier transforms:

$$\begin{eqnarray} \delta \mathbf {v}_{l,m,n,p}&=&\frac{1}{N_x N_y N_z N_t} \sum _{\lambda, \mu, \nu, \phi } \widehat{\delta \mathbf {v}}_{\lambda, \mu, \nu, \phi } e^{2\pi i (\frac{\lambda l}{N_x} + \frac{\mu m}{N_y} + \frac{\nu n}{N_z} + \frac{\phi p}{N_t})}\nonumber,\\ {\mathbf {v}_0}_{l,m,n,p}&=&\frac{1}{N_x N_y N_z N_t} \sum _{\lambda, \mu, \nu, \phi } \widehat{\mathbf {v}}_{0\lambda, \mu, \nu, \phi } e^{2\pi i (\frac{\lambda l}{N_x} + \frac{\mu m}{N_y} + \frac{\nu n}{N_z} + \frac{\phi p}{N_t})},\\ \mathbf {A}(t)&=&\frac{1}{2}[\widehat{\mathbf {A}}(\Omega) e^{-i\Omega t}+\widehat{\mathbf {A}}(-\Omega)e^{i\Omega t} ]\nonumber, \end{eqnarray} \tag{ 5 }$$

where 2π/Ω is the period of the external forcing.

For more details on how exactly the Fourier transform is applied in the radial and time directions see Section 2.4.

In Fourier space to first order in A (the strength of the perturbation) and Ωτ_c (the ratio of perturbation timescale to convective turnover time), Equation (4) is written as

$$\begin{eqnarray} \delta \hat{v}_{\lambda,\mu,\nu,\phi } &=& {-}\frac{i}{\omega _\phi } \mathbf {P}_{\lambda,\mu,\nu } [ \mathbf {\widehat{A}}(\Omega)\cdot \mathbf {\hat{v}}(\omega _\phi -\Omega, \mathbf {k}_{\lambda,\mu,\nu }) \nonumber\\ &&+\mathbf {\widehat{A}}(-\Omega)\cdot \mathbf {\hat{v}}(\omega _\phi +\Omega, \mathbf {k}_{\lambda,\mu,\nu }) ], \end{eqnarray} \tag{ 6 }$$

where P_λ,μ,ν ≡ I − k'_λ,μ,νk'_λ,μ,ν/k'²_λ,μ,ν, with $\mathbf {k^{\prime }}_{\lambda,\mu,\nu }\equiv k_{\lambda,\mu,\nu }+i\hat{z}/H_\rho$, is the discrete version of the projection operator GO defined that imposes compressibility only due to a constant density scale height.

We can then express the average rate of work done (per unit volume) on the turbulent velocities by the tide to lowest non-zero order as

$$\begin{eqnarray} \displaystyle S_{\rho,\rho ^{\prime }}&\equiv& \frac{T}{\mathcal {N}^2 N_z} \sum _{\lambda,\mu,\nu,\nu ^{\prime }} \rho _{\nu -\nu ^{\prime }}^* v_{\lambda,\mu,\nu,\rho } P_{\lambda,\mu,\nu ^{\prime }} v^{*}_{\lambda,\mu,\nu ^{\prime },\rho ^{\prime }},\nonumber\\ \displaystyle \dot{\mathcal {E}}(\Omega &=&2\pi R/T) = \hbox{Re}\lbrace S_{R,-R}+S_{R,R}\rbrace\nonumber\\ && + \sum _{r\ne 0} {\displaystyle {1 \over \pi r}} \hbox{Im}\lbrace S_{r+R, r-R} + S_{r+R, r+R}\rbrace, \end{eqnarray} \tag{ 7 }$$

where $\mathcal {N}\equiv N_x N_y N_z N_t$. In keeping with GO, we have assumed that all frequencies have an infinitesimal imaginary part, which gives rise to the first term above. The second term is entirely due to the density stratification in the box: in the case of $\rho (z)=\rm const$, it is zero. This term is the most important difference between this calculation and GO.

In deriving Equation (7), we have assumed that the density is only a function of depth. Keeping the radial dependence is necessary because the simulation box encompasses several density scale heights, while the horizontal and temporal dependence of the density is a much smaller effect, entirely due to the turbulent fluctuations in the box; and if we could average over different realizations of the turbulence, they would not be present.

In order to extract an effective viscosity, we need to express the energy dissipation rate that would occur in the presence of actual anisotropic viscosity.

Most generally, viscosity is a fourth-order tensor relating the strain, given by A(t) in this case, and the viscous stress:

$$\begin{equation} \sigma _{ij}^{\rm visc} = K_{ijmn} A_{mn}(t). \end{equation} \tag{ 8 }$$

With this definition the time-averaged dissipated power is given by:

$$\begin{equation} \dot{\mathcal {E}}_{\rm visc}(\Omega)=\frac{1}{2}\int _0^{L_z} dz K_{ijmn}(z)\, \mbox{Re}\lbrace A_{ij}(\Omega)A^*_{mn}(\Omega)\rbrace, \end{equation} \tag{ 9 }$$

where, in order to remain consistent with the Fourier transform conventions, we simply replace the integral with a sum. To get the different components of K_ijmn, we evaluate Equation (7) with A having non-zero elements at different locations, and use the above equation to find the respective viscosity coefficients.

The viscosity tensor (K_ijmn) obeys a set of symmetries that dramatically reduce the number of independent components. Since the strain rate is symmetric by definition, and the stress must be symmetric in order to keep the viscous torque on infinitesimal fluid elements finite, we must have

$$\begin{equation} K_{ijmn}=K_{jimn}=K_{ijnm}. \end{equation} \tag{ 10 }$$

In addition, the only distinct direction in the problem is that of gravity ($\hat{z}$), so we expect the viscosity tensor to be symmetric with respect to rotation around the vertical axis.

With all these symmetries we are left with only six independent components of K_ijmn: K₁₁₁₁, K₃₃₃₃, K₁₂₁₂, K₁₃₁₃, K₁₁₃₃, and K₃₃₁₁. Since the last two of these always appear together in the expression, for the energy dissipation we will assume them to be equal. The remaining non-zero components can be found from those as follows:

$$\begin{equation} \begin{array}{l} K_{2222}=K_{1111}, \\ K_{1122}=K_{2211}=K_{1111}-2K_{1212}, \\ K_{1221}=K_{2121}=K_{2112}=K_{1212}, \\ \left.\begin{array}{l@{}} \hspace{-7pt}K_{3131}=K_{3113}=K_{1331}=K_{2323}\\ \hspace{-7pt}K_{3232}=K_{3223}=K_{2332} \end{array}\right\rbrace =K_{1313},\\ K_{2233}=K_{1133},\\ K_{3322}=K_{3311}. \end{array} \end{equation} \tag{ 11 }$$

A more physically meaningful set of five viscosity components can be found by noting that under these symmetries the strain rate has only four distinct components:

$$\begin{equation} \eqalign{\mathcal {A}_0\equiv A_{11}+A_{22};&\qquad \mathcal {A}_{0^{\prime }}\equiv A_{33};\\ \ms \mathcal {A}_{1}\equiv A_{13}+iA_{23};&\qquad \mathcal {A}_{2}\equiv A_{11}-A_{22}+2iA_{12},} \end{equation} \tag{ 12 }$$

along with their complex conjugates $\mathcal {A}_{-m}=\mathcal {A}_{m}^*$, which transform under rotation by angle θ around the $\hat{z}$-axis as $\mathcal {A}_m\rightarrow e^{i\theta m}\mathcal {A}_{m}$.

Clearly then, if $\dot{\mathcal {E}}_{\rm visc}(\Omega)$ is to be invariant under such rotations, it must be of the form:

$$\begin{eqnarray} \dot{\mathcal {E}}_{\rm visc}(\Omega)&=&\frac{1}{2}\int _0^{L_z} dz \big[ 4K_1|\mathcal {A}_1|^2 + K_2|\mathcal {A}_2|^2 + K_0 \mathcal {A}_0^2 + K_{0^{\prime }} \mathcal {A}_{0^{\prime }}^2\nonumber\\ && + 2K_{00^{\prime }} \mathcal {A}_0\mathcal {A}_{0^{\prime }} \big], \end{eqnarray} \tag{ 13 }$$

where the five new viscosity coefficients can be expressed in terms of K_ijmn as follows:

$$\begin{equation} \begin{array}{c} K_0={ {1 \over 2}}(K_{1111}+K_{1122}),\\ \ms K_{0^{\prime }}=K_{3333},\\ \ms K_{00^{\prime }}=K_{1133},\\ \ms K_1=K_{1313},\\ \ms K_2={ {1 \over 2}}\left(K_{1111}-K_{1122}\right). \end{array} \end{equation} \tag{ 14 }$$

Since in Equation (9), we allow the viscosity to depend on depth, and there is no way to constrain this dependence, we have to choose some radial profile a priori. Our choice is motivated by mixing length theory:

$$\begin{equation} K_m(z)=K_m^0(\Omega) \rho \langle v^2\rangle^{1/2} H_p, \end{equation} \tag{ 15 }$$

where K⁰_m are dimensionless constants that depend on the frequency of the external shear (Ω), and H_p is the local pressure scale height. This is reasonable, since the turbulent viscosity should scale as some length scale times some velocity scale. Clearly the relevant velocity scale is that of convection, and in accordance with mixing length theory, we use the mixing length as the length scale, which is assumed proportional to the pressure scale height. If the mixing length is really the relevant quantity, we expect that the value of K⁰_m will be proportional to the mixing length parameter for the particular simulation. This same scaling has been assumed for all previous effective viscosity prescriptions (Zahn 1966, 1989; Goldreich & Keeley 1977; Goldreich & Kumar 1988; Goldreich et al. 1994).

In deriving the expression for $\dot{\mathcal {E}}$ (Equation (7)), we assumed that the perturbation to the convective velocity field due to the tide will be anelastic: ∇ · δv + v_z/H_ρ = 0, with $H_\rho =\rm const$. In this section, we define two diagnostics which measure how important the ignored compressibility is.

There are two sources of compressibility in the convective flow:

• 1.  
  The convective flow carrying parcels of matter through layers of different hydrostatic pressure, or in other words due to the stratification.
• 2.  
  Localized compression due to a possibly supersonic flow, e.g., shocks.

Ignoring the second one is justified, as long as the flow velocity is much less than the local speed of sound. In the simulations we use, that condition is met by the unperturbed flow for most of the box, with the exception of the supersonic driving region near the top. If the unperturbed flow is subsonic and hence incompressible, the perturbations due to a "small" external field can safely be assumed incompressible as well. To measure the compressibility in the simulation box, we introduce the parameter:

$$\begin{equation} \xi \equiv \tau _c\left[\mathbf {\nabla }+\hat{z}\frac{d\ln \rho }{dz}\right] \cdot \mathbf {v_0}, \end{equation} \tag{ 16 }$$

where, ρ is the density averaged over horizontal slices and time.

This quantity deviates from zero due to localized, transient compressions (e.g., shocks). Since those are unlikely to live longer than a convective turnover time, this quantity is a suitable diagnostic for the importance of such effects.

Because we are measuring the mass-averaged dissipation, in order for the perturbative treatment discussed above to be valid, we can only have ξ ≳ 1 for a negligible fraction of the mass. In Figure 1, we plot the time-averaged fraction of mass that resides in regions of the convective box which have ξ greater than some value. The value of the convective turnover time, necessary to evaluate ξ, was calculated as τ_c ≡ FWHM(v_z)/max(v_z), where FWHM(v_z) is the thickness of the layer over which v_z>max(v_z)/2.

Zoom In Zoom Out Reset image size

As we can see, in all cases, ξ>1 for less than 1% of the mass. This compression is concentrated near the top of the box, where the density, and hence the dynamic viscosity, is small. So, even though compressibility and shocks are important in determining the flow that develops, no appreciable dissipation occurs in strongly compressible regions. The situation is further improved by the fact that we apply a window in the vertical direction that significantly reduces the importance, and completely ignores part of the compressible driving region near the top of the box in determining the effective viscosity (see Section 2.4).

We partially treat the first source of compressibility discussed above by, imposing $H_\rho =\rm const$ in the continuity equation. However, this is not valid for most astrophysically interesting convective zones. In fact in all cases considered here the density scale height varies by a factor of a few between the top and the bottom of the convective layer. In some sense, assuming $H_\rho =\rm const$ is not any better than assuming H_ρ = ∞. We argue that ignoring the stratification from the continuity equation is a reasonable approximation.

The simplest way to justify this is to repeat the evaluation of the viscosity with different values of H_ρ within the range encountered in the convective layer of interest.

We can also gauge the importance of the stratification by comparing dln ρ/dz to ∂ln δv_z/∂z. We evaluate dln ρ/dz directly, and estimate:

$$\begin{equation} \frac{d\ln \delta v_z}{dz}=\frac{1}{\langle \delta {v_z}^2\rangle^\frac{1}{2}} \left<\left(\frac{\partial \delta v_z}{\partial z}\right)^2 \right>^\frac{1}{2}=\frac{1}{\big\langle v_z^2\big\rangle^\frac{1}{2}} \left<\left(\frac{\partial v_z}{\partial z}\right)^2 \right>^\frac{1}{2}. \end{equation} \tag{ 17 }$$

The last expression comes from Equation (6), and is correct when the last row of A(t) contains only a single non-zero entry. Since those are the only cases we use, this expression is sufficient for us.

In Figure 2, we compare dln ρ/dz to ∂ln δv_z/∂z (estimated as in the above expressions). We see that the logarithmic gradient of the density is approximately 2 orders of magnitude smaller than the typical logarithmic velocity gradient, and hence we are justified in ignoring it.

Zoom In Zoom Out Reset image size

Using discrete Fourier transforms to represent a data set, forces the assumption that the data is periodic in all dimensions. While this holds for each horizontal slice, it is violated for the radial and time dimensions. Ignoring the problem leads to artificially introducing spectral power at the highest frequencies because of the jumps at the boundaries. To avoid this, we need a special way to deal with the non-periodic directions.

The usual solution is to window the data so that it goes smoothly to zero at the edges of the domain. This has the effect that it makes the values near the center of the domain relatively more important than those near the boundaries. Incidentally, this is exactly what we would like in the radial direction, since the flow near the top and the bottom is affected by the artificial boundary conditions and is not representative of the actual flow that would occur in a star.

Further, as discussed in Section 2.3, we expect that the compressibility of the flow that we neglect might be significant in the upper end of the box, where the density is small and the flow is supersonic. So making this region's contribution to the overall dissipation small is exactly what we would like. In fact, in the radial direction we go a step further and limit the window to completely exclude some part of the box near the top and bottom boundaries (see Equations (18) and (19)).

In the time direction, as long as the time interval we have simulated is "representative" of the actual convection that occurs in a star, weighing the center of the interval more than the edges should not be a problem.

To confirm that the chosen window function is not affecting our final result, we derive the effective viscosity coefficient using two common windows:

$$\begin{eqnarray} \hbox{Welch}: & \mathcal {N}\left[1-\left(\frac{2t-T}{T}\right)^2\right] \hbox{max}\left[0,1-\left(\frac{2z-L_z}{\alpha L_z} \right)^2\right],\nonumber\\ \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} \hbox{Bartlett}: & \mathcal {N}\left[1-\left|\frac{2t-T}{T}\right|\right] \hbox{max}\left[0,1-\left|\frac{2z-L_z}{\alpha L_z}\right|\right],\nonumber\\ \end{eqnarray} \tag{ 19 }$$

where $\mathcal {N}$ is a normalization factor numerically equal to the inverse of the average of the squares of the window function at all the grid points, and α is a parameter determining what fraction of the radial span of the box we include in the analysis; that is, we exclude (1 − α) fraction of the linear size of the box, half from the top and half from the bottom.

The 1D stellar models used to initialize each run were computed with the YREC stellar evolution code (Guenther et al. 1992). They were calibrated to the Sun and evolved from the ZAMS. Both the 1D and 3D codes use the same realistic physics as described by Guenther & Demarque (1997), most notably the Alexander & Ferguson (1994) opacities at low temperatures, the OPAL opacities and equation of state (Iglesias & Rogers 1996), hydrogen and helium ionization and helium and heavy element diffusion.

Some details of the three models are given in Table 1. The fractional radius is given as R/R_☉, where R is radius of the stellar body and R_☉ is the radius of the Sun. Both are defined at the point where T = T_eff. The surface gravity and effective temperature are in cgs units.

The horizontal dimensions of each computational box (Row 5 in Table 1) were estimated by assuming that the granule size will roughly scale inversely with g. The final row gives the number of grid points in the two horizontal and vertical directions in the square based box.

The above splitting of the viscosity in five components was done in order to allow for anisotropic dissipation. It is interesting to see how anisotropic the derived effective viscosity really is. The general isotropic case has only two viscosity components: a bulk viscosity (ζ) and a shear viscosity (η) (see Landau & Lifshitz 1987). In terms of these components, the isotropic K_ijmn tensor is

$$\begin{equation} K_{ijmn}=\eta \left(\delta _{im}\delta _{jn} + \delta _{in}\delta _{jm} - \case{2}{3}\delta _{ij}\delta _{mn}\right) + \zeta \delta _{ij}\delta _{mn}. \end{equation} \tag{ 22 }$$

From this it can be seen that the five components of the viscosity we calculate must obey the relations:

$$\begin{eqnarray} K_1&=&K_2=\eta,\nonumber\\ K_0&=&K_{00^{\prime }}+K_1,\\ K_{0^{\prime }}&=&K_0+K_1.\nonumber \end{eqnarray} \tag{ 23 }$$

We see that the first two of these are clearly satisfied by the viscosity coefficients of Equation (21) to within the quoted uncertainties. The degree to which the last equation is not satisfied is

$$\begin{equation} \frac{K_{0^{\prime }}}{K_0+K_1}-1=0.74\pm 0.25. \end{equation} \tag{ 24 }$$

Considering the fact that the flows in our simulation boxes are not exactly like those inside the stars, and the loosely estimated errors in Equation (21); we can conclude that the effective viscosity we find is only mildly anisotropic, and it is perhaps reasonable to approximate it as completely isotropic bulk and shear viscosities with the following values:

$$\begin{equation} \begin{array}{l} \eta =0.020\pm 0.003,\\ \zeta \in (0.018,0.056). \end{array} \end{equation} \tag{ 25 }$$

The reason for ζ not being well determined by the viscosity coefficients of Equation (21) is that those coefficients do not exactly correspond to an isotropic viscosity (see Equation (24)).