However, in the absence of gradient terms, minimisation over the
(unconstrained) fields $J_{\rho},J_{m}$ becomes a purely local (algebraic)
problem (it is important here that the IEPR constraint (23) does not involve
gradient terms). The computation proceeds along similar lines to those taken
in the homogeneous case of Sec. V.

First, for stationary profiles the IEPR constraint in (28) becomes
$a=\int_{0}^{1}dx\alpha(x)$ with $\alpha$ defined in (24). Then, minimizing
the action over $J_{\rho}(x)$ and $m(x)$, one arrives at the algebraic
relations (32) and (39) but which are now space dependent. That is

$J_{\rho}(x)=\frac{m(x)\alpha(x)\lambda}{\rho(x)(1-\rho(x))+m(x)^{2}}\quad,\quad
J_{m}(x)=\frac{\rho(x)\alpha(x)\lambda}{\rho(x)(1-\rho(x))+m(x)^{2}}.$ (88)

where $m(x)=m\left[\alpha(x),\rho(x)\right]$ is given by $m=0$ in regions
where $\alpha(x)<a_{1}(\rho(x),\text{Pe})$ with $a_{1}$ as in (36); otherwise
$m$ is obtained by solving

$\alpha(x)=\left[\bar{a}(\rho(x))+(1-\rho(x))m(x)^{2}\right]\sqrt{1+\frac{2}{\text{Pe}^{2}}\frac{1}{(1-\rho(x))\sqrt{\rho(x)^{2}-m(x)^{2}}}}\;.$
(89)

Since this analysis is directly analogous to Sec. V, one arrives at (84) with
$I_{H}$ given by (38).

###  B.2 Deriving the sub-leading gradient terms

We now derive the sub-leading order correction to the expression (84), at
small $\ell$: that is (40). As noted in Sec. VI.2, this is derived under the
condition that reflection symmetry is unbroken, so it allows analysis of SM
and PS states but not TB and CM. Specifically, we assume that $J_{\rho}(x)=0$
(exactly) while $m(x)$ vanishes as $\ell\to 0$

$m(x)=\ell\tilde{m}(x)+\mathcal{O}\left(\ell^{2}\right).$ (90)

This implies that

$\mathbb{C}^{-1}=\frac{1}{2D}\left\\{\begin{bmatrix}\frac{1}{\rho\left(1-\rho\right)}&0\\\
0&\frac{1}{\rho}\end{bmatrix}-\frac{\ell\tilde{m}}{\rho^{2}}\begin{bmatrix}0&1\\\
1&0\end{bmatrix}+\frac{\ell^{2}\tilde{m}^{2}}{\rho^{3}}\begin{bmatrix}1&0\\\
0&1\end{bmatrix}\right\\}+\mathcal{O}\left(\ell^{3}\right)$ (91)

and the mean currents appearing in the action functional (12) are

$\displaystyle\frac{\bar{J}_{\rho}}{\lambda}$ $\displaystyle=$
$\displaystyle\ell\left[-\frac{\partial_{x}\rho}{\text{Pe}}+\tilde{m}\left(1-\rho\right)\right],$
$\displaystyle\frac{\bar{J}_{m}}{\lambda}$ $\displaystyle=$
$\displaystyle\rho\left(1-\rho\right)-\ell^{2}\frac{\partial_{x}\tilde{m}}{\text{Pe}}.$
(92)

As $J_{\rho}=0$ we have that the local IEPR field is

$\alpha(x)=\frac{J_{m}\left(1-\rho\right)}{\lambda},$ (93)

and the constraint to be satisfied by the IEPR in (28) is

$a=\int_{0}^{1}dx\,\alpha(x).$ (94)

####  B.2.1 Action functional sub-leading correction

We now compute the correction at $O(\ell^{2})$ to the action, which is to be
minimised in (28). First, from the dynamical equation for $m$ (5) we have

$K=-\frac{\partial_{x}J_{m}}{2\gamma}=-\frac{\text{Pe}\,\ell}{2}\partial_{x}\left(\frac{\alpha}{1-\rho}\right).$
(95)

where the second equality is (93). Thus, to leading order in $\ell$ we find
that

$\displaystyle\mathcal{L}_{K}$ $\displaystyle=$
$\displaystyle\gamma\ell^{2}\frac{\left[\frac{\text{Pe}}{2}\partial_{x}\left(\frac{\alpha}{1-\rho}\right)+\tilde{m}\right]^{2}}{2\rho}+\mathcal{O}\left(\ell^{4}\right).$
(96)

Using (92), (93), and that $J_{\rho}=0$, then the other contribution to the
action (12) gives

$\displaystyle\frac{\mathcal{L}_{J}}{\gamma}$ $\displaystyle=$ $\displaystyle
I_{H}\left[\rho(x),\alpha(x)\right]+\ell^{2}\frac{\text{Pe}^{2}}{4}\left\\{\left[\frac{\left[\partial_{x}\rho/\text{Pe}-\tilde{m}\left(1-\rho\right)\right]^{2}}{\rho\left(1-\rho\right)}+\frac{2}{\text{Pe}}\left(\frac{\alpha-\bar{a}}{1-\rho}\right)\left[\partial_{x}\left(\frac{\tilde{m}}{\rho}\right)+\text{Pe}\frac{\tilde{m}^{2}\left(1-\rho\right)}{\rho^{2}}+\frac{\text{Pe}}{2}\frac{\tilde{m}^{2}}{\rho^{3}}\left(\frac{\alpha-\bar{a}}{1-\rho}\right)\right]\right]\right\\}$
(97) $\displaystyle+$ $\displaystyle\mathcal{O}(\ell^{3}).$

Using integration by parts one can eliminate the derivatives of $\tilde{m}$ to
arrive at

$\displaystyle\frac{\mathcal{L}_{J}}{\gamma}$ $\displaystyle=$ $\displaystyle
I_{H}\left[\rho(x),\alpha(x)\right]+\ell^{2}\frac{\text{Pe}^{2}}{4}\left\\{\left[\frac{\left[\partial_{x}\rho/\text{Pe}-\tilde{m}\left(1-\rho\right)\right]^{2}}{\rho\left(1-\rho\right)}+\left(\frac{\alpha-\bar{a}}{1-\rho}\right)\left[2\frac{\tilde{m}^{2}\left(1-\rho\right)}{\rho^{2}}+\frac{\tilde{m}^{2}}{\rho^{3}}\left(\frac{\alpha-\bar{a}}{1-\rho}\right)\right]-\frac{2}{\text{Pe}}\partial_{x}\left(\frac{\alpha-\bar{a}}{1-\rho}\right)\frac{\tilde{m}}{\rho}\right]\right\\}$
(98) $\displaystyle+$ $\displaystyle\mathcal{O}(\ell^{3}).$

Putting everything together, we have the total action is given by

$\displaystyle\frac{\mathcal{A}}{\gamma T}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}dx\left\\{I_{H}\left[\rho(x),\alpha(x)\right]+\ell^{2}\left\\{\frac{1}{4}\frac{\left[\partial_{x}\rho-\text{Pe}\,\tilde{m}\left(1-\rho\right)\right]^{2}}{\rho\left(1-\rho\right)}+\frac{\left[\text{Pe}\,\partial_{x}\left(\frac{\alpha}{1-\rho}\right)+2\tilde{m}\right]^{2}}{8\rho}\right\\}\right\\}$
(99) $\displaystyle+$
$\displaystyle\ell^{2}\frac{\text{Pe}^{2}}{4}\int_{0}^{1}dx\left\\{\left(\frac{\alpha-\bar{a}}{1-\rho}\right)\left[2\frac{\tilde{m}^{2}\left(1-\rho\right)}{\rho^{2}}+\frac{\tilde{m}^{2}}{\rho^{3}}\left(\frac{\alpha-\bar{a}}{1-\rho}\right)\right]-\frac{2}{\text{Pe}}\partial_{x}\left(\frac{\alpha-\bar{a}}{1-\rho}\right)\frac{\tilde{m}}{\rho}\right\\}+\mathcal{O}(\ell^{3}).$

This action is still to be minimized over the fields $\rho,\alpha,\tilde{m}$
subject to the constraints (94) and $\rho_{0}=\int_{0}^{1}dx\rho$. Note that
(99) does not depend on gradients of $\tilde{m}$, and this field does not
feature in either constraint. Hence, its optimal value is obtained by direct
minimization of (99) with respect to $\tilde{m}$. Performing the algebraic
(although cumbersome) minimization we arrive at the announced result (40) with
the symmetric matrix $\mathbb{A}(\rho,\alpha)$ whose coefficients are given by

$\displaystyle\mathbb{A}_{1,1}$ $\displaystyle=$
$\displaystyle\frac{\alpha^{2}\text{Pe}^{2}+2(1-\rho)^{3}}{8\rho(1-\rho)^{4}}-\frac{\text{Pe}^{2}}{\frac{2}{\rho}+\text{Pe}^{2}(1-\rho+\frac{\alpha^{2}}{\rho^{3}(1-\rho)^{2}})}$
$\displaystyle\mathbb{A}_{1,2}=\mathbb{A}_{2,1}$ $\displaystyle=$
$\displaystyle\text{Pe}^{2}\frac{\alpha}{8\rho(1-\rho)^{3}}\quad,\quad\mathbb{A}_{2,2}=\frac{\text{Pe}^{2}}{8\rho(1-\rho)^{2}}$

####  B.2.2 Square gradient theory in the SM state

In the SM state, the square gradient minimisation (40) reduces to the single-
field minimisation (48). To see this, note that the IEPR of the SM state obeys
(47): on plugging that result into (40), one may eventually derive (48).
However, a simpler route to the same result is to plug (47) directly into (99)
to find that

$\frac{\mathcal{A}}{\gamma
T}=\ell^{2}\left\\{\frac{1}{4}\frac{\left[\partial_{x}\rho-\text{Pe}\,\tilde{m}\left(1-\rho\right)\right]^{2}}{\rho\left(1-\rho\right)}+\frac{\left[\text{Pe}\,\left(1-2\rho\right)\partial_{x}\rho+2\tilde{m}\right]^{2}}{8\rho}\right\\}+\mathcal{O}\left(\ell^{4}\right)$
(100)

Performing the minimization with respect to $\tilde{m}$ now we arrive at

$\tilde{m}=\frac{2\rho\text{Pe}}{2+\text{Pe}^{2}\,\left(1-\rho\right)}\partial_{x}\rho,$
(101)

Plugging this expression back into (100) and combining with (40) finally
yields (48).

##  Appendix C Finding the non-convex regions in $I_{H}$ (38)

This appendix derives the regions where $I_{H}$ is not convex, as illustrated
in Fig. 5.

###  C.1 Local non-convexity

The region of local non-convexity is defined by having a negative Hessian
determinant:

$\text{Det}\left[\text{Hess}(I_{H})\right]=\frac{\partial^{2}I_{H}}{\partial{\rho_{0}}^{2}}\frac{\partial^{2}I_{H}}{\partial{a}^{2}}-\frac{\partial^{2}I_{H}}{\partial{\rho_{0}}\partial
a}<0$ (102)

In the region $a<a_{1}$, the expression (38) for $I_{H}$ is explicit. Plugging
it in (102) one finds

$\left(\frac{2}{3}-\rho_{0}\right)\left[a-\bar{a}(\rho_{0})\right]<0.$ (103)

That is, in the suppressed IEPR region $a<\bar{a}$, the rate function $I_{H}$
is locally non-convex for all densities $\rho_{0}<2/3$, while for the enhanced
region bounded by the critical value $\bar{a}<a<a_{1}$, it is locally non-
convex for $\rho_{0}>2/3$.

In the region $a>a_{1}$ the rate function $I_{H}$ is only given in a
parametric form. To find its second derivatives we use the following
identities

