Scaling of Small-scale Dynamo Properties in the Rayleigh–Taylor Instability

The Rayleigh–Taylor instability (RTI) is ubiquitous in astrophysical contexts due to the generality of the conditions needed for its onset. The RTI operates in regions where the density gradient is misaligned with the direction of the local gravitational field or acceleration (Chandrasekhar 1961). The instability characteristically evolves with rising bubbles of lighter fluid and sinking spikes of heavier fluid that propagate away from the unstable region, leading to mixed fluids and relaxation of the unstable density gradient. In this manner, fluid mixing and transport is enhanced by the RTI in many astrophysical scenarios such as the solar corona (Isobe et al. 2005; Berger et al. 2011; Hillier 2018), gamma-ray burst scenarios (Gull & Longair 1973; Levinson 2009; Duffell & Macfadyen 2013; Duffell & MacFadyen 2014), and supernova explosions (Hillebrandt & Niemeyer 2000; Cabot & Cook 2006; Duffell & Kasen 2016). The fluid is typically a highly conducting plasma whose ability to generate and sustain magnetic fields can make the magnetohydrodynamic (MHD) version of the RTI differ from its hydrodynamic counterpart.

The role of the magnetic field can be categorized by whether the initial field strength is dynamically strong or weak. In the strong field limit, large-scale magnetic fields are able to stabilize a wide range of wavenumbers due to the restoring force of magnetic tension (Chandrasekhar 1961; Ruderman et al. 2014). This affects the nonlinear saturation and mixing rates of the RTI with a strong dependence on details of the initial geometry of the magnetic fields. Typically, magnetic tension forces suppress development of secondary shear instabilities and can cause bubble/spike structures to rise/fall more rapidly than in the hydrodynamic case (Stone & Gardiner 2007). This limit has application in the contexts of solar prominences (Hillier 2018) and pulsar wind nebulae (Porth et al. 2014), for example.

On the other hand, the weak field limit corresponds to an initial magnetic field with dynamically insignificant strength, which allows the RTI to evolve purely hydrodynamically, at least at early times. The turbulence in the nonlinear phases of the RTI has the potential to amplify the magnetic field through dynamo action. In the dynamo literature, this weak field limit is also known as the kinematic dynamo regime. If the turbulence is sufficiently vigorous and the initial magnetic energy is not too small, then the magnetic energy could grow to and saturate in equipartition with the kinetic energy of the turbulence. This results in decaying MHD turbulence post-saturation of the RTI. Otherwise, the magnetic energy is amplified, but remains at sub-equipartition energies and results in decaying hydrodynamic turbulence post-saturation of the RTI.

While early, low-resolution simulations of the RTI confirmed operation of the small-scale dynamo (SSD; Jun et al. 1995), the quantitative scaling of dynamo properties has been largely unstudied despite applications in a variety of scenarios. Here, we discuss two example applications at opposite extremes: the fireball and laser plasma experiments. Hydrodynamical simulations of the fireball propose that the high-Reynolds-number RTI-driven turbulence generates equipartition magnetic fields that can explain the observed levels of synchrotron radiation emission from these sources (Duffell & Macfadyen 2013). On the other hand, several laser plasma experiments of the RTI at modest Reynolds numbers have found magnetic energy growth and explain the observations by using the Biermann effect (Gao et al. 2012; Manuel et al. 2012; Nilson et al. 2015; Matteucci et al. 2018). However, RTI turbulence itself could be a significant contributor as well, in particular as the Reynolds number and duration of future experiments increases (Galmiche & Gauthier 1996; Bott et al. 2021). A quantitative framework for either justifying the use of the equipartition argument or otherwise predicting the level of sub-equipartition magnetic energy amplification would be useful in the general case. This motivates a careful study of the RTI-turbulence-driven dynamo.

Paper Outline—Section 2 presents a model for the SSD in the weak field limit based on assumptions of isotropic turbulence and freefall scaling for the outer velocity and length scale of the turbulence. The model makes predictions for scaling laws between the SSD growth rate, amplification factor, and parameters of the RTI (Atwood number, gravitational acceleration, viscosity, and length scale) in each phase of evolution of the RTI. Determination of the correct scaling laws is important for quantitatively extending results from practical simulation parameters to realistic experimental or astrophysical parameters, which can be separated by orders of magnitude.

Section 3 tests the model using sets of three-dimensional (3D) visco-resistive MHD direct numerical simulations that resolve the turbulent viscous scales, allowing for results independent of numerical resolution. Discrepancies between the model prediction and the numerical results are analyzed and the assumptions that are likely breaking down are identified. Section 4 applies the model to the case of possible dynamo action in RTI-unstable regions of the outer envelope of the post-merger of a binary neutron star collision. In Section 5, we summarize our findings, identify directions for future work, and conclude.

We present scaling arguments to determine the kinematic small-scale dynamo growth rate and amplification factor of the magnetic energy in a Rayleigh–Taylor unstable, conducting, collisional fluid over the course of the growth, saturation, and decay of the instability.

Setup—Without loss of generality, we consider the standard, idealized RT setup in three-dimensions with a discontinuous density jump from ρ(z) = ρ_b for z < 0 to ρ(z) = ρ_t > ρ_b for z > 0 in a bounded domain of size ∼ L and initial velocity perturbations near z = 0. Any similar setup (e.g., with a continuous positive density gradient instead) can be easily related to the idealized setup, and the same scaling arguments will apply. The initial seed magnetic field is taken to be random at all scales and arbitrarily weak so that the magnetic energy always remains much lower than the kinetic energy at all hydrodynamic scales, resulting in purely hydrodynamic evolution of the velocity field. A standard model of the RTI in fully ionized plasmas is the set of visco-resistive MHD equations given by:

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{u}}\right)=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\partial }_{t}(\rho {\boldsymbol{u}})+{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{u}}{\boldsymbol{u}}-{\boldsymbol{B}}{\boldsymbol{B}}+P\right)=-\rho g\hat{z}+\rho \nu {{\rm{\nabla }}}^{2}{\boldsymbol{u}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{B}}-{\rm{\nabla }}\times \left({\boldsymbol{u}}\times {\boldsymbol{B}}\right)=\eta {{\rm{\nabla }}}^{2}{\boldsymbol{B}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\partial }_{t}E+{\rm{\nabla }}\cdot \left[(E+P){\boldsymbol{u}}-{\boldsymbol{B}}({\boldsymbol{B}}\cdot {\boldsymbol{u}})\right]=-\rho {{gu}}_{z}+\rho \kappa {{\rm{\nabla }}}^{2}T,\end{eqnarray} \tag{ 4 }$$

where u is the velocity, B is the magnetic field, E is the total energy density, P is the sum of the gas and magnetic pressure, T is the temperature, g is the gravitational acceleration, ν is the viscosity, η is the resistivity, and κ is the thermal diffusivity.

We consider the subsonic limit where the fluid flows are nearly incompressible. Any initial vertical velocity perturbation with horizontal wavenumber k triggers the RTI and grows exponentially with rate $n={\left({gkA}+{\nu }^{2}{k}^{4}\right)}^{\tfrac{1}{2}}-\nu {k}^{2}$, where $A=\tfrac{{\rho }_{t}-{\rho }_{b}}{{\rho }_{t}+{\rho }_{b}}\leqslant 1$ is the Atwood number and effects of thermal diffusivity are neglected. If all wavenumbers are present, the fastest-growing wavenumber will be ${k}_{c}\sim {\left({Ag}/{\nu }^{2}\right)}^{\tfrac{1}{3}}$ and smaller wavenumbers k < k_c will approximately have the inviscid growth rate $n\approx {\left({Agk}\right)}^{\tfrac{1}{2}}$.

RTI evolution and definitions—The RTI can be split into four distinct phases: linear growth, mixing (or nonlinear), saturation, and decay. The linear phase begins with exponential growth of any initial vertical velocity perturbations and ends once the fluid has displaced a vertical distance comparable to the wavelength of the initial mode. The fluid then rapidly becomes turbulent with outer velocity scale u(t), integral length scale l_i (t), and Reynolds number ${Re}(t)=u(t){l}_{i}(t)/(2\pi \nu )$ evolving in time throughout the remaining three phases. It is convenient in what follows to define a characteristic velocity ${u}_{\mathrm{sat}}=\sqrt{{AgL}}$, dynamical time t_dyn = L/u_sat, and Reynolds number ${{Re}}_{\mathrm{sat}}\,=\,{u}_{\mathrm{sat}}L/(2\pi \nu )$.

