text:
# Dual theories for mixed symmetry fields. Spin-two case: (1,1) versus (2,1) Young symmetry type fields.

## I Introduction

In general there is a lot freedom to choose variables for the description of a physical system. For example, in some cases it is desirable to have a formulation where some symmetries are explicit, and this requires the use of a redundant set of variables to describe the system configuration. In other cases, as when we perform a canonical quantization, it might be more convenient to have a minimal, non redundant, set of variables. The actual proof of the equivalence between different descriptions is usually a non-trivial task. Equivalent descriptions of a given physical theory in terms of a different choice of fields are said to be dual to each other, and the relation among the fields provides the corresponding duality transformation .

One of the simplest examples is the scalar-tensor duality. It corresponds to the equivalence between a free massless scalar field $`\phi `$, with field strength $`f_\mu =_\mu \phi `$, and a massless antisymmetric field $`B_{\mu \nu }`$, the Kalb-Ramond field, with field strength $`H_{\mu \nu \sigma }=_\mu B_{\nu \sigma }+_\nu B_{\sigma \mu }+_\sigma B_{\mu \nu }`$ . The former description provides the equation of motion $`^\mu f_\mu =0`$, together with the Bianchi identities $`^\mu ϵ_{\mu \nu \sigma }{}_{}{}^{\tau }f_{\tau }^{}=0`$. In turn, the dual description provides the equations of motion $`^\mu H_{\mu \nu \sigma }=0`$ and the Bianchi identities $`^\mu ϵ_\mu ^{\nu \sigma \tau }H_{\nu \sigma \tau }=0`$. The duality transformation

$$f_\mu \frac{1}{3}ϵ_\mu ^{\nu \sigma \tau }H_{\nu \sigma \tau }$$
(1)
interchanges dynamical equations with Bianchi identities, giving a full correspondence between both descriptions.

An important predecessor of the modern approach to duality is the electric-magnetic symmetry $`(\stackrel{}{E}+i\stackrel{}{B})e^{i\varphi }(\stackrel{}{E}+i\stackrel{}{B})`$ of the free Maxwell equations. When there are charged sources this symmetry can be maintained by introducing magnetic monopoles . This transformation provides a connection between weak and strong couplings via the Dirac quantization condition. At the level of Yang-Mills theories with spontaneous symmetry breaking this kind of duality is expected, due to the existence of topological dyon-type solitons . The extension of electromagnetic duality to $`SL(2,Z)`$ is usually referred to as S-duality, and plays an important role in the non-perturbative study of field and string theories .

These basic ideas have been subsequently generalized to arbitrary forms in arbitrary dimensions. Well known dualities are the ones between massless $`p`$ -forms and $`(dp2)`$-forms fields and between massive $`p`$ and $`(dp1)`$-forms in $`d`$ dimensional space-time. These dualities among free fields have been proved using parent Lagrangians  as well as the canonical formalism . They can be extended to include source interactions .

The above duality among forms can be understood as a relation between fields in different representations of the Lorentz group. The origin of this equivalence can be traced using the little group technique for constructing the representations of the Poincare group in $`d`$ dimensions. Given a standard momentum for a massive or a massless particle, the actual degrees of freedom are determined by its spin components, which are given by the irreducible representations of $`SO(d1)`$ and $`SO(d2)`$ respectively. These are expressed by traceless tensors with a definite permutation symmetry characterized by the Young diagrams. For the orthogonal groups $`O(n)`$, the sum of the lengths of the first two columns of the Young diagrams is constrained to be less than or equal to $`n`$. Two Young diagrams having their first column of length $`l`$ and $`nl`$ respectively, with the rest of the diagrams being identical, are called associated. For the group $`SO(n)`$ associated Young diagrams correspond to the same representation . The existence of two different tensorial realizations for one irreducible spin representation suggests that we can construct two theories with fields of different tensorial character for the same physical system. The relation between two equivalent theories expressed in terms of fields corresponding to associated Young diagrams should be interpreted as a duality transformation.

In particular, the usual duality among $`p`$-forms can be interpreted in this way. A (traceless) $`p`$-form belongs to a Young diagram with only one column and $`p`$ sites, and a $`(dp1)`$-form belongs to a Young diagram with one column and $`(dp1)`$ sites. Thus they give the same representation of $`SO(d1)`$ and have the same spin; therefore, if the masses are equal there is a possibility of constructing alternative theories for the same physical system. The same can be said regarding the massless case if we consider $`p`$ and $`(dp2)`$ forms.

The preceding discussion suggests the possibility of generalizing the duality transformations among $`p`$-forms to tensorial fields with arbitrary Young symmetry types. The simplest generalization in four dimensions is provided by a massive second-rank symmetric tensor $`h_{\mu \nu }`$, with Young symmetry $`(1,1)`$, and a three index tensor $`T_{\mu \nu \rho }`$ with Young symmetry $`(2,1)`$. Consistent massless free , interacting , and massive  theories of mixed Young symmetry tensors have been constructed in the past, but the attempts to prove a dual relation between these descriptions did not lead to a positive answer . Additional interest in this type of theories arises from the recent formulation of $`d=11`$ dimensional supergravity as a gauge theory for the $`\mathrm{osp}(32|1)`$ superalgebra. It includes a totally antisymmetric fifth-rank Lorentz tensor one form $`b_\mu ^{abcde}`$, whose mixed symmetry piece does not have any related counterpart in the standard $`d=11`$ supergravity theory .

In this work we show that such dual descriptions can be constructed. We use a generalization of the parent Lagrangian method to construct a family of Lagrangians for the massive field $`T_{(\mu \nu )\rho }`$, which are dual to the standard Fierz-Pauli Lagrangian for a massive spin two field $`h_{\mu \nu }`$. In contrast with the $`p`$-forms case, we obtain a multiparameter family of duality transformations. Most notably, the kinetic part of the dual $`T`$-Lagrangians is unique, and the parameters appear only in the mass term. The massless limit of these descriptions is well-behaved in the sense that it has the correct two propagating degrees of freedom, in contrast with the absence of degrees of freedom found in previous attempts.

## II Massive spin-two irreducible representations

In the following, we perform a detailed construction for the massive spin-two fields in four dimensions. There are two associated Young diagrams, which correspond to a traceless symmetric rank-two tensor, and to a traceless rank-three tensor, antisymmetric in two indices, and satisfying a Jacobi identity. The most usual known description is given in terms of a symmetric tensor $`h_{\mu \nu }=h_{\nu \mu }`$, with the dynamics defined by the Fierz-Pauli Lagrangian

$$=_\mu h^{\mu \nu }_\alpha h_\nu ^\alpha +\frac{1}{2}_\alpha h^{\mu \nu }^\alpha h_{\mu \nu }+_\mu h^{\mu \nu }_\nu h_\alpha ^\alpha \frac{1}{2}_\alpha h_\mu ^\mu ^\alpha h_\nu ^\nu \frac{M^2}{2}\left(h_{\mu \nu }h^{\mu \nu }h_\mu ^\mu h_\nu ^\nu \right).$$
(2)
¿From the equation of motion it follows that

$$(^2+M^2)h_{\mu \nu }=0,_\mu h^{\mu \nu }=0,h_\mu {}_{}{}^{\mu }=0,$$
(3)
which leads to the five degrees of freedom corresponding to a massive spin-two irreducible representation.

An alternative description for the free massive spin-two field corresponding to the associated $`(2,1)`$ Young diagram has been already proposed. It is based on the tensor field $`T_{(\mu \nu )\sigma }`$ that satisfies

$$T_{(\mu \nu )\sigma }=T_{(\nu \mu )\sigma };T_{(\mu \nu )\sigma }+T_{(\nu \sigma )\mu }+T_{(\sigma \mu )\nu }=0.$$
(4)
At this stage $`T_{(\mu \nu )\sigma }`$ has $`20`$ independent components. The proposed Lagrangian is

$$=\frac{1}{36}\{\left(F_{(\mu \nu \sigma )\tau }\right)^23(F_{(\mu \nu \sigma )}{}_{}{}^{\mu })^23M^2[\left(T_{(\mu \nu )\sigma }\right)^22(T_{(\mu \nu )}{}_{}{}^{\mu })^2]\}$$
(5)
with the field strength given by

$$F_{(\mu \nu \sigma )\tau }=_\mu T_{(\nu \sigma )\tau }+_\nu T_{(\sigma \mu )\tau }+_\sigma T_{(\mu \nu )\tau }.$$
(6)
In this case the equations of motion imply:

$$(^2+M^2)T_{(\mu \nu )\sigma }=0,T_{(\mu \nu )}{}_{}{}^{\mu }=0,^\mu T_{(\mu \nu )\sigma }=0,^\sigma T_{(\mu \nu )\sigma }=0.$$
(7)
The first algebraic condition gives four identities. The second derivative condition is a consequence of the first one plus the cyclic identity in (4). The first derivative condition includes four identities, $`^\mu ^\nu T_{\mu \nu \rho }=0`$, plus one more, $`^\mu T_{\mu \rho }^\rho `$, when we consider independently the zero-trace condition. This leads to $`11`$ independent derivative conditions. The final count produces $`5`$ independent degrees of freedom, as appropriate to a massive spin-two field. It is natural to suspect that both theories could be related by a duality transformation, but the attempts to construct such a transformation have had no success . As we will show below the Lagrangian (5) is only one of the possible descriptions for an irreducible massive spin two field in terms of a $`T_{(\mu \nu )\sigma }`$ tensor, and in fact is not a suitable choice to construct a dual transformation.

## III Parent Lagrangian for dual theories

A standard procedure to construct alternative descriptions for a physical system, which turn out to be related by a duality transformation, consists in finding a quadratical parent Lagrangian that contains both types of fields, from which each theory can be obtained by eliminating either one of them through the corresponding equations of motion. In this work we show that a generalization of this method, already successfully applied to $`p`$-forms, allows us to construct dual theories for the massive spin-two representation. Contrary to the approach in , we start from a tensor $`T_{(\mu \nu )\sigma }`$ which is only antisymmetric in the $`\mu \nu `$ indices, without a priori satisfying the cyclic identity, together with the standard symmetric tensor $`h_{\mu \nu }`$. Eliminating the field $`T_{(\mu \nu )\sigma }`$ from the parent Lagrangian we impose the resulting Lagrangian for $`h_{\mu \nu }`$ to be the Fierz-Pauli one describing a massive spin-two system. This provides relations among the parameters of the model. Once these relations are implemented in the parent Lagrangian, the elimination of field $`h_{\mu \nu }`$ leads to the required dual Lagrangian in terms of field $`T_{(\mu \nu )\sigma }`$. The necessary Lagrangian constraints leading to the $`(2,1)`$ Young symmetry type of the field together with the required five independent degrees of freedom come from the corresponding Euler-Lagrangian equations.

Let us recall the general structure of the parent Lagrangian used to establish duality between a massive $`dq1`$ form $`L`$ and a massive $`q`$ form $`B`$

$$_P=\frac{1}{2}L{}_{}{}^{}L+LdB+\frac{M^2}{2}B{}_{}{}^{}B,$$
(8)
which we will generalize to our case. The most general first order bilinear Lagrangian for $`T_{(\mu \nu )\sigma }`$ and $`h_{\mu \nu }`$ has seven bilinears in $`T`$