$\displaystyle\left.\frac{\partial I_{H}}{\partial\rho_{0}}\right|_{a}$
$\displaystyle=$ $\displaystyle\left.\frac{\partial
I_{H}}{\partial\rho_{0}}\right|_{m}+\left.\frac{\partial I_{H}}{\partial
m}\right|_{\rho_{0}}\left.\frac{\partial
m}{\partial\rho_{0}}\right|_{a}=\left.\frac{\partial
I_{H}}{\partial\rho_{0}}\right|_{m}-\left.\frac{\partial I_{H}}{\partial
m}\right|_{\rho_{0}}\frac{\left.\frac{\partial
a}{\partial\rho_{0}}\right|_{m}}{\left.\frac{\partial a}{\partial
m}\right|_{\rho_{0}}},$ $\displaystyle\left.\frac{\partial I_{H}}{\partial
a}\right|_{\rho_{0}}$ $\displaystyle=$ $\displaystyle\left.\frac{\partial
I_{H}}{\partial m}\right|_{\rho_{0}}\left.\frac{\partial m}{\partial
a}\right|_{\rho_{0}}=\left.\frac{\partial I_{H}}{\partial
m}\right|_{\rho_{0}}\frac{1}{\left.\frac{\partial a}{\partial
m}\right|_{\rho_{0}}}.$ (104)

Employing the same differentiation rules successively, one can compute the
Hessian determinant in the region $a>a_{1}$:

$\text{Det}\left[\text{Hess}(I_{H})\right]=\chi(\rho,m)[m^{2}-\rho_{0}^{2}-\rho_{0}+1],$
(105)

where the function $\chi(\rho,m)>0$ (we do not write this function explicitly,
for compactness). Thus, the determinant of the Hessian is negative for
$m^{2}-\rho_{0}^{2}-\rho_{0}+1<0$. Using the relation (39) we find that this
is the case for $a<a_{\text{non-convex}}$ with

$a_{\text{non-
convex}}=(2\rho_{0}-1)\sqrt{(1-\rho_{0})^{2}+\frac{2\sqrt{1-\rho_{0}}}{\text{Pe}^{2}}}.$
(106)

Thus, beyond the critical value $a_{1}$, the rate function is locally non-
convex in the region bounded by $a_{1}(\rho_{0})<a<a_{\text{non-
convex}}(\rho_{0})$. This concludes the mapping of the locally non-convex
regions in $I_{H}$ that are marked in orange shading in Fig. 5 (B).

###  C.2 Global non-convexity

The globally non-convex region is where $I_{H}$ is amenable to a two
dimensional common tangent construction in the $(\rho_{0},a)$ plane. This
includes the SM region where tangent construction is degenerate and the convex
hull is given by $I_{H}=0$ (see discussion at the end of Sec. VIII). The
tangent construction beyond the SM region is dominated by the point
$(\rho_{0}=1,a=0)$ where the function $I_{H}$ has a pinch point singularity,
see Fig. 6. As a consequence, the tangent construction is degenerate here as
well since all tangents comprising the common tangent construction share this
point as one of their ends 444We could not provide an algebraic proof. We
verified this by using the ConvexHullMesh command in Mathematica. We are left
with finding the location of the other end of all tangents which would
comprise the boundary of the globally non-convex region, see Fig. 7. Denote by
$(\rho_{h}=1,a_{h}=0)$ and $(\rho_{l},a_{l})$ the high and low density end
points of the tangents, then the tangent construction reads

$\left[I_{H}+(1-\rho_{l})\frac{\partial
I_{H}}{\partial\rho_{l}}-a_{l}\frac{\partial I_{H}}{\partial
a_{l}}\right]_{(\rho_{l},a_{l})}=0.$ (107)

For $(\rho_{l},a_{l})$ such that $a_{l}<a_{1}(\rho_{l})$, the expression for
$I_{H}$ is explicit (38). Substituting it into (107) one finds that the
boundary of the globally non-convex region is given by the vertical line
$\rho_{l}=1/2$ for $\bar{a}\leq a\leq
a_{1}(1/2)=(\sqrt{1+8/\text{Pe}^{2}})/8$.

For $a>a_{1}(1/2)$ one must use the parametric representation for $I_{H}$
(38). Plugging it into (107) together with the substitution for the
differentiation rules (104), one arrives at the relation

$m_{l}^{2}=\frac{\sqrt{2\rho_{l}-1}}{2}(1+\rho_{l}\sqrt{2\rho_{l}-1}).$ (108)

Plugging this back into (39) we arrive at the curve (45) which marks the
boundary of the globally non-convex region for $a>a_{1}(1/2)$. This concludes
the mapping of the globally non-convex regions in $I_{H}$ that are marked in
purple shading in Fig. 5 (B).

##  Appendix D Landau theory

We start with first identifying the critical bias $\Lambda_{c}$ beyond which a
constant solution $\rho_{0}$ to the minimization problem (56) becomes linearly
unstable. We then use this value to establish a Landau theory that captures
the growth of the spatial excitation amplitude close to this critical value.

###  D.1 Linear stability of homogeneous solutions

The linear stability of the constant solution $\rho_{0}$ is determined by
examining the second variation $\delta\tilde{\psi}$ of the functional (56)
with respect to small density variation $\rho=\rho_{0}+\delta\rho$

$\displaystyle-\delta{\psi}(\tilde{\Lambda})=\int_{0}^{1}dx\left[M(\rho_{0})(\partial_{x}\delta\rho)^{2}+\tilde{\Lambda}(2-3\rho_{0})(\delta\rho)^{2}\right].$
(109)

The first mode to destabilize the quadratic form (109) is the principal one
$\delta\rho\sim\cos(2\pi x)$. This happens as soon as:

$\tilde{\Lambda}\operatorname{sign}(3\rho_{0}-2)>|\tilde{\Lambda}_{c}|\quad\text{where}\quad\tilde{\Lambda}_{c}=\frac{4\pi^{2}M(\rho_{0})}{3\rho_{0}-2}.$
(110)

To find the critical value of the IEPR at which this instability takes place,
we consider the Legendre transform (53) and recall that the state before the
instability is homogeneous, so its SCGF is (54). Hence, the critical value of
$a$ corresponding to the critical bias $\tilde{\Lambda}_{c}$ is

$a_{2}(\rho)=\left(1+\frac{2}{\text{Pe}^{2}}\Lambda_{c}\right)\bar{a}(\rho).$
(111)

Plugging (110) and using $\tilde{\Lambda}=\Lambda/\ell^{2}$, one arrives at
(62).

###  D.2 Establishing the Landau theory close to $\Lambda_{c}$

As we have seen in the previous Appendix D.1, instability of the homogeneous
solution sets in via the principal mode excitation. To find the amplitude of
this mode one has to go beyond the leading order expansion when close to the
critical point. Such an expansion appeared in several previous works that
studied similar second order transitions within the MFT framework Lecomte  _et
al._ (2012); Baek  _et al._ (2018); Dolezal and Jack (2019); Jack  _et al._
(2015). The relevant small parameter is

$\epsilon=\tilde{\Lambda}_{c}-\tilde{\Lambda}\;.$ (112)

The mode expansion up to the sub-leading order turns out to be

$\rho=\rho_{0}+A\cos(2\pi x)+B\cos(4\pi x)+\mathcal{O}(\epsilon^{3/2})\;,$
(113)

where close to the critical point we have that the variational parameters
scale as $A\sim\sqrt{\epsilon}$ and $B\sim\epsilon$ Lecomte  _et al._ (2012);
Dolezal and Jack (2019). Plugging this expansion into (56) and keeping terms
up to order $\epsilon^{2}$ we find

$\Psi\simeq\min_{B}\left[(16\pi^{2}M_{0}-\frac{\tilde{\Lambda}_{c}}{2}\bar{a}_{0}^{\prime\prime})\frac{B^{2}}{2}+\epsilon\bar{a}_{0}^{\prime\prime}\frac{A^{2}}{4}+(3\pi^{2}M_{0}^{\prime}-\frac{\tilde{\Lambda}_{c}}{8}\bar{a}_{0}^{\prime\prime\prime})A^{2}B+(\frac{\pi^{2}}{4}M_{0}^{\prime\prime})A^{4}\right].$
(114)

Here primes denote derivative with respect to the argument, and the subscript
0 denotes evaluating at $\rho_{0}$. That is,
$M^{\prime\prime}_{0}=d^{2}M/d\rho^{2}|_{\rho_{0}}$. Plugging (110) we arrive
at

$\displaystyle\Psi\simeq\min_{B}\left[\epsilon\bar{a}_{0}^{\prime\prime}\frac{A^{2}}{4}+\frac{\mu_{2}}{2}B^{2}+\mu_{3}A^{2}B+\frac{\mu_{4}}{2}A^{4}\right],$
(115)

with

$\mu_{2}=12\pi^{2}M_{0}\quad;\quad\mu_{3}=3\pi^{2}M_{0}^{\prime}-\frac{\pi^{2}M_{0}\bar{a}_{0}^{\prime\prime\prime}}{\bar{a}_{0}^{\prime\prime}}\quad;\quad\mu_{4}=\frac{\pi^{2}}{2}M_{0}^{\prime\prime}.$
(116)

Minimizing with respect to $B$ we find

$B=-\frac{\mu_{3}}{\mu_{2}}A^{2}=-A^{2}\left[\frac{M_{0}^{\prime}}{4M_{0}}-\frac{\bar{a}_{0}^{\prime\prime\prime}}{12\bar{a}_{0}^{\prime\prime}}\right].$
(117)

Plugging this relation back into (115), we get a Landau type of expansion with
respect to the principal mode excitation amplitude $A$

$\Psi\simeq\epsilon\bar{a}^{\prime\prime}_{0}\frac{A^{2}}{4}+\frac{\beta(\rho_{0},\text{Pe})}{2}A^{4}$
(118)

with

$\beta=\mu_{4}-\frac{\mu_{3}^{2}}{\mu_{2}}=\frac{\pi^{2}}{2}M_{0}^{\prime\prime}-\frac{(3\pi^{2}M_{0}^{\prime}\bar{a}_{0}^{\prime\prime}-\pi^{2}M_{0}\bar{a}_{0}^{\prime\prime\prime})^{2}}{12\pi^{2}(\bar{a}_{0}^{\prime\prime})^{2}M_{0}}.$
(119)

For $\beta>0$, when crossing the critical point (110), the amplitude $A$ grows
as

$A=\sqrt{\frac{-\bar{a}_{0}^{\prime\prime}\epsilon}{4\beta}}.$ (120)

###  D.3 Tricritical point

As discussed in Sec. VIII.4, points with $\beta=0$ are tricritical. For any
value of Pe, the function $\beta=\beta(\rho_{0},\text{Pe})$ has at least two
roots $\rho_{c1}(\text{Pe})<2/3$ and $\rho_{c2}(\text{Pe})>2/3$. As an
example, for Pe=2 one finds $\rho_{c1}=0.600\dots$ and $\rho_{c2}=0.709\dots$,
see Fig. 15. This signals the onset of a first order transition.

Figure 15: The function $\beta(\rho_{0},\text{Pe}=2)$. The black points are
its roots, marking the value at which a second order phase transition is
overtaken by a first order one. The dashed vertical line at $\rho_{0}=2/3$
differentiates between regions where the homogeneous solution is linearly
unstable for enhanced ($\rho_{0}>2/3$) IEPR fluctuations or suppressed
($\rho_{0}<2/3$) fluctuations.

This is verified by numerically solving the minimization problem for the CGF
$\psi(\Lambda)$ (56). As shown in Fig. 16, once crossing the tricritical point
at $\rho_{c}(\text{Pe})$, a first order phase transition overtakes the second
order one at a value $|\tilde{\Lambda}_{c}^{1}|$ smaller than the prediction
of the linear stability analysis $|\tilde{\Lambda}_{c}|$ (110). Here the
amplitude of the density modulation grows discontinuously. The magnitude of
the amplitude discontinuity is increased when advancing along the first order
transition line away from the tricritical point, see Fig. 9 of the main text.