Following the linear phase, the initial perturbations develop in the mixing phase into the characteristic RTI bubble/spike structures that rise/fall away from the boundary, driven by buoyancy forcing. The mixing phase ends when the upward/downward propagating fronts of bubbles/spikes reach the system scale L. This begins the saturation phase of the RTI where the majority of available potential energy has been released into turbulent kinetic energy. After a dynamical timescale, the viscous dissipation in the turbulence becomes larger than the energy input from buoyant forcing and thus begins the decay phase where the total kinetic energy decreases. The fluid eventually settles into a stable state with a negative density gradient.

The small-scale dynamo will be active when the Reynolds number Re(t) in the turbulent phases is above the critical Reynolds number ${{Re}}^{c}({Pm})$ of the flows, where Pm = ν/η is the magnetic Prandtl number. We characterize the dynamo with the instantaneous, exponential dynamo growth rate $\gamma (t)\,=\tfrac{d}{{dt}}\mathrm{ln}{ME}(t)$ and the net exponential amplification factor of the magnetic energy given by:

$$\begin{eqnarray}&&{\rm{\Delta }}=\mathop{\max }\limits_{t}\mathrm{ln}\left(\displaystyle \frac{{ME}(t)}{{ME}(0)}\right)=\mathop{\max }\limits_{t}{\int }_{0}^{t}\gamma (t^{\prime} ){dt}^{\prime} ,\end{eqnarray} \tag{ 5 }$$

where $\mathrm{ME}(t)=\tfrac{1}{2}\int {\boldsymbol{B}}{\left({\boldsymbol{x}},t\right)}^{2}{d}^{3}x$ is the total magnetic energy. If u, l_i , and Re are enough to characterize the flow, then dimensional constraints force $\gamma (t)={C}_{\gamma }(u/{l}_{i}){{Re}}^{{D}_{\gamma }}$ (assuming a power-law form for the Re dependence when ${Re}\,\gg \,{{Re}}^{c}$), where C_γ and D_γ are dimensionless constants. Since we consider the case of an arbitrarily weak initial magnetic energy, the dynamo always remains in the kinematic regime and the dynamo growth rate is thus dominated by the fastest eddy turnover time in the turbulence of each phase of the RTI. In isotropic turbulence at high Pm, the growth is known to scale with the turnover time of Kolmogorov-scale eddies $\gamma \sim \delta {v}_{{l}_{\nu }}/{l}_{\nu }\sim (u/{l}_{i}){{Re}}^{1/2}$, where $\delta {v}_{l}={\left(\epsilon l\right)}^{1/3}$ is the velocity increment at length scale l, ∼ u³/l_i is the turbulent kinetic energy dissipation rate, and the Kolmogorov length scale l_ν is defined by equating the momentum diffusion time and eddy turnover time $\nu /{l}_{\nu }^{2}\sim \delta {v}_{{l}_{\nu }}/{l}_{\nu }$ for an eddy of size l ∼ l_ν (Rincon 2019). This corresponds to D_γ = 1/2 . In the following analysis, we restrict to the high Pm limit and assume that the turbulence is sufficiently isotropic at viscous scales such that D_γ = 1/2 in the mixing, saturation, and decay phases of the RTI. However, this assumption may need to be revisited in the mixing and decay phases, where the anisotropic effects of gravity are particularly important. Not only does gravity affect horizontal versus vertical fluid motions differently, there is also asymmetry between rising and falling fluid structures at increasingly higher Atwood numbers (for a review, see Zhou 2017).

Model outline—In summary, we model the evolving velocity field driven by the RTI as isotropic turbulence with time dependent outer velocity scale u(t) and integral scale l_i (t) that drives the SSD with a time dependent dynamo growth rate $\gamma (t)={C}_{\gamma }(u/{l}_{i}){{Re}}^{\tfrac{1}{2}}$. The net amplification of magnetic energy, Δ, can then be found by integrating γ(t) in time across the duration of each phase of the RTI. We emphasize that we assume that the dynamo remains in the kinematic regime while the RTI evolves through the mixing, saturation, and decay phases, which is valid as long as the magnetic energy remains smaller than the kinetic energy at viscous scales.

A main goal of the paper is to determine how Δ scales with parameters of the RTI, {A, g, L, ν}. We begin now by splitting up the time integral in Equation (5) into the three consecutive, turbulent phases and consider the positive contribution of each phase one at a time:

$$\begin{eqnarray}&&{\rm{\Delta }}={{\rm{\Delta }}}^{\mathrm{mix}}+{{\rm{\Delta }}}^{\mathrm{sat}}+{{\rm{\Delta }}}^{\mathrm{decay}}.\end{eqnarray} \tag{ 6 }$$

RTI Mixing phase—In the mixing phase, a mixing region of vertical extent h(t) propagates away from the original interface. A variety of models have been proposed to predict the time dependence of h(t), and the results of simulations and experiments are on various levels of disagreement that are primarily attributed to a sensitivity to initial conditions, system size effects, or effects of diffusivity. The general form of most models predict a freefall scaling given by

$$\begin{eqnarray}&&h(t)/L=\displaystyle \frac{1}{2}\alpha {\tau }^{2}+2\sqrt{\alpha {h}_{0}/L}\tau +{h}_{0}/L,\end{eqnarray} \tag{ 7 }$$

where τ = t/t_dyn, α is a constant, and h₀ is the length scale near when the mixing region first became nonlinear (see the review by Boffetta & Mazzino 2017 and references within). The turbulent velocity and speed of the mixing region boundary are taken to be comparable and self-consistently given by $u(t)/{u}_{\mathrm{sat}}\,=\alpha \tau +2\sqrt{\alpha {h}_{0}/L}$. These scalings have been found to break down when the size of the bubbles and spikes become comparable to the horizontal or vertical extent of the domain or if there is a presence of a dominant mode l_D , both of which could lead to terminal velocity scalings $h(t)\sim \sqrt{{{gl}}_{D}}t$ and $u(t)\sim \sqrt{{{gl}}_{D}}$ (Dimonte 2004; Banerjee & Andrews 2009; Lecoanet et al. 2012). Additionally, a larger, absolute thermal diffusivity reduces buoyancy effects and therefore slows the development of the mixing region (Abarzhi 2010). For instance, dimensional arguments including diffusivity at low Atwood number can predict an alternative time dependence $h(t)\sim {{gt}}^{2}/\mathrm{ln}({{gt}}^{2}/{h}_{0})$ (Abarzhi et al. 2005).

To model the dynamo growth rate, we assume that the turbulent integral scale is comparable to the height of the mixing region l_i (t) ≈ h(t) (Chertkov 2003; Boffetta & Mazzino 2017). However, this is an assumption that may also need to be revisited because the turbulence is not in a steady state and there may be a delay between energy generation at large scales and dissipation at small scales (Livescu et al. 2009). For simplicity, we also assume the freefall scalings hold. Then the dynamo exponentially grows with rate:

$$\begin{eqnarray}&&\begin{array}{l}{\gamma }_{\mathrm{mix}}(t)\approx \displaystyle \frac{{C}_{\gamma }{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}}{{t}_{\mathrm{dyn}}}\\ \quad {\left(\displaystyle \frac{{\left(\alpha \tau +2\sqrt{\alpha {h}_{0}/L}\right)}^{3}}{\tfrac{1}{2}\alpha {\tau }^{2}+2\sqrt{\alpha {h}_{0}/L}\tau +{h}_{0}/L}\right)}^{\tfrac{1}{2}}.\end{array}\end{eqnarray} \tag{ 8 }$$

After a transient time $\tau \approx \sqrt{{h}_{0}/(\alpha L)}$, the growth will asymptotically scale as ${\gamma }_{\mathrm{mix}}(t)\sim {t}_{\mathrm{dyn}}^{-1}{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}{\tau }^{\tfrac{1}{2}}$. Depending on the geometry, the mixing phase will last for a time interval Δt_mix ∼ t_dyn when the mixing region reaches the vertical system scale h(t) ≈ L. Thus, with freefall scalings, Δ^mix is asymptotically found to be

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{\mathrm{mix}}\approx \displaystyle \frac{2\sqrt{2}\alpha {C}_{\gamma }}{3}{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}{\left(\displaystyle \frac{{\rm{\Delta }}{t}_{\mathrm{mix}}}{{t}_{\mathrm{dyn}}}\right)}^{\tfrac{3}{2}}.\end{eqnarray} \tag{ 9 }$$

