# The Rayleigh–Taylor instability in a self-gravitating two-layer viscous sphere

The Rayleigh–Taylor instability in a self-gravitating two-layer viscous sphere Summary The dispersion relation of the Rayleigh–Taylor instability in the spherical geometry is of profound importance in the context of the Earth's core formation. Here we present a complete derivation of this dispersion relation for a self-gravitating two-layer viscous sphere. Such relation is, however, obtained through the solution of a complex transcendental equation, and it is difficult to gain physical insights directly from the transcendental equation itself. We thus also derive an empirical formula to compute the growth rate, by combining the Monte Carlo sampling of the relevant model parameter space with linear regression. Our analysis indicates that the growth rate of Rayleigh–Taylor instability is most sensitive to the viscosity of inner layer in a physical setting that is most relevant to the core formation. Core, Numerical solutions, Planetary interiors 1 INTRODUCTION The Rayleigh–Taylor instability of a self-gravitating fluid sphere is important in the context of planet formation, especially in Earth's core formation (e.g. Ida et al. 1989; Stevenson 1990; Hoogenboom & Houseman 2006; Sasaki & Abe 2007). A dynamical situation typically envisioned for an early, growing Earth can be characterized by an inner layer of highly viscous silicate materials and an outer layer of denser molten metal (e.g. Honda et al. 1993; Rubie et al. 2015). This gravitationally unstable structure is susceptible to the Rayleigh–Taylor instability at the interface of these two layers. Such instability may lead to the dripping of the metallic layer into the silicate layer, which is an important part of core formation. However, the growth timescale of this Rayleigh–Taylor instability depends on several parameters such as viscosity, density and layer thicknesses, and an accurate dispersion relation becomes desirable. The dispersion relation for the Rayleigh–Taylor instability in a two-layer self-gravitating sphere was published previously by Ida et al. (1989) for both inviscid and viscous cases. This dispersion relation was utilized to discuss the Earth's core formation (Ida et al. 1987). Whereas the inviscid case is correctly presented, however, their derivation for the viscous case is incorrect. In this research note, we derive a correct dispersion relation for the Rayleigh–Taylor instability in a two-layer self-gravitating viscous sphere. As the growth rate depends on a number of parameters, it can be challenging to understand the variation of the growth rate in the high-dimensional parameter space. Also, the estimation of the growth rate in case of finite Prandtl number is based on the solution of a non-linear transcendental equation, which can be very cumbersome in the spherical geometry. Therefore, we develop a practical approach to derive an approximate formula for the growth rate by combining Monte Carlo sampling and linear regression for l = 2, the mode most relevant to the Earth's core formation. The l = 1 mode is the translational mode and results simply in an overall shift of the centre of mass of the whole spherical system. Despite being the fast growing mode for the case of μ1 ≫ μ2, therefore, it is not relevant to the core formation (e.g. Elsasser 1963; Honda et al. 1993). Modes with l > 1 are important for core formation, but with increasing l, the viscous dissipation becomes higher in proportion to ≈l(l + 1), thereby making growth timescale longer. The l = 2 mode thus has the maximum contribution, and we restrict ourselves to this particular case when deriving an approximate formula. The structure of this note is following. We begin with the theoretical formulation of the Rayleigh–Taylor instability in a self-gravitating two-layer viscous sphere. We then present a few representative dispersion curves for different combinations of layer viscosities. We next discuss the practical approach to obtain an approximate formula for growth rate. Our result suggests that when the viscosity of the inner layer is much higher than that of the outer layer, the growth rate is controlled mostly by the properties of the inner layer as well as the density contrast of the two layers. 2 DERIVATION OF THE DISPERSION RELATION The governing equations for an incompressible viscous fluid include the equation of continuity   $$\nabla \cdot {\bf U}=0,$$ (1)the equation of motion   $$\frac{{\rm d}\mathbf {U}}{{\rm d}t}=\nabla \cdot \mathbf {\Sigma }+{\bf F},$$ (2)and the constitutive relation   $$\mathbf {\Sigma }=-P\mathbf {I}+\mu \left[\nabla \mathbf {U}+(\nabla \mathbf {U})^{T}\right],$$ (3)where U is the total velocity field, Σ is the total stress, F is the body force, P is the total pressure, I is the identity operator, and μ is the dynamic viscosity. To focus on the linear Rayleigh–Taylor instability, we express the flow-related entities as first-order perturbations from the hydrostatic reference:   \begin{eqnarray} {\bf U}&=&{\bf U}_{0}+{\bf u}, \end{eqnarray} (4)  \begin{eqnarray} \mathbf {\Sigma }&=&\mathbf {\Sigma }_{0}+\mathbf {\sigma }, \end{eqnarray} (5)  \begin{eqnarray} P&=&P_{0}+\delta p, \end{eqnarray} (6)  \begin{eqnarray} {\bf F}&=&{\bf F}_{0}+{\bf f}, \end{eqnarray} (7)where [U0, Σ0, P0, F0]T and [u, σ, δp, f]T represent the hydrostatic reference state (U0 is trivially zero) and the perturbed state, respectively. In the hydrostatic state, we may write the body force (gravity) as a gradient of potential that satisfies the Poisson equation:   \begin{eqnarray} \mathbf {F}_{0}&=&-\rho _{0}\nabla \Psi , \end{eqnarray} (8)  \begin{eqnarray} \nabla ^{2} \Psi &=& -4\pi G \rho _{0}, \end{eqnarray} (9)where ρ0 and G are the density of the reference state and the gravitational constant, respectively. We may solve the Poisson equation for the hydrostatic state (eq. 9) and obtain the unperturbed gravitational potential as   \begin{eqnarray} \Psi &=& \frac{4\pi G}{3}\left(\rho _{1}(r^{2}/2-3R_{1}^{2}/2)+3\rho _{2}(R_{1}^{2}-R_{2}^{2})/2\right), \quad \quad \forall r<R_{1}, \end{eqnarray} (10)  \begin{eqnarray} \quad&=&\frac{4\pi G}{3}\left(\rho _{1}(r^{2}/2-3R_{2}^{2}/2)+(\rho _{2}-\rho _{1})R_{1}^{3}/r\right), \quad \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (11)where ρ1, ρ2, R1 and R2 are the density of the inner layer, the density of the outer layer, the radius of the inner sphere, and the radius of the outer sphere, respectively (Fig. 1). The hydrostatic state confirms the equality of the pressure gradient and the gravitational force:   $$\nabla P_{0}=\rho _{0}\nabla \Psi ,$$ (12)which will be used in boundary conditions. Subtracting the hydrostatic reference state and neglecting higher order terms lead to   \begin{eqnarray} \nabla _{j}u_{j}&=&0, \end{eqnarray} (13)  \begin{eqnarray} \nabla _{j}\sigma _{ij}+f_{i}&=&\frac{\partial u_{i}}{\partial t}, \end{eqnarray} (14)  \begin{eqnarray} \sigma _{ij}&=&-\delta p\delta _{ij}+\mu \left(\nabla _{j}u_{i}+\nabla _{i}u_{j}\right), \end{eqnarray} (15)where fi arises solely due to density perturbation at the interfaces of two adjacent fluids, and it can also be represented as the gradient of a potential perturbation δΨ satisfying the Poisson equation:   \begin{eqnarray} f_{i}&=&-\rho _{0}\nabla \delta \Psi , \end{eqnarray} (16)  \begin{eqnarray} \nabla ^{2} \delta \Psi &=&-4\pi G \delta \rho , \end{eqnarray} (17)where δρ is the perturbation to the reference density. Together with the perturbed eqs (13)–(17), we may write the equation of motion as   $$\frac{\partial \mathbf {u}}{\partial t}=-\nabla \Phi +\nu \nabla ^{2} \mathbf {u},$$ (18)where   $$\Phi =\frac{\delta P}{\rho }+\delta \Psi .$$ (19)The system of eqs (13)–(19) describing the dynamics of the perturbed state is used to calculate the growth rate of initial density perturbations. In the 3-D spherical geometry, we expand the field variables in the spherical harmonics assuming the horizontal invariance of material properties (e.g. Chandrasekhar 1961)   \begin{eqnarray} \Xi (r,\theta ,\phi ) &=& \sum _{l,m}\Xi _{lm}(r)Y^{m}_{l}(\theta ,\phi ). \end{eqnarray} (20)Using this separation of variables, we solve for the perturbed gravitational potential generated through the deformation of the interface (assuming the outer surface of the sphere to be rigid):   \begin{eqnarray} \delta \Psi (r)&=&\left(\frac{4\pi G(\rho _{2}-\rho _{1})}{(2l+1)R_{1}^{l-1}}\delta r_{1}\right) r^{l}, \quad \quad \forall r < R_{1}, \end{eqnarray} (21)  \begin{eqnarray} \qquad&=&\left(\frac{4\pi G (\rho _{2}-\rho _{1})}{2l+1}R_{1}^{l+2}\delta r_{1}\right)\bigg/r^{l+1}, \quad \quad \forall R_{1}<r<R_{2}. \end{eqnarray} (22)We can also readily solve for the total perturbed potential Φ by taking divergence of the eq. (18) and using the continuity eq. (13):   $$\nabla ^{2} \Phi =0,$$ (23)the solution of which may further be used to solve the equation of motion (18). However, before proceeding to solve the equation of motion, it is convenient to perform poloidal–toroidal decomposition in the spherical geometry (e.g. Chandrasekhar 1961) as we will only require the poloidal part for linear stability analysis. Using such decomposition, we may express the perturbed velocity field as   $$\mathbf {u}=\sum _{l,m}\left(\nabla \times (T_{lm}(r)Y^{m}_{l}\hat{\mathbf {r}})+\nabla \times \left[\nabla \times (S_{lm}(r)Y^{m}_{l}\hat{\mathbf {r}})\right]\right),$$ (24)where Slm and Tlm denote the poloidal and toroidal coefficients, respectively. Assuming u∝exp (σt) in the linear regime with σ being the growth rate, the poloidal part of the equation of motion (eq. 18) may be written as the following second-order differential equation:   $$\sigma S_{lm}=-\nabla \Phi +\nu \left[\frac{l(l+1)}{r^{2}}-D^{2}\right]S_{lm},$$ (25)where ν and D denote the kinematic viscosity and the radial derivative, that is, d/dr, respectively. The complementary function and particular solution of the eq. (25) can be obtained for the inner and outer layers separately as linear combinations of modified Bessel functions:   \begin{eqnarray} S_{lm}(r)&=&A_{1}\sqrt{(}r)I_{l+1/2}(q_{1}r)-\frac{Q_{1}R_{1}^{2}}{\sigma }\left(\frac{r}{R_{1}}\right)^{l+1}, \quad \quad \forall r<R_{1},\nonumber \\ &=&A_{2}\sqrt{(}r)I_{l+1/2}(q_{2}r)+B_{2}\sqrt{r}K_{l+1/2}(q_{2}r)-\left[Q_{2}\left(\frac{r}{R_{1}}\right)^{l+1}+Q_{3}\left(\frac{r}{R_{1}}\right)^{-l}\right]R_{1}^{2}/\sigma , \quad \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (26)where q1 = σ/ν1, q2 = σ/ν2, ν1 and ν2 are the viscosities of the inner and outer layers, respectively, and [A1, Q1, A2, B2, Q2, Q3]T is a set of unknown constants to be obtained using the boundary conditions. The velocity field associated with the Rayleigh–Taylor instability may be obtained from the poloidal potential Slm. The radial (ur) and transverse (uθ) components are obtained as follows:   \begin{eqnarray} u_{r}&=&\sum _{l,m}\frac{l(l+1)}{r^{2}}S_{lm}(r)Y_{l}^{m}(\theta ,\phi ),\nonumber \\ &=&\sum _{l,m}\left[\frac{A_{1}l(l+1)}{r^{3/2}}I_{l+1/2}(q_{1}r)-\frac{Q_{1}l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}\right]Y_{l}^{m}(\theta ,\phi ), \quad \forall r<R_{1}, \end{eqnarray} (27)  \begin{eqnarray} &&{\quad=\sum _{l,m}\Bigg[\frac{A_{2}l(l+1)}{r^{3/2}}I_{l+1/2}(q_{2}r)+\frac{B_{2}l(l+1)}{r^{3/2}}K_{l+1/2}(q_{2}r)}\nonumber \\ &&{\quad\quad-\frac{l(l+1)}{r^{2}}\left[Q_{2}\left(\frac{r}{R_{1}}\right)^{l+1}+Q_{3}\left(\frac{r}{R_{1}}\right)^{-l}\right]R_{1}^{2}/\sigma \Bigg]Y_{l}^{m}(\theta ,\phi ), \quad \forall R_{1}<r<R_{2},} \end{eqnarray} (28)  \begin{eqnarray} u_{\theta }&=&\sum _{l,m}\frac{1}{r}\frac{{\rm d}S_{lm}(r)}{{\rm d}r}\frac{\partial Y_{l}^{m}(\theta ,\phi )}{\partial \theta },\nonumber \\ &=&\sum _{l,m}\Big[A_{1}\left((l+1)r^{-3/2}I_{l+1/2}(q_{1}r)+q_{1}r^{-1/2}I_{l+3/2}(q_{1}r)\right)-\frac{Q_{1}(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}\Big]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall r<R_{1}, \end{eqnarray} (29)  \begin{eqnarray} &&{\quad=\sum _{l,m}\Bigg[A_{2}\left((l+1)r^{-3/2}I_{l+1/2}(q_{2}r)+q_{2}r^{-1/2}I_{l+3/2}(q_{2}r)\right)+B_{2}\left((l+1)r^{-3/2}K_{l+1/2}(q_{2}r)-q_{2}r^{-1/2}K_{l+3/2}(q_{2}r)\right)}\nonumber \\ &&{\quad\quad-\frac{1}{\sigma }\left(Q_{2}(l+1)\left(\frac{r}{R_{1}}\right)^{l-1}-Q_{3}l \left(\frac{r}{R_{2}}\right)^{-(l+2)}\right)\Bigg]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall R_{1}<r<R_{2}.} \end{eqnarray} (30)The stresses may be calculated from the velocities using the following constitutive relations:   \begin{eqnarray} \sigma _{rr}&=&2\mu \frac{\partial u_{r}}{\partial r}, \end{eqnarray} (31)  \begin{eqnarray} \sigma _{r\theta }&=&\mu \left[\frac{1}{r}\frac{\partial u_{r}}{\partial \theta }-\frac{u_{\theta }}{r}+\frac{\partial u_{\theta }}{\partial r}\right]. \end{eqnarray} (32) Figure 1. View largeDownload slide The spherical system considered in the Rayleigh–Taylor instability analysis. Unstable stratification corresponds to ρ2 > ρ1. Figure 1. View largeDownload slide The spherical system considered in the Rayleigh–Taylor instability analysis. Unstable stratification corresponds to ρ2 > ρ1. It is essential to compute the stresses correctly, particularly the normal stress, as it is used in the discontinuity condition, which is directly related to the expression of the growth rate. The previously published result by Ida et al. (1989) uses an incorrect expression for the normal and tangential stresses (their eqs A-21 and A-22); a simple dimensional analysis reveals incompatibility between the both sides of these equations. To make calculations simpler, we may express the velocities and stresses in the following manner:   \begin{eqnarray} u_{r}&=&\sum _{l,m}\left[a_{1}A_{1}+b_{1}Q_{1}\right]Y_l^m(\theta ,\phi ), \quad \forall r<R_{1}, \end{eqnarray} (33)  \begin{eqnarray} \quad&=&\sum _{l,m}\left[a_{2}A_{2}+b_{2}B_{2}+c_{2}Q_{2}+d_{2}D_{2}\right]Y_l^m(\theta ,\phi ), \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (34)  \begin{eqnarray} u_{\theta }&=&\sum _{l,m}\left[e_{1}A_{1}+f_{1}Q_{1}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall r<R_{1}, \end{eqnarray} (35)  \begin{eqnarray} \quad=\sum _{l,m}\left[e_{2}A_{2}+f_{2}B_{2}+g_{2}Q_{2}+h_{2}Q_{3}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (36)  \begin{eqnarray} \tau _{rr}&=&-P+2\mu _{1}\sum _{l,m}\left[i_{1}A_{1}+j_{1}Q_{1}\right]Y_l^m(\theta ,\phi ), \quad \forall r<R_{1}, \end{eqnarray} (37)  \begin{eqnarray} \quad&=&-P+2\mu _{2}\sum _{l,m}\left[i_{2}A_{2}+j_{2}B_{2}+k_{2}Q_{2}+l_{2}Q_{3}\right]Y_l^m(\theta ,\phi ), \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (38)  \begin{eqnarray} \tau _{r\theta }&=&\mu _{1}\sum _{l,m}\left[\alpha _{1}A_{1}+\beta _{1}Q_{1}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall r<R_{1}, \end{eqnarray} (39)  \begin{eqnarray} \quad\quad&=&\mu _{2}\sum _{l,m}\left[\alpha _{2}A_{2}+\beta _{2}B_{2}+\gamma _{2}Q_{2}+\delta _{2}Q_{3}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (40)where   \begin{eqnarray} a_{1}(r)&=&\frac{l(l+1)}{r^{3/2}}I_{l+1/2}(q_{1}r), \end{eqnarray} (41)  \begin{eqnarray} b_{1}(r)&=&-\frac{l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}, \end{eqnarray} (42)  \begin{eqnarray} a_{2}(r)&=&\frac{l(l+1)}{r^{3/2}}I_{l+1/2}(q_{2}r), \end{eqnarray} (43)  \begin{eqnarray} b_{2}(r)&=&\frac{l(l+1)}{r^{3/2}}K_{l+1/2}(q_{2}r), \end{eqnarray} (44)  \begin{eqnarray} c_{2}(r)&=&-\frac{l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}, \end{eqnarray} (45)  \begin{eqnarray} d_{2}(r)&=&-\frac{l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{-(l+2)}, \end{eqnarray} (46)  \begin{eqnarray} e_{1}(r)&=&\frac{(l+1)}{r^{3/2}}I_{l+1/2}(q_{1}r)+\frac{q_{1}}{\sqrt{r}}I_{l+3/2}(q_{1}r), \end{eqnarray} (47)  \begin{eqnarray} f_{1}(r)&=&-\frac{(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}, \end{eqnarray} (48)  \begin{eqnarray} e_{2}(r)&=&(l+1)r^{-3/2}I_{l+1/2}(q_{2}r)+q_{2}r^{-1/2}I_{l+3/2}(q_{2}r), \end{eqnarray} (49)  \begin{eqnarray} f_{2}(r)&=&(l+1)r^{-3/2}K_{l+1/2}(q_{2}r)-q_{2}r^{-1/2}K_{l+3/2}(q_{2}r), \end{eqnarray} (50)  \begin{eqnarray} g_{2}(r)&=&-\frac{(l+1)}{\sigma }(r/R_{1})^{l-1}, \end{eqnarray} (51)  \begin{eqnarray} h_{2}(r)&=&\frac{l}{\sigma }\left(\frac{r}{R_{2}}\right)^{-(l+2)}, \end{eqnarray} (52)  \begin{eqnarray} i_{1}(r)&=&l(l+1)\left[(l-1)r^{-5/2}I_{l+1/2}(q_{1}r)+q_{1}r^{-3/2}I_{l+3/2}(q_{1}r)\right], \end{eqnarray} (53)  \begin{eqnarray} j_{1}(r)&=&-\frac{l(l+1)(l-1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (54)  \begin{eqnarray} i_{2}(r)&=&l(l+1)\left[(l-1)r^{-5/2}I_{l+1/2}(q_{2}r)+q_{2}r^{-3/2}I_{l+3/2}(q_{2}r)\right], \end{eqnarray} (55)  \begin{eqnarray} j_{2}(r)&=&l(l+1)\left[(l-1)r^{-5/2}K_{l+1/2}(q_{2}r)-q_{2}r^{-3/2}K_{l+3/2}(q_{2}r)\right], \end{eqnarray} (56)  \begin{eqnarray} k_{2}(r)&=&-\frac{l(l+1)(l-1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (57)  \begin{eqnarray} l_{2}(r)&=&\frac{l(l+1)(l+2)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{-(l+3)}, \end{eqnarray} (58)  \begin{eqnarray} \alpha _{1}(r)&=&2(l+1)(l-1)r^{-5/2}I_{l+1/2}(q_{1}r)+(2l+1)q_{1}r^{-3/2}I_{l+3/2}(q_{1}r)+q_{1}^{2}r^{-1/2}I_{l+5/2}(q_{1}r), \end{eqnarray} (59)  \begin{eqnarray} \beta _{1}(r)&=&-\frac{2(l-1)(l+1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (60)  \begin{eqnarray} \alpha _{2}(r)&=&2(l+1)(l-1)r^{-5/2}I_{l+1/2}(q_{2}r)+(2l+1)q_{2}r^{-3/2}I_{l+3/2}(q_{2}r)+q_{2}^{2}r^{-1/2}I_{l+5/2}(q_{2}r), \end{eqnarray} (61)  \begin{eqnarray} \beta _{2}(r)&=&2(l+1)(l-1)r^{-5/2}K_{l+1/2}(q_{2}r)-(2l+1)q_{2}r^{-3/2}K_{l+3/2}(q_{2}r)+q_{2}^{2}r^{-1/2}K_{l+5/2}(q_{2}r), \end{eqnarray} (62)  \begin{eqnarray} \gamma _{2}(r)&=&-\frac{2(l-1)(l+1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (63)  \begin{eqnarray} \delta _{2}(r)&=&-\frac{2l(l+2)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{-(l+3)}. \end{eqnarray} (64) We are now ready to use boundary conditions to obtain the unknown coefficients. We assume the outer boundary to be rigid, that is, ur = uθ = 0 at r = R2. The inner boundary is assumed to be welded and hence the field variables are continuous across the deformed boundary. The first-order Taylor expansion leads to the continuity of the velocity and transverse stress across the r = R1 boundary, that is, $$\left[u_{r}\right]_{R_{1}-}^{R_{1}+}=0$$, $$\left[u_{\theta }\right]_{R_{1}-}^{R_{1}+}=0$$, and $$\left[\sigma _{r\theta }\right]_{R_{1}-}^{R_{1}+}=0$$, while the radial stress is continuous only across the deformed boundary, $$\left[-(P_{0}+\delta P)+2\mu \partial _{r}u_{r}\right]_{R_{1}-\delta r}^{R_{1}+\delta r}=0$$, and becomes discontinuous across r = R1 due to the density jump. P0 and δP may be calculated from eqs (12), (16), and (23). The rigid and continuous boundary conditions together with the orthonormality of the spherical harmonics form the following complete set of linear equations of the unknown coefficients for each (l, m) pair:   \begin{eqnarray} a_{2}(R_{2})A_{2}+b_{2}(R_{2})B_{2}+c_{2}(R_{2})Q_{2}+d_{2}(R_{2})Q_{3}&=&0, \end{eqnarray} (65)  \begin{eqnarray} e_{2}(R_{2})A_{2}+f_{2}(R_{2})B_{2}+g_{2}(R_{2})Q_{2}+h_{2}(R_{2})Q_{3}&=&0, \end{eqnarray} (66)  \begin{eqnarray} a_{1}(R_{1})A_{1}+b_{1}(R_{1})Q_{1}&=&a_{2}(R_{1})A_{2}+b_{2}(R_{1})B_{2}+c_{2}(R_{1})Q_{2}+d_{2}(R_{1})Q_{3}, \end{eqnarray} (67)  \begin{eqnarray} e_{1}(R_{1})A_{1}+f_{1}(R_{1})Q_{1}&=&e_{2}(R_{1})A_{2}+f_{2}(R_{1})B_{2}+g_{2}(R_{1})Q_{2}+h_{2}(R_{1})Q_{3}, \end{eqnarray} (68)  \begin{eqnarray} \mu _{1}\left[\alpha _{1}(R_{1})A_{1}+\beta _{1}(R_{1})Q_{1}\right]&=&\mu _{2}[\alpha _{2}(R_{1})A_{2}+\beta _{2}(R_{1})B_{2}+\gamma _{2}(R_{1})Q_{2}+\delta _{2}(R_{1})Q_{3}], \end{eqnarray} (69)  \begin{eqnarray} \frac{4\pi G(\rho _{1}-\rho _{2})[2(l-1)\rho _{1}+3\rho _{2}]}{3(2l+1)}\delta r&=&(l+1)(\rho _{1}Q_{1}-\rho _{2}Q_{2})+l\rho _{2}Q_{3}+2\mu _{2}/R_{1}\partial _{r}u_{r2}-2\mu _{1}/R_{1}\partial _{r}u_{r1}. \end{eqnarray} (70)We may eliminate the several unknowns from these equations and obtain the following form:   \begin{eqnarray} A_{2}&=&\xi _{1}Q_{2}+\eta _{1}Q_{3}, \end{eqnarray} (71)  \begin{eqnarray} B_{2}&=&\xi _{2}Q_{2}+\eta _{2}Q_{3}, \end{eqnarray} (72)  \begin{eqnarray} A_{1}&=&X_{1}Q_{2}+Y_{1}Q_{3}, \end{eqnarray} (73)  \begin{eqnarray} Q_{1}&=&X_{2}Q_{2}+Y_{2}Q_{3}, \end{eqnarray} (74)  $$Q_{2}=\Theta Q_{3},$$ (75)  $$Q_{3}=\Sigma \delta r,$$ (76)where   \begin{eqnarray} \xi _{1}&=&-\frac{f_{2}(R_{2})c_{2}(R_{2})-b_{2}(R_{2})g_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (77)  \begin{eqnarray} \eta _{1}&=&-\frac{f_{2}(R_{2})d_{2}(R_{2})-b_{2}(R_{2})h_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (78)  \begin{eqnarray} \xi _{2}&=&-\frac{a_{2}(R_{2})g_{2}(R_{2})-e_{2}(R_{2})c_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (79)  \begin{eqnarray} \eta _{2}&=&-\frac{a_{2}(R_{2})h_{2}(R_{2})-e_{2}(R_{2})d_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (80)  \begin{eqnarray} X_{1}&=&\frac{f_{1}(R_{1})[a_{2}(R_{1})\xi _{1}+b_{2}(R_{1})\xi _{2}+c_{2}(R_{1})]-b_{1}(R_{1})[e_{2}(R_{1})\xi _{1}+f_{2}(R_{1})\xi _{2}+g_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (81)  \begin{eqnarray} Y_{1}&=&\frac{f_{1}(R_{1})[a_{2}(R_{1})\eta _{1}+b_{2}(R_{1})\eta _{2}+d_{2}(R_{1})]-b_{1}(R_{1})[e_{2}(R_{1})\eta _{1}+f_{2}(R_{1})\eta _{2}+h_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (82)  \begin{eqnarray} X_{2}&=&\frac{a_{1}(R_{1})[e_{2}(R_{1})\xi _{1}+f_{2}(R_{1})\xi _{2}+g_{2}(R_{1})]-e_{1}(R_{1})[a_{2}(R_{1})\xi _{1}+b_{2}(R_{1})\xi _{2}+c_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (83)  \begin{eqnarray} Y_{2}&=&\frac{a_{1}(R_{1})[e_{2}(R_{1})\eta _{1}+f_{2}(R_{1})\eta _{2}+h_{2}(R_{1})]-e_{1}(R_{1})[a_{2}(R_{1})\eta _{1}+b_{2}(R_{1})\eta _{2}+d_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (84)  $$\Theta =\frac{\mu _{2}[\alpha _{2}(R_{1})\eta _{1}+\beta _{2}(R_{1})\eta _{2}+\delta _{2}(R_{1})]-\mu _{1}[\alpha _{1}(R_{1})Y_{1}+\beta _{1}(R_{1})Y_{2}]}{\mu _{1}[\alpha _{1}(R_{1})X_{1}+\beta _{1}(R_{1})X_{2}]-\mu _{2}[\alpha _{2}(R_{1})\xi _{1}+\beta _{2}(R_{1})\xi _{2}+\gamma _{2}(R_{1})]},$$ (85)  $$\Sigma =\frac{4\pi G (\rho _{1}-\rho _{2})(2(l-1)\rho _{1}+3\rho _{2})}{3(2l+1)\Pi },$$ (86)and   \begin{eqnarray} \Pi &=&\left[\left((l+1)\rho _{1}X_{2}-(l+1)\rho _{2}\right)+\frac{2\mu _{2}}{R_{1}}\left(i_{2}(R_{1})\xi _{1}+j_{2}(R_{1})\xi _{2}+k_{2}(R_{1})\right)-\frac{2\mu _{1}}{R_{1}}\left(i_{1}(R_{1})X_{1}+j_{1}(R_{1})X_{2}\right)\right]\Theta \nonumber\\ &&+\left[(l+1)\rho _{1}Y_{2}+l\rho _{2}+\frac{2\mu _{2}}{R_{1}}\left(i_{2}(R_{1})\eta _{1}+j_{2}(R_{1})\eta _{2}+l_{2}(R_{1})\right)-\frac{2\mu _{1}}{R_{1}}\left(i_{1}(R_{1})Y_{1}+j_{1}(R_{1})Y_{2}\right)\right]. \end{eqnarray} (87)We may write the relation between the amplitude of perturbation and the vertical velocity at the interface as   \begin{eqnarray} \frac{{\rm d}(\delta r)}{{\rm d}t}&=&u_{r}(R_{1}). \end{eqnarray} (88)By noting that δr ∝ exp (σt) in the linear regime, eqs (34) and (88) together with the six unknowns (i.e. A1-Q3) may be combined to obtain the expression for the growth rate   $$\sigma =\left(a_{1}(R_{1})[X_{1}\Theta +Y_{1}]+b_{1}(R_{1})[X_{2}\Theta +Y_{2}]\right)\Sigma .$$ (89)This is a transcendental equation for σ, because its right-hand side depends on σ in a complex manner. 3 EXAMPLES AND APPROXIMATE FORMULA In this section, we provide examples of dispersion relation for several combinations of layer viscosities. It is convenient to present the growth rate in the dimensionless form, so we scale the growth rate with respect to the free fall rate as follows:   $$\sigma^{\prime}=\sigma /\sqrt{4\pi G (\rho _{2}-\rho _{1})/3}.$$ (90)First, we fix the viscosity of the inner layer (i.e. μ1) and compute the growth rate for a range of outer layer viscosity (i.e. μ2). Similarly, we compute the growth rates for a range of inner layer viscosity keeping the viscosity of the outer layer constant. In both of these cases, the primary contributing factor turns out to be max (μ1, μ2) (Fig. 2). For μ1 ≫ μ2, the translational mode approximately corresponds to the motion of a rigid sphere through a viscous fluid. The corresponding growth rate turns out to be the fastest; the inner rigid sphere almost instantaneously translates to a new position leading to an overall shift of the centre of mass. Viscosity variation therefore does not affect considerably the growth rate for this translational mode as long as μ1 ≫ μ2 (Figs 2a and b). The growth rates of modes with l ≥ 2 exhibit nearly linear dependence on the inner layer viscosity (Fig. 2b). In the case of μ2 ≫ μ1, little variation in growth rate is observed with μ1, while growth rates associated with modes with small angular order (i.e. l < 5) vary linearly with μ2 (Figs 2c and d). However, the case of μ2 ≫ μ1 is not relevant in the context of core formation. Figure 2. View largeDownload slide The dispersion relations for several combinations of inner and outer layer viscosities: (a) the case where μ2 is varied keeping μ1 constant and μ1 ≫ μ2, (b) the case where μ1 is varied keeping μ2 constant and μ1 ≫ μ2—most relevant to the core formation, (c) the case where μ1 is varied keeping μ2 constant and μ2 ≫ μ1, (d) the case where μ2 is varied keeping μ1 constant and μ2 ≫ μ1. Other parameters are fixed as follows: R1 = 106 m, dR = (R2 − R1) = 104.84 m, ρ1 = 103.74 kg m−3 and δρ = 103.5 kg m−3. Figure 2. View largeDownload slide The dispersion relations for several combinations of inner and outer layer viscosities: (a) the case where μ2 is varied keeping μ1 constant and μ1 ≫ μ2, (b) the case where μ1 is varied keeping μ2 constant and μ1 ≫ μ2—most relevant to the core formation, (c) the case where μ1 is varied keeping μ2 constant and μ2 ≫ μ1, (d) the case where μ2 is varied keeping μ1 constant and μ2 ≫ μ1. Other parameters are fixed as follows: R1 = 106 m, dR = (R2 − R1) = 104.84 m, ρ1 = 103.74 kg m−3 and δρ = 103.5 kg m−3. In order to study the Rayleigh–Taylor instability associated with core formation in a comprehensive manner, it would be convenient to derive an approximate formula for the growth rate for the most important mode associated with the deformation of the interface. To derive such a formula, we explore the following parameter ranges: R1 = 105 to 107 m, ΔR = (R2 − R1) = 103 to 5 × 106 m, ρ1 = 3 × 103 to 104 kg m−3, δρ = (ρ2 − ρ1) = 103 to 104 kg m−3, μ1 = 1015 to 1025 Pa s, and μ2 = 106 to 1010 Pa s. First, we consider a reference case with R1 = 106 m, ΔR = 104.84 m, ρ1 = 103.74 kg m−3, δρ = 103.5 kg m−3, μ1 = 1020 Pa s, and μ2 = 108 Pa s. The model space to be explored can be defined around this reference case with the following standard deviation:   $$\zeta (\mathbf {x})=[0.5, 0.5,0.5,1,3,1.5]^{T},$$ (91)where   $$\mathbf {x}=[\log _{10}(R_{1}),\log _{10}(\rho _{1}),\log _{10}(\delta \rho ),\log _{10}(\Delta R),\log _{10}(\mu _{1}),\log _{10}(\mu _{2})]^{T}.$$ (92)This model space is subjected to random sampling using the normal distribution, and a model vector (eq. 92) is constructed for each sample. This obtained set of random model vectors is then used to compute the growth rate corresponding to the l = 2 mode. An ensemble of growth rates is obtained through 200 such calculations. Despite the complexity of the transcendental equation (eq. 89), the growth rate is found to be predicted with reasonable accuracy using the following function:   $$\log _{10}(\sigma )=a_{0}+\sum _{j=1}^{6}a_{j}X_{j},$$ (93)where Xj is the mean subtracted input parameter that is normalized by the corresponding standard deviation (eq. 91). Constants a0 through a6 are estimated by least square regression: a0 = −10.43, a1 = 1.86, a2 = 0.87, a3 = 1.3, a4 = 0.047, a5 = −0.9752, and a6 = −0.054. Fig. 3 demonstrates a reasonable correspondence between the predicted and actual values of growth rates for the aforementioned 200 computations. This way of summarizing enables us not only to see the sensitivity of the growth rate on the model parameters, but also to quickly compute the growth rate without doing the whole calculation. We suggest the following expression for the growth rate corresponding to the l = 2 mode:   $$\sigma =10^{-10.13}R_{1}^{1.86}(\rho _{2}-\rho _{1})^{1.3}\rho _{1}^{0.87}(R_{2}-R_{1})^{0.047}\mu _{1}^{-0.98}\mu _{2}^{-0.054},$$ (94)where R, ρ, μ and σ are in m, kg m−3, Pa s and s−1, respectively. Figure 3. View largeDownload slide Correspondence between the actual and predicted values of growth rate for l = 2. Figure 3. View largeDownload slide Correspondence between the actual and predicted values of growth rate for l = 2. ACKNOWLEDGEMENTS This study is based on work supported by the NASA through the NASA Astrobiology Institute under cooperative agreement no. NNA15BB03A issued through the Science Mission Directorate. We thank Gael Choblet and Neil Ribe for their careful reviews and constructive comments. REFERENCES Chandrasekhar S., 1961. Hydrodynamic and Hydromagnetic Stability  Clarendon. Elsasser W., 1963. Early history of the Earth, in Earth Science and Meteoritics  pp. 1– 30, eds Geiss J., Goldberg E., North-Holland. Honda R., Mizutani H., Yamamoto T., 1993. Numerical simulation of Earth's core formation, J. geophys. Res.  98 2075– 2089. Google Scholar CrossRef Search ADS   Hoogenboom T., Houseman G.A., 2006. Rayleigh–Taylor instability as a mechanism for corona formation on Venus, Icarus  180 292– 307. Google Scholar CrossRef Search ADS   Ida S., Nakagawa Y., Nakazawa K., 1987. The Earth's core formation due to the Rayleigh-Taylor instability, Icarus  69 239– 248. Google Scholar CrossRef Search ADS   Ida S., Nakagawa Y., Nakazawa K., 1989. The Rayleigh–Taylor instability in a self-gravitating two-layer fluid sphere, Earth Moon Planets  44 149– 174. Google Scholar CrossRef Search ADS   Rubie D., Nimmo F., Melosh H., 2015. Formation of the Earth's core, in Treatise On Geophysics  2nd edn, vol. 9, pp. 51– 90, ed. Stevenson D., Elsevier. Sasaki T., Abe Y., 2007. Rayleigh–Taylor instability after giant impacts: imperfect equilibration of the Hf-W system and its effect on the core formation age, Earth Planets Space  59 1035– 1045. Google Scholar CrossRef Search ADS   Stevenson D., 1990. Fluid dynamics of core formation, in Origin of the Earth  pp. 231– 249, eds Newsom H., Jones J.E., Oxford Univ. Press. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# The Rayleigh–Taylor instability in a self-gravitating two-layer viscous sphere

, Volume 212 (3) – Mar 1, 2018
9 pages

/lp/ou_press/the-rayleigh-taylor-instability-in-a-self-gravitating-two-layer-e0sR3YSVjd
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx507
Publisher site
See Article on Publisher Site

### Abstract

Summary The dispersion relation of the Rayleigh–Taylor instability in the spherical geometry is of profound importance in the context of the Earth's core formation. Here we present a complete derivation of this dispersion relation for a self-gravitating two-layer viscous sphere. Such relation is, however, obtained through the solution of a complex transcendental equation, and it is difficult to gain physical insights directly from the transcendental equation itself. We thus also derive an empirical formula to compute the growth rate, by combining the Monte Carlo sampling of the relevant model parameter space with linear regression. Our analysis indicates that the growth rate of Rayleigh–Taylor instability is most sensitive to the viscosity of inner layer in a physical setting that is most relevant to the core formation. Core, Numerical solutions, Planetary interiors 1 INTRODUCTION The Rayleigh–Taylor instability of a self-gravitating fluid sphere is important in the context of planet formation, especially in Earth's core formation (e.g. Ida et al. 1989; Stevenson 1990; Hoogenboom & Houseman 2006; Sasaki & Abe 2007). A dynamical situation typically envisioned for an early, growing Earth can be characterized by an inner layer of highly viscous silicate materials and an outer layer of denser molten metal (e.g. Honda et al. 1993; Rubie et al. 2015). This gravitationally unstable structure is susceptible to the Rayleigh–Taylor instability at the interface of these two layers. Such instability may lead to the dripping of the metallic layer into the silicate layer, which is an important part of core formation. However, the growth timescale of this Rayleigh–Taylor instability depends on several parameters such as viscosity, density and layer thicknesses, and an accurate dispersion relation becomes desirable. The dispersion relation for the Rayleigh–Taylor instability in a two-layer self-gravitating sphere was published previously by Ida et al. (1989) for both inviscid and viscous cases. This dispersion relation was utilized to discuss the Earth's core formation (Ida et al. 1987). Whereas the inviscid case is correctly presented, however, their derivation for the viscous case is incorrect. In this research note, we derive a correct dispersion relation for the Rayleigh–Taylor instability in a two-layer self-gravitating viscous sphere. As the growth rate depends on a number of parameters, it can be challenging to understand the variation of the growth rate in the high-dimensional parameter space. Also, the estimation of the growth rate in case of finite Prandtl number is based on the solution of a non-linear transcendental equation, which can be very cumbersome in the spherical geometry. Therefore, we develop a practical approach to derive an approximate formula for the growth rate by combining Monte Carlo sampling and linear regression for l = 2, the mode most relevant to the Earth's core formation. The l = 1 mode is the translational mode and results simply in an overall shift of the centre of mass of the whole spherical system. Despite being the fast growing mode for the case of μ1 ≫ μ2, therefore, it is not relevant to the core formation (e.g. Elsasser 1963; Honda et al. 1993). Modes with l > 1 are important for core formation, but with increasing l, the viscous dissipation becomes higher in proportion to ≈l(l + 1), thereby making growth timescale longer. The l = 2 mode thus has the maximum contribution, and we restrict ourselves to this particular case when deriving an approximate formula. The structure of this note is following. We begin with the theoretical formulation of the Rayleigh–Taylor instability in a self-gravitating two-layer viscous sphere. We then present a few representative dispersion curves for different combinations of layer viscosities. We next discuss the practical approach to obtain an approximate formula for growth rate. Our result suggests that when the viscosity of the inner layer is much higher than that of the outer layer, the growth rate is controlled mostly by the properties of the inner layer as well as the density contrast of the two layers. 2 DERIVATION OF THE DISPERSION RELATION The governing equations for an incompressible viscous fluid include the equation of continuity   $$\nabla \cdot {\bf U}=0,$$ (1)the equation of motion   $$\frac{{\rm d}\mathbf {U}}{{\rm d}t}=\nabla \cdot \mathbf {\Sigma }+{\bf F},$$ (2)and the constitutive relation   $$\mathbf {\Sigma }=-P\mathbf {I}+\mu \left[\nabla \mathbf {U}+(\nabla \mathbf {U})^{T}\right],$$ (3)where U is the total velocity field, Σ is the total stress, F is the body force, P is the total pressure, I is the identity operator, and μ is the dynamic viscosity. To focus on the linear Rayleigh–Taylor instability, we express the flow-related entities as first-order perturbations from the hydrostatic reference:   \begin{eqnarray} {\bf U}&=&{\bf U}_{0}+{\bf u}, \end{eqnarray} (4)  \begin{eqnarray} \mathbf {\Sigma }&=&\mathbf {\Sigma }_{0}+\mathbf {\sigma }, \end{eqnarray} (5)  \begin{eqnarray} P&=&P_{0}+\delta p, \end{eqnarray} (6)  \begin{eqnarray} {\bf F}&=&{\bf F}_{0}+{\bf f}, \end{eqnarray} (7)where [U0, Σ0, P0, F0]T and [u, σ, δp, f]T represent the hydrostatic reference state (U0 is trivially zero) and the perturbed state, respectively. In the hydrostatic state, we may write the body force (gravity) as a gradient of potential that satisfies the Poisson equation:   \begin{eqnarray} \mathbf {F}_{0}&=&-\rho _{0}\nabla \Psi , \end{eqnarray} (8)  \begin{eqnarray} \nabla ^{2} \Psi &=& -4\pi G \rho _{0}, \end{eqnarray} (9)where ρ0 and G are the density of the reference state and the gravitational constant, respectively. We may solve the Poisson equation for the hydrostatic state (eq. 9) and obtain the unperturbed gravitational potential as   \begin{eqnarray} \Psi &=& \frac{4\pi G}{3}\left(\rho _{1}(r^{2}/2-3R_{1}^{2}/2)+3\rho _{2}(R_{1}^{2}-R_{2}^{2})/2\right), \quad \quad \forall r<R_{1}, \end{eqnarray} (10)  \begin{eqnarray} \quad&=&\frac{4\pi G}{3}\left(\rho _{1}(r^{2}/2-3R_{2}^{2}/2)+(\rho _{2}-\rho _{1})R_{1}^{3}/r\right), \quad \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (11)where ρ1, ρ2, R1 and R2 are the density of the inner layer, the density of the outer layer, the radius of the inner sphere, and the radius of the outer sphere, respectively (Fig. 1). The hydrostatic state confirms the equality of the pressure gradient and the gravitational force:   $$\nabla P_{0}=\rho _{0}\nabla \Psi ,$$ (12)which will be used in boundary conditions. Subtracting the hydrostatic reference state and neglecting higher order terms lead to   \begin{eqnarray} \nabla _{j}u_{j}&=&0, \end{eqnarray} (13)  \begin{eqnarray} \nabla _{j}\sigma _{ij}+f_{i}&=&\frac{\partial u_{i}}{\partial t}, \end{eqnarray} (14)  \begin{eqnarray} \sigma _{ij}&=&-\delta p\delta _{ij}+\mu \left(\nabla _{j}u_{i}+\nabla _{i}u_{j}\right), \end{eqnarray} (15)where fi arises solely due to density perturbation at the interfaces of two adjacent fluids, and it can also be represented as the gradient of a potential perturbation δΨ satisfying the Poisson equation:   \begin{eqnarray} f_{i}&=&-\rho _{0}\nabla \delta \Psi , \end{eqnarray} (16)  \begin{eqnarray} \nabla ^{2} \delta \Psi &=&-4\pi G \delta \rho , \end{eqnarray} (17)where δρ is the perturbation to the reference density. Together with the perturbed eqs (13)–(17), we may write the equation of motion as   $$\frac{\partial \mathbf {u}}{\partial t}=-\nabla \Phi +\nu \nabla ^{2} \mathbf {u},$$ (18)where   $$\Phi =\frac{\delta P}{\rho }+\delta \Psi .$$ (19)The system of eqs (13)–(19) describing the dynamics of the perturbed state is used to calculate the growth rate of initial density perturbations. In the 3-D spherical geometry, we expand the field variables in the spherical harmonics assuming the horizontal invariance of material properties (e.g. Chandrasekhar 1961)   \begin{eqnarray} \Xi (r,\theta ,\phi ) &=& \sum _{l,m}\Xi _{lm}(r)Y^{m}_{l}(\theta ,\phi ). \end{eqnarray} (20)Using this separation of variables, we solve for the perturbed gravitational potential generated through the deformation of the interface (assuming the outer surface of the sphere to be rigid):   \begin{eqnarray} \delta \Psi (r)&=&\left(\frac{4\pi G(\rho _{2}-\rho _{1})}{(2l+1)R_{1}^{l-1}}\delta r_{1}\right) r^{l}, \quad \quad \forall r < R_{1}, \end{eqnarray} (21)  \begin{eqnarray} \qquad&=&\left(\frac{4\pi G (\rho _{2}-\rho _{1})}{2l+1}R_{1}^{l+2}\delta r_{1}\right)\bigg/r^{l+1}, \quad \quad \forall R_{1}<r<R_{2}. \end{eqnarray} (22)We can also readily solve for the total perturbed potential Φ by taking divergence of the eq. (18) and using the continuity eq. (13):   $$\nabla ^{2} \Phi =0,$$ (23)the solution of which may further be used to solve the equation of motion (18). However, before proceeding to solve the equation of motion, it is convenient to perform poloidal–toroidal decomposition in the spherical geometry (e.g. Chandrasekhar 1961) as we will only require the poloidal part for linear stability analysis. Using such decomposition, we may express the perturbed velocity field as   $$\mathbf {u}=\sum _{l,m}\left(\nabla \times (T_{lm}(r)Y^{m}_{l}\hat{\mathbf {r}})+\nabla \times \left[\nabla \times (S_{lm}(r)Y^{m}_{l}\hat{\mathbf {r}})\right]\right),$$ (24)where Slm and Tlm denote the poloidal and toroidal coefficients, respectively. Assuming u∝exp (σt) in the linear regime with σ being the growth rate, the poloidal part of the equation of motion (eq. 18) may be written as the following second-order differential equation:   $$\sigma S_{lm}=-\nabla \Phi +\nu \left[\frac{l(l+1)}{r^{2}}-D^{2}\right]S_{lm},$$ (25)where ν and D denote the kinematic viscosity and the radial derivative, that is, d/dr, respectively. The complementary function and particular solution of the eq. (25) can be obtained for the inner and outer layers separately as linear combinations of modified Bessel functions:   \begin{eqnarray} S_{lm}(r)&=&A_{1}\sqrt{(}r)I_{l+1/2}(q_{1}r)-\frac{Q_{1}R_{1}^{2}}{\sigma }\left(\frac{r}{R_{1}}\right)^{l+1}, \quad \quad \forall r<R_{1},\nonumber \\ &=&A_{2}\sqrt{(}r)I_{l+1/2}(q_{2}r)+B_{2}\sqrt{r}K_{l+1/2}(q_{2}r)-\left[Q_{2}\left(\frac{r}{R_{1}}\right)^{l+1}+Q_{3}\left(\frac{r}{R_{1}}\right)^{-l}\right]R_{1}^{2}/\sigma , \quad \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (26)where q1 = σ/ν1, q2 = σ/ν2, ν1 and ν2 are the viscosities of the inner and outer layers, respectively, and [A1, Q1, A2, B2, Q2, Q3]T is a set of unknown constants to be obtained using the boundary conditions. The velocity field associated with the Rayleigh–Taylor instability may be obtained from the poloidal potential Slm. The radial (ur) and transverse (uθ) components are obtained as follows:   \begin{eqnarray} u_{r}&=&\sum _{l,m}\frac{l(l+1)}{r^{2}}S_{lm}(r)Y_{l}^{m}(\theta ,\phi ),\nonumber \\ &=&\sum _{l,m}\left[\frac{A_{1}l(l+1)}{r^{3/2}}I_{l+1/2}(q_{1}r)-\frac{Q_{1}l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}\right]Y_{l}^{m}(\theta ,\phi ), \quad \forall r<R_{1}, \end{eqnarray} (27)  \begin{eqnarray} &&{\quad=\sum _{l,m}\Bigg[\frac{A_{2}l(l+1)}{r^{3/2}}I_{l+1/2}(q_{2}r)+\frac{B_{2}l(l+1)}{r^{3/2}}K_{l+1/2}(q_{2}r)}\nonumber \\ &&{\quad\quad-\frac{l(l+1)}{r^{2}}\left[Q_{2}\left(\frac{r}{R_{1}}\right)^{l+1}+Q_{3}\left(\frac{r}{R_{1}}\right)^{-l}\right]R_{1}^{2}/\sigma \Bigg]Y_{l}^{m}(\theta ,\phi ), \quad \forall R_{1}<r<R_{2},} \end{eqnarray} (28)  \begin{eqnarray} u_{\theta }&=&\sum _{l,m}\frac{1}{r}\frac{{\rm d}S_{lm}(r)}{{\rm d}r}\frac{\partial Y_{l}^{m}(\theta ,\phi )}{\partial \theta },\nonumber \\ &=&\sum _{l,m}\Big[A_{1}\left((l+1)r^{-3/2}I_{l+1/2}(q_{1}r)+q_{1}r^{-1/2}I_{l+3/2}(q_{1}r)\right)-\frac{Q_{1}(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}\Big]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall r<R_{1}, \end{eqnarray} (29)  \begin{eqnarray} &&{\quad=\sum _{l,m}\Bigg[A_{2}\left((l+1)r^{-3/2}I_{l+1/2}(q_{2}r)+q_{2}r^{-1/2}I_{l+3/2}(q_{2}r)\right)+B_{2}\left((l+1)r^{-3/2}K_{l+1/2}(q_{2}r)-q_{2}r^{-1/2}K_{l+3/2}(q_{2}r)\right)}\nonumber \\ &&{\quad\quad-\frac{1}{\sigma }\left(Q_{2}(l+1)\left(\frac{r}{R_{1}}\right)^{l-1}-Q_{3}l \left(\frac{r}{R_{2}}\right)^{-(l+2)}\right)\Bigg]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall R_{1}<r<R_{2}.} \end{eqnarray} (30)The stresses may be calculated from the velocities using the following constitutive relations:   \begin{eqnarray} \sigma _{rr}&=&2\mu \frac{\partial u_{r}}{\partial r}, \end{eqnarray} (31)  \begin{eqnarray} \sigma _{r\theta }&=&\mu \left[\frac{1}{r}\frac{\partial u_{r}}{\partial \theta }-\frac{u_{\theta }}{r}+\frac{\partial u_{\theta }}{\partial r}\right]. \end{eqnarray} (32) Figure 1. View largeDownload slide The spherical system considered in the Rayleigh–Taylor instability analysis. Unstable stratification corresponds to ρ2 > ρ1. Figure 1. View largeDownload slide The spherical system considered in the Rayleigh–Taylor instability analysis. Unstable stratification corresponds to ρ2 > ρ1. It is essential to compute the stresses correctly, particularly the normal stress, as it is used in the discontinuity condition, which is directly related to the expression of the growth rate. The previously published result by Ida et al. (1989) uses an incorrect expression for the normal and tangential stresses (their eqs A-21 and A-22); a simple dimensional analysis reveals incompatibility between the both sides of these equations. To make calculations simpler, we may express the velocities and stresses in the following manner:   \begin{eqnarray} u_{r}&=&\sum _{l,m}\left[a_{1}A_{1}+b_{1}Q_{1}\right]Y_l^m(\theta ,\phi ), \quad \forall r<R_{1}, \end{eqnarray} (33)  \begin{eqnarray} \quad&=&\sum _{l,m}\left[a_{2}A_{2}+b_{2}B_{2}+c_{2}Q_{2}+d_{2}D_{2}\right]Y_l^m(\theta ,\phi ), \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (34)  \begin{eqnarray} u_{\theta }&=&\sum _{l,m}\left[e_{1}A_{1}+f_{1}Q_{1}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall r<R_{1}, \end{eqnarray} (35)  \begin{eqnarray} \quad=\sum _{l,m}\left[e_{2}A_{2}+f_{2}B_{2}+g_{2}Q_{2}+h_{2}Q_{3}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (36)  \begin{eqnarray} \tau _{rr}&=&-P+2\mu _{1}\sum _{l,m}\left[i_{1}A_{1}+j_{1}Q_{1}\right]Y_l^m(\theta ,\phi ), \quad \forall r<R_{1}, \end{eqnarray} (37)  \begin{eqnarray} \quad&=&-P+2\mu _{2}\sum _{l,m}\left[i_{2}A_{2}+j_{2}B_{2}+k_{2}Q_{2}+l_{2}Q_{3}\right]Y_l^m(\theta ,\phi ), \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (38)  \begin{eqnarray} \tau _{r\theta }&=&\mu _{1}\sum _{l,m}\left[\alpha _{1}A_{1}+\beta _{1}Q_{1}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall r<R_{1}, \end{eqnarray} (39)  \begin{eqnarray} \quad\quad&=&\mu _{2}\sum _{l,m}\left[\alpha _{2}A_{2}+\beta _{2}B_{2}+\gamma _{2}Q_{2}+\delta _{2}Q_{3}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (40)where   \begin{eqnarray} a_{1}(r)&=&\frac{l(l+1)}{r^{3/2}}I_{l+1/2}(q_{1}r), \end{eqnarray} (41)  \begin{eqnarray} b_{1}(r)&=&-\frac{l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}, \end{eqnarray} (42)  \begin{eqnarray} a_{2}(r)&=&\frac{l(l+1)}{r^{3/2}}I_{l+1/2}(q_{2}r), \end{eqnarray} (43)  \begin{eqnarray} b_{2}(r)&=&\frac{l(l+1)}{r^{3/2}}K_{l+1/2}(q_{2}r), \end{eqnarray} (44)  \begin{eqnarray} c_{2}(r)&=&-\frac{l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}, \end{eqnarray} (45)  \begin{eqnarray} d_{2}(r)&=&-\frac{l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{-(l+2)}, \end{eqnarray} (46)  \begin{eqnarray} e_{1}(r)&=&\frac{(l+1)}{r^{3/2}}I_{l+1/2}(q_{1}r)+\frac{q_{1}}{\sqrt{r}}I_{l+3/2}(q_{1}r), \end{eqnarray} (47)  \begin{eqnarray} f_{1}(r)&=&-\frac{(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}, \end{eqnarray} (48)  \begin{eqnarray} e_{2}(r)&=&(l+1)r^{-3/2}I_{l+1/2}(q_{2}r)+q_{2}r^{-1/2}I_{l+3/2}(q_{2}r), \end{eqnarray} (49)  \begin{eqnarray} f_{2}(r)&=&(l+1)r^{-3/2}K_{l+1/2}(q_{2}r)-q_{2}r^{-1/2}K_{l+3/2}(q_{2}r), \end{eqnarray} (50)  \begin{eqnarray} g_{2}(r)&=&-\frac{(l+1)}{\sigma }(r/R_{1})^{l-1}, \end{eqnarray} (51)  \begin{eqnarray} h_{2}(r)&=&\frac{l}{\sigma }\left(\frac{r}{R_{2}}\right)^{-(l+2)}, \end{eqnarray} (52)  \begin{eqnarray} i_{1}(r)&=&l(l+1)\left[(l-1)r^{-5/2}I_{l+1/2}(q_{1}r)+q_{1}r^{-3/2}I_{l+3/2}(q_{1}r)\right], \end{eqnarray} (53)  \begin{eqnarray} j_{1}(r)&=&-\frac{l(l+1)(l-1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (54)  \begin{eqnarray} i_{2}(r)&=&l(l+1)\left[(l-1)r^{-5/2}I_{l+1/2}(q_{2}r)+q_{2}r^{-3/2}I_{l+3/2}(q_{2}r)\right], \end{eqnarray} (55)  \begin{eqnarray} j_{2}(r)&=&l(l+1)\left[(l-1)r^{-5/2}K_{l+1/2}(q_{2}r)-q_{2}r^{-3/2}K_{l+3/2}(q_{2}r)\right], \end{eqnarray} (56)  \begin{eqnarray} k_{2}(r)&=&-\frac{l(l+1)(l-1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (57)  \begin{eqnarray} l_{2}(r)&=&\frac{l(l+1)(l+2)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{-(l+3)}, \end{eqnarray} (58)  \begin{eqnarray} \alpha _{1}(r)&=&2(l+1)(l-1)r^{-5/2}I_{l+1/2}(q_{1}r)+(2l+1)q_{1}r^{-3/2}I_{l+3/2}(q_{1}r)+q_{1}^{2}r^{-1/2}I_{l+5/2}(q_{1}r), \end{eqnarray} (59)  \begin{eqnarray} \beta _{1}(r)&=&-\frac{2(l-1)(l+1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (60)  \begin{eqnarray} \alpha _{2}(r)&=&2(l+1)(l-1)r^{-5/2}I_{l+1/2}(q_{2}r)+(2l+1)q_{2}r^{-3/2}I_{l+3/2}(q_{2}r)+q_{2}^{2}r^{-1/2}I_{l+5/2}(q_{2}r), \end{eqnarray} (61)  \begin{eqnarray} \beta _{2}(r)&=&2(l+1)(l-1)r^{-5/2}K_{l+1/2}(q_{2}r)-(2l+1)q_{2}r^{-3/2}K_{l+3/2}(q_{2}r)+q_{2}^{2}r^{-1/2}K_{l+5/2}(q_{2}r), \end{eqnarray} (62)  \begin{eqnarray} \gamma _{2}(r)&=&-\frac{2(l-1)(l+1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (63)  \begin{eqnarray} \delta _{2}(r)&=&-\frac{2l(l+2)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{-(l+3)}. \end{eqnarray} (64) We are now ready to use boundary conditions to obtain the unknown coefficients. We assume the outer boundary to be rigid, that is, ur = uθ = 0 at r = R2. The inner boundary is assumed to be welded and hence the field variables are continuous across the deformed boundary. The first-order Taylor expansion leads to the continuity of the velocity and transverse stress across the r = R1 boundary, that is, $$\left[u_{r}\right]_{R_{1}-}^{R_{1}+}=0$$, $$\left[u_{\theta }\right]_{R_{1}-}^{R_{1}+}=0$$, and $$\left[\sigma _{r\theta }\right]_{R_{1}-}^{R_{1}+}=0$$, while the radial stress is continuous only across the deformed boundary, $$\left[-(P_{0}+\delta P)+2\mu \partial _{r}u_{r}\right]_{R_{1}-\delta r}^{R_{1}+\delta r}=0$$, and becomes discontinuous across r = R1 due to the density jump. P0 and δP may be calculated from eqs (12), (16), and (23). The rigid and continuous boundary conditions together with the orthonormality of the spherical harmonics form the following complete set of linear equations of the unknown coefficients for each (l, m) pair:   \begin{eqnarray} a_{2}(R_{2})A_{2}+b_{2}(R_{2})B_{2}+c_{2}(R_{2})Q_{2}+d_{2}(R_{2})Q_{3}&=&0, \end{eqnarray} (65)  \begin{eqnarray} e_{2}(R_{2})A_{2}+f_{2}(R_{2})B_{2}+g_{2}(R_{2})Q_{2}+h_{2}(R_{2})Q_{3}&=&0, \end{eqnarray} (66)  \begin{eqnarray} a_{1}(R_{1})A_{1}+b_{1}(R_{1})Q_{1}&=&a_{2}(R_{1})A_{2}+b_{2}(R_{1})B_{2}+c_{2}(R_{1})Q_{2}+d_{2}(R_{1})Q_{3}, \end{eqnarray} (67)  \begin{eqnarray} e_{1}(R_{1})A_{1}+f_{1}(R_{1})Q_{1}&=&e_{2}(R_{1})A_{2}+f_{2}(R_{1})B_{2}+g_{2}(R_{1})Q_{2}+h_{2}(R_{1})Q_{3}, \end{eqnarray} (68)  \begin{eqnarray} \mu _{1}\left[\alpha _{1}(R_{1})A_{1}+\beta _{1}(R_{1})Q_{1}\right]&=&\mu _{2}[\alpha _{2}(R_{1})A_{2}+\beta _{2}(R_{1})B_{2}+\gamma _{2}(R_{1})Q_{2}+\delta _{2}(R_{1})Q_{3}], \end{eqnarray} (69)  \begin{eqnarray} \frac{4\pi G(\rho _{1}-\rho _{2})[2(l-1)\rho _{1}+3\rho _{2}]}{3(2l+1)}\delta r&=&(l+1)(\rho _{1}Q_{1}-\rho _{2}Q_{2})+l\rho _{2}Q_{3}+2\mu _{2}/R_{1}\partial _{r}u_{r2}-2\mu _{1}/R_{1}\partial _{r}u_{r1}. \end{eqnarray} (70)We may eliminate the several unknowns from these equations and obtain the following form:   \begin{eqnarray} A_{2}&=&\xi _{1}Q_{2}+\eta _{1}Q_{3}, \end{eqnarray} (71)  \begin{eqnarray} B_{2}&=&\xi _{2}Q_{2}+\eta _{2}Q_{3}, \end{eqnarray} (72)  \begin{eqnarray} A_{1}&=&X_{1}Q_{2}+Y_{1}Q_{3}, \end{eqnarray} (73)  \begin{eqnarray} Q_{1}&=&X_{2}Q_{2}+Y_{2}Q_{3}, \end{eqnarray} (74)  $$Q_{2}=\Theta Q_{3},$$ (75)  $$Q_{3}=\Sigma \delta r,$$ (76)where   \begin{eqnarray} \xi _{1}&=&-\frac{f_{2}(R_{2})c_{2}(R_{2})-b_{2}(R_{2})g_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (77)  \begin{eqnarray} \eta _{1}&=&-\frac{f_{2}(R_{2})d_{2}(R_{2})-b_{2}(R_{2})h_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (78)  \begin{eqnarray} \xi _{2}&=&-\frac{a_{2}(R_{2})g_{2}(R_{2})-e_{2}(R_{2})c_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (79)  \begin{eqnarray} \eta _{2}&=&-\frac{a_{2}(R_{2})h_{2}(R_{2})-e_{2}(R_{2})d_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (80)  \begin{eqnarray} X_{1}&=&\frac{f_{1}(R_{1})[a_{2}(R_{1})\xi _{1}+b_{2}(R_{1})\xi _{2}+c_{2}(R_{1})]-b_{1}(R_{1})[e_{2}(R_{1})\xi _{1}+f_{2}(R_{1})\xi _{2}+g_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (81)  \begin{eqnarray} Y_{1}&=&\frac{f_{1}(R_{1})[a_{2}(R_{1})\eta _{1}+b_{2}(R_{1})\eta _{2}+d_{2}(R_{1})]-b_{1}(R_{1})[e_{2}(R_{1})\eta _{1}+f_{2}(R_{1})\eta _{2}+h_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (82)  \begin{eqnarray} X_{2}&=&\frac{a_{1}(R_{1})[e_{2}(R_{1})\xi _{1}+f_{2}(R_{1})\xi _{2}+g_{2}(R_{1})]-e_{1}(R_{1})[a_{2}(R_{1})\xi _{1}+b_{2}(R_{1})\xi _{2}+c_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (83)  \begin{eqnarray} Y_{2}&=&\frac{a_{1}(R_{1})[e_{2}(R_{1})\eta _{1}+f_{2}(R_{1})\eta _{2}+h_{2}(R_{1})]-e_{1}(R_{1})[a_{2}(R_{1})\eta _{1}+b_{2}(R_{1})\eta _{2}+d_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (84)  $$\Theta =\frac{\mu _{2}[\alpha _{2}(R_{1})\eta _{1}+\beta _{2}(R_{1})\eta _{2}+\delta _{2}(R_{1})]-\mu _{1}[\alpha _{1}(R_{1})Y_{1}+\beta _{1}(R_{1})Y_{2}]}{\mu _{1}[\alpha _{1}(R_{1})X_{1}+\beta _{1}(R_{1})X_{2}]-\mu _{2}[\alpha _{2}(R_{1})\xi _{1}+\beta _{2}(R_{1})\xi _{2}+\gamma _{2}(R_{1})]},$$ (85)  $$\Sigma =\frac{4\pi G (\rho _{1}-\rho _{2})(2(l-1)\rho _{1}+3\rho _{2})}{3(2l+1)\Pi },$$ (86)and   \begin{eqnarray} \Pi &=&\left[\left((l+1)\rho _{1}X_{2}-(l+1)\rho _{2}\right)+\frac{2\mu _{2}}{R_{1}}\left(i_{2}(R_{1})\xi _{1}+j_{2}(R_{1})\xi _{2}+k_{2}(R_{1})\right)-\frac{2\mu _{1}}{R_{1}}\left(i_{1}(R_{1})X_{1}+j_{1}(R_{1})X_{2}\right)\right]\Theta \nonumber\\ &&+\left[(l+1)\rho _{1}Y_{2}+l\rho _{2}+\frac{2\mu _{2}}{R_{1}}\left(i_{2}(R_{1})\eta _{1}+j_{2}(R_{1})\eta _{2}+l_{2}(R_{1})\right)-\frac{2\mu _{1}}{R_{1}}\left(i_{1}(R_{1})Y_{1}+j_{1}(R_{1})Y_{2}\right)\right]. \end{eqnarray} (87)We may write the relation between the amplitude of perturbation and the vertical velocity at the interface as   \begin{eqnarray} \frac{{\rm d}(\delta r)}{{\rm d}t}&=&u_{r}(R_{1}). \end{eqnarray} (88)By noting that δr ∝ exp (σt) in the linear regime, eqs (34) and (88) together with the six unknowns (i.e. A1-Q3) may be combined to obtain the expression for the growth rate   $$\sigma =\left(a_{1}(R_{1})[X_{1}\Theta +Y_{1}]+b_{1}(R_{1})[X_{2}\Theta +Y_{2}]\right)\Sigma .$$ (89)This is a transcendental equation for σ, because its right-hand side depends on σ in a complex manner. 3 EXAMPLES AND APPROXIMATE FORMULA In this section, we provide examples of dispersion relation for several combinations of layer viscosities. It is convenient to present the growth rate in the dimensionless form, so we scale the growth rate with respect to the free fall rate as follows:   $$\sigma^{\prime}=\sigma /\sqrt{4\pi G (\rho _{2}-\rho _{1})/3}.$$ (90)First, we fix the viscosity of the inner layer (i.e. μ1) and compute the growth rate for a range of outer layer viscosity (i.e. μ2). Similarly, we compute the growth rates for a range of inner layer viscosity keeping the viscosity of the outer layer constant. In both of these cases, the primary contributing factor turns out to be max (μ1, μ2) (Fig. 2). For μ1 ≫ μ2, the translational mode approximately corresponds to the motion of a rigid sphere through a viscous fluid. The corresponding growth rate turns out to be the fastest; the inner rigid sphere almost instantaneously translates to a new position leading to an overall shift of the centre of mass. Viscosity variation therefore does not affect considerably the growth rate for this translational mode as long as μ1 ≫ μ2 (Figs 2a and b). The growth rates of modes with l ≥ 2 exhibit nearly linear dependence on the inner layer viscosity (Fig. 2b). In the case of μ2 ≫ μ1, little variation in growth rate is observed with μ1, while growth rates associated with modes with small angular order (i.e. l < 5) vary linearly with μ2 (Figs 2c and d). However, the case of μ2 ≫ μ1 is not relevant in the context of core formation. Figure 2. View largeDownload slide The dispersion relations for several combinations of inner and outer layer viscosities: (a) the case where μ2 is varied keeping μ1 constant and μ1 ≫ μ2, (b) the case where μ1 is varied keeping μ2 constant and μ1 ≫ μ2—most relevant to the core formation, (c) the case where μ1 is varied keeping μ2 constant and μ2 ≫ μ1, (d) the case where μ2 is varied keeping μ1 constant and μ2 ≫ μ1. Other parameters are fixed as follows: R1 = 106 m, dR = (R2 − R1) = 104.84 m, ρ1 = 103.74 kg m−3 and δρ = 103.5 kg m−3. Figure 2. View largeDownload slide The dispersion relations for several combinations of inner and outer layer viscosities: (a) the case where μ2 is varied keeping μ1 constant and μ1 ≫ μ2, (b) the case where μ1 is varied keeping μ2 constant and μ1 ≫ μ2—most relevant to the core formation, (c) the case where μ1 is varied keeping μ2 constant and μ2 ≫ μ1, (d) the case where μ2 is varied keeping μ1 constant and μ2 ≫ μ1. Other parameters are fixed as follows: R1 = 106 m, dR = (R2 − R1) = 104.84 m, ρ1 = 103.74 kg m−3 and δρ = 103.5 kg m−3. In order to study the Rayleigh–Taylor instability associated with core formation in a comprehensive manner, it would be convenient to derive an approximate formula for the growth rate for the most important mode associated with the deformation of the interface. To derive such a formula, we explore the following parameter ranges: R1 = 105 to 107 m, ΔR = (R2 − R1) = 103 to 5 × 106 m, ρ1 = 3 × 103 to 104 kg m−3, δρ = (ρ2 − ρ1) = 103 to 104 kg m−3, μ1 = 1015 to 1025 Pa s, and μ2 = 106 to 1010 Pa s. First, we consider a reference case with R1 = 106 m, ΔR = 104.84 m, ρ1 = 103.74 kg m−3, δρ = 103.5 kg m−3, μ1 = 1020 Pa s, and μ2 = 108 Pa s. The model space to be explored can be defined around this reference case with the following standard deviation:   $$\zeta (\mathbf {x})=[0.5, 0.5,0.5,1,3,1.5]^{T},$$ (91)where   $$\mathbf {x}=[\log _{10}(R_{1}),\log _{10}(\rho _{1}),\log _{10}(\delta \rho ),\log _{10}(\Delta R),\log _{10}(\mu _{1}),\log _{10}(\mu _{2})]^{T}.$$ (92)This model space is subjected to random sampling using the normal distribution, and a model vector (eq. 92) is constructed for each sample. This obtained set of random model vectors is then used to compute the growth rate corresponding to the l = 2 mode. An ensemble of growth rates is obtained through 200 such calculations. Despite the complexity of the transcendental equation (eq. 89), the growth rate is found to be predicted with reasonable accuracy using the following function:   $$\log _{10}(\sigma )=a_{0}+\sum _{j=1}^{6}a_{j}X_{j},$$ (93)where Xj is the mean subtracted input parameter that is normalized by the corresponding standard deviation (eq. 91). Constants a0 through a6 are estimated by least square regression: a0 = −10.43, a1 = 1.86, a2 = 0.87, a3 = 1.3, a4 = 0.047, a5 = −0.9752, and a6 = −0.054. Fig. 3 demonstrates a reasonable correspondence between the predicted and actual values of growth rates for the aforementioned 200 computations. This way of summarizing enables us not only to see the sensitivity of the growth rate on the model parameters, but also to quickly compute the growth rate without doing the whole calculation. We suggest the following expression for the growth rate corresponding to the l = 2 mode:   $$\sigma =10^{-10.13}R_{1}^{1.86}(\rho _{2}-\rho _{1})^{1.3}\rho _{1}^{0.87}(R_{2}-R_{1})^{0.047}\mu _{1}^{-0.98}\mu _{2}^{-0.054},$$ (94)where R, ρ, μ and σ are in m, kg m−3, Pa s and s−1, respectively. Figure 3. View largeDownload slide Correspondence between the actual and predicted values of growth rate for l = 2. Figure 3. View largeDownload slide Correspondence between the actual and predicted values of growth rate for l = 2. ACKNOWLEDGEMENTS This study is based on work supported by the NASA through the NASA Astrobiology Institute under cooperative agreement no. NNA15BB03A issued through the Science Mission Directorate. We thank Gael Choblet and Neil Ribe for their careful reviews and constructive comments. REFERENCES Chandrasekhar S., 1961. Hydrodynamic and Hydromagnetic Stability  Clarendon. Elsasser W., 1963. Early history of the Earth, in Earth Science and Meteoritics  pp. 1– 30, eds Geiss J., Goldberg E., North-Holland. Honda R., Mizutani H., Yamamoto T., 1993. Numerical simulation of Earth's core formation, J. geophys. Res.  98 2075– 2089. Google Scholar CrossRef Search ADS   Hoogenboom T., Houseman G.A., 2006. Rayleigh–Taylor instability as a mechanism for corona formation on Venus, Icarus  180 292– 307. Google Scholar CrossRef Search ADS   Ida S., Nakagawa Y., Nakazawa K., 1987. The Earth's core formation due to the Rayleigh-Taylor instability, Icarus  69 239– 248. Google Scholar CrossRef Search ADS   Ida S., Nakagawa Y., Nakazawa K., 1989. The Rayleigh–Taylor instability in a self-gravitating two-layer fluid sphere, Earth Moon Planets  44 149– 174. Google Scholar CrossRef Search ADS   Rubie D., Nimmo F., Melosh H., 2015. Formation of the Earth's core, in Treatise On Geophysics  2nd edn, vol. 9, pp. 51– 90, ed. Stevenson D., Elsevier. Sasaki T., Abe Y., 2007. Rayleigh–Taylor instability after giant impacts: imperfect equilibration of the Hf-W system and its effect on the core formation age, Earth Planets Space  59 1035– 1045. Google Scholar CrossRef Search ADS   Stevenson D., 1990. Fluid dynamics of core formation, in Origin of the Earth  pp. 231– 249, eds Newsom H., Jones J.E., Oxford Univ. Press. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations