2.1 Definitions

In this paper, we shall be concerned with a fluid dynamics defined by the momentum equation

(2.1)\begin{equation} \rho \frac{\mathrm{d} \boldsymbol{u} }{\mathrm{d} t} ={-}\boldsymbol{\nabla} P - \rho g \hat{\boldsymbol{z}},\end{equation}

where $\rho$ is density, $\boldsymbol {u}$ fluid velocity, $t$ time, $P$ total pressure, $g$ the constant gravitational acceleration and $\mathrm {d}/\mathrm {d} t\equiv \partial /\partial t + \boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla }$ the material derivative. Density evolves according to the continuity equation

(2.2)\begin{equation} \frac{\mathrm{d}\rho}{\mathrm{d} t}={-}\rho \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}, \end{equation}

and is related to $P$ by the equation of state

(2.3)\begin{equation} \rho = \rho(P,\boldsymbol{Q})\implies P=P(\rho,\boldsymbol{Q}). \end{equation}

The vector $\boldsymbol {Q}$ encodes the conserved material properties of the fluid (i.e. its Lagrangian invariants), which we allow to vary spatially. We shall primarily be concerned with the case of 2D MHD with out-of-plane magnetic field, for which the relevant components of $\boldsymbol {Q}$ are the specific entropy $s$ and specific magnetic flux $\chi \equiv B/\rho$, where $B$ is the magnetic-field strength; we shall specialise to this case in § 2.5. In different contexts, the components of $\boldsymbol {Q}$ might instead (or additionally) include mixing ratios (e.g. salinity or specific humidity) or the entropies of a coupled radiation field or cosmic rays. Over sufficiently small time scales and at sufficiently large spatial scales that diffusion can be neglected, $\boldsymbol {Q}$ is conserved following fluid particles, i.e.

(2.4)\begin{equation} \frac{\mathrm{d} \boldsymbol{Q}}{\mathrm{d} t}=0.\end{equation}

By neglecting diffusion, we exclude double-diffusive buoyancy instabilities (see, e.g. Garaud Reference Garaud2018, Hughes & Brummell Reference Hughes and Brummell2021 and references therein) from our analysis. We likewise exclude instabilities relating to anisotropic thermal conductivity (Balbus Reference Balbus2000; Quataert Reference Quataert2008). The application of our methods to the saturation of such instabilities is a topic to which we plan to return in future work (see the discussion in § 7.2).

2.2 Stability analysis

We now consider the stability under (2.1)–(2.4) of a 1D static equilibrium state, i.e. one for which $\boldsymbol {u} = 0$ and $P$, $\rho$ and $\boldsymbol {Q}$ depend on $z$ only, with $\mathrm {d} P/\mathrm {d} z = -\rho g$. The net force $\boldsymbol {F}$ on a small parcel of fluid moved in pressure balance with its surroundings and without diffusion (i.e. satisfying (2.4)) from initial height $z_1$ to new height $z_2$ is

(2.5)\begin{equation} \boldsymbol{F} ={-}gV_2\hat{\boldsymbol{z}}[\rho(P_2,\boldsymbol{Q}_1) - \rho(P_2,\boldsymbol{Q}_2)],\end{equation}

where $V_2$ is the volume of the parcel at $z_2$ and ${P_2 = P(z_2)}$, etc. Writing the difference in densities as an integral in $z$, we obtain

(2.6)\begin{equation} \boldsymbol{F} = gV_2\hat{\boldsymbol{z}}\int^{z_2}_{z_1}\mathrm{d} z \frac{\mathrm{d} \boldsymbol{Q}}{\mathrm{d} z} \boldsymbol{\cdot} \frac{\partial \rho(P_2,\boldsymbol{Q})}{\partial\boldsymbol{Q}}. \end{equation}

2.2.1 Linear stability

The criterion for linear stability follows from taking ${\delta z \equiv z_2 -z_1\to 0}$ in (2.6)

(2.7)\begin{equation} \mathcal{L}\equiv{-}\frac{\mathrm{d} \boldsymbol{Q}}{\mathrm{d} z} \boldsymbol{\cdot} \frac{\partial \ln \rho(P,\boldsymbol{Q})}{\partial\boldsymbol{Q}}>0,\quad \forall z.\end{equation}

The function $\mathcal {L}$ may alternatively be written as

(2.8)\begin{equation} \mathcal{L} ={-}\frac{\mathrm{d} \ln \rho}{\mathrm{d} z} + \frac{\partial \ln \rho (P,\boldsymbol{Q})}{\partial \ln P}\frac{\mathrm{d} \ln P}{\mathrm{d} z},\end{equation}

from which $\mathcal {L}>0$ is readily interpreted as the condition for the density of a fluid parcel displaced upwards (downwards) infinitesimally while conserving $\boldsymbol {Q}$ to be greater (less) than the density of the background fluid in the new position of the parcel.

2.2.2 Nonlinear stability

We can use (2.6) to write the criterion for nonlinear stability as

(2.9)\begin{equation} {-}\delta z\int^{z_2}_{z_1}\mathrm{d} z \frac{\mathrm{d} \boldsymbol{Q}}{\mathrm{d} z} \boldsymbol{\cdot} \frac{\partial \rho(P_2,\boldsymbol{Q})}{\partial\boldsymbol{Q}} > 0, \quad \forall\ z_1, z_2.\end{equation}

If (2.9) is satisfied, then the buoyancy force is always in the opposite direction to $\delta z$. In many cases of interest (including the examples listed above), the components $Q_i$ of $\boldsymbol {Q}$ can each be chosen such that an increase in $Q_i$ at fixed pressure causes expansion

(2.10)\begin{equation} \frac{\partial \ln\rho(P,\boldsymbol{Q})}{\partial Q_i}<0, \quad \forall\ i, P.\end{equation}

In the case of $Q_i = s$, for example, (2.10) holds for any fluid that expands under heating at fixed pressure. We shall assume in what follows that (2.10) holds. In that case, a sufficient condition for (2.9) to hold is

(2.11)\begin{equation} \frac{\mathrm{d} Q_i}{\mathrm{d} z}>0, \quad \forall\ i,z.\end{equation}

Importantly, (2.7) and (2.11) are equivalent if the vector $\boldsymbol {Q}$ has only one component (and (2.10) holds). A fluid of this kind cannot exist in a metastable equilibrium, because it is always nonlinearly stable if linearly stable. The archetypical example is ordinary hydrodynamics, for which the only component of $\boldsymbol {Q}$ is $s$, and therefore an equilibrium is nonlinearly stable to convection if it satisfies the Schwarzschild criterion (Schwarzschild Reference Schwarzschild1906)

(2.12)\begin{equation} \frac{\mathrm{d} s}{\mathrm{d} z}>0, \quad \forall\ z.\end{equation}

2.2.3 Metastability

For fluids for which $\boldsymbol {Q}$ has more than one component, (2.11) guarantees nonlinear stability, but is not necessary for the weaker condition of linear stability: if $\mathrm {d} Q_i / \mathrm {d} z<0$ for some $i$, other components of $\boldsymbol {Q}$ can compensate so that the linear-stability criterion (2.7) remains satisfied. In such cases, the equilibrium may be nonlinearly unstable despite being linearly stable, i.e. it may be metastable. In such cases, the condition for metastability may be deduced by recasting (2.5) as the upwards force per unit mass of moved fluid, ${\boldsymbol {F}\boldsymbol {\cdot } \hat {\boldsymbol {z}}/\rho (P_2,\boldsymbol {Q}_1)V_2 = ({\rm e}^\mathcal {R}-1)g}$, where

(2.13)\begin{equation} \mathcal{R} \equiv \ln \frac{\rho (P_2, \boldsymbol{Q}_2)}{\rho (P_2, \boldsymbol{Q}_1)} ={-}\int^{z_2}_{z_1} \mathrm{d} z \mathcal{L} +\int^{z_2}_{z_1} \mathrm{d} z \frac{\mathrm{d} \ln P}{\mathrm{d} z}[\kappa(P,\boldsymbol{Q})-\kappa(P,\boldsymbol{Q}_1)],\end{equation}

and

(2.14)\begin{equation} \kappa(P,\boldsymbol{Q})\equiv \frac{\partial \ln \rho(P,\boldsymbol{Q})}{\partial \ln P},\end{equation}

is the dimensionless compressibility. To obtain the second equality in (2.13), we have used the fundamental theorem of calculus, i.e. we have differentiated $\mathcal {R}$ with respect to $z_2$ and integrated the result from $z_1$ to $z_2$. Equation (2.13) separates the integrated linear buoyancy response (the first integral on the second line) with the nonlinear response (the second integral), revealing the latter to be determined by the path-integrated difference in $\kappa$ between the moving parcel and its surroundings.Footnote ² If the displaced fluid is more compressible than the fluid through which it moves, i.e. $\kappa (P,\boldsymbol {Q}_1)>\kappa (P,\boldsymbol {Q})$, then the second integral on the right-hand side of (2.13) has the same sign as $\delta z$, so its contribution to the buoyancy force is destabilising. If the stabilising effect of the integral involving $\mathcal {L}$ is sufficiently small, then the second integral dominates in (2.13) for $\delta z$ larger than some critical value, $\delta z_c$. Such an equilibrium is metastable.

2.3 Direction and size of displacement required for nonlinear instability

Because metastability requires more compressible fluid to be moved through less compressible fluid [see (2.13)], equilibria can be metastable only to perturbations in the direction in which the compressibility of the background fluid decreases.Footnote ³ This direction is given by the sign of the compressibility scale height

(2.15)\begin{equation} H_{\kappa} \equiv{-}\left(\frac{\partial \ln\kappa }{\partial \boldsymbol{Q}} \boldsymbol{\cdot} \frac{\mathrm{d} \boldsymbol{Q} }{\mathrm{d} z}\right)^{{-}1}.\end{equation}

We can estimate the local value of $\delta z_c$ in the limit of $\mathcal {L}\to 0$ by balancing the sizes of the two integrals in (2.13) to find that

(2.16)\begin{equation} \delta z_c \simeq \frac{ H_P H_{\kappa}{\mathcal{L}}}{\kappa},\end{equation}

where we have defined the pressure scale height $H_P \equiv |\mathrm {d} \ln P/ \mathrm {d} z|^{-1}$, neglected variation in $\mathcal {L}$ on the scale $\delta z_c$ and used the fact that $-\kappa H_{\kappa }^{-1}$ is the coefficient of $\delta z$ in the small-$\delta z$ expansion of the term in square brackets in (2.13). According to (2.16), $\delta z_c$ becomes arbitrarily small compared with any stratification scale height as $\mathcal {L}\to 0$. Nonetheless, we note that, because metastability of stratified equilibria is a nonlinear effect, it is not captured by Boussinesq-like equations that employ linear approximations of equilibrium gradients.