RTI Saturation phase—In the saturation phase, the mixing region encompasses the entire domain h(t) = L ≈ l_i (t) and the turbulent kinetic energy is maximal. The velocity scale $u(t)\approx \sqrt{2\alpha }{u}_{\mathrm{sat}}$ can be solved either from the freefall scaling in the mixing phase or by equating converted initial potential energy PE ∼ (ρ_t − ρ_n )α gL⁴ to turbulent kinetic energy ${KE}\sim \tfrac{1}{2}({\rho }_{t}+{\rho }_{b}){u}^{2}{L}^{3}$. The dynamo growth rate is then ${\gamma }_{\mathrm{sat}}(t)\approx {C}_{\gamma }{\left(2\alpha \right)}^{\tfrac{3}{4}}{t}_{\mathrm{dyn}}^{-1}{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}$. We expect the velocity scale and integral length scale to remain roughly constant for a time interval Δt_sat ∼ t_dyn comparable to the dynamical time before the turbulence enters the decay phase. The exponential amplification from this phase should scale as

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{\mathrm{sat}}\approx {C}_{\gamma }{\left(2\alpha \right)}^{\tfrac{3}{4}}{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}\left(\displaystyle \frac{{\rm{\Delta }}{t}_{\mathrm{sat}}}{{t}_{\mathrm{dyn}}}\right).\end{eqnarray} \tag{ 10 }$$

RTI Decay phase—From this point, the turbulence decays and the dynamo growth rate decreases. Again for simplicity, we assume that the free decaying turbulence is isotropic. Although the isotropy assumption is not expected to hold as the system relaxes into a stably stratified state, the simple model will be useful as a point of comparison. The total energy in isotropic turbulence decays as

$$\begin{eqnarray}&&\displaystyle \frac{{dE}(t^{\prime} )}{{dt}^{\prime} }=-\displaystyle \frac{1}{2}\overline{\rho }{C}_{E}\displaystyle \frac{u{\left(t^{\prime} \right)}^{3}}{{l}_{i}(t^{\prime} )},\end{eqnarray} \tag{ 11 }$$

where $u(t^{\prime} )={\left(2E(t^{\prime} )/\overline{\rho }\right)}^{\tfrac{1}{2}}$, $\overline{\rho }$ is the average density, $t^{\prime} =0$ is the beginning of the decay phase, and C_E is a dimensionless constant (Frisch 1995; Subramanian et al. 2006). The problem lies in determining the relation between the evolving integral scale and energy (${l}_{i}(t^{\prime} )\sim {E}^{-s}$). Typical choices for s lie in the range $\tfrac{1}{5}\leqslant s\leqslant \tfrac{1}{3}$, with s = 1/5 corresponding to Batchelor-type and s = 1/3 to Saffman-type turbulence, depending on properties and initial conditions of the turbulence (Ishida et al. 2006). However, if the integral scale is already at the system scale and the system is bounded (i.e., in a numerical simulation), then l_i cannot grow and thus remains at the system scale (${l}_{i}(t^{\prime} )\approx L$) independent of E (corresponding to s → 0) (Skrbek & Stalp 2000; Touil et al. 2002). Solving this system for $u(t^{\prime} )$ and ${l}_{i}(t^{\prime} )$, the growth rate is given by:

$$\begin{eqnarray}&&{\gamma }_{\mathrm{dec}}(t^{\prime} )\approx \displaystyle \frac{{C}_{\gamma }{\left(2\alpha \right)}^{\tfrac{3}{4}}{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}}{{t}_{\mathrm{dyn}}}{\left(1+(s+\displaystyle \frac{1}{2}){C}_{E}\tau ^{\prime} \right)}^{\tfrac{-3-2s}{2(1+2s)}}.\end{eqnarray} \tag{ 12 }$$