Figure 16: The optimal density profile $\rho(x)$ just after crossing the phase
transition. The left panel at $\rho_{0}=0.58$ is below the tricritical point
$\rho_{0}<\rho_{c1}$. Here the phase transition is second order and
corresponds to the theoretical prediction given by (57) together with (120).
This prediction is plotted in dashed blue line (indistinguishable from the
numerical black line). The right panel at $\rho_{0}=0.62$ is above the
tricritical point $\rho_{0}>\rho_{c1}$. Here the phase transition is first
order, and it overtakes the second order transition
$|\tilde{\Lambda}_{c}^{1}|<|\tilde{\Lambda}_{c}|$. Here the amplitude of the
density modulation grows discontinuously. In both panels $\text{Pe}=2$ and
$\ell=0.02$.

This behavior is consistent with the following model. Assuming the next order
is accounted by a cubic term with a positive pre-factor $\gamma$

$\Psi\simeq\epsilon
a^{\prime\prime}(\rho_{0})\frac{A^{2}}{4}+\frac{\beta(\rho_{0},\text{Pe})}{2}A^{4}+\frac{\gamma}{3}A^{6}.$
(121)

Then at the tricritical point $\beta=0$ the amplitude grows with the modified
power

$A=\left(-\frac{\epsilon a^{\prime\prime}_{0}}{4\gamma}\right)^{1/4}.$ (122)

We found this to be consistent with our numerical computations. Beyond the
tricritical point where $\beta<0$ we have that a first order transition occurs
before linear instability sets in at

$\tilde{\Lambda}_{c}^{1}=\tilde{\Lambda}_{c}-\frac{3}{4}\frac{\beta^{2}}{4a^{\prime\prime}_{0}\gamma}$
(123)

where the amplitude has a jump given by

$\Delta A=\sqrt{-\frac{3\beta}{4\gamma}}.$ (124)

###  D.4 Miscibility gap at $\ell\to 0$

As argued in Sec.VIII.4, the miscibility gap for the first order transition is
finite at vanishing $\ell$ as is clearly visible in the region of suppressed
IEPR in Fig. 10 (A). The miscibility gap is determined in general by
$\psi_{\text{AP}}$ which is the CGF corresponding to time independent
solutions for (52) (which are consistent with the additivity principle). In
our study we solved for its leading order expansion
$\psi(\Lambda)_{\text{AP}}=\ell^{2}\tilde{\psi}(\Lambda/\ell^{2})+o(\ell^{2})$.
Then denote by $a_{2}$ and $\tilde{a}_{c2}$ the two endpoints of the
miscibility gap. That is, $a_{2}$ is the onset of the first order transition,
marked by the orange dashed curve in Fig. 10 (A), and $\tilde{a}_{c2}$ is its
other end, denoted by the magenta curve. Then for suppressed IEPR deviations
($a<\bar{a}$), these values are given by
$a_{2}=\partial\psi_{\text{AP}}/\partial\Lambda|_{\Lambda_{c}^{+}}$ and
$\tilde{a}_{c2}=\partial\psi_{\text{AP}}/\partial\Lambda|_{\Lambda_{c}^{-}}$
with $\Lambda_{c}$ the point of the kink singularity in $\psi_{\text{AP}}$
Touchette (2009). Then the scaling for $\psi_{\text{AP}}$ implies that these
two values indeed have a finite separation at vanishing $\ell$. An equivalent
statement holds for enhanced IEPR fluctuations with $\Lambda_{c}^{+}$ and
$\Lambda_{c}^{-}$ interchanged between the two endpoints. Note however, that
the “slender” shape of the SM region at enhanced IEPR fluctuations (see Fig.
4) is such that this separation is bound to be small (yet finite).

###  D.5 Tricritical point at MIPS criticality

As discussed in Sec.X, the MIPS criticality of the unbiased system is also a
tricritical point that occurs in the biased ensemble when $\Lambda=0$. That
is, two conditions hold. First,
$\tilde{\Lambda}_{c}\left(\rho_{0}=3/4,\text{Pe}=4\right)=0$ so that the
critical bias vanishes, and also
$\beta\left(\rho_{0}=3/4,\text{Pe}=4\right)=0$ so that the MIPS criticality is
also a tricritical point for the biased system.

That these conditions are met can be read from the definition of $M(\rho)$ in
(49). It is convenient to define

$h(\rho_{0},\text{Pe})=2+\text{Pe}^{2}\,(1-\rho)(1-2\rho)$ (125)

such that $M\propto h^{2}$. From Kourbane-Houssene  _et al._ (2018), the
spinodal instability leading to MIPS is $h(\rho_{0},\text{Pe})=0$. At the
critical point one has additionally that $\partial h/\partial\rho_{0}=0$ (this
is the point of minimal density on the spinodal). Hence $M_{0}(\rho)\propto
h^{2}=0$ at the MIPS critical point and similarly $M_{0}^{\prime}(\rho)\propto
h(\partial h/\partial\rho)=0$ and $M_{0}^{\prime\prime}(\rho)\propto(\partial
h/\partial\rho)^{2}+h(\partial^{2}h/\partial\rho^{2})=0$ there. (Recall that
$M_{0}$ depends implicitly on Pe as well as on $\rho_{0}$.) Since that $M_{0}$
and its first two derivatives are all vanishing at the MIPS critical point,
one sees from (119) that $\beta=0$ there too.