2.4 Explosive instability

The equation of motion of a small fluid parcel displaced by a distance $\delta z\sim \delta z_c$ that is much smaller than any scale height $H$ of the stratification is

(2.17)\begin{equation} \frac{\mathrm{d}^2 \delta z}{\mathrm{d} t^2} = \left(-\delta z+\frac{\delta z^2}{\delta z_c}\right)\mathcal{L}.\end{equation}

It follows from (2.17) that, for ${0<\delta z_c \ll \delta z \ll H}$ (the case of $\delta z_c<0$ is analogous), the motion of a fluid parcel is explosive, viz.,

(2.18)\begin{equation} \delta z \propto \frac{1}{(C-t)^2}, \end{equation}

where the constant $C$ is determined by initial conditions. Explosive growth of $\delta z$ persists until either ${\delta z \sim H}$, whereupon (2.17) is no longer valid, the rising fluid element is shredded by Kelvin–Helmholtz instability or its speed approaches that of sound and thus (2.5) no longer applies.

2.5 Case of 2D MHD

In the remainder of this paper (with the exception of Appendix B, where we consider moist hydrodynamics), we focus on the case of an equilibrium supported against gravity both by thermal pressure $p$ and by the magnetic pressure associated with a straight magnetic field with spatially dependent strength $B$. Thus, ${P = p + B^2/2}$. We restrict attention to 2D dynamics in the plane perpendicular to the magnetic field (i.e. to 2D interchanges of flux tubes). Then the specific magnetic flux

(2.19)\begin{equation} \chi = \frac{B}{\rho},\end{equation}

is conserved in a Lagrangian sense. For symmetry with this definition of $\chi$, we take advantage of the fact that any monotonically increasing function of specific entropy constitutes a good choice for the conserved quantity $s$, and thus let

(2.20)\begin{equation} s = \frac{p^{1/\gamma}}{\rho}, \end{equation}

where $\gamma$ is the adiabatic index. To avoid confusion with the true specific entropy, which is proportional to $\exp (s/C_v)$ with $C_v$ the specific heat capacity, we hereafter refer to $s$ as the ‘entropy function’ (akin to potential temperature in atmospheric science, see Appendix B). Thus, $\boldsymbol {Q} = (s, \chi )$ and

(2.21)\begin{equation} P(\rho,\boldsymbol{Q})=\rho^\gamma s^\gamma + \rho^2 \chi^2/2. \end{equation}

With these choices, (2.7) yields the linear-stability condition

(2.22)\begin{equation} \mathcal{L}= \frac{c_s^2}{c^2} \frac{\mathrm{d} \ln s}{\mathrm{d} z} +\frac{v_A^2}{c^2}\frac{\mathrm{d} \ln \chi}{\mathrm{d} z}>0,\end{equation}

where $c_s\equiv \sqrt {\gamma p/\rho }$ is the sound speed, $v_A\equiv B/\sqrt {\rho }$ the Alfvén speed and $c = \sqrt {c_s^2 + v_A^2}$ is the velocity of compressive waves [(2.22) is sometimes called the modified Schwarzschild criterion]. The compressibility $\kappa$ (2.14) is

(2.23)\begin{equation} \kappa(P,s,\chi)= \frac{1+\beta}{2+\gamma \beta},\end{equation}

where $\beta \equiv 2p/B^2$ is the plasma beta (the ratio of thermal to magnetic pressures) which is determined from $P$, $s$ and $\chi$ via

(2.24)\begin{equation} \frac{1}{\beta}\left(1+\frac{1}{\beta}\right)^{{2}/{\gamma}-1} = \frac{1}{2}\left(\frac{\chi}{s}\right)^2 P^{{2}/{\gamma}-1}.\end{equation}

Equation (2.23) reveals that $\kappa$ increases monotonically with $\beta$ for $\gamma <2$, from $\kappa = 1/2$ at $\beta = 0$ to $\kappa \to 1/\gamma$ as ${\beta \to \infty }$. It follows that the nonlinear buoyancy response in (2.13) is destabilising when fluid with large $\beta$ moves through ambient fluid with smaller $\beta$. According to (2.24), the $\beta$ of a flux tube with given $s$ and $\chi$ depends on pressure and therefore is not constant during its motion.Footnote ⁴ At any given pressure, however, $\beta$ is a monotonically increasing function of the ratio $s/\chi$. Therefore, an equilibrium that is sufficiently close to marginal linear stability is always nonlinearly unstable in the direction in which $s/\chi$ decreases [the same conclusion can be obtained from (2.15) or by expanding (2.5) to quadratic order in $\delta z$ directly, see Appendix C].

2.6 Examples of metastable equilibria

Explicit examples of metastable equilibria may be obtained as follows. First, we change variables from height $z$ to the total mass supported at height $z$