The dynamo growth rate will eventually become zero when ${Re}(t^{\prime} )\approx {{Re}}^{c}({Pm})$. Assuming the dynamo is active in the decay phase from $\tau ^{\prime} =0$ to $\tau ^{\prime} ={\rm{\Delta }}{t}_{\mathrm{dec}}/{t}_{\mathrm{dyn}}\gg 1/{C}_{E}$, Δ^decay is given by:

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{\mathrm{decay}}\approx {{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}\displaystyle \frac{4{\left(2\alpha \right)}^{\tfrac{3}{4}}{C}_{\gamma }}{{C}_{E}(1-2s)}.\end{eqnarray} \tag{ 13 }$$

Result—Combining our results, each phase in Equation (6) contributes a term proportional to ${{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}$, giving

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}\approx {{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}{C}_{\gamma }\left(\displaystyle \frac{2\sqrt{2}\alpha }{3}{\left(\displaystyle \frac{{\rm{\Delta }}{t}_{\mathrm{mix}}}{{t}_{\mathrm{dyn}}}\right)}^{\tfrac{3}{2}}\right.\\ \qquad \left.+\,{\left(2\alpha \right)}^{\tfrac{3}{4}}\displaystyle \frac{{\rm{\Delta }}{t}_{\mathrm{sat}}}{{t}_{\mathrm{dyn}}}+\displaystyle \frac{4{\left(2\alpha \right)}^{\tfrac{3}{4}}{C}_{\gamma }}{{C}_{E}(1-2s)}\right),\end{array}\end{eqnarray} \tag{ 14 }$$

Since both Δt_mix ∼ t_dyn from freefall scaling and Δt_sat ∼ t_dyn from dimensional constraints, all the terms in the parenthesis are constants and thus:

$$\begin{eqnarray}&&{\rm{\Delta }}\sim {{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}\sim \displaystyle \frac{{A}^{\tfrac{1}{4}}{g}^{\tfrac{1}{4}}{L}^{\tfrac{3}{4}}}{{\nu }^{\tfrac{1}{2}}},\end{eqnarray} \tag{ 15 }$$

independent of constants in the models.

In summary, the clean result comes from the choice of freefall velocity (u ∼ u_sat) and length (l_i ∼ L) scales leading to a dynamo with growth rate $\gamma \sim {t}_{\mathrm{dyn}}^{-1}{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}$ acting over dynamical timescales t_dyn and resulting in amplification of ${\rm{\Delta }}\sim \gamma {t}_{\mathrm{dyn}}\,\sim {{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}$.

Low Pm case—The RTI is also applicable to astrophysical scenarios where the magnetic Prandtl number Pm = Rm/Re can be less than one, such as in stellar interiors during supernova explosions. The only modification is the estimate of the growth rate. The resistive scale at low Pm is inside the inertial range (l_η ∼ Rm^−3/4 L > l_ν ) and the dynamo is thought to be driven by resistive-scale eddies whose turnover times are ${t}_{\eta }\sim L/({{URm}}^{\tfrac{1}{2}})$ (Iskakov et al. 2007). For ${Re}\,\gg \,{{Re}}^{c}({Pm})$, where ${{Re}}^{c}({Pm}\,\ll 1)={ \mathcal O }({10}^{2})$, this corresponds to a reduced growth rate given by $\gamma \sim {{URm}}^{\tfrac{1}{2}}/L$ and an exponential amplification factor that instead scales as

$$\begin{eqnarray}&&{\rm{\Delta }}\sim {{Rm}}_{\mathrm{sat}}^{\tfrac{1}{2}}\sim \displaystyle \frac{{A}^{\tfrac{1}{4}}{g}^{\tfrac{1}{4}}{L}^{\tfrac{3}{4}}}{{\eta }^{\tfrac{1}{2}}}.\end{eqnarray} \tag{ 16 }$$

Saturation of the SSD—We discuss the extension of our model of the dynamo in the kinematic regime to the dynamical regime and saturation of the SSD at high Pm. Our results above are only applicable in the kinematic regime of the dynamo where Lorentz forces are unimportant and the induction equation is a linear function of the velocity field, which leads to exponential growth of magnetic energy. However, when the magnetic energy becomes comparable to the kinetic energy at viscous scales ($\langle | {\boldsymbol{B}}{| }^{2}\rangle \sim {{Re}}^{-\tfrac{1}{2}}\langle | {\boldsymbol{u}}{| }^{2}\rangle$), the dynamo enters the dynamical regime where feedback on fluid motions from Lorentz forces become important and the dynamo switches to polynomial-order growth (see Rincon 2019 and references therein). The dynamical regime is also known as the nonlinear growth phase, for which there are several models in the literature. In one possible model, the growth rate decreases due to sequential suppression of dynamo-generating motions by the Lorentz forces starting from the viscous scales and ending at the forcing scales (Schekochihin et al. 2002). Assuming the growth rate is set by the turnover time of the smallest, unquenched velocity scale and that the magnetic energy is in equipartition with all smaller velocity scales, one can easily show that the magnetic energy will grow linearly in time ME(t) ≈ ζ t, where ∼ u³/l_i is the transfer rate of kinetic energy and ζ is a dimensionless constant. When the Lorentz forces become strong at the forcing scales, the dynamo will fully saturate in near-equipartition between the total magnetic and kinetic energy with a ratio ${ME}/{KE}=f={ \mathcal O }({10}^{-1})$, where f is a dimensionless constant .

Thus, our results for Δ are an upper bound because the SSD will leave the kinematic regime during some phase of the RTI if Δ is too large. As a rough estimate, any sub-equipartition, seed magnetic field configuration with magnetic energy greater than ${ME}(t=0)\gtrsim {{Re}}_{\mathrm{sat}}^{-\tfrac{1}{2}}{{KE}}_{\mathrm{sat}}{e}^{-{\rm{\Delta }}}$ will reach the dynamical regime at some point in the RTI evolution and possibly saturate before the end of the decay phase of the RTI, which would result in MHD turbulence in the remaining RTI evolution. We briefly study this regime in Section 3.3.

While in principle one can build on our kinematic dynamo model and include the nonlinear dynamo growth phase by using a time dependent (t), we do not pursue this idea in further detail because (1) the nonlinear growth phase in simulations is short and difficult to test and (2) it is unclear how the feedback from the strong Lorentz forces across a growing range of velocity scales will affect the time evolution of the RTI velocity field. However, we do note a useful estimate for the time it takes for the dynamo to saturate during the dynamical regime Δt_d.r.. If the magnetic energy at the start of the dynamical regime is ${ME}\sim {{Re}}^{-\tfrac{1}{2}}{KE}$, at saturation is ME = f · KE, and the transfer rate of kinetic energy is roughly constant ∼ KE/t_dyn, then we simply find that Δt_d.r. is comparable to the dynamical time ${\rm{\Delta }}{t}_{d.r.}\sim {t}_{\mathrm{dyn}}(f-{{Re}}^{-\tfrac{1}{2}})/\zeta$.

In this section, we use three-dimensional direct numerical simulations of the RTI to test the theoretical scaling predictions presented in Section 2. We show and discuss the results of parameter scans with the Atwood number, gravitational acceleration, and viscosity.

3.1. Numerical Setup

We use the Athena++ code (Stone et al. 2019) to solve the visco-resistive MHD equations in a 3D Cartesian domain with a similar setup to previous works (Stone & Gardiner 2007). The horizontal −L/2 ≤ x, y ≤ L/2 directions have periodic boundary conditions while the vertical direction −L ≤ z ≤ L has reflecting boundary conditions, where L=0.1. All simulations use a resolution of N_x N_y N_z = 512² × 1024, RK3 for the timestepper, and HLLD for the Riemann solver (Miyoshi & Kusano 2005).

The initial hydrodynamic conditions set the fluid density ρ(z) = ρ_t in the top half of the domain, ρ(z) = ρ_b in the lower half of the domain, and a seed vertical velocity perturbation given by:

$$\begin{eqnarray}&&\begin{array}{l}{u}_{z}(t=0)={\boldsymbol{ \mathcal R }}\\ \quad \left[\sum _{{n}_{x},{n}_{y}}\displaystyle \frac{{\widetilde{a}}_{{n}_{x},{n}_{y}}}{\sqrt{n}}{e}^{i\tfrac{2\pi }{L}\left({n}_{x}x+{n}_{y}y\right)}\cos \left(\displaystyle \frac{\pi z}{2L}\right)\right],\end{array}\end{eqnarray} \tag{ 17 }$$

where ${\widetilde{a}}_{{n}_{x},{n}_{y}}$ is a random complex number, $n=\sqrt{{n}_{x}^{2}+{n}_{y}^{2}}$, and $0\lt n\leqslant {n}_{\max }$ with ${n}_{\max }=32$. The magnetic field is initialized with an isotropic spectrum $M(k)\approx \mathrm{const}.$ in the wavenumber range $0\lt k\leqslant {k}_{\max }=2\pi {n}_{\max }/L$. The total initial magnetic energy E_M (t = 0) = ∫M(k)dk is set to a dynamically insignificant value E_M (t = 0) ≈ 10⁻¹⁹ when studying the kinematic regime in Section 3.2, and set higher when studying of the saturation regime. Last, we fix the magnetic Prandtl number to Pm = 3 and the thermal Prandtl number to Pr = ν/κ = 1 for all runs.

Asymptotic scaling laws of hydrodynamic quantities in the mixing phase of the RTI are known to be sensitive to initial conditions (IC) and the box aspect ratio (Dimonte 2004; Boffetta & Mazzino 2017). Thus, parameter scans with different choices of IC (e.g., varying n_max) could lead to slightly different scaling laws. It is not clear which choice of IC is generally most physically applicable, since the RTI setup is already quite idealized. We discuss throughout this article how our dynamo results may vary for choices of IC different from ours. One common alternative option is mode perturbations with wavelengths that go down to the grid scale, while another is restricting to a shell of modes in Fourier space (Dimonte 2004). Our choice of IC was motivated to allow us to study the effect of changing viscosity on the dynamo while retaining a similar hydrodynamic RTI evolution by having our fastest-growing mode always be fixed at ${k}_{\max }\lt {k}_{c}$ for our chosen range of viscosities.

In order to expect our results to be generalizable, it is important to at least check that minor variations of ICs for a fixed set of parameters have a minimal effect on Δ. We have found this to generally be the case for a fiducial set of parameters where we widely varied ${n}_{\max }\geqslant 8$ for both the initial velocity and magnetic field spectra. The two exceptions we found are the edge cases of a single mode RTI (${n}_{\max }=1$) or a purely uniform weak initial magnetic field. The single mode RTI becomes nonlinear near the system scale, unlike the multi-mode RTI, leading to negligible contribution from Δ^mix and a reduced Δ. In the uniform magnetic field case, the dynamo additionally has a short and intense transient growth at the beginning of the mixing phase, due to coherent alignment of the field with the secondary shear layers that undergo the Kelvin–Helmholtz instability. We expect neither edge case to be applicable in a realistic astrophysical system.

The parameter scans are based on a fiducial simulation with parameters {A, g, ν⁻¹} = {0.67, 0.65, 3 × 10⁵}, which we use below to describe a typical simulation and explain our method of analysis. Figure 1 shows the time evolution of the kinetic energy, magnetic energy, dynamo growth rate, and approximate Reynolds number of the fiducial simulation. When the Reynolds number rises above the critical Reynolds number ${Re}(t)\gtrsim {{Re}}^{c}({Pm}\gt 1)\approx 60$, the dynamo growth rate sharply increases. The dynamo is triggered near the beginning of the mixing phase (green shaded region in Figure 1), which we define as the time when KE(t) = 0.01PE^mix, where PE^mix = 0.5gL⁴(ρ_t − ρ_b ) is the potential energy released from complete mixing of the fluids. The kinetic energy rises until a maximum in the middle of the saturation phase at the saturation time ${t}_{\mathrm{sat}}^{(n)}$ given by ${KE}({t}_{\mathrm{sat}}^{(n)})=\max ({KE}(t))$, at which point we also have the numerical values for ${u}_{\mathrm{sat}}^{(n)}=u({t}_{\mathrm{sat}}^{(n)})$ and ${t}_{\mathrm{dyn}}^{(n)}=L/{u}_{\mathrm{sat}}^{(n)}$.

Zoom In Zoom Out Reset image size

We define the mixing phase to end and the saturation phase (orange shaded region) to begin at time $t={t}_{\mathrm{sat}}^{(n)}-{t}_{\mathrm{dyn}}^{(n)}$. Similarly, the saturation phase we define to end and the decay phase (red shaded region) to begin at time $t={t}_{\mathrm{sat}}^{(n)}+{t}_{\mathrm{dyn}}^{(n)}$. At some point in the decay phase, the magnetic energy will reach its maximum and we define the decay phase to end for the purposes of the dynamo. Note that, while the exact definitions of the start and end time of each phase are somewhat arbitrary, we find that our results are not sensitive to reasonable variations of the definitions.

With the above definitions, we obtain the numerical values of the mean dynamo growth rate and amplification factor in each phase (${\overline{\gamma }}_{\mathrm{phase}}$ and Δ^phase) with the following expressions:

$$\begin{eqnarray}&&{\overline{\gamma }}_{\mathrm{phase}}=\displaystyle \frac{1}{{\rm{\Delta }}{t}_{\mathrm{phase}}}{\int }_{{t}_{\mathrm{beg}}}^{{t}_{\mathrm{end}}}\gamma (t){dt},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{\mathrm{phase}}=\mathrm{ln}\left(\displaystyle \frac{{ME}({t}_{\mathrm{end}})}{{ME}({t}_{\mathrm{beg}})}\right)={\overline{\gamma }}_{\mathrm{phase}}{\rm{\Delta }}{t}_{\mathrm{phase}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{\mathrm{phase}}={t}_{\mathrm{end}}-{t}_{\mathrm{beg}},\end{eqnarray} \tag{ 20 }$$

where t_beg and t_end are the beginning and ending times of each phase. The total amplification factor is then very closely given by:

$$\begin{eqnarray}&&{\rm{\Delta }}={{\rm{\Delta }}}^{\mathrm{mix}}+{{\rm{\Delta }}}^{\mathrm{sat}}+{{\rm{\Delta }}}^{\mathrm{dec}}\approx \mathrm{ln}\left(\displaystyle \frac{\max ({ME}(t))}{\min ({ME}(t))}\right).\end{eqnarray} \tag{ 21 }$$

For completeness, two-dimensional slices from the fiducial simulation of the density ρ and vertical magnetic field B^z during each of the three turbulent phases are shown in Figure 2. The amplified magnetic field in the mixing phase (left column) closely follows the rising bubbles and falling spike structures at intermediate scales visible in the density field. The saturation phase (middle column) has large eddies on the system scale with the magnetic field developing on small scales, as would be qualitatively expected of isotropic turbulence. Last, the decay phase (right column) shows a negative mean density gradient in which residual kinetic energy sloshes fluid around while the magnetic fields still appear on scales that are small but slightly larger than in the saturation phase.

Zoom In Zoom Out Reset image size

3.2. Scaling Results

To test the scaling relations of the dynamo model proposed in Section 2, we perform a series of parameter scans of the Atwood number, gravitational acceleration, and viscosity. We choose our fiducial simulation with values of {A, g, ν⁻¹} = {0.67, 0.65, 3e5} (corresponding to $\max ({Re}(t))\approx 300$) and then run a parameter scan across a maximum range of each parameter, while keeping all others fixed. For the gravitation acceleration parameter, we scan g ∈ {0.25, 0.35, 0.5, 0.65, 0.8, 1.0, 1.2, 1.4}, for the Atwood number, we scan A ∈ {0.2, 0.33, 0.43, 0.5, 0.6, 0.67, 0.71, 0.76, 0.82, 0.88, 0.9}, and for the viscosity, we scan ν⁻¹ ∈ {2, 2.5, 3, 3.5, 4, 4.5, 5} × 10⁵. Time evolution of the magnetic energies and dynamo growth rates for representative subsets of each parameter scan are shown in Figure 3. The lower bound for each parameter is constrained by needing ${{Re}}_{\mathrm{sat}}\,\gg \,{{Re}}^{c}$ in order for the dynamo growth rate to be approximated by its asymptotic power law $\gamma \sim {{Re}}^{\tfrac{1}{2}}$, while the upper bound is constrained by requiring the Kolmogorov scale to be larger than the grid scale. In particular, the lowest value of the explicit viscosity is chosen to be a factor of two above the numerical viscosity for the numerical setup, which is estimated based on the decay rate of Alfvén waves as described in Appendix. Note that it is critical to at least marginally resolve the Kolmogorov scale, because otherwise the arbitrary simulation grid scale introduces another length scale that breaks the scaling arguments.

Zoom In Zoom Out Reset image size

First, we examine the dynamo amplification factor across the entire RTI instability. A linear fit in log-log space between Δ and each RTI parameter gives:

$$\begin{eqnarray}&&{\rm{\Delta }}\sim {A}^{0.40\pm 0.03}{g}^{0.35\pm 0.03}{\nu }^{-0.51\pm 0.09},\end{eqnarray} \tag{ 22 }$$

For reference, the predicted freefall scaling relations are given by:

$$\begin{eqnarray}&&{\rm{\Delta }}\sim {A}^{0.25}{g}^{0.25}{\nu }^{-0.5}.\end{eqnarray} \tag{ 23 }$$

While the viscosity exponent is in good agreement with the model (within one standard deviation), the exponents of A and g are close in magnitude to the predicted value of 0.25, but are different by a statistically significantly amount (roughly five and three standard deviations). An alternative way to view the result is shown in the top panel of Figure 4 where the Δ is plotted versus ${{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}$ on a log-log scale, so a linear fit $\mathrm{ln}{\rm{\Delta }}=m\mathrm{ln}{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}+b$ with a slope of m = 1 means perfect agreement. The data do appear to qualitatively follow a power law; however, the fit also shows a minor but statistically significant disagreement with the model quantified by the measured value of m = 1.29 ± 0.08. In other words, with a measured value of b = −0.90, the simulations have a scaling ${\rm{\Delta }}\approx 0.4{{Re}}_{\mathrm{sat}}^{0.65\pm 0.04}$ instead of the predicted ${\rm{\Delta }}\sim {{Re}}_{\mathrm{sat}}^{0.5}$. Breaking up Δ further into the contribution of each phase Δ^phase in the bottom panel of Figure 4 clearly shows that the mixing phase is in disagreement while the saturation and decay phases are in good agreement with the model (although the decay phase data are highly scattered).

Zoom In Zoom Out Reset image size

To better understand the discrepancy, it is necessary to examine the scaling relations of the mean growth rate and duration of each phase with the RTI parameters. Following a linear fit in log-log space between $\{{{\rm{\Delta }}}^{\mathrm{phase}},{\overline{\gamma }}_{\mathrm{phase}},{\rm{\Delta }}{t}_{\mathrm{phase}}\}$ and {A, g, ν} in each phase separately, Table 1 shows the resulting scaling exponents. Overall, the freefall model predictions are again in good agreement for the saturation phase (all roughly within two standard deviations), while they show several disagreements for the mixing and decay phases. The deviation of the m = 1.29 ± 0.08 result for Δ from the model prediction (m = 1) can thus be explained by a likely failure of model assumptions in the mixing and decay phases. The value of m = 1.29 is still close to 1, however, because Δ_sat provides the dominant contribution to Δ. We now present a more detailed analysis of the mixing and decay phases.

Table 1. Table of the Scaling Exponents between the Dynamo Amplification Factor, Mean Dynamo Growth Rate, and Duration of Each Phase (rows) and the RTI Parameters (columns)

	A	g	ν
Δff	0.25	0.25	−0.5
Δmix	0.68 ± 0.16	0.43 ± 0.10	−0.79 ± 0.23
Δsat	0.25 ± 0.03	0.24 ± 0.02	−0.38 ± 0.04
Δdec	0.28 ± 0.23	0.37 ± 0.13	−0.35 ± 0.33
γff	0.75	0.75	−0.5
${\overline{\gamma }}_{\mathrm{mix}}$	1.01 ± 0.05	0.77 ± 0.03	−0.54 ± 0.07
${\overline{\gamma }}_{\mathrm{sat}}$	0.81 ± 0.07	0.81 ± 0.05	−0.60 ± 0.10
${\overline{\gamma }}_{\mathrm{dec}}$	0.17 ± 0.19	1.0 ± 0.27	−1.3 ± 0.36
Δtff	−0.5	−0.5	0
Δtmix	−0.33 ± 0.12	−0.34 ± 0.08	−0.24 ± 0.17
Δtsat	−0.56 ± 0.06	−0.57 ± 0.04	0.21 ± 0.10
Δtdec	0.1 ± 0.29	−0.62 ± 0.38	0.96 ± 0.41

Note. For example, the entry for Δ_mix and A means Δ_mix ∼ A^0.68±0.16. For reference, Δ_ff, γ_ff, Δt_ff denote the prediction for the exponent based on the freefall model.

Download table as:  ASCIITypeset image

Mixing phase—The goal is to find which assumptions in the model of Section 2 are violated in the mixing phase. We begin by comparing the time dependence of KE(t) and γ_mix(t) between freefall model predictions and rescaled simulation time series averaged over all runs, shown in Figure 5. The freefall model with the isotropy assumption asymptotically predicts a time dependence ${\gamma }_{\mathrm{mix}}\sim \overline{u}{\left(t\right)}^{1.5}/{l}_{i}{\left(t\right)}^{0.5}\sim {t}^{0.5}$ for the dynamo growth rate and ${KE}(t)\sim \overline{u}{\left(t\right)}^{2}h(t)\sim {t}^{4}$ for the total kinetic energy with $\overline{u}\sim t$ and l_i (t) ∼ h(t) ∼ t², where

$$\begin{eqnarray}&&\overline{u}{\left(t\right)}^{2}=\displaystyle \frac{{\int }_{-h(t)}^{h(t)}{\int }_{-L/2}^{L/2}{\int }_{-L/2}^{L/2}u{\left({\boldsymbol{x}},t\right)}^{2}{d}^{3}x}{2h(t)}.\end{eqnarray} \tag{ 24 }$$

Zoom In Zoom Out Reset image size

The numerical dynamo growth rate in the simulation data is slightly steeper, ${\gamma }_{\mathrm{mix}}^{(n)}\sim {t}^{0.71\pm 0.19}$ (within two standard deviations from the prediction), and the kinetic energy growth significantly shallower, KE⁽ⁿ⁾ ∼ t^3.06±0.19 (four standard deviations from the prediction). The discrepancy suggests a slightly alternative scaling of either $\overline{u}(t)$, h(t), or l_i (t) in our simulations. Fitting the time power laws to $\overline{u}(t)$ and h(t) independently for the fiducial simulation (not shown), we find the numerical time dependence ${\overline{u}}^{(n)}\sim {t}^{0.8\pm 0.025}$ and h⁽ⁿ⁾ ∼ t^1.55±0.06, which suggests an intermediate scaling between the freefall and terminal velocity scaling ($\overline{u}\sim {t}^{0}$, h ∼ t¹) and has also been observed in previous simulations (Lecoanet et al. 2012). One possible reason we do not obtain freefall scaling may be our choice of initial conditions, which do not set ${k}_{\max }={k}_{c}$, while another reason may be that the vertical dimension of the box is not asymptotically large enough to reduce the effect of the linear term in Equation (7). The effects of thermal diffusion may also play a role, despite the moderately high resolution of our simulations. The deviation of Δt_mix from freefall predictions in Table 1 is likely tied to the above reasons as well.

Nonetheless, the time dependence in the simulations predicts ${\gamma }_{\mathrm{mix}}\sim {{\overline{u}}^{(n)}}^{1.5}/{{h}^{(n)}}^{0.5}\sim {t}^{0.45\pm 0.05}$, which is even less steep than the freefall prediction. This motivates checking which of the remaining major assumptions break down: (1) isotropy of the dynamo-generating scales or (2) whether the integral scale and mixing height are directly proportional. The assumption of isotropy of turbulence at dynamo-generating scales appears to be supported by simulation data since ${D}_{\gamma }=d\mathrm{ln}{\overline{\gamma }}_{\mathrm{mix}}/(d\mathrm{ln}\nu )\approx 0.5$ in Table 1. Additionally, while the anisotropy of the total kinetic energy (defined by 2KE_z (t)/KE_h (t)) is large in the mixing phase, as shown in the top panel of Figure 6, the anisotropy in the total magnetic energy (defined by 2ME_z (t)/ME_h (t)) is much smaller and approaches unity by the end of the mixing phase (bottom panel of Figure 6). We expect the anisotropy of the magnetic field in the mixing phase to be even lower for larger Reynolds numbers because higher-resolution simulations of the RTI than ours find strong evidence for isotropy at Kolmogorov scales (Zingale et al. 2005; Cabot & Cook 2006).

Zoom In Zoom Out Reset image size

We check the second assumption by computing l_i (t) = (∫k⁻¹ E(k, t)dk)/(∫E(k, t)dk) for the fiducial simulation and find ${l}_{i}^{(n)}\sim {t}^{1.1\pm 0.04}$, which is statistically significantly different than the time dependence for h⁽ⁿ⁾ ∼ t^1.55±0.06. Notably, the direct substitution ${\gamma }_{\mathrm{mix}}\sim {{\overline{u}}^{(n)}}^{\tfrac{3}{2}}/{{l}_{i}^{(n)}}^{\tfrac{1}{2}}\sim {t}^{0.65\pm 0.04}$ has a better agreement with the observed ${\gamma }_{\mathrm{mix}}^{(n)}(t)\sim {t}^{0.71\pm 0.19}$. However, the large uncertainty of the time dependence of ${\gamma }_{\mathrm{mix}}^{(n)}$ makes any conclusive determination difficult.

Overall, a detailed understanding of the dynamo in the mixing phase of the RTI requires a further systematic study of the effects of initial conditions, box aspect ratio, and diffusivities, which we leave for future work.

We speculate that there are two relatively important modifications to the model not examined in this study: (1) accounting for the non-steady-state nature of the mixing phase turbulence, and (2) incorporating effects of large-scale anisotropy in the velocity field into the dynamo growth rate. The first modification requires including the time delay between the buoyancy forcing at the scale l_i (t) and the dissipation rate (t), which is dissipating cascading energy due to forcing from earlier times (Livescu et al. 2009). This is important because the dynamo operates at the dissipation scale and can more accurately be expressed as $\gamma (t)\sim {\left(\epsilon (t)/\nu \right)}^{\tfrac{1}{2}}$ (Beresnyak 2012), with $\epsilon (t)\sim u{\left(t\right)}^{3}/{l}_{i}(t)$ a good approximation only if the cascade rate is faster than the time rate of change of l_i (t) and u(t).

The second modification may require incorporating the effect of large-scale velocity anisotropy in the dynamo growth rate. As shown in Figure 7, the velocity spectra in the mixing phase are highly anisotropic at large scales with a dominant vertical component and become quasi-isotropic at intermediate and smaller scales (for kL/2π ≳ 10) while the magnetic field component spectra maintain roughly the same level of quasi-isotropy at all scales. Both of these patterns are observed in stably stratified turbulence, with the only difference being that the horizontal velocity components dominate at large scales instead (Skoutnev et al. 2021). Thus, the effect of the large-scale anisotropy on the dynamo in the mixing phase can perhaps be modeled as a reduction in the effective Reynolds number through an effective Froude number, Fr < 1, similar to the way the buoyancy Reynolds number, ${Rb}={{Fr}}^{2}{Re}$, controls the dynamo growth rate in stably stratified turbulence (Skoutnev et al. 2021). This extension may help explain the significant discrepancy of the exponent relating ${\overline{\gamma }}_{\mathrm{mix}}$ and the Atwood number in Table 1, for instance.

Zoom In Zoom Out Reset image size

Decay Phase—The scatter in the data and uncertainty in the scaling exponents are unfortunately large in the decay phase. The uncertainty for some exponents in Table 1 is comparable to the mean values, making it impossible to meaningfully compare with asymptotic predictions of the model. We attribute this to intermittency in the RTI turbulence, where occasionally an intermediate-scale, heavy parcel remains suspended for longer than usual and causes residual large-scale forcing, which can be seen as brief increases of KE(t) and 2KE_z (t)/KE_h (t) in the decay phase of some runs in Figures 5 and 6, respectively. The main assumption that the turbulence in the decay phase is freely decaying may perhaps be violated. It is not obvious whether residual forcing effects will persist in the decay phase at even higher Reynolds numbers. Additionally, the large uncertainty in ${\overline{\gamma }}_{\mathrm{dec}}$ may also be a finite Reynolds number effect, since the rapid decay of Re(t) from the moderate value of Re in our simulations may violate the assumptions that ${Re}\,\gg \,{{Re}}^{c}$ and Δt_dec ≫ t_dyn/C_E .

The asymptotic time dependence of the dynamo growth rate and kinetic energy similarly has large uncertainties. We find that l_i (t) ≈ t⁰ with l_i ≈ 0.3L throughout the entire decay phase, suggesting s = 0 as expected since the turbulence has reached the box scale. This predicts an asymptotic scaling $\overline{u}(t)\sim t{{\prime} }^{-1}$, ${KE}(t)\sim t{{\prime} }^{-2}$, and $\gamma (t)\sim t{{\prime} }^{-1.5}$, where $t^{\prime} =0$ is the beginning of the decay phase. The numerical power-law fit (Figure 5) gives ${{KE}}^{(n)}\sim t{{\prime} }^{-4.73\pm 1.27}$ and ${\gamma }^{(n)}\sim t{{\prime} }^{-4.86\pm 1.6}$, which are both significantly steeper than the model prediction. The behavior of decaying RTI turbulence and associated dynamo is clearly poorly approximated by the model of freely decaying, isotropic turbulence.

Fortunately, the decay phase has a smaller contribution to the total dynamo amplification factor (see bottom panel of Figure 4) than either the mixing or saturation phases, which mitigates the effects of the large uncertainties. Reducing the uncertainties to better understand the decay phase will likely require simulations with much higher Reynolds numbers.

3.3. Saturation

In this section, we study the case where the dynamo is able to saturate before the RTI fully relaxes. We simulate this regime by initializing the magnetic energy with ME(t = 0) ≈ 10⁻⁵ PE^mix for the fiducial run, which satisfies the ${ME}(t=0)\geqslant {{Re}}_{\mathrm{sat}}^{-\tfrac{1}{2}}{{KE}}_{\mathrm{sat}}{e}^{-{\rm{\Delta }}}={ \mathcal O }({10}^{-7})\cdot {{PE}}^{\mathrm{mix}}$ criterion for the dynamo to reach the dynamical regime as discussed in Section 2. The energy evolution is shown in the inset of Figure 8, where the dynamical regime is observed to begin around t ≈ 4 when the slope of the magnetic energy growth sharply drops. The following nonlinear growth phase appears to last only briefly until t ≈ 5 (t_dyn ≈ 0.5 for the fiducial simulation). This is expected because the ${{Re}}_{\mathrm{sat}}^{-\tfrac{1}{2}}\sim 0.03$ for the fiducial simulation and the final ratio of magnetic to kinetic energy is f ≈ 0.1, leading to a short predicted time to saturation in the dynamical regime of ${\rm{\Delta }}{t}_{d.r.}\sim {t}_{\mathrm{dyn}}\left(f-{{Re}}_{\mathrm{sat}}^{-\tfrac{1}{2}}\right)$.

Zoom In Zoom Out Reset image size

To qualitatively test the phenomenological model of the dynamo in the dynamical regime discussed in Section 2, we plot the isotropic energy spectra of a cubic volume in the region −L/2 ≤ z ≤ L/2 at two representative times t ∈ {4, 10} in Figure 8. At t = 4 (red curves), just after the dynamo has enter the nonlinear growth phase, the magnetic energy is larger on a scale-by-scale basis below an intermediate scale kL/2π ∼ 30. In the fully saturated phase of the dynamo at t = 10 (blue curves), the magnetic energy is larger than the kinetic energy at basically all but the largest scales. Both of these observations are in accord with the predictions from the current understanding of the dynamical dynamo regime in steady-state turbulence as discussed in Section 2.

Overall, these results show that the RTI is capable of generating small-scale magnetic fields in near-equipartition (${ME}/{{KE}}^{\mathrm{dec}}=f={ \mathcal O }({10}^{-1})$) with the decaying turbulent kinetic energy, KE^dec , which is a fraction on the order of ${{KE}}^{\mathrm{dec}}/{{PE}}^{\mathrm{mix}}={ \mathcal O }({10}^{-2})$ of the initially available potential energy. These magnetic fields then can act as seed fields for further processes driven by larger scales in the astrophysical object.

One open question in recent studies of binary neutron star mergers is the source, efficiency, and time scale of magnetic field amplification in the post-merger star (Price & Rosswog 2006; Giacomazzo et al. 2015; Kiuchi et al. 2015). Several studies have proposed the role of the Kelvin–Helmholtz instability (KHI) active at the contact region of the two stars in contributing to the magnetic field amplification observed in the core of the post-merger (Kiuchi et al. 2015; Aguilera-Miret et al. 2020). It is generally assumed that the Reynolds number of the KHI turbulence is high enough that the dynamo will saturate in near-equipartition with the kinetic turbulence, although it is unclear how fast that will happen (Kiuchi et al. 2018). The KHI, however, cannot explain the large amount of dynamo action that is observed in the surface layer during the merger (Kiuchi et al. 2015). We propose that the RTI is likely responsible for dynamo action in the outer parts of the merger where centrifuged denser matter can be seen falling down onto less dense matter at adjacent longitudes. Rapid saturation of the dynamo and the associated amplification of the magnetic field to near-equipartition with the RTI-driven turbulence might critically affect the prospects of long-term mass ejection (Metzger et al. 2018; Ciolfi & Kalinani 2020), as well as jet launching (Ciolfi 2020; Mösta et al. 2020) from a stable magnetar remnant. Using our model, we verify that the RTI-driven dynamo does saturate, and we provide a quantitative estimate for the time scale for saturation in the conditions of the post-merger neutron star envelope.

We first need an estimate for the Reynolds number Re_sat. The viscosity of neutron star matter is strongly dependent on the temperature and density conditions. Assuming low temperatures of T ≃ 1 MeV and densities around nuclear saturation, i.e., 2 × 10¹⁴ g cm⁻³, the kinematic shear viscosity ν ∼ 3 × 10⁻³ m² s⁻¹ in the electron-scattering-dominated regime (Shternin & Yakovlev 2008). On the other hand, for temperatures T ≃ 10 MeV, neutrino-emitting Urca processes (e.g., $n\to p\,{e}^{-}\,{\bar{\nu }}_{e}$ and p e⁻ → n ν_e ) might be the dominant contribution, leading to ν ≃ 10⁴ m² s⁻¹ at densities around saturation (Alford et al. 2018). The resistivity for a warm (T ≃ 1 MeV) neutron star crust is set by electron scattering of correlated nuclei (Harutyunyan & Sedrakian 2016) and has a value of η ≃ 5 × 10⁻⁷m² s⁻¹ (Harutyunyan et al. 2018), which puts neutron star matter well into the high Pm regime with ${Pm}={ \mathcal O }({10}^{4-11})$. The extreme density gradients at the surface of a neutron star, imply that the Atwood number A ≈ 1. We consider a thin layer exhibiting these gradients and becoming RTI-unstable on a length scale that is a fraction of a scale height L ∼ 0.1H (H ∼ 1 km). We approximate the gravitational acceleration by means of the Newtonian expression for surface gravity, i.e., $g\simeq \tfrac{{GM}}{{R}^{2}}\approx 3\times {10}^{16}\,{\rm{m}}\,{{\rm{s}}}^{-2}$. Although u_sat as obtained from the expressions above technically would be superluminal (owing to the simplified Newtonian assumptions made here), we take this as an indication that u_sat ≃ c. Finally, all these parameters correspond to an estimate ${{Re}}_{\mathrm{sat}}\sim { \mathcal O }({10}^{6}-{10}^{13})$. If we assume our model holds (${\rm{\Delta }}\approx {C}_{{\rm{\Delta }}}{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}$ with ${C}_{{\rm{\Delta }}}={ \mathcal O }(1)$), then in an RTI-unstable region, the expected amplification would be ${\rm{\Delta }}\sim { \mathcal O }({10}^{3}-{10}^{6})$. The assumption that the dynamo will reach equipartition with the kinetic energy of the viscous scales, ${KE}/{{Re}}^{\tfrac{1}{2}}$, is easily justified since ${e}^{{\rm{\Delta }}}\,\gg {KE}/({ME}(0){{Re}}^{\tfrac{1}{2}})$ even for generous estimates of order ${KE}/{ME}(0)\sim {10}^{{ \mathcal O }(10)}$, where ME(0) is the initial magnetic energy density near the surface of one of the initial neutron stars and KE is the characteristic turbulent kinetic energy in post-merger envelope.

We can now estimate how long it will take the dynamo to leave the kinematic regime (Δt_k.r.) using the dynamo growth rate $\gamma \sim {t}_{\mathrm{dyn}}^{-1}{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}\,\sim \,{ \mathcal O }({10}^{3}-{10}^{7})\mu {s}^{-1}$, where we have used t_dyn = L/u_sat ∼ 0.3 μ s for the dynamical time. The estimate is as follows:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}{t}_{k.r.}\sim \displaystyle \frac{1}{\gamma }\mathrm{ln}\left(\displaystyle \frac{{KE}}{{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}{ME}(0)}\right)\\ \,\sim \,2{\mathrm{log}}_{10}\left(\displaystyle \frac{{KE}}{{{Re}}_{\mathrm{sat}}^{\tfrac{1}{2}}{ME}(0)}\right)\\ \,\times \,\left[{10}^{-7}-{10}^{-3}\right]\mu {\rm{s}}\ll {t}_{\mathrm{dyn}}.\end{array}\end{eqnarray} \tag{ 25 }$$

After the kinematic regime, the dynamo will continue to grow in the dynamical regime and fully saturate on a timescale of

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{t}_{d.r.} & \sim & {t}_{\mathrm{dyn}}(f-{{Re}}_{\mathrm{sat}}^{-\tfrac{1}{2}})/\zeta \\ & \sim & {t}_{\mathrm{dyn}}f/\zeta \\ & \sim & { \mathcal O }({10}^{-1})\mu s\end{array}\end{eqnarray} \tag{ 26 }$$

using the model of the nonlinear dynamo growth phase described in Section 2 and assuming $f={ \mathcal O }({10}^{-1})\gg {{Re}}_{\mathrm{sat}}^{-\tfrac{1}{2}}$.

Since the duration of the kinematic dynamo regime (upper bound of nanoseconds) and duration of the dynamical dynamo regime (upper bound of microseconds) of the RTI-driven turbulence in the envelope are both much smaller than the relaxation time of the merger (order of milliseconds), we expect magnetic energies to be in near-equipartition with the kinetic energy of the RTI-driven turbulence in the envelope across essentially the entire duration of the post-merger evolution.

We present a model for the kinematic small-scale dynamo in the mixing, saturation, and decay phases of the Rayleigh–Taylor instability of an ionized, collisional plasma with freefall and isotropy assumptions for the turbulence in each of the phases. The model quantitatively predicts scaling relations between the properties of the dynamo (growth rate and total magnetic energy exponential amplification factor) and parameters of the RTI (Atwood number, gravitational acceleration, length scale, and viscosity). The model predictions are tested with sets of three-dimensional direct numerical simulations that solve the visco-resistive MHD equations using the Athena++ code. The main results are itemized below:

• 1.  
  We find that the total magnetic energy exponential amplification factor, Δ, based on simulation data scales as

  $$\begin{eqnarray}&&{\rm{\Delta }}\approx {C}_{{\rm{\Delta }}}{{Re}}_{\mathrm{sat}}^{{D}_{\gamma }},\end{eqnarray} \tag{ 27 }$$

  with constants C_Δ ≈ 0.4 and D_γ ≈ 0.65 ± 0.04. This is in fairly close agreement with the model prediction of D_γ = 0.5, but the difference is statistically significant. An analysis of the dynamo in each phase reveals that the model correctly predicts scaling relations in the saturation phase, while having several discrepancies in the mixing and decay phases.
• 2.  
  An analysis of the dynamo scaling relations and time dependence in the mixing phase shows that the freefall and isotropy assumptions are fairly well-supported, with a few minor discrepancies, which are attributed to deviations of the evolution of the hydrodynamic turbulence from freefall predictions. For example, we do not find strong agreement with the freefall predictions of quadratic scaling of mixing height with time, linear scaling of the root-mean-square velocity with time, nor constant proportionality between mixing height and the instantaneous integral scale. These discrepancies are well-known in the literature and are primarily attributed to the choice of initial conditions and the effects of finite diffusivities in simulations. We leave a more detailed analysis of the SSD in the mixing phase for future studies.
• 3.  
  In the decay phase, the dynamo scaling relations and time dependencies have large uncertainties and our model assumption of freely decaying isotropic turbulence is not a good fit. We attribute this to residual buoyant forcing and possibly finite Reynolds number effects. Fortunately, the magnetic amplification in the decay phase is small compared to the contributions from the mixing and saturation phases.
• 4.  
  In the saturation regime of the small-scale dynamo, the magnetic field is found to reach near-equipartition with the large scales and super-equipartition with the intermediate and small scales of the decaying velocity field of the RTI. We study this regime by running a single simulation with a moderate initial magnetic field so that the dynamo reaches saturation before the RTI fully relaxes.
• 5.  
  We propose that the small-scale dynamo driven by RTI turbulence helps explain observations of magnetic energy amplification in the outer regions of the post-merger in global simulations, complementary to amplification by saturation of the Kelvin–Helmholtz instability observed in the core. Applying the scaling relations to the parameter regime of RTI-unstable regions of the outer envelopes of binary neutron star mergers, the model predicts that the kinematic regime of the small-scale dynamo will end on the timescale of nanoseconds and then reach saturation on a timescale of microseconds, which are both fast in comparison to the millisecond relaxation timescale of the merger.

The flexibility of the model allows for easy extensions in future studies of the dynamo in RTI turbulence. The model primarily prescribes a time dependence for the instantaneous integral length and outer velocity scale in each phase of the RTI, which are then substituted into an equation for the dynamo growth rate. Using alternative time dependencies based on the choice of initial conditions, or allowing a time delay between forcing and dissipation, could improve understanding of the SSD in the mixing phase, for example. We leave such extensions for future work.

We acknowledge the Flatirons Center for Computational Astrophysics (CCA) and the Princeton Plasma Physics Laboratory (PPPL) for the support of collaborative CCA-PPPL meetings on plasma astrophysics where insightful comments and discussions contributed to this work. Research at the Flatiron Institute is supported by the Simons Foundation. Simulations were carried out on the Frontera cluster with NSF Frontera grant number AST20008. A.B. was supported by the DOE Grant for the Max Planck Princeton Center (MPPC). A.P. acknowledges support by the National Science Foundation under grant No. AST-1909458. E.R.M. gratefully acknowledges support from a joint fellowship at the Princeton Center for Theoretical Science, the Princeton Gravity Initiative, and the Institute for Advanced Study. V.S. was supported by the Max-Planck/Princeton Center for Plasma Physics (NSF grant PHY-1804048).

Software: Athena++ (Stone et al. 2019).

In a grid-based code like Athena++, numerical diffusivity acts as a resolution-dependent and algorithm-dependent effective viscosity that regularizes the turbulent cascade in a hydrodynamical simulation without an explicit viscosity. An explicit viscosity will only be meaningful if it is larger than the effective viscosity. We estimate the effective numerical viscosity of our numerical setup (RK3 for the timestepper and HLLD for the Riemann solver) by launching an Alfvén wave along the main diagonal in a cubical domain with zero explicit viscosity and measuring the decay rate. The Alfvén wave has wavenumber k = (1, 1, 1) · 2π/L in a background field ${{\boldsymbol{B}}}_{0}=(1,\sqrt{2},0.5)$ and the domain has triply periodic boundary conditions with a resolution N_grid ³. The kinetic energy of the wave will decay approximately exponentially ${KE}(t)\sim {e}^{-2{\nu }_{\mathrm{eff}}{{\boldsymbol{k}}}^{2}t}$. We measure ν_eff with a linear fit on a plot of $\mathrm{log}({KE}(t))$ versus time for four grid resolutions as shown in the inset of Figure 9. The main plot in Figure 9 shows a clean power-law fit ${\nu }_{\mathrm{eff}}\sim {N}_{\mathrm{grid}}^{-1.51}$. The value of the numerical viscosity is approximately ν_eff ≈ 10⁻⁶ for the resolution N_grid = 512 used in our simulations of the RTI. This informs our choice of the explicit viscosity ν ≈ 3 · 10⁻⁶ for the fiducial simulation and ν = 2 · 10⁻⁶ for our lowest choice of viscosity in the viscosity parameter scan in Section 3.

Zoom In Zoom Out Reset image size