$`T_{(\mu \nu )\sigma }T^{(\mu \nu )\sigma },T_{(\mu \nu )}^\nu T_\sigma ^{(\mu \sigma )},T_{(\mu \nu )\sigma }T^{(\mu \sigma )\nu },`$ (9)
$`ϵ^{\mu \nu \alpha \beta }T_{(\mu \nu )\sigma }T_{(\alpha \beta )}^\sigma ,`$ $`ϵ^{\mu \nu \alpha \beta }T_{(\mu \nu )\sigma }T_{\alpha )\beta }^{(\sigma },ϵ^{\mu \nu \alpha \beta }T_{(\mu \sigma )\nu }T_{\alpha )\beta }^{(\sigma },ϵ^{\mu \nu \alpha \beta }T_{(\mu \nu )\alpha }T_{(\beta \lambda )}^\lambda ,`$ (10)
two in $`h`$,

$$h_{\mu \nu }h^{\mu \nu },h_\mu ^\mu h_\nu ^\nu ,$$
(11)
and seven mixing (duality generating) terms which contain both fields,

$`T_{(\mu \nu )\sigma }^\mu h^{\nu \sigma },T_{(\mu \nu )}^\nu ^\mu h_\sigma ^\sigma ,T_{(\mu \nu )}^\nu ^\sigma h_\sigma ^\mu ,`$ (12)
$`T_{(\mu \nu )\sigma }ϵ^{\mu \nu \alpha \beta }_\alpha h_\beta ^\sigma ,`$ $`T_{(\mu \nu )\sigma }ϵ^{\mu \nu \sigma \beta }_\alpha h_\beta ^\alpha ,T_{(\mu \nu )\sigma }ϵ^{\mu \nu \sigma \beta }_\beta h_\alpha ^\alpha ,T_{(\mu \nu )\sigma }ϵ^{\mu \sigma \alpha \beta }_\beta h_\alpha ^\nu .`$ (13)
In the following, and for the sake of simplicity, we will explore in detail only the Lagrangian with one mixing term, $`T_{(\mu \nu )\sigma }ϵ^{\mu \nu \alpha \beta }_\alpha h_\beta ^\sigma `$ , because this is the most natural generalization of the term $`LdB`$ in (8). Thus, we take our parent Lagrangian to be

$`L`$ $`=`$ $`aT_{(\mu \nu )\sigma }T^{(\mu \nu )\sigma }+bT_{(\mu \nu )}{}_{}{}^{\nu }T_{}^{(\mu \sigma )}{}_{\sigma }{}^{}+cT_{(\mu \nu )\sigma }T^{(\mu \sigma )\nu }+dϵ^{\mu \nu \alpha \beta }T_{(\mu \nu )\sigma }T_{(\alpha \beta )}^\sigma `$ (17)
$`+mϵ^{\mu \nu \alpha \beta }T_{(\mu \nu )}{}_{}{}^{\sigma }T_{(\sigma \alpha )\beta }^{}+nϵ^{\mu \nu \alpha \beta }T_{(\mu \sigma )\nu }T^{(\sigma }{}_{\alpha )\beta }{}^{}+pϵ^{\mu \nu \alpha \beta }T_{(\mu \nu )\alpha }T_{(\beta \lambda )}^\lambda `$
$`+eT_{(\mu \nu )\sigma }ϵ^{\mu \nu \alpha \beta }_\alpha h_\beta ^\sigma `$
$`+fh_{\mu \nu }h^{\mu \nu }+kh_\mu {}_{}{}^{\mu }h_{\nu }^{}{}_{}{}^{\nu }.`$
The equations of motion for $`h_{\mu \nu }`$ allow us to solve algebraically this field in terms of $`T_{(\mu \nu )\sigma }`$

$$h^{\mu \nu }=\frac{e}{4f}(ϵ^{\alpha \beta \sigma \nu }_\sigma T_{(\alpha \beta )}{}_{}{}^{\mu }+ϵ^{\alpha \beta \sigma \mu }_\sigma T_{(\alpha \beta )}{}_{}{}^{\nu })g^{\mu \nu }\frac{ke}{2\left(f+4k\right)f}ϵ^{\alpha \beta \rho \kappa }_\rho T_{(\alpha \beta )\kappa }.$$
(18)
In a similar way, from the Euler-Lagrange equation for $`T^{(\mu \nu )\sigma }`$ we can algebraically solve this field in terms of $`h^{\mu \nu }`$

$$CT^{(\kappa \lambda )\sigma }=e\left(2m+n4d\right)D^{\lambda \sigma \kappa }+\frac{e}{2}\left(2a+c\right)ϵ_{\mu \nu }{}_{}{}^{\kappa \lambda }D_{}^{\nu \sigma \mu }\frac{e}{2}Eϵ_\mu {}_{}{}^{\sigma \kappa \lambda }D_{}^{\mu }\frac{e}{2}FG^{\sigma \kappa \lambda }.$$
(19)
where

$`D^{\lambda \sigma \kappa }`$ $``$ $`^\lambda h^{\sigma \kappa }^\kappa h^{\sigma \lambda },D^\mu _\beta h^{\beta \mu }^\mu h_\beta ^\beta =D_\sigma ^{\sigma \mu },`$ (20)
$`G^{\sigma \kappa \lambda }`$ $``$ $`g^{\sigma \kappa }\left(_\alpha h^{\alpha \lambda }^\lambda h_\alpha ^\alpha \right)g^{\sigma \lambda }\left(_\alpha h^{\alpha \kappa }^\kappa h_\alpha ^\alpha \right)=g^{\sigma \kappa }D^\lambda g^{\sigma \lambda }D^\kappa `$ (21)
and

$`A`$ $`=`$ $`\left(4d+m3p+2n\right)\left[\left(4d+m3p+2n\right)^2+2(2a+c+3b)\left(ac\right)\right]^1,`$ (22)
$`B`$ $`=`$ $`2(2a+c+3b)\left[\left(4d+m3p+2n\right)^2+2(2a+c+3b)\left(ac\right)\right]^1,`$ (23)
$`C`$ $`=`$ $`\left(2m+n4d\right)^2+\left(2a+c\right)^2,`$ (24)
$`E`$ $`=`$ $`\left(2m+n4d\right)\left(2bA+(m+np)B\right)+\left(2a+c\right)\left(2(m+np)A+cB\right),`$ (25)
$`F`$ $`=`$ $`\left(2m+n4d\right)\left(2(m+np)A+cB\right)+\left(2a+c\right)\left(2bA+(m+np)B\right).`$ (26)
The parent Lagrangian we are considering contains seven bilinears in $`T`$, two bilinears in $`h`$, and one mixing term, which means we have started with ten parameters. Nevertheless, only eight combinations appear in the dual transformations (18) and (19): $`a`$, $`b`$, $`c`$, $`(m+np)`$, $`(2m+n4d)`$, $`e`$, $`f`$, and $`k`$. Note that $`(4d+m3p+2n)=3(m+np)(2m+n4d)`$. We further fix appropriate combinations by using Eq. (19) to rewrite the parent Lagrangian (17) in terms of $`h^{\mu \nu }`$ only, and subsequently demanding it to be the Fierz-Pauli Lagrangian for a massive spin two field. The result is

$$C=2e^2\left(2a+c\right),E=2\left(2a+c\right),f=k=\frac{M^2}{2},$$
(27)
The parameters $`d`$, $`m`$, $`n`$ and $`p`$ only appear through the combinations $`\left(2m+n4d\right)`$ and $`(m+np)`$, each of which can be written in terms of $`a`$, $`b`$, $`c`$ and $`e`$ using Eqs. (27). Considering that the same equations of motion are obtained up to an overall factor in the Lagrangian, we have constructed a three-parameter family of parent Lagrangians for a spin-two field with mass $`M`$ leading to the Fierz-Pauli Lagrangian for $`h_{\mu \nu }`$.

## IV Spin-two mixed symmetry description

Now we consider the dual description of the massive spin two system in terms of the field $`T^{(\mu \nu )\rho }`$. To this end we use Eq. (18) and rewrite the parent Lagrangian as

$`L`$ $`=`$ $`{\displaystyle \frac{e^2}{3M^2}}(\left(F_{(\mu \nu \sigma )\tau }\right)^2{\displaystyle \frac{3}{2}}(F_{(\mu \nu \sigma )}{}_{}{}^{\mu })^2+{\displaystyle \frac{1}{2}}\left(ϵ^{\mu \nu \alpha \beta }_\beta T_{(\mu \nu )\alpha }\right)^2)`$ (30)
$`+aT_{(\mu \nu )\sigma }T^{(\mu \nu )\sigma }+bT_{(\mu \nu )}{}_{}{}^{\nu }T_{}^{(\mu \sigma )}{}_{\sigma }{}^{}+cT_{(\mu \nu )\sigma }T^{(\mu \sigma )\nu }+dϵ^{\mu \nu \alpha \beta }T_{(\mu \nu )\sigma }T_{(\alpha \beta )}^\sigma `$
$`+mϵ^{\mu \nu \alpha \beta }T_{(\mu \nu )}{}_{}{}^{\sigma }T_{(\sigma \alpha )\beta }^{}+nϵ^{\mu \nu \alpha \beta }T_{(\mu \sigma )\nu }T^{(\sigma }{}_{\alpha )\beta }{}^{}+pϵ^{\mu \nu \alpha \beta }T_{(\mu \nu )\alpha }T_{(\beta \lambda )}{}_{}{}^{\lambda }.`$
The second term of the kinetic part of this Lagrangian has a coefficient $`\left(\frac{3}{2}\right)`$ instead of the corresponding $`\left(3\right)`$ in the Lagrangian (5). This fact was recognized in  as a potential problem for dualization. However, here we do not impose the cyclic property (4) for $`T_{(\mu \nu )\sigma }`$ from the beginning, and this is manifested in the existence of the third term in the kinetic Lagrangian, which is crucial for dualization.

It turns out that the trace $`T_\mu T_{(\mu \beta )}^\beta `$ is an auxiliary field in the above Lagrangian. This is made explicit by defining

$$T_{(\mu \nu )\sigma }=\widehat{T}_{(\mu \nu )\sigma }\frac{1}{3}\left(g_{\mu \sigma }T_\nu g_{\nu \sigma }T_\mu \right),$$
(31)
where $`\widehat{T}_{(\mu \nu )\sigma }`$ is a traceless field, $`\widehat{T}_{(\mu \nu )}{}_{}{}^{\nu }=0`$. Rewriting the Lagrangian (30) in terms of $`\widehat{T}^{(\chi \psi )\sigma }`$ and $`T_\mu `$ and using the Euler-Lagrangian equation for $`T_\mu `$ we obtain

$$T^\beta =\frac{4d+m+2n3p}{2\left(2a+3b+c\right)}ϵ^{\mu \nu \alpha \beta }\widehat{T}_{(\mu \nu )\alpha },$$
(32)
which allows us to algebraically eliminate the trace in the Lagrangian (30). We further introduce the field strength $`\widehat{F}_{(\alpha \beta \gamma )\nu }`$ appropriate to $`\widehat{T}^{(\mu \nu )\sigma }`$, defined in (6). In this way, up to a global factor, (30) can be written as

$`L`$ $`=`$ $`{\displaystyle \frac{4}{9}}\widehat{F}_{(\alpha \beta \gamma )\nu }\widehat{F}^{(\alpha \beta \gamma )\nu }+{\displaystyle \frac{2}{3}}\widehat{F}_{(\alpha \beta \gamma )\nu }\widehat{F}^{(\alpha \beta \nu )\gamma }\widehat{F}_{(\alpha \beta \mu )}{}_{}{}^{\mu }\widehat{F}_{}^{(\alpha \beta \nu )}_\nu `$ (34)
$`+P\widehat{T}_{(\nu \rho )\sigma }\widehat{T}^{(\nu \rho )\sigma }+Q\widehat{T}_{(\nu \rho )\sigma }\widehat{T}^{(\nu \sigma )\rho }+Rϵ^{\mu \nu \alpha \beta }\widehat{T}_{(\mu \nu )\sigma }\widehat{T}_{(\alpha \beta )}{}_{}{}^{\sigma }+\lambda ^\mu \widehat{T}_{(\mu \beta )}{}_{}{}^{\beta },`$
where we have introduced the Lagrange multiplier $`\lambda ^\mu `$ to enforce the traceless condition upon $`\widehat{T}_{(\mu \sigma )\nu }`$. Here it is