## References

  * Marchetti  _et al._ (2013) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, Madan Rao, and R. Aditi Simha, “Hydrodynamics of soft active matter,” Reviews of Modern Physics 85, 1143–1189 (2013).
  * Bechinger  _et al._ (2016) Clemens Bechinger, Roberto Di Leonardo, Hartmut Löwen, Charles Reichhardt, Giorgio Volpe, and Giovanni Volpe, “Active Particles in Complex and Crowded Environments,” Reviews of Modern Physics 88, 045006 (2016).
  * Fodor and Cristina Marchetti (2018) Étienne Fodor and M. Cristina Marchetti, “The statistical physics of active matter: From self-catalytic colloids to living cells,” Physica A: Statistical Mechanics and its Applications 504, 106–120 (2018).
  * Cates and Tailleur (2015) Michael E. Cates and Julien Tailleur, “Motility-Induced Phase Separation,” Annual Review of Condensed Matter Physics 6, 219–244 (2015).
  * Caprini  _et al._ (2019) Lorenzo Caprini, Umberto Marini Bettolo Marconi, Andrea Puglisi, and Angelo Vulpiani, “The entropy production of Ornstein–Uhlenbeck active particles: A path integral method for correlations,” Journal of Statistical Mechanics: Theory and Experiment 2019, 053203 (2019).
  * Dabelow  _et al._ (2019) Lennart Dabelow, Stefano Bo, and Ralf Eichhorn, “Irreversibility in Active Matter Systems: Fluctuation Theorem and Mutual Information,” Physical Review X 9, 021009 (2019).
  * Ganguly and Chaudhuri (2013) Chandrima Ganguly and Debasish Chaudhuri, “Stochastic thermodynamics of active Brownian particles,” Physical Review E 88, 032102 (2013).
  * Shankar and Marchetti (2018) Suraj Shankar and M. Cristina Marchetti, “Hidden entropy production and work fluctuations in an ideal active gas,” Physical Review E 98, 020604 (2018).
  * Fodor  _et al._ (2016) Étienne Fodor, Cesare Nardini, Michael E. Cates, Julien Tailleur, Paolo Visco, and Frédéric van Wijland, “How Far from Equilibrium Is Active Matter?” Physical Review Letters 117, 038103 (2016).
  * Mandal  _et al._ (2017) Dibyendu Mandal, Katherine Klymko, and Michael R. DeWeese, “Entropy Production and Fluctuation Theorems for Active Matter,” Physical Review Letters 119, 258001 (2017).
  * Li and Cates (2021) Yuting I Li and Michael E Cates, “Steady state entropy production rate for scalar Langevin field theories,” Journal of Statistical Mechanics: Theory and Experiment 2021, 013211 (2021).
  * Caballero and Cates (2020) Fernando Caballero and Michael E. Cates, “Stealth Entropy Production in Active Field Theories near Ising Critical Points,” Physical Review Letters 124, 240604 (2020).
  * Fodor  _et al._ (2022) Étienne Fodor, Robert L. Jack, and Michael E. Cates, “Irreversibility and Biased Ensembles in Active Matter: Insights from Stochastic Thermodynamics,” Annual Review of Condensed Matter Physics 13, 215–238 (2022).
  * O’Byrne  _et al._ (2022) J. O’Byrne, Y. Kafri, J. Tailleur, and F. van Wijland, “Time irreversibility in active matter, from micro to macro,” Nature Reviews Physics 4, 167–183 (2022).
  * Nardini  _et al._ (2017) Cesare Nardini, Étienne Fodor, Elsen Tjhung, Frédéric van Wijland, Julien Tailleur, and Michael E. Cates, “Entropy Production in Field Theories without Time-Reversal Symmetry: Quantifying the Non-Equilibrium Character of Active Matter,” Physical Review X 7, 021007 (2017).
  * Borthne  _et al._ (2020) Øyvind L Borthne, Étienne Fodor, and Michael E Cates, “Time-reversal symmetry violations and entropy production in field theories of polar active matter,” New Journal of Physics 22, 123012 (2020).
  * Solon  _et al._ (2018a) Alexandre P Solon, Joakim Stenhammar, Michael E Cates, Yariv Kafri, and Julien Tailleur, “Generalized thermodynamics of motility-induced phase separation: Phase equilibria, Laplace pressure, and change of ensembles,” New Journal of Physics 20, 075001 (2018a).
  * Cagnetta  _et al._ (2017) F. Cagnetta, F. Corberi, G. Gonnella, and A. Suma, “Large Fluctuations and Dynamic Phase Transition in a System of Self-Propelled Particles,” Physical Review Letters 119, 158002 (2017).
  * Whitelam  _et al._ (2018) Stephen Whitelam, Katherine Klymko, and Dibyendu Mandal, “Phase separation and large deviations of lattice active matter,” The Journal of Chemical Physics 148, 154902 (2018).
  * Tociu  _et al._ (2019) Laura Tociu, Étienne Fodor, Takahiro Nemoto, and Suriyanarayanan Vaikuntanathan, “How Dissipation Constrains Fluctuations in Nonequilibrium Liquids: Diffusion, Structure, and Biased Interactions,” Physical Review X 9, 041026 (2019).
  * Gradenigo and Majumdar (2019) Giacomo Gradenigo and Satya N Majumdar, “A first-order dynamical transition in the displacement distribution of a driven run-and-tumble particle,” Journal of Statistical Mechanics: Theory and Experiment 2019, 053206 (2019).
  * Nemoto  _et al._ (2019) Takahiro Nemoto, Étienne Fodor, Michael E. Cates, Robert L. Jack, and Julien Tailleur, “Optimizing active work: Dynamical phase transitions, collective motion, and jamming,” Physical Review E 99, 022605 (2019).
  * Cagnetta and Mallmin (2020) F. Cagnetta and E. Mallmin, “Efficiency of one-dimensional active transport conditioned on motility,” Physical Review E 101, 022130 (2020).
  * Chiarantoni  _et al._ (2020) P Chiarantoni, F Cagnetta, F Corberi, G Gonnella, and A Suma, “Work fluctuations of self-propelled particles in the phase separated state,” Journal of Physics A: Mathematical and Theoretical 53, 36LT02 (2020).
  * Fodor  _et al._ (2020) Étienne Fodor, Takahiro Nemoto, and Suriyanarayanan Vaikuntanathan, “Dissipation controls transport and phase transitions in active fluids: Mobility, diffusion and biased ensembles,” New Journal of Physics 22, 013052 (2020).
  * Yan  _et al._ (2022) Jiawei Yan, Hugo Touchette, and Grant M. Rotskoff, “Learning nonequilibrium control forces to characterize dynamical phase transitions,” Physical Review E 105, 024115 (2022).
  * GrandPre  _et al._ (2021) Trevor GrandPre, Katherine Klymko, Kranthi K. Mandadapu, and David T. Limmer, “Entropy production fluctuations encode collective behavior in active matter,” Physical Review E 103, 012613 (2021).
  * Keta  _et al._ (2021) Yann-Edwin Keta, Étienne Fodor, Frédéric van Wijland, Michael E. Cates, and Robert L. Jack, “Collective motion in large deviations of active particles,” Physical Review E 103, 022603 (2021).
  * Agranov  _et al._ (2021) Tal Agranov, Sunghan Ro, Yariv Kafri, and Vivien Lecomte, “Exact fluctuating hydrodynamics of active lattice gases—typical fluctuations,” Journal of Statistical Mechanics: Theory and Experiment 2021, 083208 (2021).
  * Agranov  _et al._ (2022) Tal Agranov, Sunghan Ro, Yariv Kafri, and Vivien Lecomte, “Macroscopic Fluctuation Theory and Current Fluctuations in Active Lattice Gases,” (2022), arXiv:2208.02124 [cond-mat] .
  * Peshkov  _et al._ (2012) Anton Peshkov, Igor S. Aranson, Eric Bertin, Hugues Chaté, and Francesco Ginelli, “Nonlinear Field Equations for Aligning Self-Propelled Rods,” Physical Review Letters 109, 268701 (2012).
  * Bertin  _et al._ (2013) Eric Bertin, Hugues Chaté, Francesco Ginelli, Shradha Mishra, Anton Peshkov, and Sriram Ramaswamy, “Mesoscopic theory for fluctuating active nematics,” New Journal of Physics 15, 085032 (2013).
  * Wittkowski  _et al._ (2014) Raphael Wittkowski, Adriano Tiribocchi, Joakim Stenhammar, Rosalind J. Allen, Davide Marenduzzo, and Michael E. Cates, “Scalar $\Phi$4 field theory for active-particle phase separation,” Nature Communications 5, 4351 (2014).
  * Solon  _et al._ (2018b) Alexandre P. Solon, Joakim Stenhammar, Michael E. Cates, Yariv Kafri, and Julien Tailleur, “Generalized thermodynamics of phase equilibria in scalar active matter,” Physical Review E 97, 020602 (2018b).
  * Tjhung  _et al._ (2018) Elsen Tjhung, Cesare Nardini, and Michael E. Cates, “Cluster Phases and Bubbly Phase Separation in Active Fluids: Reversal of the Ostwald Process,” Physical Review X 8, 031080 (2018).
  * Pietzonka and Seifert (2018) Patrick Pietzonka and Udo Seifert, “Entropy production of active particles and for particles in active baths,” Journal of Physics A: Mathematical and Theoretical 51, 01LT01 (2018).
  * Seifert (2012) Udo Seifert, “Stochastic thermodynamics, fluctuation theorems and molecular machines,” Reports on Progress in Physics 75, 126001 (2012).
  * Kourbane-Houssene  _et al._ (2018) Mourtaza Kourbane-Houssene, Clément Erignoux, Thierry Bodineau, and Julien Tailleur, “Exact Hydrodynamic Description of Active Lattice Gases,” Physical Review Letters 120, 268003 (2018).
  * Bertini  _et al._ (2015) Lorenzo Bertini, Alberto De Sole, Davide Gabrielli, Giovanni Jona-Lasinio, and Claudio Landim, “Macroscopic fluctuation theory,” Reviews of Modern Physics 87, 593–636 (2015).
  * Cates  _et al._ (2022) Michael E. Cates, Étienne Fodor, Tomer Markovich, Cesare Nardini, and Elsen Tjhung, “Stochastic Hydrodynamics of Complex Fluids: Discretisation and Entropy Production,” Entropy 24, 254 (2022).
  * Jack  _et al._ (2015) Robert L. Jack, Ian R. Thompson, and Peter Sollich, “Hyperuniformity and Phase Separation in Biased Ensembles of Trajectories for Diffusive Systems,” Physical Review Letters 114, 060601 (2015).
  * Dolezal and Jack (2019) Jakub Dolezal and Robert L Jack, “Large deviations and optimal control forces for hard particles in one dimension,” Journal of Statistical Mechanics: Theory and Experiment 2019, 123208 (2019).
  * Lecomte  _et al._ (2012) Vivien Lecomte, Juan P Garrahan, and Frédéric van Wijland, “Inactive dynamical phase of a symmetric exclusion process on a ring,” Journal of Physics A: Mathematical and Theoretical 45, 175001 (2012).
  * (44) Tal Agranov, Michael E. Cates, and Robert L. Jack, in preparation .
  * Touchette (2009) Hugo Touchette, “The large deviation approach to statistical mechanics,” Physics Reports 478, 1–69 (2009).
  * Jack (2020) Robert L. Jack, “Ergodicity and large deviations in physical systems with stochastic dynamics,” The European Physical Journal B 93, 74 (2020).
  * Pietzonka  _et al._ (2016) Patrick Pietzonka, Andre C. Barato, and Udo Seifert, “Universal bounds on current fluctuations,” Physical Review E 93, 052145 (2016).
  * Horowitz and Gingrich (2017) Jordan M. Horowitz and Todd R. Gingrich, “Proof of the finite-time thermodynamic uncertainty relation for steady-state currents,” Physical Review E 96, 020103 (2017).
  * Bodineau and Lagouge (2010) T. Bodineau and M. Lagouge, “Current Large Deviations in a Driven Dissipative Model,” Journal of Statistical Physics 139, 201–218 (2010).
  * Baek  _et al._ (2017) Yongjoo Baek, Yariv Kafri, and Vivien Lecomte, “Dynamical Symmetry Breaking and Phase Transitions in Driven Diffusive Systems,” Physical Review Letters 118, 030604 (2017).
  * Baek  _et al._ (2018) Yongjoo Baek, Yariv Kafri, and Vivien Lecomte, “Dynamical phase transitions in the current distribution of driven diffusive channels,” Journal of Physics A: Mathematical and Theoretical 51, 105001 (2018).
  * Bodineau and Derrida (2004) T. Bodineau and B. Derrida, “Current Fluctuations in Nonequilibrium Diffusive Systems: An Additivity Principle,” Physical Review Letters 92, 180601 (2004).
  * Hurtado and Garrido (2009) Pablo I. Hurtado and Pedro L. Garrido, “Test of the Additivity Principle for Current Fluctuations in a Model of Heat Conduction,” Physical Review Letters 102, 250601 (2009).
  * Hurtado and Garrido (2010) Pablo I. Hurtado and Pedro L. Garrido, “Large fluctuations of the macroscopic current in diffusive systems: A numerical test of the additivity principle,” Physical Review E 81, 041102 (2010).
  * Hurtado  _et al._ (2011) P. I. Hurtado, C. Perez-Espigares, J. J. del Pozo, and P. L. Garrido, “Symmetries in fluctuations far from equilibrium,” Proceedings of the National Academy of Sciences 108, 7704–7709 (2011).
  * Prados  _et al._ (2011) A. Prados, A. Lasanta, and Pablo I. Hurtado, “Large Fluctuations in Driven Dissipative Media,” Physical Review Letters 107, 140601 (2011).
  * Meerson  _et al._ (2014) Baruch Meerson, Arkady Vilenkin, and P. L. Krapivsky, “Survival of a static target in a gas of diffusing particles with exclusion,” Physical Review E 90, 022120 (2014).
  * Meerson (2015) Baruch Meerson, “Full absorption statistics of diffusing particles with exclusion,” Journal of Statistical Mechanics: Theory and Experiment 2015, P04009 (2015).
  * Agranov  _et al._ (2016) Tal Agranov, Baruch Meerson, and Arkady Vilenkin, “Survival of interacting diffusing particles inside a domain with absorbing boundary,” Physical Review E 93, 012136 (2016).
  * Agranov and Meerson (2017) Tal Agranov and Baruch Meerson, “Fluctuations of absorption of interacting diffusing particles by multiple absorbers,” Physical Review E 95, 062124 (2017).
  * Agranov and Meerson (2018) Tal Agranov and Baruch Meerson, “Narrow Escape of Interacting Diffusing Particles,” Physical Review Letters 120, 120601 (2018).
  * Agranov  _et al._ (2019) Tal Agranov, P. L. Krapivsky, and Baruch Meerson, “Occupation-time statistics of a gas of interacting diffusing particles,” Physical Review E 99, 052102 (2019).
  * Bertini  _et al._ (2005) L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, “Current Fluctuations in Stochastic Lattice Gases,” Physical Review Letters 94, 030601 (2005).
  * Bodineau and Derrida (2005) T. Bodineau and B. Derrida, “Distribution of current in nonequilibrium diffusive systems and phase transitions,” Physical Review E 72, 066110 (2005).
  * Bertini  _et al._ (2006) L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, “Non Equilibrium Current Fluctuations in Stochastic Lattice Gases,” Journal of Statistical Physics 123, 237–276 (2006).
  * Shpielberg and Akkermans (2016) O. Shpielberg and E. Akkermans, “Le Chatelier Principle for Out-of-Equilibrium and Boundary-Driven Systems: Application to Dynamical Phase Transitions,” Physical Review Letters 116, 240603 (2016).
  * Zarfaty and Meerson (2016) Lior Zarfaty and Baruch Meerson, “Statistics of large currents in the Kipnis–Marchioro–Presutti model in a ring geometry,” Journal of Statistical Mechanics: Theory and Experiment 2016, 033304 (2016).
  * Andelman  _et al._ (1987) David Andelman, M. E. Cates, D. Roux, and S. A. Safran, “Structure and phase equilibria of microemulsions,” The Journal of Chemical Physics 87, 7229–7241 (1987).
  * Bodineau and Derrida (2007) Thierry Bodineau and Bernard Derrida, “Cumulants and large deviations of the current through non-equilibrium steady states,” Comtes Rendus Physique 8, 540–555 (2007).
  * Kaiser  _et al._ (2018) Marcus Kaiser, Robert L. Jack, and Johannes Zimmer, “Canonical Structure and Orthogonality of Forces and Currents in Irreversible Markov Chains,” Journal of Statistical Physics 170, 1019–1050 (2018).


# Clouds on partial atmospheres of lava planets and where to find them

T. Giang Nguyen Department of Physics, McGill University, 3600 rue University,
Montréal, QC H3A 2T8, Canada  Nicolas B. Cowan Department of Physics, McGill
University, 3600 rue University, Montréal, QC H3A 2T8, Canada  Department of
Earth & Planetary Sciences, 3450 rue University, Montréal, QC H3A 0E8, Canada
Lisa Dang Trottier Institute for Research on Exoplanets and Department of
Physics, Université de Montréal, 1375 Ave.Thérèse-Lavoie-Roux Montréal, QC H2V
0B3, Canada

###### Abstract

With dayside temperatures hot enough to sustain a magma ocean and a silicate
atmosphere, lava planets are the best targets to study the atmosphere of a
rocky world. In the absence of nightside heating, the entire atmosphere
collapses near the day-night terminator, so condensation seems inevitable, but
the impact of clouds on radiative transfer, dynamics, and observables has not
yet been studied in the non-global atmospheric regime. Therefore, we simulate
cloud formation and determine which lava planets should be most affected by
clouds. We find that despite the scattering of visible light by clouds, heat
advection compensates for the cooling effect of clouds in the atmosphere. On
the other hand, surface temperatures are significantly affected and can drop
100-200 K under a cloudy sky. We find that among our targets, HD213885b and
HD20329b are most affected by cloud formation: there is a discernable
difference between having clouds and not having them, but the precision
required to make such an inference is at the limit of current instruments.

exoplanet, terrestrial, atmosphere, modelling

##  1 Introduction

###  1.1 Background

Hundreds of Earth-sized planets have been discovered (Adams et al., 2021) but
inferring their surface and atmospheric properties remains a significant
challenge. Our best options for in-depth characterization of rocky planets are
currently ultra-short period planets. Phase curves obtained with the Spitzer
Space Telescope have enabled tentative inferences of atmospheres on highly
irradiated super-Earths K2-141b (Zieba et al., 2022) and 55 Cnc-e (Demory et
al., 2016b, a; Mercier et al., 2022).

A more powerful instrument is needed for definitive atmospheric detection and
characterization of thin silicate atmospheres. Indeed, the recently launched
James Webb Space Telescope (JWST) has observed ultra-short period planets
K2-141b (Dang et al., 2021; Espinoza et al., 2021) and 55 Cnc-e (Hu et al.,
2024), and will observe TOI-561b (Teske et al., 2023) and K2-22b (Wright et
al., 2023). Such observations can provide unique clues about terrestrial
planet evolution.

Ultra-short period rocky planets with dayside temperatures $>$1800K are
theorized to harbour a permanent dayside magma oceans and are often referred
to a lava planets. Since the silicate melt is directly exposed to the
atmosphere, atmospheric dynamics are intrinsically dependent on interior
dynamics (Dorn & Lichtenberg, 2021; Kite et al., 2016; Boukaré et al., 2022).
Magma oceans are a common result of planetary formation (Schaefer & Elkins-
Tanton, 2018) and play a critical role in sequestering volatiles (Moore et
al., 2023), so lava planets provide a unique window into the early solar
system and Earth (Hirschmann, 2012; Zahnle et al., 2020). For a complete
overview of silicate outgassing and the corresponding chemistry and radiative-
convective dynamics, see Wordsworth & Kreidberg (2022).

Tidally locked planets risk having their atmosphere condense out entirely on
the permanent nightside (Wordsworth, 2015). However, tidally-locked
terrestrial planets may sustain a global atmosphere if the initial outgassed
atmosphere is thick enough (Joshi, 2003). Therefore, general circulation
models (GCMs) are adequate for simulating a thick-atmosphere regime for lava
planets such as 55 Cnc-e Hammond & Pierrehumbert (2017). A general source of
atmospheric circulation on tidally locked planets can be found in
Pierrehumbert & Hammond (2019). Despite such robust work on lava planet
atmospheres, the non-global atmospheric regime is less studied and cannot be
simulated with a GCM. Two approaches are therefore used to simulate a thin
lava planet atmosphere: either focusing on the dayside radiative transfer or
the hydrodynamics.

Radiative transfer simulations of lava planets show a silicate atmosphere with
a strong vertical temperature inversion (e.g., Ito et al., 2015; Zilinskas et
al., 2022; Piette et al., 2023). These studies often find the atmosphere
hotter than the surface as heat cannot be advected to cooler regions. On the
other hand, hydrodynamical models of non-global atmospheres often simplify, if
not outright ignore radiative transfer (e.g., Castan & Menou, 2011; Nguyen et
al., 2020). However, this method can trace the horizontal heat and mass
transport due to day-night winds. Both approaches ultimately lack self-
consistent cloud formation which affects both the fluid dynamics and radiative
transfer.

###  1.2 Objective

In anticipation of JWST observations of lava planets, we simulate cloud
formation and quantify their effect on thin-atmospheric dynamics and
observations. Although there are many published works specific to exoplanet
clouds (for a review see Helling, 2019), it is impossible to use a
sophisticated cloud model within an existing non-global atmosphere modelling
framework. We therefore use a simple parameterization of cloud formation
consistent with the turbulent boundary layer model of Nguyen et al. (2022)
which is the only published work that have coupled radiative transfer with
hydrodynamical processes for non-global atmospheres.

We use this model to simulate the atmospheres of five lava planets: K2-141b
(Malavolta et al., 2018), TOI-431b (Osborn et al., 2021), HD20329b (Murgas et
al., 2022), HD213885b (Espinoza et al., 2020), and 55-cnc e (Winn et al.,
2011, but we note that its recently-published eclipse spectrum suggests a
volatile-full planet: Hu et al., 2024). The selected candidates have some of
the best signal-to-noise ratios for rocky exoplanets and they will likely be
observed extensively. Therefore, we compare these lava planets to predict
where cloud formation might be expected.

##  2 Methods

###  2.1 The hydrodynamical model

Thin and tenuous atmospheres are common on moons and dwarf planets of the
outer solar system. Huge temperature fluctuations on these bodies induce
periods and locations of atmospheric collapse and massive surface outgassing,
e.g., on Io (Tsang et al., 2016) or Pluto (Gladstone & Young, 2019). Due to
the stark changes in fluid regime, it is difficult to implement GCMs and
hydrodynamics are best approximated using boundary layer analysis (Ingersoll
et al., 1985). Because the atmospheres of lava planets are expected to behave
similarly, the same calculations are useful, if not necessary.

The expectation of synchronous rotation due to tidal locking makes lava
planets easier to deal with than solar system objects. Having a permanent
dayside and nightside allows for steady-state flow. Due to the complexities in
handling non-linear terms around evaporation and condensation, we model a pure
atmosphere. We are also forced to neglect atmospheric escape which should
happen for lighter volatiles (Ito & Ikoma, 2021) but is justified if we track
heavy molecules like SiO. Finally, we neglect Coriolis force as to preserve
the axial symmetry of tidally locked objects. Generally, for planets with
rotational period of multiple days, the km/s wind speeds are strong enough to
make Coriolis force insignificant (Castan & Menou, 2011). However, our targets
have a period of $\leq$ 1 day and so we acknowledge that Coriolis force may
play a large role in dynamics but implementing it is beyond our current
capability.

By applying the approximations described above, we can construct a framework
that describes the flow along the planetary boundary layer. The formulation of
the base evaporation-driven flow is akin to the shallow-water equations of
Ingersoll et al. (1985). The equations essentially calculate pressure, winds,
and temperature of the atmosphere at half the scale height. Detailed
formulation can be found in the Appendix.

###  2.2 Radiative transfer

Our implementation of radiative transfer largely follows Nguyen et al. (2022)
and a detailed formulation can be found in the appendix. Nguyen et al. (2022)
split up the radiative transfer into three bands (UV, visible, IR) with an
averaged absorption cross-section within each band. To improve on this scheme,
we now calculate radiative absorption and emission line by line in the
spectrum, doing away with averaging spectral bands. We validate our radiative
transfer scheme by directly comparing our wind-free silicate atmosphere with
Zilinskas et al. (2022). We find that our atmospheric temperatures match well
with theirs (within 5%) whereas the old scheme of Nguyen et al. (2022) can be
off by up to 25%, especially for hotter and brighter stars.

We adopt the absorption cross-section of SiO from in Yurchenko et al. (2022).
As for the stellar spectra, we found proxy analogues in the MUSCLES treasury
surveys for K2-141, TOI-431 and 55-Cnc: their closest analogues are HD85512,
HD40307, and HD97658, respectively (France et al., 2016). For our targets of
HD213885 and HD20329, their properties are closest to the Sun and so we use
the solar spectrum; the spectra are shown on the bottom panel of Fig. 1.
Absolute magnitude of each analogue’s spectrum is scaled to the target’s
effective stellar temperature and orbital distance.

Figure 1: _Top:_ the absorption cross-section of SiO (Yurchenko et al., 2022).
The blue, green, and red show the UV, visible, and IR range. _Bottom:_ the
normalized flux of each star using close analogues.

###  2.3 Cloud implementation

Clouds are abundant on exoplanets but their modelling and detection have
largely been confined to hot Jupiters (Demory et al., 2013) and sub-Neptunes
(Kreidberg et al., 2014). Some of these clouds likely consist of silicate
material like SiO and MgSiO3 (Powell et al., 2018; Gao et al., 2021). Although
lava planets have exotic atmospheric dynamics, their cloud composition is not,
at least by the standards of exoplanets. Cloud particles on lava planets may
be analogous to dust found on catastrophically evaporating rocky worlds Curry
et al. (2024).

Rapid mass exchange between the surface and the atmosphere of a lava planet
may lead to meteorological processes driven by the recycling of rocky material
between the surface and atmosphere (Boukaré et al., 2022). Volatile components
like Na and K may form clouds due to interactions with stellar winds but these
light material are subjected to atmospheric loss (Ito & Ikoma, 2021) or freeze
out in nightside alkali glaciers (Nguyen et al., 2020).

Because we are focusing on lava planets with Earth-like bulk composition
(excluding 55 Cnc-e), we can reasonably assume that volatiles such as water or
CO2 are absent from the atmosphere (cf. Graham et al., 2021). Silicate
condensates remain the dominant constituent but if the dew point is too close
to the surface, the silicate vapour may not condense midair but simply deposit
directly on the surface of the magma ocean (Schaefer & Fegley, 2009). We do
not know how much of the SiO vapour will condense into cloud particles nor do
we know the specific properties of these particles. Since clouds might have
significant impact on the atmosphere of a lava planet, ignoring them
altogether is not a viable option. We therefore opt to parameterize cloud
effects into our hydrodynamical model in the same spirit as how radiative
transfer was implemented, i.e., as simply as possible.

As cloud formation is inherently tied to the condensation rate $C$, we define
a cloud optical depth proportional to $C$: $\tau_{c}=\text{f}_{c}\ C$. The
dimensionless fc parameterizes the fraction of condensible material that
actually turns into cloud particles as well as the particle properties that
affects optical depth such as particle size and density. Although much nuance
is lost in this formulation, current observations only allow us to probe the
temperature of the atmosphere and surface. Therefore, our main aim is to
estimate how much clouds affect the overall thermodynamics and see if clouds
can produce a detectable signal.

By creating and testing a range of fc, we are essentially inferring how the
planet behaves under varying strength of cloud coverage. This helps to
identify important aspects such as minimum detectable cloud formation, which
planets are most affected by clouds, and direct inference of cloud properties.

We calculate the resultant opacity ($\epsilon_{\text{c}}$) following Eq. 22.
Here we make another simplifying cloud assumption: cloud particles scatter
visible light isotropically so that half the incident visible stellar light
($F_{*,\text{vis}}$) is reflected upwards when the cloud opacity is 1. This
assumes that the clouds are transparent in the IR and that most of the UV is
absorbed by the upper levels of the atmosphere. Note that ignoring cloud
absorption in the IR neglects their greenhouse effect, removing significant
complexities to the system of equations. This yields $F_{\text{c}}$, the
stellar flux that is reflected because of clouds:

$F_{\text{c}}=\frac{1}{2}\ \epsilon_{\text{c}}\ F_{*,\text{vis}}.$ (1)

Because SiO is nearly transparent at the visible wavelengths (Fig. 1), the
atmospheric energy balance is unaffected by $F_{c}$. However, the surface
energy balance is significantly affected as visible light is absorbed entirely
by the surface. This stems from treating the surface as a blackbody but is
justified in the expected low geometric albedo of magma surfaces (Essack et
al., 2020): the surface energy is calculated from the bulk bolometric flux
minus stellar radiation absorbed by the atmosphere plus the emissions from the
atmosphere minus visible stellar flux that is reflected by clouds minus latent
heat when evaporating. The surface energy balance then becomes:

$\sigma T_{s}=F_{*}-F_{\text{stel}}+F_{\text{RC}}-F_{c}-Q_{\text{lat}},$ (2)

where $\sigma$ is the Stefan-Boltzmann constant. The surface temperature,
$T_{s}$ depends on the overlying atmospheric properties and vice versa (i.e.
Eq. 9, Eq. 14), coupling the surface and atmosphere. $F_{*}$ is the bulk
stellar flux, $F_{stel}$ is the stellar flux that has been absorbed by the
atmosphere, $F_{\text{RC}}$ is the thermal emission from the atmosphere,
$F_{c}$ is the stellar flux that was reflected by clouds, and $Q_{\text{lat}}$
is the latent heat of vaporization.

###  2.4 Simulated observations

This subsection will describe how state variables $P,V$ and $T$, computed as
in the Appendix, are used to simulate observations. Because a lava planet’s
atmosphere is so thin and does not extend much past the terminator,
transmission spectroscopy of the gravitationaly bound atmosphere is
challenging (Nguyen et al., 2022). We therefore focus on planetary emission.
We first calculate the top-of-atmosphere outgoing radiative flux
$F_{\text{TOA}}$:

$F_{\text{TOA}}(\lambda)=\epsilon_{\lambda}B(\lambda,T)+(1-\epsilon_{\lambda})B(\lambda,T_{s}),$
(3)

where the first term is atmospheric emission and the second surface emission.
We take into account the adiabatic profile where $T(P)$ is calculated at the
pressure where optical depth is 2/3. For any location $\theta$, the flux can
be calculated and so we convert temperature maps into flux maps.

We convert $F_{\text{TOA}}(\theta)$ into a map of co-latitude ($\alpha$) and
longitude ($\phi$) assuming polar symmetry about the sub-stellar point. We
then integrate the flux across the visible hemisphere following Cowan & Agol
(2008):

$F_{e}(\lambda,\zeta)=\int_{0}^{\pi}\int_{\zeta-\pi/2}^{\zeta+\pi/2}F(\alpha,\phi,\lambda)\sin^{2}(\alpha)\cos(\phi-\zeta)d\phi
d\alpha,$ (4)

where $F_{e}$ is the emission signal, $\zeta$ is the position of the planet in
its orbit ($\zeta=0$ at superior conjunction and $\zeta=\pi$ at inferior
conjunction). Keeping $\lambda$ constant and varying $\zeta$ between 0 and
$2\pi$ yields a phase curve at wavelength $\lambda$. Keeping $\zeta$ at 0
while varying $\lambda$ yields the emission spectrum at eclipse.

##  3 Results

###  3.1 Output from model

All the equations from the previous section together form a self-consistent
model that couples the atmosphere and surface along with rudimentary cloud
formation. To simulate our specific targets, we apply the planetary parameters
listed in Table 1.

Table 1: System parameters used for the hydrodynamical model. Although these values do not show up in the equations mentioned in the previous section with the exception of $g$, the orbital properties are needed to calculate the incident stellar flux pattern, $F_{*}(\theta)$, following Kopal (1954). Planet Name | Radius (RE) | Stellar Radius (R∗) | Distance (a/R∗) | Stellar Temp. (K) | g (m/s2) | Source
---|---|---|---|---|---|---
K2-141b | 1.54 | 0.681 | 2.30 | 4373 | 21.8 | Barragán et al. (2018)
TOI-431b | 1.28 | 0.731 | 3.32 | 4850 | 18.4 | Osborn et al. (2021)
HD20329b | 1.72 | 1.13 | 3.42 | 5596 | 24.6 | Murgas et al. (2022)
HD213885b | 1.74 | 1.10 | 3.93 | 5978 | 28.5 | Espinoza et al. (2020)
55 cnc e | 1.87 | 0.964 | 3.52 | 5317 | 22.3 | Bourrier et al. (2018)

We show the initial results without any cloud formation in Fig. 2. We see that
atmospheric pressure is directly related to surface temperature: the hotter
the surface, the thicker the atmosphere. However, despite having cooler
surface temperature, HD213885b and HD20329b have much hotter atmospheres; they
also have stronger winds.

Figure 2: State variables half a scale height above the surface for selected
lava planets. From top to bottom, the panels show pressure, wind speed,
atmospheric temperature, and surface temperature. The solutions end when the
flow no longer behaves as a continuous fluid and the calculations break down;
this is essentially where horizontal mass advection reach near 0 and the
atmospheric collapse is complete. K2-141b has the thickest atmosphere as it
receives the most stellar flux and therefore has the hottest surface
temperature. However, K2-141b’s atmosphere is significantly cooler than
HD213885b and HD20329b as atmospheric temperature is largely dependent on
stellar UV flux. The temperature jump for 55 Cnc-e at the edge of the
atmosphere can be attributed to a dynamical effect described by Ingersoll et
al. (1985) where the air viscosity begin to induce significant heating at
certain pressure and temperature levels; the same phenomenon can be seen in
Castan & Menou (2011).

We plot the contributions of all radiative terms in Fig. 3. With the exception
of K2-141b, most of the atmospheric absorption is in the UV owing the
absorption features of SiO in that range. Most of the radiative cooling occurs
in the IR although HD213885b and HD20329b are hot enough to emit significant
UV radiation.

Despite the enormous magnitude of the radiative fluxes, incoming and outgoing
radiation tend to reach equilibrium: the net radiative balance is 1–2 order of
magnitude smaller than each individual radiative terms.

Figure 3: Radiative budget of atmosphere of our five lava planets. The first
column shows absorption of stellar flux, the second shows absorption of
blackbody radiation from the surface, and the third shows radiative cooling.
K2-141 is less UV bright than the other stars, so the stellar absorption of
its planet is mostly in the IR. Because of the high UV stellar absorption on
HD213885b and HD20329b, their atmospheres become hot enough to emit
significant UV radiation.

###  3.2 Results with clouds

To review, we parameterized cloud optical effects using the parameter fc in
the Methods section. Note that our simple implementation of clouds neglect
certain feedbacks such as the greenhouse effect of clouds. With the assumption
that cloud formation must be proportional to the condensation rate $C$,
optical depth and opacity must subsequently also depend on $C$. Therefore, fc
acts as a substitute to estimate how much of the condensible material becomes
cloud, as well as cloud properties such as droplet size; we explore these
properties in the Discussion section.

We restrict fc between 102–104; any lower and cloud effects are negligible
across all planets while anything higher will make opacity uniformly 1.
Because the SiO vapour absorption is weak in the visible, clouds do not
directly affect the radiative budget of the atmosphere but do so via
atmosphere–surface interactions. The biggest impact, therefore, is on the
surface energy budget as expected from our formulation (see Eq. 2). We show
the condensation rate, optical depth, opacity, and surface temperature of
cloudy simulations in Fig. 4.

Figure 4:  The impact of clouds on lava planet surface temperature. The first
column is the condensation rate. We zoom in on $\theta$ where condensation
occurs. We impose the factors fc onto the condensation rate to create the
optical depth plots shown in the second column, and the opacity is plotted in
the third column. The final column is the surface temperature below the cloud
deck. Clouds cool the surface significantly (up to 200 K) but the changes in
atmospheric temperature are an order of magnitude smaller.

###  3.3 Simulated observations

Given a global map for surface and atmospheric temperature, we calculate the
outgoing radiative flux at the top of the atmosphere using Eq. 3. We fix the
time at eclipse and integrate the flux over the visible hemisphere using Eq. 4
at different wavelength $\lambda$ to yield Fig. 5.

Figure 5: Simulated eclipse spectra of selected planets. The spectrum is
calculated from the top-of-the-atmosphere emission which ultimately depends on
the temperatures of both the atmosphere and the surface. In strong SiO
spectral features, emission originates from upper layers of the atmosphere.
For an adiabatic T-P profile, this leads to extremely cold temperature reading
and thus the eclipse spectra show deep absorption features. Although clouds do
cool the atmosphere slightly, the changes in the eclipse depth are barely
visible. Since clouds cool the surface more, the overall blackbody curve of
the spectra is lowered and the difference between cloudy and cloudless is seen
relatively most clearly on HD20329b and HD213885.

##  4 Discussion

###  4.1 On radiative energy balance

We see that the stellar spectrum predictably dictates how hot the atmosphere
gets. Despite the surface of K2-141b receiving more bolometric stellar flux
than HD20329b and HD213885b, its atmosphere absorbs less stellar flux as the
star K2-141 emits mostly in the visible and IR. Meanwhile, HD213885b and
HD20329b absorbs more stellar UV and visible radiation, thus providing more
heat to their atmopsheres.

The coarser radiative scheme of Nguyen et al. (2022) tended to overestimate UV
radiative balance while underestimating the IR. With a finer spectral
resolution, we see that UV emission never dominates and that SiO spectral
features in the IR are strong enough to emit flux even when the pressure is
low.

###  4.2 Vertical temperature profile

One-dimensional radiative-convective models predict that while the lower
silicate atmosphere of lava planets does follow an adiabatic profile, there is
a sharp temperature inversion at the top of the atmosphere due to shallow UV
penetration in the atmosphere Zilinskas et al. (2022); Piette et al. (2023).
Hydrodynamical studies, on the other hand, predict that advection of heat
weakens this inversion Nguyen et al. (2022).

Although we adopted an adiabatic vertical temperature profile, we can guess
what happens when the vertical temperature lapse rate is inverted. This can be
done by increasing the $\beta$ value (this parameter is calculated from
vertical pressure integration in Eq. 17). Nguyen et al. (2022) showed that
arbitrarily increasing $\beta$ leads to stronger horizontal flow but an
overall cooler atmosphere. Within the calculations, raising $\beta$ increases
the horizontal pressure gradient force, thus increasing wind velocity, $V$.
However, the energy balance does not change so any gains in the kinetic term,
$V^{2}/2$, translate to a loss in the internal heat term, $CpT$ (see Eq. 8).
The physical reasoning is that because the atmospheric inversion inhibits
mixing, horizontal winds are significantly increased, but since energy must be
conserved, the air’s internal heat energy must be lower.

When balancing the incoming radiation in the atmosphere, we find that stellar
absorption is comparable to surface blackbody absorption for K2-141b,
TOI-431b, and 55 Cnce; this suggests a rather isothermal T-P profile as there
is as much heat from the top as there is at the bottom. However, for HD213885b
and HD20329b, there is significantly more stellar absorption than surface
absorption. With significantly more heating from the top, most of which is UV
flux, we would expect a strong temperature inversion as predicted by Zilinskas
et al. (2022). While we have kept a consistent vertical temperature profile in
this study to isolate the effects of clouds, we note that future simulations
of lava planets should aspire to use self-consistent vertical profiles given
the spectrum of the host star.

###  4.3 Latent and chemical heat

As the radiative terms tend to cancel out, the net energy budget of the
atmosphere is controlled by dynamics terms such as latent heat and sensible
heat. We show in Fig. 6 that the latent heat is negligible, as predicted by
Kite et al. (2016). Even when coupling the latent heat with the atmosphere
where radiative transfer is hampered by the absorption cross-section, latent
heating changes the atmosphere by less than 5%. However, K2-141b has a
relatively larger condensation rate and a relatively lower radiative
absorptivity due to its cooler star; temperatures can increase up to 10% when
latent heating is accounted for. Sensible heating is still the main driver of
the flow and it is composed of heat flux from two directions: horizontally via
flow-induced heat advection and vertically via temperature mixing between the
surface and atmosphere.

Figure 6:  We show the changes in atmospheric and surface temperature when
latent heat is implemented self-consistently in our code. Latent heat cools
the surface when evaporating and warms the atmosphere when condensing, because
we assume that condensation occurs aloft. As predicted by Kite et al. (2016),
temperature changes due to latent heating are negligible at the surface but
warms the atmosphere a little. However, the changes are minor with the
exception of K2-141b where atmospheric temperatures can warm by 10%. This is
due to the relatively large condensation rate but relatively low radiative
heating and cooling caused by orbiting around a cooler star.

Although latent heat is relatively small, the vaporization of the surface and
subsequent condensation is not a simple phase change. The magma ocean is
principally composed of SiO2 (assuming a Bulk Silicate Earth composition;
Schaefer & Fegley, 2009), while the atmosphere contains mostly SiO, hence
evaporation and condensation involve the dissociation and recombination of
oxygen atoms. In the context of hot gaseous planets, Bell & Cowan (2018)
showed that the energy of dissociation and recombination of hydrogen molecules
can be significant. We estimate the energy of the reaction SiO2 $\Rightarrow$
SiO + (1/2)O2 based on the difference in enthalpy of formation between SiO and
SiO2 (enthalpy of O2 is 0). With this difference being 811 kJ/mol (Chase &,
US), the energy per unit mass of SiO is 2$\times 10^{7}$ J/kg, an order of
magnitude magnitude greater than the latent heat of vaporization ($10^{6}$
J/kg, Kite et al., 2016). Unlike the latent heat, this chemical energy may
dominate the energy budget, but we leave this problem to future studies, since
our current model can only accommodate a single gas species.

###  4.4 The effects of clouds on lava planet climate

The formation of clouds has a weaker effect on the atmosphere than on the
surface, due to the fact that the atmosphere is optically thin at visible
wavelengths even in the absence of clouds. Visible light reflected by clouds
would have mostly passed through the atmosphere in any case; it is the surface
that loses that energy and thus can cool by 200 K.

Smaller fc correspond to larger cloud particles, which are radiatively less
important. Assuming spherical droplets, we can decompose fc into cloud
particle properties via:

$\text{f}_{c}=\frac{3}{8}\ \frac{t}{\rho\ d}\ \text{p}_{e}.$ (5)

The density, $\rho$, of the condensed SiO2 is 2650 kg/m3, pe is the fraction
of condensible material that forms cloud droplets instead of condensing
directly on the surface of the magma ocean—this depends on the altitude at
which clouds form. The particle diameter, $d$, is typically sub-micron to 10s
of microns for silicate clouds in hot Jupiter atmospheres (Powell et al.,
2018) while dust grain analysis of catastrophically evaporating rocky world
yield 0.1–1 micron Budaj (2013). Finally, $t$ is the timescale over which the
cloud droplet stays in the atmosphere. Two timescales could be at play: the
time it takes for the cloud droplet to fall to the surface or the time it
takes for the cloud droplet to be advected from the point of nucleation to
where the atmosphere collapses.

We estimate the terminal velocity of free-falling cloud particles falling
through 2500 K gas under an acceleration $g$ = 25 m/s2 using Stokes’ law:
terminal velocity, $V=(2/9)(\rho_{\text{c}}-\rho_{\text{a}})gr^{2}/\eta$,
where $\rho_{\text{c}}$ is the particle density and $\rho_{\text{a}}$ is the
air density. With micron-size particles, the terminal velocity is a very slow
2$\times 10^{-4}$ m/s. We surmise that the horizontal advection timescale is
likely the limiting factor. It takes $\sim$ 5000 s to transport a parcel of
gas 45∘ of angular distance for a planetary radius of 1.5 R${}_{\text{E}}$ and
a wind speed of 1.5 km/s. Therefore we get f${}_{c}\sim 10^{4}$ given
$d=10^{-6}$ m, and p${}_{e}=0.1$.

The effect of clouds on the atmospheric energy budget is minimal because we
have assumed clouds are transparent in the IR, but the surface is influenced
significantly. Because we have treated the magma ocean as a still surface, it
does not have any mechanism to advect heat (cf. Boukaré et al., 2022).
Therefore, the surface is entirely dependent on instellation and losing up to
half of all visible light due to clouds results in significant temperature
decreases that might be detectable due via astronomical observations.

In the simulated eclipse spectra shown in Fig. 4. We see that the greater the
cloud opacity, the lower the overall emission. Clouds have the biggest impact
on the dynamics of HD213885b and HD20329b due to the high incident stellar
radiation at visible wavelengths. However, this difference in overall signal
is too small for current instruments to precisely characterize. Ultimately,
K2-141b remains the best target lava planet for atmospheric detection because
of signal-to-noise but is the worst target to differentiate between cloudy and
clear sky conditions.

##  5 Conclusions

Clouds are a fundamental part of a planet’s meteorology. Lava planets are most
likely to have, at minimum, non-global atmospheres. In this regime, clouds
should appear as the atmosphere starts to condense out. Therefore, we
simulated the atmosphere of lava planets K2-141b, TOI-431b, HD20329b,
HD213885b, and 55-cnc e to see which of these planets are best for cloud
detection and characterization.

With or without clouds, our calculations show that atmospheric temperature is
more dependent on the stellar spectrum than bolometric flux. HD213885b and
HD20329b, whose stars emit more UV light, have a hotter atmosphere than
K2-141b despite receiving less instellation.

We implement clouds by scaling the optical depth to the condensation rate
using the fc parameter to roughly approximate weak, medium, to strong cloud
effects. We find that atmospheric properties are not significantly affected by
visible light scattering of clouds: the existing heat advection from the
evaporating region through sensible heat as well as latent is enough to offset
clouds. The surface is more affected by clouds and the effect is predictably
stronger for greater optical instellation. As such, the overall effect of
clouds is best seen on HD213885b and HD20329b but is still challenging to
distinguish.

Implementing an adiabatic profile up to the top of atmosphere create deep
absorption features which render the spectroscopic signature of clouds to be
indistinguishable; adopting the silicate T-P profiles of Zilinskas et al.
(2022) and Piette et al. (2023) should yield more accurate observables which
may reveal an opportunity to detect clouds but the results of such 1D models
must be reconciled with atmospheric dynamics.

## Acknowledgements

This research was partially funded by the Canadian Space Agency’s James Webb
Space Telescope observer’s program and the Heising-Simons Foundation. NBC
acknowledges support from an NSERC Discovery Grant, a Tier 2 Canada Research
Chair, and an Arthur B. McDonald Fellowship. The authors also thank the
Trottier Space Institute and l’Institut de recherche sur les exoplanètes for
their financial support and dynamic intellectual environment.

## References

  * Adams et al. (2021) Adams, E. R., Jackson, B., Johnson, S., et al. 2021, The Planetary Science Journal, 2, 152
  * Barragán et al. (2018) Barragán, O., Gandolfi, D., Dai, F., et al. 2018, Astronomy & Astrophysics, 612, A95
  * Bell & Cowan (2018) Bell, T. J., & Cowan, N. B. 2018, The Astrophysical Journal Letters, 857, L20
  * Boukaré et al. (2022) Boukaré, C.-É., Cowan, N. B., & Badro, J. 2022, The Astrophysical Journal, 936, 148
  * Bourrier et al. (2018) Bourrier, V., Dumusque, X., Dorn, C., et al. 2018, Astronomy & Astrophysics, 619, A1
  * Budaj (2013) Budaj, J. 2013, Astronomy & Astrophysics, 557, A72
  * Castan & Menou (2011) Castan, T., & Menou, K. 2011, The Astrophysical Journal Letters, 743, L36
  * Chase & (US) Chase, M. W., & (US), N. I. S. O. 1998, NIST-JANAF thermochemical tables, Vol. 9 (American Chemical Society Washington, DC)
  * Cowan & Agol (2008) Cowan, N. B., & Agol, E. 2008, ApJ, 678, L129, doi: 10.1086/588553
  * Curry et al. (2024) Curry, A., Booth, R., Owen, J. E., & Mohanty, S. 2024, Monthly Notices of the Royal Astronomical Society, 528, 4314
  * Dang et al. (2021) Dang, L., Cowan, N. B., Hammond, M., et al. 2021, JWST Proposal. Cycle 1, 2347
  * Demory et al. (2016a) Demory, B.-O., Gillon, M., Madhusudhan, N., & Queloz, D. 2016a, Monthly Notices of the Royal Astronomical Society, 455, 2018
  * Demory et al. (2013) Demory, B.-O., De Wit, J., Lewis, N., et al. 2013, The Astrophysical Journal Letters, 776, L25
  * Demory et al. (2016b) Demory, B.-O., Gillon, M., De Wit, J., et al. 2016b, Nature, 532, 207
  * Dorn & Lichtenberg (2021) Dorn, C., & Lichtenberg, T. 2021, The Astrophysical Journal Letters, 922, L4
  * Espinoza et al. (2020) Espinoza, N., Brahm, R., Henning, T., et al. 2020, Monthly Notices of the Royal Astronomical Society, 491, 2982
  * Espinoza et al. (2021) Espinoza, N., Bello-Arufe, A., Buchhave, L. A., et al. 2021, JWST Proposal. Cycle 1, 2159
  * Essack et al. (2020) Essack, Z., Seager, S., & Pajusalu, M. 2020, The Astrophysical Journal, 898, 160
  * France et al. (2016) France, K., Loyd, R. P., Youngblood, A., et al. 2016, The Astrophysical Journal, 820, 89
  * Gao et al. (2021) Gao, P., Wakeford, H. R., Moran, S. E., & Parmentier, V. 2021, Aerosols in exoplanet atmospheres, Wiley Online Library
  * Gladstone & Young (2019) Gladstone, G. R., & Young, L. A. 2019, Annual Review of Earth and Planetary Sciences, 47, 119
  * Graham et al. (2021) Graham, R., Lichtenberg, T., Boukrouche, R., & Pierrehumbert, R. T. 2021, The Planetary Science Journal, 2, 207
  * Hammond & Pierrehumbert (2017) Hammond, M., & Pierrehumbert, R. T. 2017, The Astrophysical Journal, 849, 152
  * Helling (2019) Helling, C. 2019, Annual Review of Earth and Planetary Sciences, 47, 583
  * Hirschmann (2012) Hirschmann, M. M. 2012, Earth and Planetary Science Letters, 341, 48
  * Hu et al. (2024) Hu, R., Bello-Arufe, A., Zhang, M., et al. 2024, Nature, 1
  * Ingersoll et al. (1985) Ingersoll, A. P., Summers, M. E., & Schlipf, S. G. 1985, Icarus, 64, 375
  * Ito & Ikoma (2021) Ito, Y., & Ikoma, M. 2021, Monthly Notices of the Royal Astronomical Society, 502, 750
  * Ito et al. (2015) Ito, Y., Ikoma, M., Kawahara, H., et al. 2015, The Astrophysical Journal, 801, 144
  * Joshi (2003) Joshi, M. 2003, Astrobiology, 3, 415
  * Kang et al. (2023) Kang, W., Nimmo, F., & Ding, F. 2023, The Astrophysical Journal Letters, 949, L20
  * Kite et al. (2016) Kite, E. S., Fegley Jr, B., Schaefer, L., & Gaidos, E. 2016, The Astrophysical Journal, 828, 80
  * Kopal (1954) Kopal, Z. 1954, Monthly Notices of the Royal Astronomical Society, 114, 101
  * Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Nature, 505, 69
  * Malavolta et al. (2018) Malavolta, L., Mayo, A. W., Louden, T., et al. 2018, The Astronomical Journal, 155, 107
  * Mercier et al. (2022) Mercier, S. J., Dang, L., Gass, A., Cowan, N. B., & Bell, T. J. 2022, The Astronomical Journal, 164, 204
  * Miguel et al. (2011) Miguel, Y., Kaltenegger, L., Fegley, B., & Schaefer, L. 2011, The Astrophysical Journal Letters, 742, L19
  * Moore et al. (2023) Moore, K., Cowan, N. B., & Boukaré, C.-É. 2023, Monthly Notices of the Royal Astronomical Society, 526, 6235
  * Murgas et al. (2022) Murgas, F., Nowak, G., Masseron, T., et al. 2022, Astronomy & Astrophysics, 668, A158
  * Nguyen et al. (2020) Nguyen, T. G., Cowan, N. B., Banerjee, A., & Moores, J. E. 2020, Monthly Notices of the Royal Astronomical Society, 499, 4605
  * Nguyen et al. (2022) Nguyen, T. G., Cowan, N. B., Pierrehumbert, R. T., Lupu, R. E., & Moores, J. E. 2022, Monthly Notices of the Royal Astronomical Society, 513, 6125
  * Osborn et al. (2021) Osborn, A., Armstrong, D. J., Cale, B., et al. 2021, Monthly Notices of the Royal Astronomical Society, 507, 2782
  * Pierrehumbert & Hammond (2019) Pierrehumbert, R. T., & Hammond, M. 2019, Annual Review of Fluid Mechanics, 51, 275
  * Piette et al. (2023) Piette, A. A., Gao, P., Brugman, K., et al. 2023, arXiv preprint arXiv:2306.10100
  * Powell et al. (2018) Powell, D., Zhang, X., Gao, P., & Parmentier, V. 2018, The Astrophysical Journal, 860, 18
  * Schaefer & Elkins-Tanton (2018) Schaefer, L., & Elkins-Tanton, L. T. 2018, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 376, 20180109
  * Schaefer & Fegley (2009) Schaefer, L., & Fegley, B. 2009, The Astrophysical Journal, 703, L113
  * Teske et al. (2023) Teske, J., Dang, L., Grant, D., et al. 2023, JWST Proposal. Cycle 2, 3860
  * Tsang et al. (2016) Tsang, C. C., Spencer, J. R., Lellouch, E., Lopez-Valverde, M. A., & Richter, M. J. 2016, Journal of Geophysical Research: Planets, 121, 1400
  * Winn et al. (2011) Winn, J. N., Matthews, J. M., Dawson, R. I., et al. 2011, The Astrophysical Journal Letters, 737, L18
  * Wordsworth (2015) Wordsworth, R. 2015, The Astrophysical Journal, 806, 180
  * Wordsworth & Kreidberg (2022) Wordsworth, R., & Kreidberg, L. 2022, Annual Review of Astronomy and Astrophysics, 60, 159
  * Wright et al. (2023) Wright, J. T., Tusay, N., Beatty, T. G., et al. 2023, JWST Proposal. Cycle 2, 3315
  * Yurchenko et al. (2022) Yurchenko, S. N., Tennyson, J., Syme, A.-M., et al. 2022, Monthly Notices of the Royal Astronomical Society, 510, 903
  * Zahnle et al. (2020) Zahnle, K. J., Lupu, R., Catling, D. C., & Wogan, N. 2020, The Planetary Science Journal, 1, 11
  * Zieba et al. (2022) Zieba, S., Zilinskas, M., Kreidberg, L., et al. 2022, Astronomy & Astrophysics, 664, A79
  * Zilinskas et al. (2022) Zilinskas, M., Van Buchem, C., Miguel, Y., et al. 2022, Astronomy & Astrophysics, 661, A126

## Appendix

###  5.1 Base hydrodynamical formulation

The framework for the hydrodynamical calculation is a system akin to the
shallow-water equations describing the conservation of mass (Eq. 6), momentum
(Eq. 7), and energy (Eq. 8) following Ingersoll et al. (1985):

$\partial_{x}(\rho hV)=m\ E,$ (6) $\partial_{x}(\rho
hV^{2})=-\partial_{x}\int_{z}Pdz+\tau,$ (7) $\partial_{x}\bigg{(}\rho
hV\left(\frac{V^{2}}{2}+C_{p}T\right)\bigg{)}=Q_{\text{sens}}+Q_{\text{lat}}+Q_{\text{RC}},$
(8)

where $x$ and $z$ denote the horizontal and vertical direction coordinates.
$P$, $V$, and $T$ are state variables pressure, velocity, and temperature.
Variables $\rho$ and $h$ are air density and column thickness. Parameters $m$
and $C_{p}$ are the molecular mass and heat capacity respectively. Eq. 6
describes how the change of material moving along the flow is controlled by
the mass exchange between the surface and the atmosphere via evaporation, $E$.
Eq. 7 describes how the momentum of the air parcel is controlled by the
horizontal pressure gradient force and surface drag, $\tau$. Eq. 8 describes
how the kinetic energy and internal heat energy of the parcel is controlled by
the energy fluxes which includes sensible (Q${}_{\text{sens}}$), latent
(Q${}_{\text{lat}}$), and radiative heating terms (Q${}_{\text{RC}}$) and is
expanded upon later in the text.

The mass flux, $E$ [molecules/s/m2], is the term that describes what is added
and what is removed from the atmosphere. This is controlled by evaporation:

$E=\frac{P-P_{v}(T_{s})}{m\sqrt{2\pi RT_{s}}},$ (9)

where $P_{v}$ is the saturated vapour pressure taken from Miguel et al. (2011)
and $T_{s}$ is the surface temperature. $R$ is the gas constant. To further
calculate the fluxes, such as surface drag $\tau$ and components of $Q$, we
define a mean flow velocity $V_{e}=mE/\rho$ and an eddy velocity
$V_{d}=V_{*}^{2}/V$. The frictional velocity $V_{*}$ is calculated iteratively
from the implicit relation with $V$:

$V=2.5V_{*}\log\bigg{(}\frac{9.0\ V_{*}\ H\rho}{2\eta}\bigg{)},$ (10)

where $H$ is the scale height, $\rho$ is the density, and $\eta$ is the
temperature dependent viscosity 111$\eta=1.8\times
10^{-5}(T/291)^{3/2}(411)/(T+120)$ following Sutherland’s ideal gas formula
(Castan & Menou, 2011). With these velocities, we further define advection
coefficients that describe how heat and momentum is dynamically transferred
between the surface and atmosphere, $w_{a}$ and $w_{s}$:

$\displaystyle w_{a}$ $\displaystyle=\frac{2V_{d}^{2}}{V_{e}+2V_{d}},$
$\displaystyle(V_{e}>0)$ (11a) $\displaystyle w_{s}$
$\displaystyle=\frac{V_{e}^{2}+2V_{d}V_{e}+2V_{d}^{2}}{V_{e}+2V_{d}},$ (11b)

$\displaystyle w_{a}$
$\displaystyle=\frac{V_{e}^{2}-2V_{d}V_{e}+2V_{d}^{2}}{-V_{e}+2V_{d}},$
$\displaystyle(V_{e}\leq 0)$ (12a) $\displaystyle w_{s}$
$\displaystyle=\frac{2V_{d}^{2}}{-V_{e}+2V_{d}}.$ (12b)

The above equations yield the dimensionless advection factors when there is
net evaporation ($V_{e}>0$) and deposition ($V_{e}\leq 0$). These parameters
are then used to calculate the frictional drag, $\tau$ and sensible heat,
$Q_{\text{sens}}$:

$\tau=-\rho\ V\ w_{a},$ (13)

$Q_{\text{sens}}=\rho\bigg{(}w_{s}C_{p}T_{s}+w_{a}(\frac{V^{2}}{2}+C_{p}T)\bigg{)}.$
(14)

In Eq. 14, the first term is the heat contribution from the surface and the
second term is the atmospheric advection. The latent heat term
$Q_{\text{lat}}$, is simply:

$Q_{\text{lat}}=L\ m\ C,$ (15)

where $L$ is the latent heat of vaporization and $C$ is the condensation rate.
This term is defined such that latent heat is added to the atmosphere when
condensing but is taken away from the surface when evaporating; this term will
subsequently show up in the surface energy balance equation. The last
remaining flux term is the pressure gradient force where we must integrate the
vertical pressure with respect to height from a specified altitude to
infinity:

$\int_{z=z_{0}}^{\infty}Pdz=\int_{0}^{P_{0}}\frac{\Phi}{g}dP=\frac{C_{p}}{g}\int_{0}^{P_{0}}(T_{0}-T)dP.$
(16)

Eq. 16 switches from the local $z$ altitude coordinates to $P$ levels by
substituting the geopotential variable $\Phi$. By expresing $\Phi$ as work
done by moving the air parcel vertically under the gravitational acceleration
$g$, we can relate it to enthalpy which yields a function dependent on the
pressure-temperature profile. For an adiabatic profile, completing the
vertical integration yields:

$\frac{C_{p}}{g}\int_{0}^{P_{0}}(T_{0}-T)\ dP=\frac{C_{p}\ P_{0}\ T_{0}\
\beta}{g},$ (17)

where $\beta=R/(R+C_{p})$. Eq. 17 essentially provides an analytical solution
to the pressure integral of a specific profile but one can generate another
temperature profile and solve Eq. 16 numerically.

With the exception of $Q_{\text{RC}}$, equations 9, 13, 14, and 17 show how
all of the fluxes are calculated in equations 6-8. The next step is to convert
the local coordinates to global coordinates in terms of $\theta$, the angular
distance from the substellar point. Steps detailing the full conversion can be
found in Ingersoll et al. (1985) but the net result is shown:

$\displaystyle\frac{1}{r\ \sin(\theta)}\frac{d}{d\theta}\bigg{(}\frac{V_{*}\
P_{*}\ \sin(\theta)}{g}\bigg{)}=m\ E,$ (18)

$\displaystyle\frac{1}{r\
\sin(\theta)}\frac{d}{d\theta}\bigg{(}\frac{(V_{*}^{2}+\beta\ C_{p}\
T_{*})P_{*}\ \sin(\theta)}{g}\bigg{)}=\frac{\beta\ C_{p}\ T_{*}\ P_{*}}{g\ r\
\tan(\theta)}+\tau,$ (19)

$\displaystyle\frac{1}{r\
\sin(\theta)}\frac{d}{d\theta}\bigg{(}\frac{(V_{*}^{2}/2+\ C_{p}\ T_{*})V_{*}\
P_{*}\ \sin(\theta)}{g}\bigg{)}=Q_{\text{sens}}+Q_{\text{lat}}+Q_{\text{RC}},$
(20)

The asterisk subscript on the state variables denotes its value at the of the
turbulent boundary layer. Because surface drag induces a logarithmic vertical
wind profile, the wind speed at half the scale height is closest to the
weighted mean flow of the air column. Therefore, the turbulent boundary layer
is taken to be at half the scale height to harmonize our treatment of mass,
momentum, and energy advection.

The equations are solved as a set of ordinary differential equation with
respect to $\theta$. We use the “shooting” method where we guess the boundary
conditions at the substellar point. Because we are dealing with non-linear
mass, momentum, and energy fluxes, a steady-state solution is particularly
sensitive to the initial conditions. The only definitive condition is
$V(0)=0$; for $T(0)$ and $P(0)$, we iteratively have to solve for the
conditions that is most numerically stable as described in Nguyen et al.
(2022). Once a stable solution is established, the flow inevitably reaches
supersonic speed and we linearly interpolate the transition from subsonic to
supersonic flow.

As the atmosphere gradually condenses out, there comes a point when the
atmosphere transitions into an exosphere where molecules behave less like a
continuous fluid and more like independent particles in ballistic motion. As
the flow regime changes, our calculations fail. The result is a spatially
bounded atmosphere spanning more or less a hemisphere with steady state
evaporation and condensation.

###  5.2 radiative transfer formulation

To calculate the radiative contributions to the atmosphere and surface, we
first define the vertical optical depth, $\tau_{d}$ as a function of pressure:

$\tau_{d}=x_{a}P_{*}/mg,$ (21)

where $x_{a}$ is the wavelength-dependent absorption cross-section (see Fig.
1). The absorptivity/emissivity $\epsilon$ of the overlying atmospheric column
is related to the optical depth:

$\epsilon=1-e^{-\tau_{d}}.$ (22)

The emissivity $\epsilon$ helps to calculate the individual radiative
contribution from stellar radiation absortion ($F_{\text{stel}}$), surface
blackbody radiation ($F_{\text{surf}}$), and radiative cooling
($F_{\text{RC}}$). The total heat contribution from radiative transfer,
$Q_{\text{RT}}$ becomes:

$\displaystyle Q_{\text{RT}}=$ $\displaystyle
F_{\text{stel}}+F_{\text{surf}}-2F_{\text{RC}}-F_{c}$ (23) $\displaystyle=$
$\displaystyle\int_{\lambda}\epsilon(\lambda)F_{*}(\theta,\lambda)d\lambda+\int_{\lambda}\epsilon(\lambda)B(T_{s},\lambda)d\lambda-2\int_{\lambda}\epsilon(\lambda)B(T,\lambda)d\lambda-\epsilon$

where $F_{*}$ is the stellar flux and $B$ is the Planck blackbody function
which depends on both wavelength and temperature. Since lava planets are so
close to their star, parallel ray approximations no longer hold and $F_{*}$ is
not easily solvable; we use the formulation of Kang et al. (2023) to account
for finite light source effect. Note that the radiative cooling term is
doubled to signify that the atmosphere radiates both down into the ground and
outwards into space.

Although Eq. 23 is identical to what we had in Nguyen et al. (2022) barring
cloud reflection, we have significantly improved the numerical integration of
the stellar flux and blackbody radiation. Previously, we averaged the
emissivity across the entire ultraviolet and infrared wavelength bands, while
neglecting wavelengths in the optical completely. The problem is that the
absorption cross-section of SiO varies enormously with changing temperature
and pressure. We therefore calculate absorption and emission line by line with
varying pressure and temperature ahead of the simulation. Stellar absorption
depends mostly on pressure and a function of $F_{*}(P)$ is extracted for easy
implementation. Radiative cooling is much more tricky as it depends on both
pressure and temperature, therefore a lookup table of flux is generated (see
Fig. 7 to quickly determine the radiative transfer terms.

Figure 7: Radiative cooling as a function of pressure and temperature: We
calculate the radiative cooling and seperate the fluxes into three channel: UV
(left), visibile (middle), and IR (right).