(2.25)\begin{equation} m = \int_z^{\infty} \rho(z')\,\mathrm{d} z',\end{equation}

so that the equilibrium condition becomes ${P = mg}$, which may be expressed as

(2.26)\begin{equation} \rho^{\gamma}s^{\gamma}+\frac{1}{2}\rho^2 \chi^2 = mg,\end{equation}

or, using (2.24), as

(2.27)\begin{equation} \frac{1}{\beta}\left(1+\frac{1}{\beta}\right)^{{2}/{\gamma}-1} = \frac{1}{2}\left( \frac{\chi}{s} \right)^2 (mg)^{{2}/{\gamma}-1}.\end{equation}

Throughout this work, we shall find it convenient to use the supported mass $m$ as a proxy for height because the mass of a flux tube is preserved as it moves, while its cross-sectional area is not (note that $m$ decreases with increasing $z$).

We seek an equilibrium that is close to marginal linear stability [so that the second integral in (2.13) stands a chance of dominating the first]. We therefore demand that

(2.28)\begin{equation} \mathcal{L}=\frac{\epsilon}{H_P},\end{equation}

where $\mathcal {L}$ is given by (2.22), $H_P$ is the total-pressure scale height $[\mathrm {d} \ln P /\mathrm {d} z]^{-1}=m/\rho$ and $\epsilon \ll 1$ is a small number. Equation (2.28) constitutes a differential equation involving $s$, $\chi$ and their gradients in $m$ – to solve it, we specify a relationship between $s$, $\chi$ and $m$ that defines the particular equilibrium under consideration. We choose to specify as a function of $m$ the ratio $s/\chi$, which controls the compressibility of the fluid (§ 2.5). We therefore recast (2.28) [with $\mathcal {L}$ given by (2.22)] as

(2.29)\begin{equation} \frac{\mathrm{d} \ln s}{\mathrm{d} m} = \frac{2}{\gamma \beta + 2}\frac{\mathrm{d}}{\mathrm{d} m}\ln\left(\frac{s}{\chi}\right)-\frac{\epsilon}{m},\end{equation}

which we integrate numerically for $s$ as a function of $m$. We determine the dependence of all other quantities on $m$ (or $z$) via their definitions and (2.25) and (2.27).

Some example cases are visualised in figure 3, where we show $s$, $\chi$, $s/\chi$ and $\beta$ as functions of $z$ (lower row of panels) together with the force (2.5) per unit mass on a small parcel of fluid moved from height $z_1$ to $z_2$ as a function of $z_1$ and $z_2$ (upper two rows of panels). Figure 3a shows the case of

(2.30)\begin{equation} \frac{s}{\chi} = (m_{\mathrm{tot}}g)^{1/\gamma - 1/2}\left(3\frac{m}{m_\mathrm{tot}}+0.3\right),\end{equation}

for which $s/\chi$ increases with $m$. The most compressible material is therefore at the bottom of the atmosphere and the equilibrium is unstable to upwards displacements (towards larger $z$). This is the equilibrium whose relaxation is visualised in figures 1 and 2. Figure 3b shows the case of

(2.31)\begin{equation} \frac{s}{\chi} = (m_{\mathrm{tot}}g)^{1/\gamma - 1/2} \left(3\frac{m}{m_\mathrm{tot}}+0.3\right)^{{-}1},\end{equation}

for which $s/\chi$ decreases with $m$, so produces an equilibrium unstable to downward displacements (towards smaller $z$). Figure 3c shows the case of

(2.32)\begin{equation} \frac{s}{\chi} = (m_{\mathrm{tot}}g)^{1/\gamma - 1/2} \left[3\exp\left(-\frac{(m/m_{\mathrm{tot}} -0.2)^2}{0.1^2}\right)+0.3\right],\end{equation}

for which $s/\chi$ has a maximum at $m=0.2m_{\mathrm {tot}}$; the equilibrium is therefore unstable to both upwards and downwards displacements in the vicinity of the maximum. Finally, figure 3d shows the case of

(2.33)\begin{equation} \frac{s}{\chi} = (m_{\mathrm{tot}}g)^{1/\gamma - 1/2} \left[3\left(1-\exp\left(-\frac{(m/m_{\mathrm{tot}} -0.2)^2}{0.1^2}\right)\right)+0.3\right],\end{equation}

for which $s/\chi$ has a minimum at $m=0.2m_{\mathrm {tot}}$. The most compressible fluid is therefore situated at the top and bottom of the atmosphere, so it is unstable to both downwards and upwards displacements.

Figure 3. Top row: the upward force (2.5) per unit mass on a small fluid parcel moved in pressure balance and without diffusion from height $z_1$ to $z_2$ for each of the profiles described in § 2.6. Values of $\epsilon _0$ in (2.34) are indicated on each panel. Middle row: the same as the top row, but with $\epsilon _0=0$ (marginal linear stability in the bulk). Bottom row: profiles of the entropy function $s=p^{1/\gamma }/\rho$, specific magnetic flux $\chi =B/\rho$, their ratio and the plasma $\beta$ (2.27) as a function of height $z$ (these profiles correspond to the top row specifically, but the profiles that correspond to the middle row look essentially the same because the differences in $\epsilon _0$ are very small). Panel (a) corresponds to (2.30), (b) to (2.31), (c) to (2.32) and (d) to (2.33).

In figures 1–3 and in the rest of the paper, we choose units of mass and length such that $\rho =1$ and ${m = m_{\mathrm {tot}}=1}$ at the bottom of the atmosphere (${z=0}$), so the total-pressure scale height $H_P=1$ there. We choose units of time such that $g=1$; it follows that the time taken for a compressive wave to traverse a total-pressure scale height ${H_P/c\sim \sqrt {H_P/g}}$ is $1$ at the bottom of the atmosphere. We limit the total height of the atmospheres to ${z_{\mathrm {max}}=3.5}$ in these units, where $\rho \lesssim 10^{-2}$ in each of the four cases described above. We stabilise the equilibria near $z=0$ and $z=z_{\mathrm {max}}$ by choosing

(2.34)\begin{equation} \epsilon = \epsilon_0 + \tanh\left(\frac{z-z_u}{\varDelta_u}\right) - \tanh\left(\frac{z-z_l}{\varDelta_l}\right) + 2,\end{equation}

in (2.28), with $\varDelta _l=\varDelta _u = 0.25$, $z_l = 0.25$, $z_u = 3.25$. This ensures that artificial boundary effects do not impact our numerical or theoretical analyses. We choose ${\epsilon _0 = 10^{-2}}$ for the simulations visualised in figures 1 and 2.

The top row of panels in figure 3 correspond to equilibria with small positive values of $\epsilon _0$; these are linearly stable, and exhibit the various sorts of metastability described above. Unless stated otherwise, we shall in the rest of the paper focus attention on the case of marginal linear stability, i.e. $\epsilon _0 = 0$ in (2.34). This is because we expect nonlinear instability to be triggered for equilibria close to marginal linear instability in practice. The middle row of figure 3 corresponds to the $\epsilon _0 = 0$ case.

3.1 Estimating the available energy

Before proceeding to solve the LSA problem, which can only be achieved numerically in most cases, let us try to estimate the outcome analytically: What is the typical available potential energy of a metastable atmosphere? This turns out to be a small fraction of the total potential energy, a fact we that we utilise in § 5 to predict the relaxation of destabilised equilibria.

The change in potential energy $\delta E$ that results from moving a slice of atmosphere upwards from supported mass $m_a$ to new supported mass $m_b$, shuffling downwards the slices that it passes on the way, is

(3.8)\begin{align} \frac{\delta E}{\Delta m} & = \mathcal{E}(m_b, m_a)-\mathcal{E}(m_a, m_a)+\sum_{i=b}^{a-1}[\mathcal{E}(m_{i+1},m_i)-\mathcal{E}(m_{i},m_i)]\nonumber\\ & \simeq \int^{m_b}_{m_a}\mathrm{d} m \left[\frac{g}{\rho(m,m_a)}-\frac{g}{\rho(m,m)}\right] \nonumber\\ & ={-} g\int^{z_a}_{z_b}\mathrm{d} z \left[\frac{\rho(m(z),m(z))}{\rho(m(z),m_a)}-1\right], \end{align}

where $m(z_a) = m_a$ and $m(z_b) = m_b$. We have used (3.7) in moving from the first to the second line of (3.8). The integrand that appears in the last line of (3.8) is straightforwardly recognised as the net buoyancy force (2.5) per unit mass of fluid, so, sensibly, $\delta E$ is just the work done by this force on the moving slice.Footnote ⁵ Evidently, the energy that can be liberated by this process is maximal when the ratio that appears inside the integrand in the last line of (3.8) – i.e., of the density of the ambient fluid to the density of the moving slice – is maximal. We can evaluate this ratio using (2.13), which yields, in the most optimistic case of (i) marginal linear stability, i.e. $\mathcal {L}=0$, (ii) $\kappa =1/\gamma$ inside the slice (maximally compressible large-$\beta$ fluid) and (iii) $\kappa = 1/2$ for the ambient fluid (minimally compressible small-$\beta$ fluid)

(3.9)\begin{equation} \frac{\rho(m,m_a)}{\rho(m,m)} = \left(\frac{m}{m_a}\right)^{1/\gamma - 1/2}.\end{equation}

Substituting this into (3.8), choosing $m_a=m_{\mathrm {tot}}$ (the slice originates from the bottom of the atmosphere) and evaluating integrals, we find that

(3.10)\begin{equation} \frac{\delta E}{E_0} = \frac{3}{4}\frac{\Delta m}{m_{\mathrm{tot}}}\left[\frac{\gamma}{\gamma - 1}\left(\left(\frac{m_b}{m_{\mathrm{tot}}}\right)^{1-1/\gamma}- 1\right)-2\left(\left(\frac{m_b}{m_{\mathrm{tot}}} \right)^{1/2}-1\right)\right],\end{equation}

where $E_0$ is the initial total energy of the atmosphere (assuming that it is mostly populated with small-$\beta$ fluid, as is consistent with assumption (iii) above).

We plot (3.10) in figure 4 for the case of $\gamma = 5/3$. We observe that the fraction of energy liberated per mass fraction of fluid moved decreases sharply as a function of $m_b/m_{\mathrm {tot}}$ (even though its limiting value of $3/8$ as $m_b/m_{\mathrm {tot}}\to 0$ is finite). This is significant as, because fluid parcels exclude each other (i.e. different slices cannot all be assigned to the same supported mass), the vast majority of slices that are interchanged in a global rearrangement experience an order-unity change in supported mass. Such reassignments are far less profitable that ones that take a slice all the way to the top of the atmosphere (i.e. $m_b/m_{\mathrm {tot}}\to 0$). For example, figure 4 shows that the fraction of energy liberated per mass fraction of fluid moved is approximately $10^{-2}$ for $m_b = m_{\mathrm {tot}}/2$. The reason is that an order-unity change in pressure only results in a small change in density, owing to the smallness of the exponent $1/\gamma - 1/2 = 0.1$ that appears on the right-hand side of (3.9). A more detailed calculation (see Appendix D) reveals that $10^{-2}$ is indeed a good estimate for the fractional available energy: the maximum fractional available energy of a (marginally) stable atmosphere consisting of small-$\beta$ fluid above a layer of large-$\beta$ fluid is $1.75\,\%$.

Figure 4. The fraction (3.10) of energy liberated when a slice of fluid with mass $\Delta m$ is moved from the bottom of an atmosphere at marginal linear stability to a new position where the supported mass is $m_b$, under the most optimistic assumptions about the compressibility of the slice and that of the fluid through which it moves.

The above estimates correspond to the most optimistic assumptions: we shall find that the fractional available energies of the equilibria described in § 2.6 and represented in figure 3 are less than $1\,\%$ (by roughly an order of magnitude), because (i) these equilibria have finite $\beta$, and thus do not have the extremal values of the fluid compressibility considered here; and (ii) the stable buffer regions contribute to the potential energy of the equilibrium (and prevent the relaxing fluid from accessing large pressure differences) but do not participate in the reassignment.

3.2 The Hungarian algorithm

The discrete energy-minimisation problem (3.6) can be solved without recourse to numerical optimisation only in the limits $\beta \to \infty$ and $\beta \to 0$, i.e. the cases where only one of thermal or magnetic pressure support the atmosphere against gravity. For ${\beta \to \infty }$, ${\mathcal {E} \to \gamma s(mg)^{1-1/\gamma }/(\gamma - 1)}$, which increases monotonically with both $m$ and $s$. The arrangement with least total energy is therefore the one for which the slice with largest $s$ has smallest $m$, the slice with the next largest $s$ has the next smallest $m$, and so on. It follows that the profile with the smallest energy is the unique rearrangement for which $s$ is a monotonically increasing function of height. Indeed, we already know this state to be nonlinearly stable by (2.12). A similar conclusion is obtained for ${\beta \to 0}$, for which ${\mathcal {E} \to 2(mg)^{1/2} \chi }$ and so the nonlinearly stable atmosphere is the one with $\chi$ increasing monotonically with height. Analogous simple constructions for the minimum-energy configuration do not exist in the finite-$\beta$ case: the optimal solution that balances the competing imperatives of ‘entropy should increase upwards’ and ‘flux should increase upwards’ is non-trivial.

Solution of the LSA problem in general relies on the observation that the modified cost matrix $\tilde {\mathcal {E}}(m_{j}, m_i)$, where

(3.11)\begin{equation} \tilde{\mathcal{E}}(m_{j}, m_i) = \mathcal{E}(m_{j}, m_i) - a(m_{j}) - b(m_{i}),\end{equation}

has the same optimal assignment of agents to tasks as does the original cost matrix ${\mathcal {E}(m_{j}, m_{i})}$. This is intuitive: if the cost of assigning a given agent (slice) to each task (position) decreases by the same amount, their optimal assignment will not change (although the total cost of the solution will decrease). Likewise, if a particular task becomes more expensive by the same amount for all agents, the optimal choice of agent for that task will remain the same. $\tilde {\mathcal {E}}$ is called the ‘normal form’ of the cost matrix $\mathcal {E}$ if it satisfies the following properties (Burkard et al. Reference Burkard, Dell'Amico and Martello2012):

1. (i) $\tilde {\mathcal {E}}(m_{j}, m_i)\geq 0$,

2. (ii) There exists at least one bijection $\sigma$ such that $\tilde {\mathcal {E}}(m_{\sigma (i)}, m_i)=0, \ \forall \, i$.

A proof that an $\tilde {\mathcal {E}}(m_j,m_i)$ with these properties can always be found with suitable choices of $a(m_j)$ and $b(m_i)$ is provided by the Hungarian algorithm (Kuhn Reference Kuhn1955; Munkres Reference Munkres1957), which consists of a series of row and column operations on the cost matrix that are guaranteed to reduce it to normal form in polynomial (in $N$) time. The utility of the normal form $\tilde {\mathcal {E}}(m_j,m_i)$ is that each of the bijections $\sigma$ constitutes an optimal assignment.Footnote ⁶

3.3 Energy minimisation: numerical results and interpretation

In this section, we present the outcomes of energy minimisation using the Hungarian algorithm (§ 3.2) for each of the equilibria described in § 2.6 and represented in figure 3.

3.3.1 Metastable upwards, (2.30)

We first consider the equilibrium defined by (2.30), which is nonlinearly unstable to upward displacements (figure 3a). The solution of the LSA problem is visualised in figure 5, which compares the initial and minimum-energy assignments of the slices. In the minimum-energy assignment, material initially from the bottom of the metastable part of the atmosphere (i.e. $z\gtrsim z_l$; see (2.34)) is reassigned to the top $z\lesssim z_u$, as is intuitive. The order in which the slices that are re-assigned are stacked also reverses. This happens because, as the material moves to smaller $m$, its $\beta$ increases (2.27), and therefore the contribution of $s$ to the linear-stability criterion becomes more important relative to $\chi$. As a consequence, the stacking order reverses so that $s$ increases upwards.

Figure 5. Visualisation of the assignment of 1D slices that minimises the total energy, (3.6), with $\Delta m = 5\times 10^{-4}m_{\mathrm {tot}}$, for the upwards-unstable profile defined by (2.30). Panels on the left show the initial profiles of $s$ and $\chi$ as functions of height $z$, while panels on the right show the minimum-energy assignment. The slices are coloured by their height $z$ in the initial state to aid comparison. Blue slices from $0.8< z<1.0$ are moved to $2.1< z<3.1$, reversing order and foliating with red slices originally from $2.7< z<3.1$. Each slice has vertical extent ${\Delta z = \Delta m/\rho}$ with $\rho$ given by (2.26).

A striking feature of the minimum-energy state is that slices that were adjacent at the bottom of the atmosphere do not remain adjacent in the minimum-energy state. Instead, slices from the bottom become foliated with those at the top. The scale of the foliation is set by $\Delta m$, and therefore is arbitrarily small as $\Delta m \to 0$. This is illustrated by figure 6, which shows, for three different values of $\Delta m$, the normal-form cost matrix $\tilde {\mathcal {E}}(m_j,m_i)$ [defined by (3.11)], with its zeros (which indicate the optimal assignment $j=\sigma (i)$) marked with white circles. Because the scale of foliation depends on $\Delta m$, the optimal assignment of slices does not converge as $\Delta m \to 0$, although it does converge in a coarse-grained sense: the proportion of slices assigned to any small finite range of $m$ that originated from any similarly small given range of $\mu$ converges as $\Delta m \to 0$ (provided that $\tilde {\mathcal {E}}$ converges as $\Delta m\to 0$).Footnote ⁷

Figure 6. The normal form of the cost matrix $\tilde {\mathcal {E}}(m,\mu )$ [see (3.11) for its definition] for the unstable-upwards profile defined by (2.30), for three different choices of the discretisation scale $\Delta m$. White dots show the optimal assignment.

Foliation occurs when the material properties of the fluid vary with $m$ at a different rate in the part of the atmosphere from which a slice originates than in the part to which it is reassigned. We demonstrate this fact by considering the motions of two slices from the bottom of the atmosphere to the top, each displacing downwards the slices through which they pass (as in § 3.1). Let the first slice have initial assignment $m_a$ and new assignment $m_b$. Because the new assignment is optimal, the total energy has a local minimum as a function of displacement of this slice when its density in its new location, $\rho (m_b,m_a)$, is equal to the density $\rho (m_b,m_b)$ of the slices that surround it there (neutral buoyancy); this makes the integrand in the second line of (3.13) zero at $m=m_b$.Footnote ⁸ Now let us consider a second slice of fluid that initially neighbours the first, i.e. that originates from supported mass $m_a + \Delta m$. This slice reaches neutral buoyancy at a different supported mass $m_b+\delta m$. If the density of the background equilibrium changes more slowly with supported mass at $m_b$ than at $m_a$, then $\delta m > \Delta m$. Setting the density of the second slice, $\rho (m_b+\delta m,m_a+\Delta m)$, equal to that of the ambient fluid at supported mass $m_b+\delta m$, i.e. ${\rho (m_b+\delta m,m_b+\delta m)}$ (because we are concerned with the motion of only two slices of infinitesimal thickness, we neglect the fact that the reassignment of the first slice might have changed the identity of the fluid at $m_b$), we discover that

(3.12)\begin{equation} \left.\frac{\delta m}{\Delta m}=\left.\frac{\partial \rho(m_b,\mu)}{\partial \mu}\right|_{\mu=m_a}\right/ \left.\frac{\partial \rho(m_b,\mu)}{\partial \mu}\right|_{\mu=m_b}.\end{equation}

Thus, the adjacency of the slices if not preserved ($\delta m\neq \Delta m$) if the rate of change with supported mass of the density that a slice would have if moved to $m_b$ is different for slices originating at $m_a$ than at $m_b$.

Although (3.12) reveals the physical reason for foliation, it fails if the fraction on the right-hand side is small ($\delta m \lesssim \Delta m$ is not allowed for the discrete problem because slices exclude each other). It also does not apply in the case where a substantial mass of fluid is reassigned, in which case the background equilibrium through which the first slice moves is different from that through which the last slice does. A general treatment of such cases is as follows. Let us suppose that we are somehow given all the optimal assignments $\sigma (i)$, except those to some small range of $m$, i.e. those with $m_{\sigma (i)}=m_b+\delta m_{\sigma (i)}$, and seek the condition under which the optimal choice of the remaining assignments will be a foliated state. The contribution to the total energy of the slices that remain to be assigned is

(3.13)\begin{align} \delta E & = \Delta m \sum_{i} \mathcal{E}(m_b+\delta m_{\sigma(i)}, m_i) \nonumber\\ & = \Delta m \sum_{i} \left[ \mathcal{E}(m_b, m_i)+ \frac{\delta m_{\sigma(i)}g}{ \rho(m_{b}, m_{i})}+\mathcal{O}(\delta m_{\sigma(i)}^2)\right], \end{align}

where sums are over all indices $i$ of the slices that remain to be assigned, and we have used (3.7) to obtain the second equality. The first term inside the square bracket in the second line of (3.13) is independent of the assignment $\sigma$. The second term takes the form of the differential supported mass $\delta m_{\sigma (i)}$ multiplied by a quantity that does not depend on $\sigma (i)$, viz., $1/{\rho (m_b,m_i)}$. This yields a local stacking rule: $\delta E$ is minimised when the slice with largest $\rho (m_{b}, m_i)$ is assigned to the largest supported mass, the next largest $\rho (m_b, m_i)$ is assigned to the next largest supported mass, and so on. This is physically intuitive: if more dense slices were situated above less dense ones, the equilibrium would be Rayleigh–Taylor unstable. We deduce that in order for slices from different initial locations to become foliated in the final state they must have

(3.14)\begin{equation} \rho(m_b,m_{i'})-\rho(m_b,m_i)=\mathcal{O}(\Delta m).\end{equation}

As $\Delta m\to 0$, slices can be foliated in the vicinity of some $m$ only if they have the same density at that $m$.

A foliated minimum-energy state may be considered a 1D representation of a minimum-energy state that is actually 2D, as follows. Because the scale of foliation is set by $\Delta m$, we can reconfigure the foliated state (figure 7a) into a 2D state (figure 7b) by rearranging fluid parcels locally in $z$, i.e. over a small vertical distance $\delta z \sim \Delta m / \rho$, and then sorting globally in $x$ (with $P$, $s$ and $\chi$ fixed for each parcel) to remove small-scale variation (figure 7c,d). Because the foliated state is a stable equilibrium with respect to local rearrangements in $z$, the force acting on each fluid parcel in the new state will be proportional to $\delta z$, and therefore vanishes as $\Delta m \to 0$. Thus, the 2D state is also an equilibrium state with the same energy as the foliated one.Footnote ⁹ In § 4, we shall show with direct numerical simulations that minimum-energy 2D states are the result of the nonlinear relaxation of a destabilised metastable state in certain regimes.

Figure 7. 2D minimum-energy states. Panel (a) shows the 1D stacking with minimum energy in the $x-z$-plane. Panel (b) shows a discrete approximation to the 2D ground state, obtained by arranging vertical slices sequentially into fixed vertical bins of fixed extent $\Delta z$, while preserving the $P$, $s$ and $\chi$ of each slice. Each slice occupies the same area as before (because its mass and density each remain constant) but now the slices are arranged horizontally within the bin, with left to right corresponding to increasing height in panel (a). Panel (c) shows the state in panel (b) but sorted horizontally, which does not change the energy and removes the imprint of the foliation. Panel (d) shows the expected equivalent of Panel (c) at very large resolution, but is obtained differently, by taking the small-thermodynamic-temperature limit of Lynden-Bell statistical mechanics (§ 5.2).

To conclude our discussion of minimum-energy states of equilibria defined by (2.30), we present in figure 8 a comparison of energy-minimising assignments for different values of the parameter $\epsilon _0$, which controls the distance from marginal stability [see (2.34)]. First, we note that $\epsilon _0 = 0$, i.e. marginal linear stability, is not a special point as far as the minimum-energy assignments are concerned: the assignments with $\epsilon _0 = 0.0075$ (linearly stable), $\epsilon _0 = 0.0$ (marginal) and $\epsilon _0 = -0.015$ (linearly unstable) are qualitatively similar, although, in the unstable cases of $\epsilon _0 = -0.015$ and $\epsilon _0 = -0.075$, foliation occurs over a much wider range of supported masses (and material from three different initial locations are foliated together near the top of the atmosphere).

Figure 8. The energy-minimising assignments of slices from initial supported mass $\mu$ to new supported mass $m$ for the initial profile (2.30) with different values of the parameter $\epsilon _0$, which controls linear stability via (2.34).

Figure 9 shows the ratio of the available energy ${E_{\mathrm {avail}}=E_0-E_{\mathrm {min}}}$ to the original potential energy $E_0$ as a function of $\epsilon _c-\epsilon _0$, where $\epsilon _c\simeq 1.7\times 10^{-2}$ is the largest value of $\epsilon _0$ for which the atmosphere has a restacking with smaller energy. We see that $E_{\mathrm {avail}}/ E_0$ is small: it is around $10^{-3}$ for $\epsilon _0 = 0$, which is the value that corresponds to figures 5 and 6. This is despite the fact that the minimum-energy assignment involves significant rearrangement of the atmosphere (around $10\,\%$ by mass of the atmosphere is reassigned upwards for $\epsilon _0 = 0$, see inset to figure 9). As explained in § 3.1, the reason for the smallness of $E_{\mathrm {avail}}/ E_0$ is the fact that fluid slices exclude each other and so only very few of them can experience a significant change in total pressure as a result of reassignment.

Figure 9. The available energy $E_{\mathrm {avail}}=E_0-E_{\mathrm {min}}$ as a fraction of the initial potential energy $E_0$, plotted as a function of $\epsilon _c-\epsilon _0$, where $\epsilon _c\simeq 1.7\times 10^{-2}$ is the largest value of $\epsilon _0$ in (2.34) for which the initial state is metastable. The inset shows the fraction of fluid that is assigned to a smaller supported mass than its initial one under optimal reassignment. Red lines correspond to $\epsilon _0=0$.

We observe from figure 9 that

(3.15)\begin{equation} \frac{E_{\mathrm{avail}}}{E_0}\propto (\epsilon_c-\epsilon_0)^2 \ \mathrm{as}\ \epsilon_c-\epsilon_0\to 0;\end{equation}

this scaling is readily interpreted as a result of the fact that both (i) the typical amount of energy liberated when a slice is reassigned from the bottom of the atmosphere to the top, and (ii) the number of slices that are reassigned in this way, are proportional to $\epsilon _c - \epsilon _0$ when the latter is small [see inset to figure 9; $s$, $\chi$ and, therefore, the buoyancy force on a displaced fluid element, (2.5), depend linearly on $\epsilon _0$, by (2.29)]. This argument being independent of the particular profile under consideration, we expect a quadratic dependence of $E_{\mathrm {avail}}/E_0$ on $\epsilon _c - \epsilon _0$ for any metastable profile equilibrium with $\epsilon _0$ sufficiently close to $\epsilon _c$. We shall find in § 3.3.2 that a quadratic scaling is indeed reproduced for the profile represented by (2.31).

3.3.2 Metastable downwards, (2.31)

We now turn to the second example case introduced in § 2, (2.31), which describes an initial state that is metastable to downwards displacements. The 1D minimum-energy state associated with this profile [with $\epsilon _0=0$ in (2.34)] is shown figure 10, which is the analogue for (2.31) of figure 5. The minimum-energy assignment is similar qualitatively to the one examined in § 3.3.1: in this case, material from the top of the atmosphere moves to the bottom, reverses stacking order, and becomes foliated with the material already there (in fact, material from three different initial heights becomes foliated).

Figure 10. Visualisation of the minimum-energy assignment for the equilibrium defined by (2.31), with $\epsilon _0=0$ in (2.34). Details are the same as for figure 5.

Figure 12 shows the dependence of the optimal assignment on the value of $\epsilon _0$ in (2.34). A qualitatively new feature appears in the cases of ${\epsilon _0=-0.025}$ and ${\epsilon _0=-0.069}$: for these unstable equilibria, there exists a range of the initial supported mass coordinate (in the vicinity of $\mu \simeq 0.25$ for the former case and $\mu \simeq 0.4$ for the latter) for which slices that are neighbouring in the initial state are alternately assigned to two different final locations. Like foliation, this phenomenon can be interpreted as a consequence of our seeking a 1D optimisation when, in fact, the true optimal assignment is higher-dimensional in the continuous limit $\Delta m \to 0$. In this case, the optimal assignment involves splitting the fluid at given height in a horizontal sense, and reassigning it to multiple new locations. In Appendix E, we derive a necessary condition for this kind of one-to-many assignment and prove that two is, in fact, the largest possible value of ‘many’.

Figure 11. Minimum-energy assignments for the profile (2.31). We observe that the optimal assignment is one to two over certain ranges of $m_1$ in the cases with $\epsilon _0<0$.

Figure 12. Fractional available energy as a function of $\epsilon _c - \epsilon _0$ for the profile (2.31), where $\epsilon _c=1.2\times 10^{-2}$ is the largest value of $\epsilon _0$ in (2.34) for which the equilibrium is metastable.

3.3.3 Bi-directional metastability, (2.32)

We visualise in figure 13 the optimal assignment for the profile defined by (2.32), which has a local maximum in its profile of $s/\chi$ at $z\simeq 1.75$ (see figure 3c) and therefore the fluid is metastable both to upward and downward displacements there. It is intuitive, therefore, that the minimum-energy state should be obtained by reassignment of slices from the vicinity of $z\simeq 1.75$ to both the top and the bottom of the region that is at marginal linear stability. This is indeed the case in figure 13 (with foliation between moved and ambient fluid).

Figure 13. Visualisation of the minimum-energy assignment for the equilibrium defined by (2.32), with $\epsilon _0=0$ in (2.34). Details are the same as for figure 5.

Figure 14 shows the optimal assignments for different values of $\epsilon _0$. The available energy associated with this profile is a fraction $6\times 10^{-4}$ of the initial total energy for $\epsilon _0 = 0$, which is comparable to the equilibria considered in §§ 3.3.1 and 3.3.2.

Figure 14. Minimum-energy assignments for the profile (2.32) for different values of $\epsilon _0$ in (2.34).

3.3.4 Overturning metastability, (2.33)

In figure 15, we visualise the optimal assignment for the profile defined by (2.33), which, in contrast to (2.32), has a local minimum in its profile of $s/\chi$ at $z\simeq 1.75$ (see figure 3d). The fluid at the top of the atmosphere is therefore nonlinearly unstable to downwards motions, while the fluid at the bottom is nonlinearly unstable to upwards motions – we expect therefore that the minimum-energy state will be reached by an ‘overturning’ of the atmosphere. This is roughly what we observe in figure 15, although the precise optimal assignment is remarkably complex (see figure 16 for a visualisation of the optimal assignments for different values of $\epsilon _0$). The available energy associated with this profile at $\epsilon _0 = 0$ is $2\times 10^{-3}$ of the initial total, which is somewhat more than that of the equilibria considered in §§ 3.3.1–3.3.3.

Figure 15. Visualisation of the minimum-energy assignment for the equilibrium defined by (2.33), with $\epsilon _0=0$ in (2.34). Details are the same as for figure 5.

Figure 16. Minimum-energy assignments for the profile (2.33) for different values of $\epsilon _0$ in (2.34).

4.1 Conditions for diffusion to be absent

In order for diffusion of $s$ and $\chi$ to be negligible during the whole period of relaxation, the flow generated during relaxation must decay before it is able to mix $s$ and $\chi$ to scales $l$ such that their diffusion time scales ($l^2/K$ and $l^2/\eta$, respectively, where $K$ and $\eta$ are the thermal and magnetic diffusivities) become comparable to the flow's decay time scale. The latter is $H^2/\nu$ ($\nu$ is the kinematic viscosity) independently of whether the flow is laminar or turbulent, because 2D turbulence does not cascade energy to smaller scales. A simple regime in which diffusion may be negligible is the one in which the flow is laminar. This requires a Reynolds number $\mathrm {Re}\equiv UH/\nu \lesssim 1$, where $U$ is the characteristic velocity developed during relaxation and $H$ the stratification height, which we take to be the characteristic outer scale of the flow (this being necessary for complete relaxation). By (3.8), $\delta \rho \sim \rho E_{\mathrm {avail}}/E_0$, so the net gravitational force on a fluid element is $g\delta \rho \sim \rho c^2 E_{\mathrm {avail}}/HE_0$, where $c$ is the characteristic speed of compressive waves (assumed constant here for simplicity). Balancing this with the viscous force $\rho \nu U/H^2$ to estimate $U\sim Hc^2 E_{\mathrm {avail}}/\nu E_0$, we find that relaxation is laminar if

(4.1)\begin{equation} \nu \gtrsim Hc\sqrt{\frac{E_{\mathrm{avail}}}{E_0}}.\end{equation}

If (4.1) is satisfied, the flow turns over at most once before it decays, so $s$ and $\chi$ are not mixed to smaller scales. They are therefore well conserved provided that

(4.2)\begin{equation} \mathrm{Pr}_t,\ \mathrm{Pr}_m\gg 1,\end{equation}

where $\mathrm {Pr}_{t}\equiv \nu /K$ and $\mathrm {Pr}_{m}\equiv \nu /\eta$ are the thermal and magnetic Prandtl numbers, respectively.

For $\mathrm {Re}\gg 1$, the outer-scale flow has velocity ${U\sim (E_{\mathrm {avail}}/E_0)^{1/2}c}$ and is turbulent. It turns over $\mathrm {Re}$ times before decaying (in two dimensions), so $s$ and $\chi$ are mixed to the scale ${l\sim H\exp (-\mathrm {Re})}$. Diffusion at this scale can be neglected if its time scale is longer than the decay time of the turbulence, i.e. if

(4.3)\begin{equation} \ln \mathrm{Pr}_{t},\ \ln \mathrm{Pr}_{m} \gg \mathrm{Re}\sim \frac{cH}{\nu} \sqrt{\frac{E_{\mathrm{avail}}}{E_0}}.\end{equation}

In addition to mixing by the flow at scale $H$, Rayleigh–Taylor instability at interfaces of fluid with different densities may generate small-scale vortices that mix the fluid (this effect was evident in the ${t=70}$ panel of figure 2). The fastest-growing Rayleigh–Taylor mode, which develops at scale $L_{\mathrm {RT}}$, is limited by viscosity: its growth rate is ${\gamma _{\mathrm {RT}}\sim (gL_{\mathrm {RT}}\delta \rho /\rho )^{1/2}\sim \nu / L_{\mathrm {RT}}^2}$, whence ${L_{\mathrm {RT}}\sim \mathrm {Re}^{-2/3} H}$ and ${\gamma _{\mathrm {RT}}\sim \mathrm {Re}^{1/3} U/H}$. Taking the nonlinear turnover rate of the developed vortex to be equal to $\gamma _{\mathrm {RT}}$, this mode turns over $\mathrm {Re}^{4/3}$ times before the outer-scale turbulence decays (at which point we assume that no Rayleigh–Taylor-unstable interfaces remain). An argument similar to the one that led to (4.1) indicates that the Rayleigh–Taylor vortices will not establish diffusion-scale structure in $s$ and $\chi$ provided that

(4.4)\begin{equation} \ln \mathrm{Pr}_{t},\ \ln \mathrm{Pr}_{m} \gg \mathrm{Re}^{4/3}\sim \left(\frac{cH}{\nu}\right)^{4/3} \left(\frac{E_{\mathrm{avail}}}{E_0}\right)^{2/3}.\end{equation}

Equation (4.4) is a stricter criterion than (4.3); because of their faster turnover rate than the flow at scale $H$, the Rayleigh–Taylor-generated vortices are more effective at mixing.

If (4.4) is satisfied [or if (4.1) and (4.2) are, for the case of $\mathrm {Re}\lesssim 1$], the relaxation flow decays before $s$ and $\chi$ diffuse via thermal conduction or ohmic heating. Provided that the initial destabilisation was sufficiently thorough (so that the system does not become trapped in a new metastable state), we expect the final static state of the system to be the one with smallest potential energy subject to fluid-element-wise conservation of $s$ and $\chi$. Because the amount of heat per unit mass generated by viscosity is $U^2\sim c^2 E_{\mathrm {avail}}/E_0\ll c^2$, the fractional change in $s$ of each fluid element due to viscous heating during relaxation is small. Thus, the final state reached by the system ought to be, to first approximation, the states of minimum energy calculated in § 3. In the next section, we verify with numerical simulations that this is indeed the case.

4.2 Numerical results

Figure 17 visualises the relaxation of the equilibrium defined by (2.30) after the application of an impulsive force that accelerates the fluid to a velocity

(4.5)\begin{equation} \boldsymbol{u} = u_0 \boldsymbol{\hat{z}}\sin\left(\frac{2{\rm \pi} x}{L_x}\right)\exp\left(-\frac{(z-z_0)^2}{\Delta z^2}\right),\end{equation}

where $L_x = 2$ is the size of the simulation domain in the $x$ direction, $u_0=0.1$, $z_0=1.0$ and $\Delta z = 0.5$ (see § 2.6 for an explanation of our system of units). This corresponds to an initial kinetic energy $E_{\mathrm {kin,}0} \simeq 0.3 E_{\mathrm {avail}}$, where the available energy $E_{\mathrm {avail}}\sim 10^{-3}E_0$. The kinematic viscosity is $1.6\times 10^{-3}$ in these units, so the Reynolds number of the flow at the initial time is $\mathrm {Re}\sim u_0 L_x/\nu \sim 10^{2}$. The magnetic and thermal Prandtl numbers are $\mathrm {Pr}_m = 400$ and $\mathrm {Pr}_t = 670$, respectively. We provide further details of the numerical set-up in Appendix A.

Figure 17. The relaxation of the equilibrium defined by (2.30) at $\mathrm {Re}\sim 10^2$. The initial velocity field is given by (4.5). Upper panels show the evolution of $\ln (s/\chi )$ in $x$–$z$ space. Lower panels show the same quantity but sorted horizontally at each $z$, with contours of the theoretical minimum-energy state (figure 7d) overlaid.

Although $\mathrm {Re}\sim 100>1$, figure 17 shows that $\mathrm {Re}$ is insufficiently large for turbulence to develop, either at the outer scale or driven by Rayleigh–Taylor instability (which does nonetheless lead to the development of structure at scales smaller than $L_x\sim H$). Thus, relaxation takes place without significant diffusion of $s$ and $\chi$ because the Prandtl numbers are large (4.2). We observe that the upwards plume generated by the initial impulse forms a long-lived 2D state (upper panels of figure 17), which is indeed consistent with the minimum-energy state obtained in § 3.3.1 (see lower panels of figure 17).

In order to assess the sensitivity of the final state to the initial perturbation, we visualise in figure 18 the late-time state developed by simulations identical to the one shown in figure 17, but with different values of the initial kinetic energy. Specifically, we choose $u_0=0.05$ (centre panel) and $u_0=0.025$ (left panel), so that $E_{\mathrm {kin,0}}\simeq 0.1 E_{\mathrm {avail}}$ and $E_{\mathrm {kin,0}}\simeq 0.02 E_{\mathrm {avail}}$, respectively. We observe that there is some sensitivity to the amplitude of the initial perturbation: somewhat more material is displaced upwards at $E_{\mathrm {kin,0}}\simeq 0.3 E_{\mathrm {avail}}$ than $E_{\mathrm {kin,0}}\simeq 0.02 E_{\mathrm {avail}}$. This weak, but measurable, dependence of the final state on initial conditions is despite the fact that the equilibrium is initially at marginal linear stability, so there is no potential barrier to be overcome in order to trigger instability. However, partial relaxation stabilises the atmosphere (see Appendix D), so that, while the first magnetic-flux tube to move upwards experiences no potential barrier, later ones do.

Figure 18. The distribution of $\ln (s/\chi )$ at $t=300$ for simulations analogous to the one visualised in figure 17 but for three different values of $u_0$. As in figure 17, upper panels visualise the state of the simulation, while lower panels are sorted horizontally and overlaid with the 2D minimum-energy state (figure 7d).

In Appendix F, we present analogous simulations of viscous relaxation for the equilibrium defined by (2.31) (i.e. metastable to downwards perturbations). The results are qualitatively similar to those presented in this section.

5.1 Equivalence of 1D and 2D microstates

Even through the relaxing system explores 2D non-equilibrium microstates, maximising the number of them that are consistent with a given macrostate (i.e. the multiplicity of the macrostate) is formally equivalent to doing so for the equilibria obtainable by 1D rearrangements (similar to those we used to find minimum-energy states in § 3). This is because, for every 2D non-equilibrium state in horizontal pressure balance, there exists a distinct 1D equilibrium state with the same total energy and mass of fluid with each value of $s$ and $\chi$, as we now demonstrate.

We define the set of 2D microstates as follows. We partition the initial equilibrium into horizontal slices of fixed width $\Delta z$, which we further subdivide into parcels of equal mass $\Delta m$. The number of parcels in each horizontal slice is not fixed, and parcels may ‘spill’ over into adjacent slices if the total mass in the slice is not an integer multiple of $\Delta m$, although the fraction of parcels that do vanishes as $\Delta m \to 0$. We consider labelling each parcel by an index that increases along each slice from left to right, starting from the lowest slice and then moving upwards. Then, the full set of microstates is the set of possible permutations of these indices, i.e. the set of 1D shuffles of parcels with fixed $s$ and $\chi$, assembled in 2D space as described (see figure 19a,b).

Figure 19. At $\mathrm {Re}\gg 1$, turbulence mixes the advected scalars $s$ and $\chi$ [(a); this is a subsection of the state visualised in the $t=70$ panel of figure 23]. Panel (b) is a copy of (a) but with the fluid discretised into parcels of equal mass. The 2D non-equilibrium states obtainable by shuffling these parcels while maintaining horizontal pressure balance constitute microstates in our theory. As explained in the main text, there exists a bijection between such 2D states and 1D static equilibria with the same energy (c). We may therefore take microstates to be 1D equilibria.

A given permutation of parcel indices is not a complete description of the microstate, as the pressure $P$ in each parcel remains to be specified. As noted above, the 2D microstates are not equilibria: horizontal density variations induce baroclinic torques. On the other hand, we do not expect significant horizontal variation in pressure: because ${E_{\mathrm {avail}}/E_0\ll 1}$, the flow that develops during relaxation is subsonic, i.e. $U/c\ll 1$. Fluctuations $\delta P$ of the total pressure $P$ about its horizontal mean $P_0$ are therefore small, ${\delta P/P_0\sim U^2/c^2 \ll 1}$. We can evaluate $P_0(z)$ by integrating the $z$-component of the momentum equation (2.1) over $x$ and from the given $z$ to $z=\infty$, neglecting inertial terms (which correspond to the pressure fluctuations). This yields $P_0= m g$, where $m$ is the total mass of fluid above height $z$ per unit horizontal length, given by

(5.1)\begin{equation} m = \frac{1}{L_x}\int_0^{L_x} \mathrm{d}\kern0.07em x \int^{\infty}_z\mathrm{d} z' \rho(x,z'),\end{equation}

where $L_x$ is the horizontal extent of the system.

To leading order in $\Delta z$, (5.1) evaluated at a given parcel is the total mass of parcels that have greater indices. Consequently, if we were to rearrange the parcels to be stacked vertically in order of their indices (preserving the $P$, $s$ and $\chi$ from the 2D state), $P_0$ is the pressure they would have in equilibrium. Restacking in this way produces no change in the total energy as $\Delta z, \Delta m \to 0$, so the total energy of both states can be evaluated using (3.4). The leading-order term in $\delta P$ in the integrand is ${\sim c^2 \delta P^2/P_0^2 \sim (U^2/c^2)U^2\ll U^2}$, so the contribution of $\delta P$ to the energy of the state can be neglected. Thus, there exists a correspondence between 2D non-equilibrium (but horizontally pressure-balanced) states and 1D equilibria: for every 2D state, we can find a 1D equilibrium state with the same energy, by the process described above (visualised in figure 19c). Thus, in what follows, we consider microstates to be 1D equilibrium states, with energy given by (3.6).

5.2 Lynden-Bell statistical mechanics of MHD atmospheres

As explained in § 5.1, we may construct our statistical mechanics on the space of 1D equilibria, this being fully equivalent to doing so on the space of 2D non-equilibria in horizontal pressure balance. We follow the formulation of Reference Lynden-BellLynden-Bell's (Reference Lynden-Bell1967) statistical mechanics of distinguishable particles with an exclusion principle – originally derived for collisionless stellar systems and plasma – by Chavanis (Reference Chavanis2003).

We introduce $\mathcal {P}(m, \mu )\,\mathrm {d} \mu$ as the probability of finding the material with initial supported mass in the range ${[\mu, \mu +\mathrm {d} \mu ]}$ to have a supported mass of $m$ in the final state. We obtain $\mathcal {P}(m, \mu )$ by maximising the number of microstates with which it is consistent after coarse graining. This corresponds to maximising the mixing entropy (Robert & Sommeria Reference Robert and Sommeria1991)

(5.2)\begin{equation} S ={-} \int \mathrm{d} m \int \mathrm{d} \mu \mathcal{P}(m, \mu) \ln \mathcal{P}(m, \mu).\end{equation}

We maximise $S$ subject to the constraints of fixed total probability (i.e. the normalisation of $P$)

(5.3)\begin{equation} \int \mathrm{d} \mu \mathcal{P}(m, \mu) = 1,\quad\forall\ m;\end{equation}

fixed potential energy $E_{\mathrm {pot}}$

(5.4)\begin{equation} \int \mathrm{d} m \int \mathrm{d} \mu \mathcal{E}(m,\mu) \mathcal{P}(m, \mu) = E_{\mathrm{pot}};\end{equation}

and fixed mass of fluid with each value of $\mu$

(5.5)\begin{equation} \int \mathrm{d} m \mathcal{P}(m, \mu) = 1, \quad\forall\ \mu.\end{equation}

This constrained maximisation of $S$ is equivalent to unconstrained maximisation of

(5.6)\begin{align} & -\int \mathrm{d} m \int \mathrm{d} \mu \mathcal{P}(m, \mu) \ln \mathcal{P}(m, \mu) - \beta_T \int\mathrm{d} m\lambda(m) \left[\int \mathrm{d} \mu \mathcal{P}(m, \mu) - 1\right] \nonumber\\ & \quad - \beta_T \left[\int \mathrm{d} m\int \mathrm{d} \mu \, \mathcal{E}(m,\mu) \mathcal{P}(m, \mu) - E_{\mathrm{pot}}\right] - \beta_T \int \mathrm{d} \mu \psi(\mu) \left[\int \mathrm{d} m \mathcal{P}(m, \mu) - 1\right], \end{align}

over the probability $\mathcal {P}(m,\mu )$ and the Lagrange multipliers $\beta _{T}$, $\lambda (m)$ and $\psi (\mu )$. The solution is

(5.7)\begin{equation} \mathcal{P}(m, \mu) = {\rm e}^{-\beta_T [\mathcal{E}(m, \mu)-\psi(\mu)-\lambda(m)]},\end{equation}

where the Lagrange multipliers $\beta _T$ (the thermodynamic beta, to be identified with the inverse of the statistical mechanical temperature), $\psi (\mu )$ and $\lambda (m)$ are determined from the constraints (5.3)–(5.5).Footnote ¹¹

In Appendix H, we prove that, as $\beta _T\to \infty$, $\mathcal {P}(m,\mu )$ becomes increasingly sharply peaked around the solution to the LSA problem of § 3 (see figure 20). Thus, the minimum-energy states are the $\beta _T\to \infty$ limit of our statistical mechanics. Mathematically, this happens because the exponent of (5.7) has the form of a modified cost matrix (see (3.11)) multiplied by $\beta _T$; as $\beta _T\to \infty$, $\mathcal {P}(m,\mu )$ vanishes except in the vicinity of the zeros of the normal form of the cost matrix (see § 3.2), which represent the optimal assignment.

Figure 20. Convergence of $\mathcal {P}(m,\mu )$ ((5.7)) to the solution of the LSA problem (black circles; see § 3.2) as $\beta _T E_{\mathrm {avail}}\to \infty$ for the case of the unstable-upwards profile (2.30) with all integrals discretised at scale $\Delta m = 0.01$. Contours of the normal from of the cost matrix $\tilde {\mathcal {E}}_{ij}$ visualised in figure 6 are plotted in white.

The algorithm we use to evaluate $\mathcal {P}(m,\mu )$ numerically for given $E_{\mathrm {pot}}$ is analogous to the one proposed by Ewart, Nastac & Schekochihin (Reference Ewart, Nastac and Schekochihin2023). A brief summary is as follows. First, we choose a trial value of $\beta _T$ and, with suitable discretisation for $m$ and $\mu$, calculate

(5.8)\begin{equation} \mathcal{P}(m, \mu) = \frac{{\rm e}^{-\beta_T [\mathcal{E}(m, \mu)-\psi(\mu)]}}{\displaystyle\int \mathrm{d} \mu'\, {\rm e}^{-\beta_T [\mathcal{E}(m, \mu')-\psi(\mu')]} },\end{equation}

where we obtain $\psi (\mu )$ from the iterative formula

(5.9)\begin{equation} {\rm e}^{-\beta_T \psi_{n+1} (\mu)} = \int \mathrm{d} m \frac{{\rm e}^{-\beta_T \mathcal{E}(m, \mu)}}{\displaystyle\int \mathrm{d} \mu'\, {\rm e}^{-\beta_T [\mathcal{E}(m, \mu')-\psi_n(\mu')]}}.\end{equation}

Equation (5.9) is the result of integrating (5.8) over $m$ and using (5.5). Equation (5.8) satisfies the constraints (5.5) and (5.3), but does not necessarily correspond to the correct energy $E_{\mathrm {pot}}$; in this case, we increment $\beta _T$ and repeat the procedure described until the desired energy is obtained.

5.3 Diffusion

The function $\mathcal {P}(m,\mu )$ gives the fractional abundances of parcels from supported mass $\mu$ at new supported mass $m$ in the ‘most mixed’ state accessible by ideal rearrangements. We shall use these abundances to determine the result of diffusion: when sufficiently small scales are developed by mixing, diffusion homogenises nearby fluid parcels. Thus, diffused states may be obtained by taking suitable moments of $\mathcal {P}(m,\mu )$.

The standard method for deriving predictions from Lynden-Bell probability distribution functions is to argue that, due to the presumed stochastic nature of the underlying microstate, physically measurable quantities correspond to expectation values. In our case, these are

(5.10)\begin{gather} \langle s \rangle \equiv \int \mathrm{d} \mu s(\mu) \mathcal{P}(m, \mu), \end{gather}

(5.11)\begin{gather}\langle \chi \rangle \equiv \int \mathrm{d} \mu \chi(\mu) \mathcal{P}(m, \mu). \end{gather}

Unlike in the traditional contexts, however, expectation values – in particular, (5.10) – are unsuitable as a model of the state that develops after diffusion acts. This is because diffusive processes do not preserve the mean value of $s$; instead, thermal conduction and ohmic heating increase thermal entropy.Footnote ¹² It is readily verified that energy is not conserved under a collapse of the distributions of $s$ and $\chi$ onto their expectation values

(5.12)\begin{equation} \int \mathrm{d} m \mathcal{E}(mg, \langle s\rangle, \langle \chi \rangle) \neq E_{\mathrm{pot}},\end{equation}

because $\mathcal {E}(m,\mu )$ is nonlinear in $s$ and $\chi$. In the terminology of Chavanis (Reference Chavanis2003), energy is a ‘fragile integral’ – its value is different depending on whether it is computed using the coarse- or fine-grained distributions of $s$ and $\chi$. The difference between the left- and right-hand sides of (5.12) is typically much greater than the available energy of the original equilibrium – around ten times greater in the case of the unstable-upwards profile defined by (2.30).

Equation (5.10) being unsuitable, we instead determine the entropy function after diffusion, $\bar {s}$, from the conservation of energy, i.e. from

(5.13)\begin{equation} \mathcal{E}(mg,\bar{s},\bar{\chi}) = \int \mathrm{d} \mu \mathcal{E}(mg,s(\mu),\chi(\mu)) \mathcal{P}(m,\mu),\end{equation}

with the diffused magnetic flux $\bar {\chi }$ given by

(5.14)\begin{equation} \bar{\chi}=\langle\chi\rangle;\end{equation}

(the diffusion of straight magnetic-field lines does conserve magnetic flux). We plot $\bar {s}$ and $\bar {\chi }$ against $z$ (with density ${\rho (m,\bar {s} ,\bar {\chi } )}$) with a solid cyan line in figure 21. For comparison, we plot $\langle s \rangle$ and $\langle \chi \rangle$ against $z$ (with density ${\rho (m,\langle s \rangle,\langle \chi \rangle )}$) with a gold dashed line – these profiles are appreciably different, and we will find in § 6 that the former is indeed a better predictor of numerical simulations.

Figure 21. The predictions of the Lynden-Bell statistical mechanics (§ 5.2) for each of the profiles described in § 2.6. Panel (a) corresponds to (2.30), (b) to (2.31), (c) to (2.32) and (d) to (2.33). In each case, the black dashed line corresponds to the initial profile, the gold dashed line to $\langle s \rangle$ or $\langle \chi \rangle$ [see (5.10) and (5.11)] and the cyan solid line to the predictions $\bar {s}$ and $\bar {\chi }$ for the result of diffusion [see (5.13) and (5.14)]. Other coloured lines correspond to iterations of the statistical mechanical calculation, as described in § 5.4.

5.4 Secondary relaxation

Equations (5.13) and (5.14) need not, in general, represent the final state reached by relaxation, because the turbulent mixing flow remains present even after diffusion. This flow can cause further reorganisation if the post-diffusion state of the system has more than one accessible microstate, i.e. if it is not a nonlinearly stable minimum-energy state.Footnote ¹³ Remarkably, this often turns out to be the case: diffusion [in the sense of (5.13) and (5.14)] of the state predicted by Lynden-Bell statistical mechanics tends to produce new states that are unstable to further (ideal) dynamics. This phenomenon has no analogue in the relaxation of collisionless stellar systems and plasma, for which the coarse-grained probability-distribution function is always a state of minimum energy with respect to rearrangements of phase space. That it is possible for MHD atmospheres is consequence of the fact that diffusion produces changes in buoyancy. A well-known example of this is the phenomenon of ‘buoyancy reversal’ in the terrestrial atmosphere: the nonlinear dependence of density on the advected Lagrangian invariants (in that context, the mixing ratio of water and potential temperature; see Appendix B) means that a buoyant parcel of fluid that rises and mixes with denser ambient fluid can become denser than the ambient fluid, and sink as a result (see, e.g. Stevens Reference Stevens2005).

The 1D equilibrium states given by (5.13) and (5.14) may be linearly or nonlinearly unstable (metastable). Examples of each case are illustrated in figure 22. Panel (a) plots the force (2.5) per unit mass on a small parcel of fluid displaced from height $z_1$ to $z_2$ for the post-diffusion state of the metastable-upwards equilibrium (2.30). This state is linearly unstable between $z \simeq 1.4$ and $z \simeq 2.3$, and also close to $z=3$. Panel (b) is analogous, but for a slightly larger energy of $E = E_0+0.3 E_{\mathrm {avail}}$ in (5.4) (this is the initial energy of the numerical simulation to be presented in figures 23 and 24 in § 6.1). In this case, the new profile is linearly stable, but is unstable nonlinearly (metastable), and turns out to have states with lower energy (as can be confirmed by solving its LSA problem). Likewise, the post-diffusion states corresponding to the ‘bi-directional’ ((2.32)) and ‘overturning’ ((2.33)) profiles are not minimum-energy states. On the other hand, it turns out that diffusion does yield a minimum-energy state in the case of the unstable-downwards profile (2.31).

Figure 22. The force (2.5) per unit mass of a small fluid parcel moved in pressure balance and without diffusion from height $z_1$ to $z_2$ for: left panel, the state that corresponds to the cyan lines in figure 21a; and, right panel, the analogue of this profile for $E=E_0+0.3 E_{\mathrm {avail}}$. The left panel exhibits linear instability, the right panel nonlinear instability (metastability).

Figure 23. Numerical simulation of the relaxation of the equilibrium defined by (2.30) at $\mathrm {Re}\sim 10^5$. The quantity plotted is $\ln (s/\chi )$. The initial velocity field is given by (4.5) with $u_0 = 0.1$, $z_0 = 1.0$ and $\Delta z = 0.5$, which corresponds to ${E_{\mathrm {kin},0}\simeq 0.3 E_{\mathrm {avail}}}$. A movie version of this figure is available at https://doi.org/10.1017/S0022377824001521.

Figure 24. Horizontally averaged profiles of $s$ and $\chi$ plotted at intervals of $1$ code time unit, between $t=0$ (blue) to $t=600$ (red), for the simulation visualised in figure 23. Also plotted are the profiles that correspond to expectation values of $\mathcal {P}(m,\mu )$ (gold dashed line; which we claim is not a suitable model of diffusion), from (5.13) and (5.14) (cyan dashed line; these model energy- and flux-conserving diffusion), and after iterating the statistical mechanical prediction from the profile based on the cyan dashed line (pink dashed line). For reference, we also plot the minimum-energy state: this is as shown in figure 5.

As noted above, in cases where the new state given by (5.13) and (5.14) has more than one accessible microstate with the same energy, continued turbulent mixing can reshuffle fluid parcels and produce further diffusion until a minimum-energy state is reached. A simple model of this process is to apply the procedure described in §§ 5.2 and 5.3 iteratively. This corresponds to the diffused system exploring the full space of states that are energetically accessible under ideal arrangements before diffusing again. We show in figure 21 the profiles of $\bar {s}$ and $\bar {\chi }$ at the second, third and fourth iterations, plotted in pink, green and purple, respectively. The procedure terminates (i.e. reaches a minimum-energy state) after two iterations in the unstable-upwards (2.30) and ‘overturning’ (2.32) cases [see panels (a) and (d)], and after four iterations in the unstable-‘bi-directional’ case (2.33). [Note that, in the ‘overturning’ case, we find at the fourth iteration that the state is sufficiently close to the ground state to become sensitive to the discretisation that we employ to compute the integrals in (5.11) and (5.13): see the jagged structure of the purple line in figure 21c. At the fifth iteration, $\beta _T$ becomes so large as to preclude accurate computation of the relevant integrals, so we terminate the process at the fourth iteration.]

The difference in the profiles of $\bar {s}$ and $\bar {\chi }$ between the first and last stages of the iterative procedure turn out to be slight. Therefore, if, as in reality, secondary relaxations are incomplete or diffusion occurs concurrently with them, the final state reached ought not to be very different. We shall see in § 6 that accounting for these diffusive rearrangements is, in practice, a precision overkill – greater discrepancies between numerical experiment and the theoretical prediction arise which appear to be a result of the tendency of relaxing profiles to become ‘stuck’ in other metastable states. Nonetheless, it is interesting to note that, on a qualitative level, the chief outcome of the iterative procedure is the formation of plateaus (corresponding to thorough mixing and diffusion in regions where the first relaxed state is close to marginal linear stability). In the case of the unstable-upwards profile, for example, we see from figure 21a that, between the first and second iterations, material at $z\lesssim 1.75$ moves upwards to settle in the range $1.75 \lesssim z \lesssim 2.75$, producing a flatter region between $2.5 \lesssim z \lesssim 3.0$ and in the vicinity of $z\simeq 2$. Similar plateaus are observed in panels (c) and (d), with panel (c) resembling a staircase. The intriguing possibility of modelling the formation of staircases (which are observed in myriad diffusing systems in geo- and astro-physical contexts) statistical mechanically is a topic to which we shall return in future work (see § 7.2).

6.1 Metastable upwards, (2.30)

We first consider the case of the unstable-upwards profile defined by (2.30). Figure 23 visualises the distribution of $s/\chi$ (the advected scalar that controls the fluid compressibility – see § 2) during the relaxation that follows perturbation by a velocity field given by (4.5) with $u_0 = 0.1$ and $z_0=1.0$ (this corresponds to an initial kinetic energy $E_{\mathrm {kin},0}\simeq 0.3 E_{\mathrm {avail}}$). We observe that a substantial plume of material rises until reaching the stable region at the top of the simulation domain, where it overturns and develops Rayleigh–Taylor instabilities (upper left panel). Under advection by the increasingly chaotic flow, small-scale structures are developed (upper middle and right panels). These structures diffuse, ultimately leading to a one-dimensional but non-homogeneous state with no fine-scale structure (lower three panels).

Figure 24 compares the horizontally averaged instantaneous profiles of $s$ and $\chi$ developed in the simulation with both $\langle s\rangle$ and $\langle \chi \rangle$ [as defined in (5.10) and (5.11)] and $\bar {s}$ and $\bar {\chi }$ [as defined in (5.13) and (5.14)]. The quantities are computed from $\mathcal {P}(m,\mu )$ calculated with $E_{\mathrm {pot}}=E_0+E_{\mathrm {kin},0}$ in (5.4).Footnote ¹⁴

We make the following observations. First, the late-time profile of $s$ is almost everywhere larger than $\langle s \rangle$ ((5.10); gold dashed line in figure 24). This validates the reasoning we used to reject (5.10) as a predictor of the final state – evidently, dissipation causes the entropy of the fluid to grow. Secondly, $\bar {s}$ and $\bar {\chi }$ computed from (5.13) and (5.14) (cyan dashed line in figure 24) constitute a very reasonable prediction of the relaxed state. The chief discrepancies are in the range $2.5 \lesssim z \lesssim 3.0$, where $\bar {s}$ and $\bar {\chi }$ are, respectively, somewhat smaller and larger than in the late-time profiles developed by the simulation. We interpret this as a consequence of the system not ‘exploring’ the full surface of constant energy in configuration space, owing to the fluid that rises from the bottom of the equilibrium becoming trapped in a metastable state at the top. In support of this interpretation, we note that, in the range $0.9 \lesssim z \lesssim 1.5$, the cyan lines somewhat over- and under-predict the simulation result, respectively, indicating that ‘too much’ material rose in the initial plume (and became stuck).

A second discrepancy between the cyan line and the simulation result is that the latter exhibits a clear plateau in the range $1.5 \lesssim z \lesssim 2.3$. On the other hand, the cyan line is nonlinearly unstable [see figure 3(b) for its force diagram] – the dashed pink line in figure 24 shows the result of taking it as the initial state for a secondary relaxation (see § 5.4). The pink line features a plateau over roughly the same range of $z$ as the one that forms in the simulation, although the predicted plateau has slightly larger $s$ (smaller $\chi$). This might be interpreted as a consequence of the fact that the large-$s$ fluid that would have risen upwards under the secondary relaxation to form the plateau was, in the simulation, already displaced to the top of the atmosphere. As a consequence, the plateau forms with somewhat smaller $s$ than it would otherwise have had.

Figures 25 and 26 are analogous to figure 24 but for simulations with ${u_0 = 0.05}$ ($E_{\mathrm {kin},0}\simeq 0.1 E_{\mathrm {avail}}$) and ${u_0 = 0.025}$ ($E_{\mathrm {kin},0}\simeq 0.02 E_{\mathrm {avail}}$), respectively. In these cases, the quantitative agreement between simulation and theory is less good than in figure 24: with a smaller initial impulse, relaxation is incomplete. Figure 27, shows that, at peak, between $50\,\%$ and $60\,\%$ of the available kinetic plus potential energy is in the form of kinetic energy for all three simulations, showing that the liberation of available potential energy during relaxation is fairly efficient.

Figure 25. As in figure 24, but for a simulation initialised with $u_0=0.05$ in (4.5) ($E_{\mathrm {kin},0}\simeq 0.1 E_{\mathrm {avail}}$).

Figure 26. As in figure 24, but for a simulation initialised with $u_0=0.025$ in (4.5) ($E_{\mathrm {kin},0}\simeq 0.02 E_{\mathrm {avail}}$).

Figure 27. Evolution of the kinetic energy as a fraction of the total available energy, which is the kinetic energy plus the available potential energy of the initial state, for the simulations visualised in figures 24–26.

6.2 Metastable downwards, (2.31)

We now report the results of analogous simulations to those in § 6.2, but for the unstable-downwards profile (2.31) and with $z_0 = 2.25$ in (4.5).

Figure 28 visualises the evolution in $x$–$z$ space for ${u_0=0.14}$ ($E_{\mathrm {kin},0}\simeq 0.2E_{\mathrm {avail}}$), following the tracer $s/\chi$. Similarly to figure 23, we observe that a descending plume reaches the stable buffer region at the bottom of the simulation domain, develops small-scale structure due to Rayleigh–Taylor instability and advection by chaotic motions, and ultimately diffuses. The evolution of the horizontally averaged profiles is displayed in figures 29–31 for the cases of $u_0 = 0.14$ ($E_{\mathrm {kin},0}\simeq 0.2E_{\mathrm {avail}}$), $u_0 = 0.1$ ($E_{\mathrm {kin},0}\simeq 0.1E_{\mathrm {avail}}$), and $u_0 = 0.05$ ($E_{\mathrm {kin},0}\simeq 0.03 E_{\mathrm {avail}}$), respectively. Again, we observe reasonable agreement between the statistical mechanical prediction and the late-time limit of the simulations, although the predictions do somewhat under-predict $s$ and over-predict $\chi$ for $z\gtrsim 1.5$. Again, the reason appears to be partial relaxation: the degree of inaccuracy is greater in the simulations with smaller initial velocity fields. We also note that the theory fails to predict the plateau that forms in the vicinity of $z\simeq 0.5$ (no secondary relaxation is possible from the diffused state that corresponds to the cyan line, as it is nonlinearly stable, see § 5.4). Nonetheless, the theoretical predictions agree reasonably well in this region with all three simulation profiles, the flatness notwithstanding. Finally, we show the evolution of the kinetic energy as a fraction of the total available energy in figure 32. As in figure 27, the relaxation is fairly efficient at liberating available potential energy. Large initial perturbations are not required to do so: over a factor $\simeq 30$ difference in the initial kinetic energy, the peak ratio of kinetic energy to available energy varies only between around $30\,\%$ and $45\,\%$.

Figure 28. Numerical simulation of the relaxation of the equilibrium defined by (2.31) at ${\mathrm {Re}\sim 10^5}$. The quantity plotted is $\ln (s/\chi )$. The initial velocity field is given by (4.5) with $u_0 = 0.14$, $z_0 = 2.25$ and $\Delta z = 0.5$; this corresponds to an initial kinetic energy $\simeq 0.2$ times the available potential energy. A movie version of this figure is available at https://doi.org/10.1017/S0022377824001521.

Figure 29. Horizontally averaged profiles of $s$ and $\chi$ plotted at intervals of $1$ code time unit, between $t=0$ (blue) to $t=600$ (red), for the simulation visualised in figure 28. Also plotted are the profiles obtained by taking expectation values of $\mathcal {P}(m,\mu )$ (gold dashed line) and from (5.13) and (5.14) (cyan dashed line).

Figure 30. As in figure 29, but for a simulation initialised with $u_0=0.1$ in (4.5), which corresponds to ${E_{\mathrm {kin},0}\simeq 0.1 E_{\mathrm {avail}}}$.

Figure 31. As in figure 29, but for a simulation initialised with $u_0=0.05$ in (4.5), which corresponds to ${E_{\mathrm {kin},0}\simeq 0.03 E_{\mathrm {avail}}}$.

Figure 32. Evolution of the kinetic energy as a fraction of the total available energy, which is the kinetic energy plus the available potential energy of the initial state, for the simulations visualised in figures 29–31, as well as for one with $u_0 = 0.025$ in (4.5).