$$P=\frac{2M^2}{3e^2}\left(\frac{1}{B}(2a+c)\right),Q=\frac{2M^2}{3e^2}\left(\frac{2}{B}+(2a+c)\right),R=\frac{M^2}{2e^2}\left(2m+n4d\right).$$
(35)
The equations (27) translate in the following relations between these parameters

$`12P^23Q^2+16R^2`$ $`=`$ $`0,`$ (36)
$`{\displaystyle \frac{1}{4}}(2P+Q)+4{\displaystyle \frac{R^2}{(2P+Q)}}`$ $`=`$ $`M^2\text{.}`$ (37)
This family of Lagrangians, dual to the Fierz-Pauli Lagrangian (2) by construction, is the main result of this work.

Let us recall that in (34) the field $`\widehat{T}_{(\mu \nu )\sigma }`$ is only antisymmetric and traceless, therefore having 20 components, but the Euler-Lagrange equations give rise to the necessary Lagrangian constraints leading to the final required five independent degrees of freedom. The final equations of motion together with the independent Lagrangian constraints which completely determine the dynamics are

$`_\alpha (^\alpha \widehat{T}^{(\nu \rho )\sigma }+^\nu \widehat{T}^{(\rho \alpha )\sigma }+^\rho \widehat{T}^{(\alpha \nu )\sigma }){\displaystyle \frac{1}{4}}(Q+2P)\widehat{T}^{(\nu \rho )\sigma }{\displaystyle \frac{1}{2}}Rϵ^{\nu \rho \alpha \beta }\widehat{T}_{(\alpha \beta )}{}_{}{}^{\sigma }=0,`$ (38)
$`\widehat{T}_{(\mu \beta )}{}_{}{}^{\beta }=0,ϵ_{\alpha \beta \gamma \lambda }\widehat{T}^{(\beta \gamma )\lambda }=0,_\sigma \widehat{T}_{(\alpha \beta )}{}_{}{}^{\sigma }=0,`$ (39)
$`_\nu \widehat{T}^{(\nu \rho )\sigma }+{\displaystyle \frac{R}{2P+Q}}(ϵ^{\nu \rho \alpha \beta }_\nu \widehat{T}_{(\alpha \beta )}{}_{}{}^{\sigma }+ϵ^{\nu \sigma \alpha \beta }_\nu \widehat{T}_{(\alpha \beta )}{}_{}{}^{\rho })=0.`$ (40)
To be sure we have the correct number of degrees of freedom and to identify the mass parameter it is convenient to use plane wave solutions in a rest frame, where $`k^\mu =(\mu ,\mathrm{𝟎})`$. In this frame the equations of motion read

$$\left(\frac{1}{4}(Q+2P)+\mu ^2\right)\widehat{T}^{(\nu \rho )\sigma }+\mu ^2\left(g^{0\nu }\widehat{T}^{(\rho 0)\sigma }+g^{0\rho }\widehat{T}^{(0\nu )\sigma }\right)+\frac{1}{2}Rϵ^{\nu \rho \alpha \beta }\widehat{T}_{(\alpha \beta )}^\sigma =0$$
(41)
with the constraints

$$\widehat{T}_{(\alpha \beta )}{}_{}{}^{0}=0,\widehat{T}^{(0i)k}=\frac{R}{2P+Q}(ϵ^{imn}\widehat{T}_{(mn)}{}_{}{}^{k}+ϵ^{kmn}\widehat{T}_{(mn)}{}_{}{}^{i}).$$
(42)
Therefore, all the degrees of freedom can be expressed in terms of the nine components $`\widehat{T}_{(ij)k}`$, which satisfy the four constraints

$$\widehat{T}_{(ki)}{}_{}{}^{i}=0,ϵ_{ijk}\widehat{T}^{(ij)k}=0.$$
(43)
and hence, as expected, five degrees of freedom remain. Using these constraints the equations of motion become

$$\left(\frac{1}{4}(Q+2P)+4\frac{R^2}{Q+2P}+\mu ^2\right)\widehat{T}^{(jk)l}=0,$$
(44)
which using (37) gives the mass $`\mu =M`$ for the $`\widehat{T}^{(jk)l}`$ field.

The divergence $`_\kappa T^{(\kappa \lambda )\sigma }`$ is non null if $`\left(2m+n4d\right)0`$. From Eq. (19) the on-shell duality relation can be written

$$T^{(\kappa \lambda )\sigma }=eC^1\left(\left(2m+n4d\right)\left(^\lambda h^{\sigma \kappa }^\kappa h^{\sigma \lambda }\right)+\left(2a+c\right)ϵ_{\mu \nu }^{\kappa \lambda }^\nu h^{\sigma \mu }\right),$$
(45)
and we see that the divergence is proportional to $`h^{\sigma \lambda }`$ on the equations of motion

$$_\kappa T^{(\kappa \lambda )\sigma }=eC^1M^2\left(2m+n4d\right)h^{\sigma \lambda }.$$
(46)
The linear combination of the tensors appearing on the right hand side of (45) is the most general form that can take the tensor $`T^{(\mu \nu )\sigma }`$ of symmetry $`(2,1)`$ as a function of $`h^{\sigma \kappa }`$ on the equations of motion. Taking $`\left(2m+n4d\right)=0`$, the second term on the right hand side gives the standard form of the duality involving the Levi-Civita tensor. For null spatial momentum this becomes $`T^{(ij)k}ϵ_l^{ij}h^{lk}`$ . This is the duality, in the sense of representations of $`SO(n)`$, for the tensors of types $`(1,1)`$ and $`(2,1)`$ in three spatial dimentions. On the other hand the term proportional to $`\left(^\lambda h^{\sigma \kappa }^\kappa h^{\sigma \lambda }\right)`$ makes the symmetric part of the divergence $`_\kappa T^{(\kappa \mu )\nu }`$ non null and proportional to $`h_{\mu \nu }`$. It corresponds to giving $`T^{(\mu \nu )\sigma }`$ some longitudinal components, as can be seen in the zero spatial momentum where $`T^{(0i)j}h^{ij}`$. Thus, the physical degrees of freedom of $`h_{\mu \nu }`$ are mapped either on the null divergence part of $`T^{(\mu \nu )\sigma }`$ or on the symmetric part of the divergence $`_\kappa T^{(\kappa \mu )\nu }`$, and in the general case on a linear combination of both. The case where $`h_{\mu \nu }`$ is completely mapped on the divergence of $`T^{(\mu \nu )\sigma }`$ can not be achieved with the Lagrangian (17) because it would require $`\left(2a+c\right)=0`$, in which case the transformation becomes ill defined, but can be realized including the aditional duality term $`T_{(\mu \nu )\sigma }^\mu h^{\nu \sigma }`$. We note that although all the parent Lagrangians with three free parameters represent the same theory, they give place to different duality mappings, characterized by one parameter.

We can see that in the parity-conserving case $`2m+n4d=0`$($`R=0`$), the equations of motion and the Lagrangian constraints simplify to

$$(^2+M^2)\widehat{T}^{(\nu \rho )\sigma }=0,\widehat{T}^{(\nu \rho )}{}_{\rho }{}^{}=0,ϵ_{\kappa \nu \rho \sigma }\widehat{T}^{(\nu \rho )\sigma }=0,_\alpha \widehat{T}^{(\alpha \nu )\sigma }=0,_\alpha \widehat{T}^{(\nu \rho )\alpha }=0.$$
(47)
## V Massless spin-two dual theories

Here we will include only a few comments in relation with the zero-mass case, as a detailed discussion is postponed for a separate publication. To count the independent degrees of freedom we can use the Hamiltonian analysis of the zero-mass limit $`(P=Q=R=0)`$ of the Lagrangian (34). The $`(3+1)`$ splitting of degrees of freedom produces the coordinates $`\widehat{T}^{i00},\widehat{T}^{0ij},\widehat{T}^{ij0},\widehat{T}^{ijk}`$ together with their corresponding canonically conjugated momenta $`\mathrm{\Pi }_{i00},\mathrm{\Pi }_{0ij},\mathrm{\Pi }_{ij0},\mathrm{\Pi }_{ijk}`$, satisfying the standard equal-time non-zero Poisson brackets. The whole set of constraints, which include up to terciary constraints, is

$`\mathrm{\Gamma }^i=\mathrm{\Pi }^{i00},\mathrm{\Gamma }^{ij}=\mathrm{\Pi }^{0ij},\mathrm{\Gamma }=\mathrm{\Pi }4F^0,\mathrm{\Delta }=\widehat{T}^{0m}{}_{m}{}^{},\mathrm{\Delta }^i=\widehat{T}_{}^{i00}+\widehat{T}^{im}{}_{m}{}^{},`$ (48)
$`\mathrm{\Sigma }=ϵ^{ijk}_k\mathrm{\Pi }_{ij0},\mathrm{\Sigma }^{ij}=_m\mathrm{\Pi }^{mij}{\displaystyle \frac{1}{3}}g^{ij}_m\pi ^m,(\mathrm{\Sigma }^{ij}g_{ij}=0),`$ (49)
$`\mathrm{\Phi }_{jm}=_m^q\mathrm{\Pi }_{jq0},(^j\mathrm{\Phi }_{jm}=0,g^{jm}\mathrm{\Phi }_{jm}=0),`$ (50)
with $`\mathrm{\Pi }=\frac{1}{2}ϵ_{ijk}\mathrm{\Pi }^{ijk},F^0=\frac{1}{2}ϵ_{ijk}^i\widehat{T}^{jk0},\pi ^m=\mathrm{\Pi }^{mk}_k`$. In parenthesis we have indicated the identities satisfied by the corresponding constraints. Due to this we have $`31`$ independent constraints. The first class constraints turn out to be the $`91=8`$ traceless components of $`\mathrm{\Gamma }^{ij}`$ together with the $`931=5`$ independent components of $`\mathrm{\Phi }^{jm}`$. Thus, $`13`$ constraints are first class, while the remaining 18 are second class. This gives $`\frac{1}{2}(2\times 242\times 1318)=2`$ independent degrees of freedom in the configuration space.

Finally, a comment on the possible dual systems that could arise in the massless case. Let us recall the situation in the case of $`p`$-forms. Basically there are two ways of taking the zero-mass limit. One arises directly from the Lagrangian (8) and gives

$$_P=\frac{1}{2}L{}_{}{}^{}L+LdB.$$
(51)
The second possibility is to rescale the forms $`BM^1`$ $`B,`$ $`LML`$ in (8), which leads to the alternative massive Lagrangian

$$_P=\frac{M^2}{2}L{}_{}{}^{}L+LdB+\frac{1}{2}B{}_{}{}^{}B,$$
(52)
together with the corresponding zero mass limit

$$_P=LdB+\frac{1}{2}B{}_{}{}^{}B.$$
(53)
In the case of $`M0`$ the theories generated by the Lagrangians (8) and (52) are equivalent. Nevertheless, for $`M=0`$ the situation is different. In fact, the Lagrangian (51) gives the equation of motion $`dL=0L=dA,`$ where $`A`$ is a $`(dq2)`$-form, thus producing a duality of the type $`q(dq2)`$. In an analogous way, the Lagrangian (53) produces the equation of motion $`dB=0`$, thus introducing a $`(q1)`$-form $`C`$ such that $`B=dC.`$ Here the duality is of the type $`(q1)(dq1).`$ Both types of dualities are not equivalent. The two possible dualities mentioned above arise from the mixing term $`LdB`$, plus the choice of the auxiliary field. In our case, the corresponding mixing term is

$$eT_{(\mu \nu )\sigma }ϵ^{\mu \nu \alpha \beta }_\alpha h_\beta ^\sigma ,$$
(54)
which gives the alternative field equations

$$_\alpha (ϵ^{\mu \nu \alpha \beta }T_{(\mu \nu )}{}_{}{}^{\sigma }+ϵ^{\mu \nu \alpha \sigma }T_{(\mu \nu )}{}_{}{}^{\beta })=0,ϵ^{\mu \nu \alpha \beta }_\alpha h_{\sigma \beta }=0.$$
(55)
Let us emphasize that they do not correspond to a condition of the form $`dB=0`$, whose solution is given in terms of the Poincare lemma. Recent generalizations of the Poincare lemma for mixed symmetry tensors  could be useful to find the general solutions of the equations above, leading to an explicit construction of the corresponding dual systems.

Acknowledgments: This work was partially supported by CONICET-Argentina, CONACYT-México, and Universidad Nacional Autónoma de México under grant DGAPA- UNAM IN-100397. RM acknowledges interesting talks with C. Fosco and LFU also acknowledges useful conversations with A. García.


format: text

text:
# Layer Formation in Semiconvection

## 1 Introduction

Though semiconvection may be an important mechanism for mixing of chemicals in late stages of stellar evolution, there have been few numerical simulations of semiconvective scenarios. In part, this is due to the fact that the numerical resolution required for a simulation has, until recently, made computation prohibitive. Coupled to that, the fact that the existing mixing length theories (such as Stevenson (1979), Langer et al. (1983) and Spruit (1992)) which both draw on theoretical and experimental arguments, seemed adequate when compared to observations.

However, recent observations of SN1987A have been at odds with models of late stages of stellar evolution and one main suspect are the theories of elemental mixing used in those models (Maeder & Conti (1994)). As with all turbulence models, the more sophisticated semiconvective models (Canuto (1999), Grossman (1996), Grossman et al. (1993)) which arose in the wake of these new data require closure assumptions. Implicit in these are assumptions about the character of the flow which have yet to be tested against numerical experiments of semiconvection. Numerical simulations are now sophisticated enough that we can hope to distinguish between the assumptions made in the mixing length theories and settle the question of the heat transport, chemical mixing and velocity field associated with semiconvective zones. In this work, I shall address one particular question related to all three of these topics: does layered convection develop when a compositionally stratified fluid, with very low viscosity, is subjected to a vertical heat flux?

### 1.1 Semiconvection and Double Diffusion

Kato (1966a) performed a local analysis of the convective instability of gases with a varying compositional gradient subjected to an upward heat flux and showed three criteria for instability. If the molecular weight increases upward, $`d\mu /dz>0`$, the fluid is unstable to the Rayleigh-Taylor instability. If

$$_{rad}>_{Led}=_{ad}+\frac{\beta }{43\beta }_\mu $$
(1)
then the fluid is exponentially unstable to classical convection. Here, $`_{ad}`$ is the adiabatic temperature gradient, $`_{ad}=(\mathrm{ln}T/\mathrm{ln}p)_S`$ and $`_{rad}`$ is the corresponding temperature derivative due to radiative diffusion. $`_\mu =\mathrm{ln}\mu /\mathrm{ln}p`$, $`\mu `$ is the molecular weight and $`\beta `$ is the ratio of gas pressure to the total pressure. This criterion was discovered by Ledoux (1947) and bears his name. Finally, Kato demonstrated that when

$$_{ad}<_{rad}<_{Led}$$
(2)
convection sets in through overstable oscillations. The physical explanation of this instability requires considering blobs of fluid rising in a hydrostatically stratified medium. If the medium is stable to the Ledoux criterion, such blobs will fall back to their initial positions. If, further, we allow heat to diffuse from the blobs (by far the fastest diffusing process in astrophysical contexts) then the returning blobs will have lost a small amount of heat on their round trip and will overshoot their initial positions. This mechanism relies on the existence of thermal diffusion and the growth rate of overstable modes is in fact proportional to the thermal diffusivity. Thus, overstable convection grows much more slowly than classical thermal convection and it would be reasonable to expect that the convection, itself, is less vigorous than its classical counterpart.

The weakness of linear oscillatory convection may be a completely irrelevant fact since the bifurcation that leads to overstability is subcritical. That is to say, the fluid can go unstable to small, but finite amplitude perturbations, at values of the temperature gradient far below that required for overstability. This is a weakly non-linear effect and is due to the fact that a small rearrangement in the molecular weight (due to breaking gravity waves, shear instability, etc) can drive the fluid to classical thermal convection, since the stabilizing compositional gradient may be locally destroyed.

Weakly non-linear analysis was performed by Veronis (1968) for an incompressible, Boussinesq fluid with a stabilizing salt gradient and destabilizing temperature gradient. Though superficially appropriate for only the oceanographic situation for which it was derived, it was pointed out by E.A. Spiegel (see Kato (1966a)) that the linear stability equations of semiconvection are equivalent to the incompressible equations with the salt playing the role the molecular weight. In the oceanographic context such overstability is referred to as double diffusion.

The effect of the redistribution of the stabilizing component is most dramatically seen in the now famous experiment by Turner & Stommel (1964) and reproduced by Turner (1968), Huppert & Linden (1979) and Fernando (1987) and Fernando (1989). In the prototype and all of its successors, a tank of water with a stable salinity concentration gradient is heated from below. At first the water becomes unstable to rapidly rising convective plumes but, since they are rising into a region of less salt concentration, the plumes are restricted to a thin region at the bottom of the tank. The salt is rapidly mixed thus there is no gradient in this region and the fluid motion becomes turbulent. Subsequently, the layer proceeds to slowly entrain from the nonconvecting fluid above it while transporting heat to that region and at all times increasing in temperature. After a critical layer thickness is reached, the boundary layer separating the convecting fluid from the quiescent fluid above is seen to undergo oscillatory instability. This, in turn, redistributes the salt in the second convecting layer which is now separated from the first by a buoyancy jump. The convection becomes turbulent in this second layer and the process continues in subsequent layers. At intermediate times, the tank has formed a vertical stack of convecting layers with well mixed temperature and salt within each. Though mixing between layers is limited by the saline density jump across each, the lower layers are able to slowly erode this density jump over much longer time scales by the breaking of internal waves on the interfaces between the layers.

Fernando (1987) & Fernando (1989), provide the most accurate quantitative data for the water-salt experiment. In particular, he confirms the existing theory for the development of the first layer and gives a criterion for the stability of a layer. Defining $`l`$ as the integral length scale of the turbulence, the scale at which the turbulence is forced, $`u_{}`$ as the root mean square turbulent velocity in the convecting layer, $`\rho _{}`$, the average density of the convecting layer and $`\mathrm{\Delta }b=g(\mathrm{\Delta }\rho )`$ as the buoyancy contrast between the convecting layer and the fluid above, then

$$Ri=\frac{l\mathrm{\Delta }b}{\rho _{}u_{}^2}$$
(3)
is the Richardson number of the flow. This parameter generally arises in the study of a stable buoyancy gradient in the presence of a horizontal shear flow. In the layered convection case it amounts to the ratio of the potential energy difference over the scale of the largest eddies in the turbulent flow to the kinetic energy of the flow. Fernando showed that, in the case of water-salt experiment, there are three regimes of interest delimited by the Richardson number.

At low Richardson numbers, the turbulent eddies are able to penetrate deeply into the static fluid and the interface migrates upward. At intermediate Richardson numbers, mixing across the stable density interface is primarily due to the breaking of interfacial (internal) waves and the scouring of this surface by a weaker flow that can be maintained in the upper fluid. At higher Richardson number, the effect of interfacial waves and penetration becomes negligible as convective plumes flatten out when they hit the density interface. Transport across the layers then takes place purely by diffusion, though enhanced by the steep gradient at the interface.

Only when the Richardson number, based on the uppermost layer and the static fluid, is of intermediate or high values does there develop an oscillatory instability in the boundary between those two regions: another layer forms.

Let us now return to astrophysical semiconvection, contrasting the semiconvective theories of Spruit (1992) and Stevenson (1979). Using the intuition of the laboratory experiments, Spruit constructed a mixing length theory of layered convection separated by thin molecular weight gradients and thicker thermal gradients. Transport across the layers is solely due to diffusion through the thin thermal and helium transition regions at the interfaces of the layers. In fact, this theory is similar to the double diffusive theory constructed for confined geometries by Knobloch & Merryfield (1992). However, confined geometries cannot support surface waves like the internal waves in layered convection and these could contribute to even greater cross interface transport.

Stevenson (1979) described a theory of semiconvection in which growing overstable modes resonantly feed energy to smaller scales and eventually break, whereupon the molecular weight gradient is redistributed. By making reference to the salt-water experiments and energetic arguments for double diffusive convection, he compellingly argued that, should layers form as a result of this wave breaking, they would be unstable if $`Pr<\tau ^{1/2}`$. The Prandtl number, defined as $`Pr=\nu /\kappa _T`$, and the Lewis number, $`\tau =\kappa _{He}/\kappa _T`$, are ratios of the kinematic viscosity and molecular diffusivity of the second species to the thermal diffusivity, respectively. (The two species are helium in a hydrogen convection zone.) In semiconvection, the thermal diffusivity is dominated by radiation whereas both viscosity and helium diffusion are molecular processes and, as such, $`1Pr>\tau `$. In contrast to Spruit, Stevenson’s theory would lead to the conclusion that layers are dynamically unimportant in stellar evolution.

Further evidence of layer instability is provided by the only numerical experiment on semiconvection to date, Merryfield (1995). Using the anelastic approximation for a two dimensional fluid, Merryfield studied the onset of overstable oscillations in a low Prandtl number fluid: the thermal parameters that he used put the simulations within the regime of equation (2). Although motion set in though overstable oscillations and the compositional gradient was seen to mix through the breaking of internal waves, these tended to occur on length scales comparable to the domain, and layers did not form. In one of the experiments, an initial condition of two layers separated by a stable density interface was simulated. The interface became quickly unstable to internal waves, disrupting the layers and leading to well mixed, global thermal convection. As yet, layers have not been observed to form, much less be stable, in any semiconvective simulation at low Prandtl number.

On the other hand, layers have been observed to form in a numerical experiment tailored to explaining the observations of Huppert & Linden (1979) and Fernando (1989). Molemaker & Dykstra (1997) considered an incompressible fluid at $`Pr=7`$ (water) with vertical salinity gradient subjected to cooling from above. Due to Boussinesq symmetry, this is completely equivalent to heating from below as is done in laboratory experiments. The numerical simulation well reproduced the observations of Fernando (1989) and confirmed the measurement of a critical Richardson number delineating the transition from layer growth by entrainment to quasi-static convective layers separated by diffusive interfaces.

### 1.2 Current Work

The numerical and theoretical evidence provided thus far verifies the formation of layers in the high Prandtl number regime appropriate for water. However, the results for low Prandtl number are still inconclusive. Does a multiple layered structure form in the low Prandtl number regime? If so, how does the critical Richardson number depend on Prandtl number?

The formation of a second layer above the first requires a quiescent thermal boundary layer separating the convective zone from the static fluid. In this thermal boundary the slow oscillatory instability can incubate and grow before convective erosion from below becomes significant. A smaller Prandtl number means that the viscous boundary layers separating static from convective zones is smaller than the thermal boundary layer. Convective plumes will tend to decelerate less before reaching the buoyancy interface and be able to penetrate more energetically. Smaller Prandtl number also means smaller scales for kinetic energy dissipation. The interfacial waves which are generated by the convective plumes will tend to break at smaller wavelengths. The interfacial waves will then be more effective at mixing the composition gradient across the interfaces. Both the energetic waves and smaller scale breaking tend to suppress the growth of oscillatory instability in the static fluid.

In the present work I present results of two dimensional numerical simulations of a fully compressible fluid. The molecular weight of the medium decreases upward and the diffusivity of the second species is always small. Results for successively decreasing Prandtl number and at different values of the heat flux are compared.

In section 2 the equations, choice of parameters and numerical technique are presented. The results of the simulations are presented and described in §3 and these results are compared to previous mixing length theories and numerical simulations in §4. In §5 the results are interpreted in the context of stellar models.

## 2 Formulation of the problem

### 2.1 Equations and Boundary Conditions

This work attempts to bridge the divide between terrestrial double-diffusive convection and stellar semiconvection simulations. To this end I have, at times, neglected certain processes which are less relevant to the formation and disruption of layers, though they are obviously present in realistic stars. Radiation pressure tends to decrease the contrast between the semiconvective and convective stability thresholds (equations 1, 2) and would correspondingly mitigate layer formation. Although radiation pressure is not negligible in stellar convection, for the purposes of this work, it has been neglected. This serves the added purpose of singling out differences between astrophysical and terrestrial double diffusion, since the ideal gas equation of state at constant pressure tends to the Boussinesq approximation for incompressible fluids. Any formation of layers in the present simulations can be seen as a necessary, but not sufficient condition for them to form in stars.

The equations of mass, momentum, energy and helium concentration conservation which govern a compressible fluid of hydrogen gas with a varying concentration of helium, $`Y`$, are (Spiegel & Veronis (1960))

$$\frac{\rho }{t}=(\rho 𝐯)$$
(4)
$$\rho \left(\frac{𝐯}{t}+(𝐯)𝐯\right)=p+𝐠\rho +\nu \rho \left(^2𝐯+\frac{1}{3}(𝐯)\right)$$
(5)
$$\rho C_v\left(\frac{T}{t}+𝐯T\right)=p(𝐯)+(\kappa T)+\nu \rho ((𝐯𝐯)\frac{2}{3}(𝐯)^2)$$
(6)
$$\frac{Y}{t}+𝐯Y=\frac{\kappa _{He}}{\rho }(\rho Y)$$
(7)
$$p=\frac{\rho RT}{\mu },$$
(8)
where $`\nu ,\kappa ,\kappa _{He}`$ are the constant kinematic viscosity, thermal diffusivity and helium molecular diffusivity, respectively. $`R`$ is the gas constant, $`\mu (Y)`$ is the molecular weight for an ionized gas of hydrogen and helium

$$\mu =(2\frac{5}{4}Y)^1$$
(9)
and the specific heat at constant volume is

$$C_v=\frac{3}{2}\frac{R}{\mu }.$$
(10)
Computations are performed on a rectangular box of vertical dimension, $`d`$ and horizontal dimension $`\lambda d`$ with downward directed, constant gravity. The initial conditions for the simulations are constant, negative helium gradient throughout the domain. In the bulk of the domain the temperature is chosen such that the fluid is only marginally unstable to the Schwarzschild criterion but very stable according to Ledoux.

An additional heat flux is supplied at the bottom of the domain and the temperature gradient is continuously set to match it, exponentially, over a very thin region. This has the effect of increasing $`_{ad}`$ to the point where the fluid is significantly unstable in a thin region at the bottom. This transition region is necessary in order to avoid any numerical instabilities which arise from impulsive heating.

In all examples, the helium concentration varies from $`Y=1`$ at the bottom of the domain linearly to $`Y=0`$ at the top.

The side boundaries are periodic while the following top and bottom boundary conditions are used,

$`v_z=_zv_x=0`$ $`\mathrm{at}`$ $`z=0,d`$ (11)
$`_zY=0`$ $`\mathrm{at}`$ $`z=0,d`$ (12)
$`_zT=\alpha `$ $`\mathrm{at}`$ $`z=d`$ (13)
$`_zT=(\alpha +F)`$ $`\mathrm{at}`$ $`z=0.`$ (14)
Physically, the velocity boundary conditions correspond to impenetrable, stress free boundaries. The helium boundary conditions are chosen so that, when convection sets in, there will not be sharp boundary layers between the convective region and the boundaries. Such layers would be computationally expensive to resolve and be irrelevant to the physics. The initial temperature profile is given by $`T(z,x,t=0)=T_d\alpha (zd)+Fz_{}e^{\frac{z}{z_{}}}`$ where $`\alpha `$ and $`F`$ are both positive and $`z_{}d`$. $`z_{}`$ is chosen so that the thermal boundary layer at the bottom is very thin and the temperature at the lower boundary is dominated by the $`Fz_{}`$ term in the initial condition. That the heat flux is so low at the top boundary may seem to be cause for concern, since the driven convection will continue to accumulate energy faster than it can dissipate it from the top boundary. In practice, if the temperature changes significantly at the top boundary, it is due to the fact that a convective plume has arrived from the lower boundary. In such circumstances, there is no layer formation in the domain (only global thermal convection) and the simulations are stopped.

### 2.2 Numerical Method

The code is based on one developed by Rubini et al. (1997) and employs a mixed spatial scheme: a psuedo spectral, periodic routine in the horizontal direction and $`6^{th}`$ order compact finite difference in the vertical direction. The latter is more appropriate for the boundary conditions under consideration and the accuracy of the derivatives is very close to the spectral accuracy of the horizontal derivatives if double the resolution is used.

The time step for low Mach number compressible simulations is limited by the CFL condition for sound waves, $`\mathrm{\Delta }t<\mathrm{\Delta }x/c_s`$ where $`\mathrm{\Delta }x`$ is the mesh spacing and $`c_s`$ is the sound speed. The condition corresponds to the requirement that sound waves must be tracked across each grid element. In all cases considered here, this condition is more restrictive than the corresponding condition for thermal diffusion, $`\mathrm{\Delta }t<(\mathrm{\Delta }x)^2/\kappa `$. Since sound waves play little role in the convection, it would be beneficial to disregard them and their CFL condition as is done in the anelastic simulations of Merryfield (1995).

The simulations herein are advanced in time by a second order Adams Bashforth scheme, which is explicit, and were therefore constrained by the sound CFL condition (see ) (Canuto et al., (1988)). The code employed was actually developed with the intention of overcoming just this constraint, having the capability of semi-implicit time-stepping. Under appropriate temporal discretization, this allows to jump the sound CFL condition (see Rubini et al. (1997) for details). Some simulations reported herein were reproduced with the semi-implicit method and the two results were found to agree quite well. Furthermore, after accounting for a factor of 10 speed increase due to the sound waves and a factor of 2 decrease due to the implicit method, the semi-implicit scheme was found to have a net factor of about 5 speed increase. Future work will be carried out using this method, but all of the current results are due to the explicit scheme.

The results were obtained from the SGI Origin 2000 machines, at Los Alamos and Argonne National Laboratories. Any given simulation used from 32 to 96 processors.

Finally note that all of the simulations are fully resolved and do not suffer from the need to filter numerical instabilities at the grid scale. Though time consuming, this means that no numerical diffusion has been introduced thus the Prandtl and Lewis numbers which are quoted in table 1 are the actual ones in the simulation.

### 2.3 Physical Parameters

The dynamical equations are solved by non-dimensionalizing the variables in the standard fashion. All lengths are scaled to the height of the domain, $`d`$, and time to the thermal diffusion time across the domain $`t_{diff}=d^2/\kappa `$. This naturally leads to the introduction of the Rayleigh numbers,

$$Ra_T\frac{gd^4(_{rad}_{ad})}{\kappa \nu H_p}$$
(15)
for heat and

$$Ra_{He}\frac{gd^4_\mu }{\kappa \nu H_p}$$
(16)
for helium.

These parameters arise from the consideration of the linear instability of a fluid configuration whose initial state corresponds to a constant molecular weight gradient and constant thermal gradient. Linear and weakly non-linear theory has been extensively studied with regard to the double diffusive problem (see Veronis (1965), Veronis (1968), and the thorough review by Turner (1985)), but the numerical experiments that follow drive the fluid far from the static equilibrium making both linear and weakly non-linear discussions irrelevant.

For example, even in a thin slot geometry where layering is seen to form in laboratory experiments (such as Biello (1997)) and weakly non-linear analysis is greatly simplified, Balmforth & Biello (1998) show that the subcritical nature of the instability prohibits the solutions from saturating at small amplitude and the layering phenomena are not predicted. Thus far, the only theories which successfully predict layers are the heuristic one-dimensional turbulence models of Kerstein (1999) and the dynamical mixing-length models of Balmforth et al. (1998). Though the latter was derived for shear forced layering it may hold insights for the buoyancy forced layering considered here. In light of the futility of non-turbulent approaches to the study of layers, I shall therefore dispense with further discussion of the predictions of linear and weakly non-linear theories and concentrate only on the highly non-linear problem at hand.

The Prandtl and Lewis numbers, already defined above, are the relative measures of viscous and molecular diffusion to thermal dissipation. $`C_s`$ is the non-dimensional sound speed of the medium and $`H_p=(\mathrm{ln}p/z)^1`$, the pressure scale height, has been non-dimensionalized in Table 1 but not in equations (15)-(16).

The nondimensional parameters are given in Table 1, $`\tau =0.1`$ in all simulations and $`Y(t=0)=1z/d`$ in all except $`VII`$ (see below). Astrophysically realistic values have Rayleigh numbers in the range $`10^{12}10^{14}`$, $`\tau 10^6`$, $`Pr10^5`$ and $`C_s10^5`$. To perform a completely realistic simulation of these parameters is a computationally impossible task as we would have to resolve the smallest gradients in the the domain (in this case, due to molecular diffusivity) and track sound waves across them.

Keeping in mind that our primary aim is a study of the the robustness of layers in the limit of small Prandtl number, we have chosen fluid parameters in accord with that goal, limited by computational feasibility and with an ordering relevant to astrophysics. Therefore, $`\tau <\sigma <<1`$ preserves the relative sizes of the gradients of helium, velocity and temperature, respectively. $`C_s>>1`$ in all simulations and, correspondingly, the mach numbers of the flow were also $`<1`$ (typical maximum mach numbers are $`0.4`$).

The numerical experiments are organized as follows. $`I`$ and $`II`$ consider large $`(>1)`$ Prandtl numbers for varying forcing and domain size. In particular, $`I`$ is most similar to the Boussinesq regime and will provide the standard against which the possibility of layer formation in all of the subsequent experiments, which are run at small Prandtl number, will be assessed. Simulations $`IIIa,b`$ contrast two failed attempts at producing layers at weak and strong forcing. This analysis of this work focuses on the results of $`IV`$ and $`V`$ where a single layer is formed at low Prandtl number and migrates upward. In particular, the mixing associated with layer growth is studied in detail in $`V`$. In $`VI`$ an attempt is made to study the effect of large aspect ratio on the results of the previous two simulations. Finally, experiment $`VII`$ is set up as a very stable two layer system with convection setting in rapidly. In this setting, one can easily distinguish the mechanisms of boundary layer erosion, breaking internal waves and penetrative convection.

## 3 Results

The laboratory experiments on salt-water demonstrate that layer formation is necessarily a transient phenomenon. The time asymptotic state of any layered convective fluid would necessarily be one of well mixed composition undergoing classical convection. For convective layers to be dynamically important to stellar evolution, not only must they must be demonstrated to form by some relevant physical mechanism, they must do so on time scales short enough for their presence to be felt, and  they must persist over comparably long time scales. In late stages of stellar evolution, such time scales are of order the evolution time for a star. Irrespective, such times are always much longer than the thermal diffusion time, which is the reference time for all convection calculations.

Due to Rayleigh-Taylor instabilities of the thermal boundary layer at the bottom of the fluid the first layer is seen to form rapidly in the laboratory. The second layer forms from oscillatory instabilities of the internal boundary layer above the first. Since the oscillatory instabilities occur on thermal diffusive time scales, tracing their complete evolution would be computationally prohibitive with the current compressible simulation; to study them over stellar evolution times would be impossible. There is, however, an interesting result from Bruenn & Dinhva (1996) that double-diffusive instabilities may play a role in the supernova process, itself, which could render moot the necessity of considering large time scales. Nonetheless, all of the simulations reported hereafter cease at approximately $`10\%`$ of a thermal diffusion time; too little time to see a second well mixed layer.

To circumvent this problem, I shall use phenomena observed during layer formation in the laboratory as indication that a second layer will form. Firstly, as a necessary condition, the height of the interface separating the convecting fluid from the nearly stationary flow must be quasi-static, for experiments (in particular Fernando (1989)) it is seen to grow as $`ht^{1/2}`$. This is tantamount to the restriction that significant entrainment across the interface cannot occur, or else it would quickly migrate upward.

Secondly, overstability must been seen to set in above the convecting flow. This latter requirement cannot be underestimated, since a layer which is quasi-static over the time scales of the integration hardly can be regarded as a sufficient condition for a second layer to form. In particular, one can envision a scenario where the lower fluid develops turbulent convection and the instability in the upper fluid sets in so slowly as to make it imperceptible over these time scales. On time scales much longer than that observed in the laboratory, the convecting region could slowly entrain fluid from above, avoiding layer formation entirely.

In the descriptions below, I shall show when these phenomena are seen or not in the high Prandtl number simulations and use those results as a guide to study the low Prandtl number cases.

### 3.1 High Prandtl Number

#### 3.1.1 Simulation $`I`$: Layering in a viscous fluid.

Figure 1 shows the onset of instability at early times, for $`Pr=3.5`$. Since the molecular diffusion is small, the helium fraction, which began horizontally constant and decreasing linearly, will indicate the convective mixing by the flow. The color coding runs from lightest for $`Y=1`$ to darkest, $`Y=0`$, and is chosen because most of the convective activity will occur in regions of high helium concentration; at the bottom. At $`0.715\%`$ of a thermal diffusion time, shown in the figure, the Rayleigh-Taylor plumes are clearly evident. The symmetry of the flow betrays the symmetry of the initial conditions which were chosen in order to highlight the relevant mechanisms in this illustrative test case.

In the bottom left panel are plotted the horizontal averages of density $`(\overline{\rho })`$, temperature $`(\overline{T})`$ and helium fraction $`(\overline{Y})`$ along with their initial conditions, respectively, from top to bottom. The Rayleigh-Taylor instability is an effect of the initial density inversion seen at the bottom of the fluid. It is also clear that the flow rapidly smooths out that inversion and layering becomes evident below $`z=0.3`$ where $`\overline{Y}`$ and $`\overline{T}`$ flatten. Moreover, owing to the relative diffusivities, the transition from mixed fluid (little or no gradient) at the bottom to constant gradient at the top is much sharper for $`\overline{Y}`$ than for $`\overline{T}`$.

The final panel shows the local value of $`ϵ_S=_{ad}`$ (Schwarzschild criterion) and $`ϵ_L=_{Led}`$ (Ledoux criterion) determined using the horizontally averaged thermodynamic quantities and compared to their respective initial values. Initially $`ϵ_L`$ is everywhere negative except for the temperature transition region at the bottom where it rises to coincide with the maximum value of $`ϵ_S`$ (since $`_\mu =0`$ near the bottom boundary). Conversely, $`ϵ_S`$ is only slightly supercritical in the bulk of the domain and also rises at the lower boundary. Though the region above the transition layer is stable, the momentum of the Rayleigh-Taylor plumes allow them to penetrate deeply and consequently homogenize the temperature and helium.

Figure 2 shows a well developed, almost fully mixed lower layer at $`t=2.74\times 10^2`$. The figure does not reval the steady internal wave of wavelength equal to the domain width and height, $`h`$, $`0.25`$ which is also present. In the horizontal averages there is a flat distribution of helium and almost flat distribution of temperature below $`z=0.4`$. The internal helium boundary layer, extending $`z=[0.45,0.5]`$, is of thickness $`\delta _{He}0.05`$ whereas the internal temperature boundary layer, $`z=[0.4,0.55]`$, is of thickness $`\delta _T0.15`$. This is consistent with the predictions of Spruit (1992) that the relative thickness of diffusive internal layers of temperature and helium in layered convection should be of order $`d_{He}/d_T\tau ^{1/2}(1/3`$ for all of the results herein). Such estimates are straightforward to derive by arguing that the convection homogenizes the interior and a convective turnover time sets the scale over which which material (or heat) can diffuse out of a fluid blob at an internal layer (or boundary layer, see Shraiman (1987)).

A sign of efficient convection is that the fluid be almost adiabatic, and this is evidenced by the last panel of figure 2 where $`ϵ_Tϵ_S0`$ over $`z=[0.15,0.4]`$. Furthermore, the internal thermal transition layer and the entire region above it is unstable to oscillations yet stable to overturning convection. In fact, a weak two cell oscillatory flow is seen in the velocity vector plot and will be elaborated upon below.

Three mechanisms contribute to mixing helium across the interface, the first being molecular diffusion enhanced by the sharp gradients there. From figure 2 it is also clear that there are significant regions where low helium concentration fluid is down-welled into the convection zone. Such regions are a remnant of both the original Rayleigh-Taylor plumes and the large amplitude interfacial waves which they excite. As time proceeds the lower layer becomes more well mixed, the interface sharpens and mixing by down-welling becomes less relevant. Finally, the scouring of the helium internal layer is a major source of vertical helium transport at all times.

The layer height is determined from $`\overline{Y}`$ by finding the height where $`|\mathrm{ln}\mu (\overline{Y})/z|`$ is maximum. Plotted in the first panel of figure 3 is the interface height as a function of time. After an initial growth phase, and except for internal wave oscillations, the layer height is seen to remain relatively constant with time. This implies that the sum of the interfacial mixing mechanisms do not contribute enough to destabilize the interface. Therefore the necessary condition for the formation of a second layer, a quasi-static interface, is satisfied.

The third mixing mechanism, internal layer entrainment, has the added effect of driving oscillations in the upper fluid in the following manner. As the interfacial wave rises, it drives an upward flow directly above the wave crest. Since the flow above the interface rises on hot plumes from below, this behavior is often referred to as the buoyancy coupling of two layers. When the wave turns downward, the flow in the upper domain continues upward and, through shear drag, is able to strip some of the helium rich fluid from the interface. Locally, fluid that has been stripped has greater helium concentration than its environment, yet due to its momentum it continues upward. Before being able to overturn completely, the density excess decelerates the flow and the roll turns around: this is the oscillatory instability of the upper fluid.

The helium rich, upward flowing fluid at $`(x,z)=(1.0,0.75)`$ in figure 2 is a product of this internal layer stripping. Plotted in figure 4 is the time trace of the vertical velocity field and $`500`$ times the helium deviation from its initial value at this point. The solid line clearly shows growing oscillations with period $`1\%`$ of a diffusion time and a longer growth time. The helium deviation is slightly out of phase with the velocity, which drives the oscillations. Parenthetically note that the ripples on the velocity trace clear at $`t=2\%`$ and elsewhere are due to sound waves.

Though a layer has not yet fully developed by the end of the simulation $`I`$, one is obviously in the process of doing so. Since it contains both phenomena which are observed in laboratory layer formation, I conclude that it would likely form one if allowed to continue.

### 3.2 Simulation $`II`$: A vortex/shear layer.

The important differences of the second simulation are the lower aspect ratio, weaker forcing and initial conditions. For this run the initial velocity is set to a random spectrum of very low amplitude vortical waves with a wavenumber cut off at an intermediate length scale $`(k_{max}20)`$.

Again, the instability sets in through isolated, randomly spaced Rayleigh-Taylor plumes rising from the lower boundary. However, soon the flow becomes coherent and mixes the helium well in the lower half of the domain. At $`t=0.07`$ figure 5 shows that the interfacial layer height has saturated at $`z=0.8`$.

The flow in figure 6 is exceedingly elegant in its simplicity and merits some elaboration. The stable concentration jump is clear at $`z=0.75`$ in that $`ϵ_L<0`$, while $`ϵ_S>0`$ indicates that overstable oscillations should eventually set in. Moreover, an interfacial wave of the size of the domain propagates to the right. In the bulk of the well mixed region a strong shear layer has spontaneously arisen, reminiscent of the shear flow seen at high Rayleigh number Rayleigh-Bénard convection (though much less turbulent).

The shear at the bottom boundary aligns the flow so that heat fluxes into the fluid through one coherent thermal, located at $`x0.8`$ in the figure. Being carried by the shear, this thermal propagates leftward.

The entire shear layer is not stationary, rather it oscillates between closed flow lines and a single closed vortex, the eye of which is reforming at $`(x,z)=(0.15,0.25)`$. Merryfield (1995) produced a coherent structure in a weakly forced simulation, though in that case the vortex was a double gyre. In dynamical systems literature, an oscillating jet/vortex flow is referred to as separatrix reconnection, and takes place when a propagating wave is modulated at another frequency. In this case, the second period is provided by the leftward propagating thermal at the base of the flow which interacts with the rightward propagating interfacial wave above.

Jets of collimated flow, like the one below the interface in figure $`3.6`$ are notoriously stable to advective transport across them, and one would expect molecular diffusivity to be the primary mode of helium mixing across the interface. However, separatrix reconnection provides a mixing channel so that, even in this simple flow, it is able to entrain helium poor fluid through the trough of the interfacial wave (particularly evident in the plot of helium fraction).

Finally, take note of the small knee in the temperature average and the corresponding peak in $`ϵ_S,ϵ_L`$. This is the imprint of the coherent thermal which, owing to the vortex and shear, mixes with the rest of the flow through discreet release events.

Though performed at high Prandtl number, this simulation holds an important lesson for astrophysical layering. Even in a weakly forced scenario where the interfacial waves are very low amplitude, mixing by separatrix reconnection will provide much greater cross interface helium transport than molecular diffusion enhanced by the sharp gradient there.

### 3.3 Low Prandtl Number

#### 3.3.1 Simulations $`III`$: No layering

Simulations $`IIIa,b`$ are presented as two cases in which, according to the analytic theory described in §5.1, interfaces should form. In both cases layered convection does not occur, but for different reasons.

In $`IIIa`$ the stabilizing compositional gradient is not sufficient to impede the rising Rayleigh-Taylor plumes. Consequently, they hit the top boundary before turning back and eventually coalescing. By $`t=0.061`$ in figure 7 the helium is already well mixed and the convection is global.

A relation for the quasi-static interfacial layer height in terms of the heat flux, the initial buoyancy frequency of the domain and empirically determined factors was given by Fernando (1989) and will be discussed at length below. Suffice for the moment to say that the critical height varies as the square root of the heat flux at the lower boundary and is thus proportional to $`\sqrt{Ra_T}`$ there. Guided by this and the previous failure, a heat flux one fifth that of $`IIIa`$ is chosen, with the expectation that a layer should saturate within the height of the domain. However, by $`t=0.07`$ in $`IIIb`$ the lower thermal boundary layer has diffused to $`z=0.5`$ though convection has yet to set in. Correspondingly, the whole region $`z<0.5`$ has become Ledoux unstable. It is reasonable to conclude that the growth rate of the plumes is too low to observe convection before the Ledoux unstable region encompasses the entire computational domain. Global convection would thus ensue.

This is an interesting result, but probably less relevant for astrophysics on two counts. First, astrophysically high Rayleigh numbers yield correspondingly large growth rates. Second, the simulation was initialized with a small velocity field; more realistic in nature would be a spectrum of finite amplitude gravity waves, allowing to avoid the, somewhat long, linear growth phase.

#### 3.3.2 Simulations $`IV`$, $`V`$, $`VI`$: Layer formation and interface migration

Stars have both small Prandtl number and exceedingly large Rayleigh numbers and, to the limit of computational feasibility, these conditions are reproduced in simulations $`IV,V`$ and $`VI`$.

Figures 8 \- 10 are the $`t=0.025,`$ $`0.035`$ and $`0.06`$ snapshots of simulation $`IV`$ in which a small initial velocity spectrum is used to seed the Rayleigh-Taylor instability. By $`t=0.025`$ the plumes have reached $`z=0.5`$ and slowed significantly. Notice, though, that smaller plumes with a large helium contrast are still able to rise from the interface of the main plumes; this plume splitting is seen at the center of the $`Y`$ concentration plot. Though mostly having stopped near the midplane the plumes are still slightly buoyant, therefore the primary mechanism slowing them must be their drag.

Even though viscosity is weak, it is well known that drag can be significantly enhanced by the adiabatic expansion of rising plumes, the entrainment of fluid into the plumes and the breakup of the plume tip. Since it is clear that parts of the plumes remain buoyant after breakup and continue to penetrate up to $`z=0.7`$ then, necessarily, buoyancy and thermal diffusion out of the plumes play a smaller role in their deceleration up to the midplane.

In figure 9 the plumes have turned around and begin to resemble convecting cells mixing the bottom half of the domain. However, during the deep penetration, there is significant mixing of helium rich fluid up to $`z=0.8`$, the shadow of which are the faint plumes seen in the upper half plane. Like the scoured interfacial layer in simulation $`I`$, these faint plumes have a $`Y`$ greater than their surroundings and can be expected to descend. Also like their counterparts, these plumes are in a region which is apparently oscillatorily unstable. However, a clear interface has yet to form above the convecting region and no oscillatory instability is apparent for $`z>0.5`$ at $`t=0.035`$.

From interface formation to $`t=0.06`$ in figure 10 the flow becomes turbulent, but with a persistent large scale component. In turn, these few large vortices and plumes drive an interfacial wave of wavelength approximately equal to the domain length. At this point $`ϵ_L`$, $`ϵ_S`$ are irrelevant measures of convective stability as the interface is clearly not horizontal. Unlike the high Prandtl number case, this interfacial wave is of very large amplitude and all of the mixing mechanisms thus far discussed take place across it. Needless to say, molecular diffusion enhanced by the interfacial gradient plays a relatively small role in the total transport.

Most obvious from figure 10 is the strong down-well obliquely crossing the domain and upon which a helium rich internal wave is about to break. This is a more dramatic version of the separatrix reconnection demonstrated in $`II`$. Traces of previous mixing events are seen as helium poor filaments in the otherwise mixed region. Less obviously, the upper left of the thick helium boundary layer is experiencing scouring by a flow along the interface (around $`(x,z)=(0.5,0.9)`$). Moreover, a down-flow near the mouth of the down-well (at the center of the domain) is clearly shearing two helium lamina, though where they will mix is not clear. Irrespective, this flow is violently mixing the domain and by $`t=0.07`$ the interfacial wave collides with the upper boundary so that no second layer can be formed.

Interface height is plotted against time in figure 11. Except for the last fraction of a diffusion time ($`10^3`$), the interface compares well with that of figure 3. The spikes betray the large amplitude internal waves on the interface and the collapse of the interface near the end reflects the overturning convection which sets in at that time.

Each of the dynamic mixing mechanisms have at their source the fact that the amplitude of the interfacial wave is very large, comparable to its wavelength. This, in turn, may be a consequence of both the fact that the flow under consideration is 2-dimensional and that the aspect ratio of the domain is small, $`\lambda =2`$. In 3-d simulations this inverse cascade would not occur while, since rolls will not likely have wide aspect ratio themselves, it would be irrelevant for large aspect ratio 2-d boxes.

In simulation $`VI`$ an attempt is made to study the effect of large aspect ratio on, in particular, the amplitude of the interfacial wave. Unfortunately, the simulation lost resolution and had to be ended before an interface was well developed. Nonetheless it is evident from the temperature field shown in figure 12 that rising thermals are well spaced along the bottom boundary and there is some hope that, should many weak thermals persist, the flow would not produce such large scale structures particularly, the catastrophic internal wave.

The most promising case for layer formation is provided by simulation $`V`$. Akin to $`IV`$, it has lower heat flux at the bottom boundary and larger aspect ration ($`\lambda =3`$), the objective being that the layer height should saturate at smaller $`z`$ and that the convection should not be dominated by large structures. Figure 13 shows that this is indeed the case since by $`t=0.05`$ the interface has reached a quasistationary height of $`z=0.35`$. Notice that the time is much later than in previous simulations since the flow is much weaker. The thermodynamic criteria are favorable to onset of oscillations just above the interface and already at this early time such oscillations are seen in the fluid above. Despite weak convection in the lower layer, the interface in figure 13 is strikingly irregular. Though the interfacial waves remain small, mixing is again dominated by large structures being engulfed between the plumes below the interface.

By $`t=0.096`$ the interface has continued to migrate upward, and the density jump is below $`z=0.6`$ in figure 14. Interfacial scouring due to oscillatory convection in the upper fluid creates filaments of helium rich fluid above the interface and mixing below remains dominated by large blobs of fluid between the plumes.

Large plumes have penetrated to $`z=0.8`$ by $`t=0.136`$ (figure 15) effectively destroying the oscillatory flow which had been growing there. Two plumes dominate the flow at this time and are clearly able to entrain large, helium poor fluid elements downward into the convection as is seen at $`(0.6,0.3)`$ in the figure. Interfacial splashing, evident in a thin filament above this blob, overtakes scouring as the principal mechanism for helium transport above the interface.

Though the vestiges of a layer are evident in the mean helium field, the interfacial waves have grown large enough to nearly reach the top boundary, suggesting that a second layer will not form before the interface grows to encompass the whole fluid. Such waves, though less marked, are clear in figure 16. Nonetheless, this is the only case discussed thus far where the theoretical prediction of Fernando (1989), $`h\sqrt{t}`$, provides an excellent fit. The curve $`h=1.8\sqrt{t}`$ agrees well with the data after $`t=0.01`$ and despite the growing interfacial waves at late times. Furthermore this fit yields an effective diffusion coefficient across the interface of $`\tau _{eff}=1.62`$: it will be shown subsequently that this greatly exceeds the mixing provided in existing models.

Even in this most weakly forced example where the aspect ratio of the convecting fluid remains greater than $`6`$ for most of the simulation, there persists an inverse cascade to large scale structures. Figure 17 shows a space time plot of the temperature perturbation on the bottom boundary. Hot spots, or thermals, along this boundary are seen as bright filaments in this plot; from here plumes rise and convection is driven. It is clear that the main plume dynamic is one of coalescence for, even though approximately $`10`$ thermals arise at the onset of convection, only two large thermals are apparent by the end. Furthermore, weak plumes are continually formed and incorporated into the stronger ones since the convection tends to direct flow toward strong plumes on the boundary. Two plumes implies four large scale rolls and since at $`t=0.136`$, $`h0.7`$ the aspect ratio of each roll is $`1.1`$, less than, but not atypical of large structure in regular Rayleigh-Benard convection. One must therefore conclude that the aspect ratio of layered convection will remain $`1`$ irrespective of the layer height.

Therefore, as to the the question of whether a second layer forms, the results for low Prandtl number are inconclusive but not promising. Since the large amplitude internal wave seen in $`IV`$ and $`V`$ drive instabilities along the interface and are a source for much cross interface mixing it remains possible that any layer can continually entrain fluid from above before oscillatory convection becomes well developed. Whether this is simply due to the strong forcing or to the relatively small aspect ratio under consideration remains to be clarified. If due to the latter, then it can be argued that in a stellar context where aspect ratios are large, such a large interfacial wave will not exist.

If, in fact, the fluid is forced too strongly in $`IV`$ and $`V`$ then there still remain several possibilities. It is clear that the convection continues to inverse cascade to larger scales at each stage, so that plumes within a layer are always spaced at a distance about twice the layer height. The plumes generate the interfacial waves reaching amplitudes of significant fractions of the layer thickness (particularly clear in figure 15) which may persist, constantly entraining fluid from above in broad down-flows. In such a circumstance an oscillatory instability would be unable to set in above the interface and eventually a state of global, classical convection would exist.

On the other hand, the interface may continue to migrate upward, gradually stabilizing until the plumes flatten along the buoyancy jump at the interface and are no longer able to excite large waves there. This mechanism occurs in the laboratory and is the motivation for introducing the Richardson number (equation 3). A low value of $`Ri`$ implies that the interface should migrate upward and Fernando (1989) measured a critical value for this behavior in the laboratory. The bottom panels in figures 3, $`3.5`$, 11 and 16 trace the interfacial Richardson number with time for simulations $`I`$, $`II`$, $`IV`$ and $`V`$, respectively. That $`Ri`$ is smallest in case $`IV`$ is consistent with the observed layer migration; in fact, in case $`IV`$ it is just within the range of the critical number measured in the laboratory experiments of Fernando (1989) where $`Ri_{crit}230`$.

Throughout simulation $`V`$, however, the Richardson number is greater than $`50`$, well above the critical value measured by Fernando (1989). However, the interface continually migrates upward and interfacial waves increase with time, as opposed to decreasing. Compare this with simulation $`I`$ (figure 3) where $`Ri<50`$ for a short period of time yet the interface does not migrate appreciably. Such observations can be most readily explained if the interfacial Richardson number at late times is below the critical value and thus provide strong evidence that the critical interfacial Richardson number should increase with decreasing Prandtl number.

## 4 Low Prandtl Number Convection Impinging on a Stable Compositional Interface.

In order to isolate the mixing experienced at a stable interface from the question of the formation of that interface, simulation $`VII`$ is initialized with $`Y=0`$ in the top half of the domain, $`Y=1`$ in the bottom and a thin ($`.05`$) transition region between.

By $`t=0.014`$, (figure 18) the convection has set in throughout the lower layer. Though the density jump still stabilizes the interface, waves are apparent. Again, mixing takes place at the troughs of these waves and in the helium plot there are two significant helium rich structures penetrating into the upper layer. Some interfacial scouring has occured, but is apparently much weaker than the wave breaking. Although $`ϵ_S>0`$ just above the interface, there is yet no motion in the upper layer beyond its interaction with the interfacial wave.

By $`t=.027`$ the interfacial wave breaking is very significant and a large mixing plume is seen in figure 19. The concentration jump has barely been eroded and the interfacial wave is still of moderate height ($`0.07`$) yet its self interaction is able to generate large structures, much like the interaction of two solitary waves on the surface of water.

In the upper half of the domain, not only is flow visible around the helium plume, but also ahead of it, resembling the upstream field of a weak vortex (notice the velocity in 19). This weak vortical flow, however, is unable to fully develop into convection before it becomes disrupted by the interfacial wave below. Clearly in figure 20 the flow in the upper half plane is soon dominated by the interfacial wave, which by $`t=0.042`$ has grown quite large.

Though this simulation deals with an idealized setting, it is obvious that the structures in the flow are less simple. In particular, the distinction between “helium mixing plume”, which results from the interaction of two gravity waves, and a large amplitude gravity wave is not precise. Gravity waves below an interface are also referred to as penetrative convection and have been studied numerically (see Massaguer et al. (1984), Hurlburt et al. (1994)) and analytically (Schmitt et al. (1984), in the context of a fixed stability interface). The current results are complicated by the fact that the stable interface itself interacts with the flow, yielding a highly nonlinear wave pattern and blurring previous distinctions of plume and wave. The structure in figure 19 is classified as mixing because, for the most part, it becomes entrained in the uppper fluid. The large interfacial wave in figure 20 is described as such because, although there is significant scouring of its upper boundary and entrainment at its crest, it mostly returns to the lower layer intact.

Using the Richardson number as a criterion for buoyancy stability, one concludes that the interface is stable throughout the simulation This is in stark contrast to the fact that interfacial waves are significant and at $`t=0.05`$, figure 21, are of amplitude equal to the layer height. As in the prior simulation, this discrepancy results from the fact that at late times the up-flow at the lower boundary concentrates into one dominant rising thermal. Recall that the Richardson number, plotted in figure 22 for this case, is the ratio of the buoyancy across an eddy height to the average kinetic energy in the lower layer: a mean field argument. Since the buoyancy and resulting momentum of the thermal plume is much greater than would be anticipated from a mean field argument it is able to penetrate much deeper than expected.

From figure 22 the conclusion should be drawn that, since the average kinetic energy in the lower domain is always $`>1/40`$ the potential energy wall of the interface, convection should not penetrate deeply. However, it is obvious from figure 21 that the flow is collimated in the strong thermal up-flows, thus the average kinetic energy is not a good representation of the dynamics in the lower layer.

Therefore it is also clear from this simulation that if the lower layer develops a large scale convective roll, then the interfacial wave will have a large scale and correspondingly large amplitude. Furthermore, figure 21 dramatically shows that in such instance, mixing is dominated by entrainment below the wave trough and by scouring of the interface from above through the weaker flow generated by the interfacial wave above it.

## 5 Discussion

### 5.1 Comparison with previous work

I remind the reader of the mixing theories of Stevenson (1979), Spruit (1992) and Fernando (1989), in order to compare them with these simulations. The former two provide mixing length theories of semiconvection in the presence of layers and, therefore, effective diffusion coefficients for transport across these layers. The latter derives the equilibrium height of the first layer through energy arguments and provides a criterion for onset of convection above the first layer.

First let us consider Fernando (1989). Using conservation of energy and molecular weight and with the assumption that the kinetic energy available to convection is proportional to the heat flux from the boundary, he derived the following relation

$$h_cC\left[\frac{Ra_{T,bot}}{Ra_{He}^{\frac{3}{4}}Pr^{\frac{1}{2}}}\right]^{\frac{1}{2}}$$
(17)
for $`h_c`$, the critical interface height of the first layer. The constant, $`C`$ must be determined experimentally as it contains the parameters of convective efficiency and its dependence on $`Pr`$ is not at all known. Calibrating the numerical experiments to the theory may not be reasonable, but the theory can be used to compare the simulation results amongst themselves. By calibrating the constant to simulation $`II`$, where $`h_c=0.8`$, equation (17) yields the critical layer heights $`h_c=0.50,\mathrm{\hspace{0.17em}\hspace{0.17em}1.4},\mathrm{\hspace{0.17em}\hspace{0.17em}0.63},\mathrm{\hspace{0.17em}\hspace{0.17em}1.0},\mathrm{\hspace{0.17em}\hspace{0.17em}0.55}`$ for simulations $`I,III(a),III(b),IV`$ and $`V`$ respectively. Though not quantitatively correct, these agree qualitatively with the observation of a layer in $`I`$ and $`V`$, plumes impinging on the upper boundary in $`III(a)`$ and interfacial waves hitting the upper boundary in $`IV`$. The lack of layering in $`III(b)`$ cannot, of course, be explained by this theory.

Of the necessary criteria suggested at the outset for the formation of a second layer, only simulation $`I`$ satisfied them both. Since it would be useful to refine them, we must now turn our consideration to the dynamics of the interface above the layer. Fernando (1989) reproduced the observation of Turner (1968) that the onset of a second layer requires the Rayleigh number in the interfacial region be sub-critical to overturning convection. If the interfacial thickness is denoted $`\delta `$ then this amounts to the requirement that $`\delta h_c`$. Again, using energetic arguments, it can be shown that

$$\frac{\delta }{h_c}=\frac{Ra_{crit}^{\frac{1}{4}}}{\sqrt{2}C}\left[\frac{Ra_TPr}{Ra_{He}}\right]^{\frac{1}{4}}$$
(18)
where the Rayleigh numbers are to be measured at the interface level and $`C`$ is the same as in equation (17). The critical Rayleigh number, $`Ra_{crit}`$ must also be determined experimentally and was found to be rather higher than that of classical convection, making the combination of constants of order one.

This can be reworked into necessary condition for the onset of oscillatory convection above the first layer and, amounts to

$$\mathrm{\Sigma }_{Fern}\frac{ϵ_S}{_\mu }\frac{4C^4}{Ra_{crit}Pr}.$$
(19)
Again experimental results suggest that the combination of constants on the right hand side are order one.

The theories of Spruit (1992) and Stevenson (1979) give effective diffusion coefficients for helium transport across an interface under the respective assumptions of enhanced diffusion or internal wave breaking. If the rate of diffusion is greater than the growth rate of oscillations above the interface, then the region above will slowly be entrained into the convection below. Otherwise, there is some hope for the growth of oscillations above. The growth rate of weakly forced overstable modes of wavenumber, $`k`$, is

$$\gamma =\frac{ϵ_S}{2|ϵ_L|}k^2$$
(20)
while the effective molecular diffusivities are

$$\tau _{wave}=\left[\frac{ϵ_S}{_\mu }\right]^2$$
(21)
and

$$\tau _{diff}=\sqrt{\tau }\frac{ϵ_S}{_\mu }$$
(22)
for Stevenson’s wave model and Spruit’s enhanced diffusion model, respectively. The rate of diffusion of modes $`k`$ is proportional $`k^2\tau _{eff}`$ and the requirement that oscillations grow more rapidly than a diffusion time yields the stability criteria

$$\mathrm{\Sigma }_{wave}\frac{ϵ_S}{_\mu }\frac{|ϵ_L|}{_\mu }1$$
(23)
and

$$\mathrm{\Sigma }_{diff}\sqrt{\tau }\frac{|ϵ_L|}{_\mu }1$$
(24)
respectively for the two theories.

The effective diffusion coefficients and stability parameter predicted by each theory can be measured by evaluating $`ϵ_S`$ and $`_\mu `$ at the interface. Figure 23 compares $`\tau _{wave}`$ and $`\tau _{diff}`$ with the fit $`h=1.8\sqrt{t}`$ to simulation $`V`$ which yields $`\tau _{eff}1.62`$. While the difference of both theories with the numerical data is astonishing, the enhanced diffusion model further suffers from the fact that it is proportional to $`\sqrt{\tau }`$, a negligibly small number in stars.

The variation of the stability criteria with time for simulation $`VII`$ is plotted in figure 24 and they all notably remain less than one, but not significantly so. Again, the criterion $`\mathrm{\Sigma }_{diff}`$ is proportional to $`\sqrt{\tau }(1/3`$ here) meaning that, in the limit of small molecular diffusivities, layered convection will always occur. This conclusion simply results from the assumption that only enhanced diffusive transport acts across the interfaces, inconsistent with the wave breaking seen in the simulations.

Both $`\mathrm{\Sigma }_{Fern}`$ and $`\mathrm{\Sigma }_{wave}`$ were derived considering buoyancy effects, yet the fact that they agree well throughout most of the simulation is still striking. This agreement results from the fact that the Ledoux criterion is dominated by the compositional gradient at the interface, $`|ϵ_L|_\mu `$. Nonetheless, it is clear that $`\tau _{wave}`$ underestimates the true layer entrainment and thus the stability criterion is overestimating the stability of the layer. Therefore, it seems unlikely that a second convecting region will form above the interface in $`VII`$ even though the Richardson number criteria indicates that it should be stable.

### 5.2 Conclusions and future work

These simulations have reproduced layer formation in the context of laboratory flows and distinguish the mechanisms for transport across a compositional interface: enhanced diffusion, scouring and wave breaking. They also confirm the difficulty of forming layers in a low Prandtl number fluid. Characteristic of all of the results is the inverse cascade of velocity to large length scales resulting in strong isolated upward rising thermals which are able to excite high amplitude, long wavelength waves on the interface. Mixing across the interface is dominated by wave breaking in all of the low Prandtl number cases; even, and most importantly, when the interface is prescribed and highly stable ($`VII`$).

Though it remains to be determined whether strong, isolated thermals develop in a three dimensional fluid, it can yet be safely concluded that, at low Prandtl number, cross interfacial mixing will be dominated by wave breaking. This is especially evident at early times in $`VII`$ where even small interfacial waves drive large “splashes”, mixing helium high into the stable layer. Therefore, the criteria associated with the Stevenson (1979) or Fernando (1989) theories of interfacial mixing seem to be the most appropriate for determining whether a second layer, and thus multiple layers, will form. Even so, these theories greatly underestimate the mixing across the interface.

This underestimate is linked to the the large interfacial waves and therefore the critical Richardson number for the stability of an interface. Comparing simulation $`V`$ to the low Prandtl number case, $`I`$, it is clear that the critical Richardson number must increase with decreasing Prandtl number.

Further study of whether the estimate of Stevenson (1979) truly underestimates mixing requires determining whether or not the large interfacial waves persist in 3-dimensional simulations and at large aspect ratio. Though the former would be most valuable, a parameter space study is not computationally feasible. The best prospects for sorting out the issue of isolated thermals rests in continuing this work at larger aspect ratios.

The simulations I have presented envision a scenario where an energy flux is abruptly turned on below a region which is marginally unstable to semiconvection. This results in the formation of layered convection due to the redistribution of a compositional gradient in a thin layer and on time scales short compared to a thermal diffusion time. The convection within such a layer is concentrated in plumes spaced such that the largest rolls have aspect ratio approximately $`1`$. The interface itself is characterized by large internal waves and smaller breaking waves, making it ineffective as a barrier to transport. Since, in all cases considered, the fluid is very well mixed after a small fraction of a thermal diffusion time, it is unlikely that layering caused by buoyancy disruption would play any role in stellar convection. In fact, the very large effective diffusion measured herein suggests that such layered structures are so unstable as to make them dynamically unimportant even when a compositional gradient is redistributed by other mechanisms such as rotation or large internal waves, which set up initial conditions similar to $`VII`$.

I gratefully acknowledge support and computing resources from the ASCI/Alliances Center for Astrophysical Thermonuclear Flashes at the University of Chicago under DOE subcontract B341495. Preliminary analytic studies were also supported by an NSF Graduate Fellowship and NASA Space Theory Grant NAG5-8495. This work was motivated by suggestions of Prof. Robert Rosner and benefited greatly from continuing conversations with him. Without Dr. Francesco Rubini, who provided the numerical expertise and code, this project would have taken an extra year. Without Dr. Diana Hall, this project would have taken six months less.