∎

11institutetext:  S.-i. Sasa 22institutetext: Department of Physics, Kyoto
University, Kyoto 606-8502, Japan 22email<EMAIL_ADDRESS>33institutetext: N. Nakagawa 44institutetext: Department of Physics, Ibaraki
University, Mito 310-8512, Japan 44email<EMAIL_ADDRESS>
# Non-equilibrium phase coexistence in boundary-driven diffusive systems

Shin-ichi Sasa  Naoko Nakagawa

###### Abstract

Liquid-gas phase coexistence in a boundary-driven diffusive system is studied
by analyzing fluctuating hydrodynamics of a density field defined on a one-
dimensional lattice with a space interval $\Lambda$. When an interface width
$\ell$ is much larger than $\Lambda$, the discrete model becomes the standard
fluctuating hydrodynamics, where the phase coexistence condition is given by
the local equilibrium thermodynamics. In contrast, when $\ell<\Lambda$, the
most probable density profile is determined by a new variational principle,
where the chemical potential at the interface is found to deviate from the
equilibrium coexistence chemical potential. This means that metastable states
at equilibrium stably appear near the interface as the influence of the
particle current. The variational function derived in the theoretical analysis
is also found to be equivalent to the variational function formulated in an
extended framework of thermodynamics called global thermodynamics. Finally,
the validity of the theoretical result is confirmed by numerical simulations.

###### Keywords:

phase coexistence, out of equilibrium, fluctuating hydrodynamics

††journal: Journal of Statistical Physics

##  1 Introduction

A rich variety of phenomena exhibit non-equilibrium phase coexistence, such as
boiling heat transfer, pattern formation in crystal growth, and motility-
induced phase separation boiling ; crystal ; Cannell ; Yoshida ; Zhong ;
Ahlers ; Urban ; mips . Although many such impressive phenomena are dynamic
and complex, a non-trivial and surprising phenomenon has been predicted in
calm and simple phase coexistence out of equilibrium. One example is the
quantitative prediction that, in liquid-gas coexistence under heat conduction,
the temperature of the liquid-gas interface is lower than the equilibrium
coexistence temperature for the pressure Global-PRL ; Global-JSP , where the
equilibrium phase coexistence occurs at the first-order phase transition point
far from the critical point. This phenomenon means that metastable states at
equilibrium are stabilized by a steady current even in the linear response
regime.

The prediction was presented in an extended framework of thermodynamics that
we call global thermodynamics. This framework was first proposed as a natural
extension of the minimum principle of free energy with the key concept of
global temperature Global-PRL . Applying the framework to a van-der Waals
fluid revealed that the temperature of the liquid-gas interface is different
from the first-order transition temperature at equilibrium. Then, the
formulation was carefully arranged so that quantitative predictions could be
made for real materials Global-JSP . The equivalence among different ensembles
was discussed, and finally, the maximum entropy principle was formulated for
enthalpy-conserving heat conduction systems Global-PRR . The entropy defined
in the formulation is found to possess a non-additive term in addition to the
space integral of the local entropy density field. This formulation enables us
to apply global thermodynamics to various systems.

There have been no experimental reports on the predictions of global
thermodynamics. Nevertheless, it is worth mentioning that numerical simulation
of the Hamiltonian Potts model in heat conduction shows results consistent
with the quantitative prediction of global thermodynamics for large enough
systems KNS . The singular nature of the phase coexistence has also been
discussed by analyzing a mesoscopic model describing the order parameter
dynamics in heat conduction Global-PRE . However, this analysis involves some
phenomenological assumptions in calculating the singular part. Furthermore,
the model is too complicated to extract the microscopic mechanism of the
deviation in interface temperature from the equilibrium phase coexistence
temperature.

On the basis of this background, we study a simple system that exhibits non-
equilibrium phase coexistence. We consider a system in which the number
density field is a single dynamical variable and the density field is directly
driven by a boundary condition of the chemical potential, where the
temperature is given as a constant parameter of the system. The stochastic
time evolution of the density field is described in terms of a discrete form
of fluctuating hydrodynamics with a space interval $\Lambda$. When the width
of an interface in phase coexistence $\ell$ is much larger than $\Lambda$, the
model is equivalent to the standard fluctuating hydrodynamics Schmitz ; Eyink
; Spohn ; Bertini-RMP . Then, local fluctuation properties of thermodynamic
quantities are described by local equilibrium distribution Kuramoto . In
contrast, when $\ell\ll\Lambda$, the local distribution may be out of
equilibrium Yana . For this case, we derive the variational function that
determines the most probable density profile. We then find that the chemical
potential at the interface of the most probable profile deviates from the
equilibrium coexistence value. This means that metastable states at
equilibrium stably appear near the interface of the driven system. The formula
describing the deviation takes the same form as those phenomenologically
predicted by global thermodynamics. Indeed, we can formulate extended
thermodynamic relations for this system by using the variational function,
which is equivalent to global thermodynamics. We also confirm the validity of
the theoretical calculation by numerical simulations.

The rest of this paper is organized as follows. In Sec. 2, we introduce a
model we study in this paper. We then review phase coexistence conditions for
equilibrium systems, and summarizes basic issues for non-equilibrium phase
coexistence. In Sec. 3, we derive a variational function by analyzing the
Zubarev-McLennan representation of the stationary distribution. In Sec. 4, we
rewrite the variational equation as the form giving the chemical potential at
the interface. In Sec. 5, we formulate extended thermodynamic relations based
on the variational function and discuss the relationship with global
thermodynamics. In Sec. 6, we show results of numerical simulations. Section 7
provides some concluding remarks.

##  2 Setup

###  2.1 Model

Figure 1:  Schematics of a system we study. A one-dimensional density field is
defined on a ring with a battery at $x=0$.

We consider a collection of stochastic and diffusive particles in a closed
tube which are driven by an external battery at one surface $x=0$. See Fig. 1
for the illustration of the setup. We describe the system by an averaged
particle density $\rho(x)$ defined in a one-dimensional region $[0,L]$. More
precisely, we define

$\rho(x)\equiv\frac{1}{A}\int dydz~{}\rho(x,y,z)$ (2.1)

for the three-dimensional particle density $\rho(x,y,z)$, where the area of
the cross section of the tube is denoted by $A$. We then assume a standard
continuum description of fluctuating hydrodynamics of $\rho(x)$. The density
field $\rho(x,t)$ satisfies the continuity equation

$\partial_{t}\rho+\partial_{x}j=0.$ (2.2)

Based on the mean field picture in the cross section, we assume the current as

$j(x,t)=-\sigma(\rho(x))\left[\partial_{x}\frac{\delta{\cal
F}}{\delta\rho(x)}-\phi\delta(x)\right]+\sqrt{\frac{2\sigma(\rho(x))T}{A}}\cdot\xi(x,t).$
(2.3)

$T$ is the temperature, $\sigma(\rho)$ is a conductivity as a smooth function
of $\rho$, and $\phi$ represents the voltage of a battery located at $x=0$.
$\xi$ is Gaussian-white noise satisfying $\left\langle\xi\right\rangle=0$ and

$\left\langle\xi(x,t)\xi(x^{\prime},t^{\prime})\right\rangle=\delta(x-x^{\prime})\delta(t-t^{\prime}).$
(2.4)

The symbol “$\cdot$” in front of $\xi$ in (2.3) represents the Stratonovich
product in space and the Ito product in time. Here, because the variance of
the surface average $\Xi(x,t)\equiv\int dydz\xi(x,y,z,t)/A$ of the three-
dimensional Gaussian-white noise $\xi(x,y,z,t)$ with unit variance is $1/A$,
we set $\Xi(x,t)=\xi(x,t)/\sqrt{A}$, which determines the $A$ dependence of
(2.3). The free energy functional ${\cal F}$ of the density profile
${\bm{\rho}}=(\rho(x))_{0\leq x\leq L}$ is expressed as

${\cal
F}({\bm{\rho}})=\int_{0}^{L}dx\left[f(\rho(x))+\frac{\kappa}{2}(\partial_{x}\rho)^{2}\right].$
(2.5)

As a specific example, we consider the case that

$f(\rho)=-\frac{1}{2}(\rho-1.5)^{2}+\frac{1}{4}(\rho-1.5)^{4}$ (2.6)

by introducing dimensionless length and energy in this form. Note that our
argument below is independent of the specific form if $f(\rho)$ contains two
local minima. The length unit is assumed to be the order of particle
distances. $\kappa$ characterizes the interface free energy, which is relevant
when $\partial_{x}\rho$ is large. In particular, when $0.5<\bar{\rho}<2.5$ in
noiseless equilibrium systems with $T=0$ and $\phi=0$, phase coexistence
occurs with two interfaces. Then, $\sqrt{\kappa}$ determines the interface
width in the phase coexistence.

Now, we notice that there is a cutoff length $\Lambda$ of the continuum
description. Because the noise is assumed to be white in space, $\Lambda$
should be larger than the microscopic length, which is set to the order of
unity. On the other hand, it is obvious that $\Lambda$ should be much smaller
than the system size $L$. In many cases, the calculation result of fluctuating
hydrodynamics is independent of the cutoff $\Lambda$ even in the limit
$\Lambda\to 0$, while there is a case where a singular cutoff dependence is
observed Kado . Here, let us recall that the width of interfaces in phase
coexistence is estimated as a microscopic length. Thus, this may be smaller
than the cutoff length of fluctuating hydrodynamics. Such a case cannot be
studied by the continuum model. We thus need to propose and analyze a discrete
model in which the microscopic cutoff $\Lambda$ is explicitly introduced.

With this background, we consider a sequence of $N$-boxes in a one-dimensional
ring. Let $\rho_{i}$ be the density of particles at the $i$-th box, where
$1\leq i\leq N$. Mathematically, $\rho_{i}$ is defined on the $i$-th site in
the one-dimensional lattice $\\{i|1\leq i\leq N,\ i\in\mathbb{Z}\\}$ with the
periodic boundary condition $\rho_{0}\equiv\rho_{N}$ and
$\rho_{N+1}=\rho_{1}$. The horizontal size of the box is denoted by $\Lambda$
and the cross section area of the box is $A$. The system size $L$ is given by
$L=\Lambda N$. The free energy functional given in (2.5) is then replaced with

${\cal
F}({\bm{\rho}})=\Lambda\sum_{i=1}^{N}\left[f(\rho_{i})+\frac{\kappa}{2\Lambda^{2}}(\rho_{i+1}-\rho_{i})^{2}\right]$
(2.7)

for ${\bm{\rho}}=(\rho_{i})_{1\leq i\leq N}$. The density $\rho_{i}$ satisfies
the continuity equation

$\frac{d\rho_{i}}{dt}+\frac{j_{i}-j_{i-1}}{\Lambda}=0$ (2.8)

with $j_{0}=j_{N}$. The current $j_{i}$ is defined on the $i$-th bond
connecting the $i$-th site and the $(i+1)$-th site. Using the generalized
chemical potential $\tilde{\mu}_{i}$ defined by

$\tilde{\mu}_{i}\equiv\frac{1}{\Lambda}\frac{\partial{\cal
F}}{\partial\rho_{i}},$ (2.9)

we replace the current given in (2.3) by

$j_{i}(t)=-\frac{\sigma(\rho_{i}^{\rm
m})}{\Lambda}(\tilde{\mu}_{i+1}-\tilde{\mu}_{i}-\phi\delta_{i,N})+\sqrt{\frac{2\sigma(\rho_{i}^{\rm
m})T}{A\Lambda}}\cdot\xi_{i}(t),$ (2.10)

where we set $\rho^{\rm m}_{i}=(\rho_{i}+\rho_{i+1})/2$ to satisfy the
detailed balance condition. $\xi_{i}(t)$ is Gaussian-white noise satisfying
$\left\langle\xi_{i}\right\rangle=0$ and

$\left\langle\xi_{i}(t)\xi_{j}(t^{\prime})\right\rangle=\delta_{ij}\delta(t-t^{\prime}).$
(2.11)

The $\Lambda$ dependence of the noise intensity in (2.10) is understood from
the replacement of $\delta(x-x^{\prime})$ with $\delta_{ij}/\Lambda$.
Explicitly, $\tilde{\mu}_{i}$ is written as

$\tilde{\mu}_{i}=\mu(\rho_{i})-\frac{\kappa}{\Lambda^{2}}(\rho_{i+1}+\rho_{i-1}-2\rho_{i})$
(2.12)

with $\mu(\rho_{i})=\partial f(\rho_{i})/\partial\rho_{i}$. The total number
of particles $\sum_{i=1}^{N}\rho_{i}=\bar{\rho}N$ is conserved in the time
evolution. The average density $\bar{\rho}$ is a parameter of the system.

The steady state of this system is characterized by five parameters
$(\kappa_{\Lambda},T_{\rm eff},\phi,\bar{\rho},N)$, where

$\kappa_{\Lambda}\equiv\frac{\kappa}{\Lambda^{2}},\qquad T_{\rm
eff}\equiv\frac{T}{A}.$ (2.13)

That is, systems with the same values of $(\kappa_{\Lambda},T_{\rm
eff},\phi,\bar{\rho},N)$ exhibit the same steady state. In the argument below,
the $\Lambda$ dependence appear only through $\kappa_{\Lambda}$ dependence.
When $\kappa_{\Lambda}\gg 1$ and $N\to\infty$, the system behavior of (2.8)
and (2.10) is understood by analyzing (2.2) and (2.3) because (2.8) and (2.10)
correspond to an accurate approximation of (2.2) and (2.3). In contrast, the
system behavior for the case $\kappa_{\Lambda}<1$ cannot be described by (2.2)
and (2.3). For such cases, we have to analyze the discrete model with focusing
on the limiting case that $\kappa_{\Lambda}\ll 1$. See Fig. 2 for the
illustration of two limiting cases. Finally, we note that $A\gg 1$ because $A$
is the square of a macroscopic length. We thus consider the weak noise limit
$T_{\rm eff}\to 0$ for the steady state realized in the limit $t\to\infty$.

Figure 2:  Illustration of two limiting cases (a) $\kappa_{\Lambda}\ll 1$ and
(b) $\kappa_{\Lambda}\gg 1$.

###  2.2 Review of equilibrium phase coexistence

We review phase coexistence states for the equilibrium system with $\phi=0$.
The stationary distribution of ${\bm{\rho}}$ is given by

${\cal P}_{{\rm
eq}}({\bm{\rho}};\bar{\rho})=\frac{1}{Z}\exp\left[-\frac{1}{T_{\rm eff}}{\cal
F}({\bm{\rho}})\right]\delta\left(\sum_{i=1}^{N}\rho_{i}-\bar{\rho}N\right).$
(2.14)

See Appendix A for the derivation. The most probable profile in the weak noise
limit $T_{\rm eff}\to 0$ is determined as the density profile that minimizes
${\cal F}({\bm{\rho}})$ defined by (2.7). The variational equation is obtained
as

$\tilde{\mu}_{i}=\tilde{\mu}_{i+1},$ (2.15)

which means that $\tilde{\mu}_{i}$ is a constant independent of $i$.

We study two limiting cases $\kappa_{\Lambda}\gg 1$ and $\kappa_{\Lambda}\ll
1$. First, for the case $\kappa_{\Lambda}\gg 1$, we can derive the phase
coexistence condition from the analysis of the continuum limit of (2.15) with
$N\to\infty$, which is given in (B.2). Importantly, the solution of the
variational equation is unique under the boundary condition that
$\partial_{x}\rho(0)=\partial_{x}\rho(L)=0$ with $\rho(0)>\rho(L)$. As shown
in Appendix B, phase coexistence occurs when $\bar{\rho}$ satisfies $\rho^{\rm
G}_{\rm c}\leq\bar{\rho}\leq\rho^{\rm L}_{\rm c}$, where $\rho^{\rm L}_{\rm
c}$ and $\rho^{\rm G}_{\rm c}$ are determined by

$\mu(\rho^{\rm L}_{\rm c})=\mu(\rho^{\rm G}_{\rm c}),\qquad p(\rho^{\rm
L}_{\rm c})=p(\rho^{\rm G}_{\rm c})$ (2.16)

with the thermodynamic pressure $p(\rho)$ defined by

$p(\rho)\equiv\rho\mu(\rho)-f(\rho).$ (2.17)

Hereafter, the equilibrium value of the coexistence chemical potential is
denoted by $\mu_{c}\equiv\mu(\rho^{\rm L}_{\rm c})=\mu(\rho^{\rm G}_{\rm c})$.
We assume $\rho^{\rm G}_{\rm c}<\rho^{\rm L}_{\rm c}$ without loss of
generality. Then, $\rho^{\rm L}_{\rm c}$ and $\rho^{\rm G}_{\rm c}$ represent
the densities of the liquid and gas in the phase coexistence, respectively.
Furthermore, once $\rho^{\rm L}_{\rm c}$ and $\rho^{\rm G}_{\rm c}$ are
obtained, the fraction of the liquid region $X^{\rm eq}$ is uniquely
determined from $\Lambda\sum_{i}\rho_{i}=\bar{\rho}L$, which is expressed by

$\rho^{\rm L}_{\rm c}X^{\rm eq}+\rho^{\rm G}_{\rm c}(1-X^{\rm
eq})=\bar{\rho}.$ (2.18)

See Fig. 3(c) for the profile.

Figure 3:  Illustration of metastable profile and equilibrium profile.

For the other case $\kappa_{\Lambda}\ll 1$, which we mainly study in this
paper, the derivation method in the continuum limit cannot be applied.
However, even in this case, the variational principle determines the most
probable profile. The variational equation (2.15) is rewritten as
$\mu_{i}=\mu_{i+1}$ in the limit $\kappa_{\Lambda}\to 0$ so that the chemical
potential is uniform. In this case, there is a one-parameter family of
solutions $\rho_{i;X}^{\phi=0}$ characterized by an interface position $X$. We
here attempt to write the solution $\rho_{i;X}^{\phi=0}$ explicitly. First,
let $\mu_{X}^{\phi=0}$ denote the uniform value of the chemical potential for
the solution. Referring to Fig. 3(a), we set ${\rho}^{\rm L}_{X}$ and
${\rho}^{\rm G}_{X}$ as those satisfying

$\mu({\rho}^{\rm L}_{X})=\mu({\rho}^{\rm G}_{X})=\mu_{X}^{\phi=0}$ (2.19)

with ${\rho}^{\rm G}_{X}<{\rho}^{\rm L}_{X}$. We can then express the solution
as

$\rho_{i;X}^{\phi=0}={\rho}^{\rm L}_{X}\chi(i\in[1,NX])+{\rho}^{\rm
G}_{X}\chi(i\in[NX,N])$ (2.20)

with a half-integer $NX$ and $0\leq X\leq 1$, where $\chi(P)$ represents the
characteristic function that takes $1$ if $P$ is true and 0 otherwise. Then,
for a given $X$, the condition

$X{\rho}^{\rm L}_{X}+(1-X){\rho}^{\rm G}_{X}=\bar{\rho}$ (2.21)

determines the value of $\mu_{X}^{\phi=0}$.

Note that the value of $X$ is not determined from the variational equation
(2.15) with $\kappa_{\Lambda}\to 0$. Physically, the solutions
${\bm{\rho}}^{\phi=0}_{X}$ form a family of metastable states characterized by
$X$. The most probable value of $X$, which is denoted by $X_{*}$, is derived
from the minimum free energy principle formulated as follows. We first define
the variational function

${\cal F}_{\rm eq}(X;\bar{\rho})\equiv{\cal F}({\bm{\rho}}_{X}^{\phi=0}).$
(2.22)

Then, the minimum free energy principle means

${\cal F}_{\rm eq}(X_{*};\bar{\rho})=\min_{X}{\cal F}_{\rm eq}(X;\bar{\rho}),$
(2.23)

where, in the present case, $X_{*}$ satisfies

$\left.\frac{d{\cal F}_{\rm eq}(X;\bar{\rho})}{dX}\right|_{X=X_{*}}=0.$ (2.24)

Substituting the expression

${\cal F}({\bm{\rho}}_{X}^{\phi=0})=LXf({\rho}^{\rm
L}_{X})+L(1-X)f({\rho}^{\rm G}_{X})$ (2.25)

into (2.24), we obtain

$f({\rho}^{\rm L}_{X_{*}})-f({\rho}^{\rm
G}_{X_{*}})=\mu_{X_{*}}^{\phi=0}({\rho}^{\rm L}_{X_{*}}-{\rho}^{\rm
G}_{X_{*}})$ (2.26)

using (2.19) and (2.21). Comparing (2.19) and (2.26) with (2.16), we find

$\mu_{c}=\mu_{X_{*}}^{\phi=0},\quad\rho^{\rm L}_{\rm c}={\rho}^{\rm
L}_{X_{*}},\qquad\rho^{\rm G}_{\rm c}={\rho}^{\rm G}_{X_{*}}.$ (2.27)

Moreover, $X_{*}$ is equal to $X^{\rm eq}$. The result means that the
equilibrium states for the two limiting cases, $\kappa_{\Lambda}\gg 1$ and
$\kappa_{\Lambda}\ll 1$, are equivalent.

###  2.3 Preliminaries for non-equilibrium phase coexistence

We concentrate on the case that $\bar{\rho}$ satisfies $\rho^{\rm G}_{\rm
c}<\bar{\rho}<\rho^{\rm L}_{\rm c}$ where liquid and gas coexist. When the
voltage of the battery $\phi$ is positive, the stationary distribution is not
written as the canonical form (2.14). Therefore, we do not have a general
variational principle for determining the most probable profile for $\phi>0$.
Nevertheless, we divide the problem for determining the most probable profile
into two steps. As the first step, we consider stationary solutions of the
deterministic equations without noise. If the stationary solution is unique,
this is the most probable profile in the weak noise limit. In contrast, if
stationary solutions form a one parameter family, we proceed to the second
step, where we will formulate a variational principle for selecting the most
probable profile among the stationary solutions as we have examined in (2.23)
for the equilibrium system. In this section, we focus on the first step.

We consider stationary solutions of the deterministic equation given in (2.2)
and (2.3) with $T_{\rm eff}=0$. Let $J$ be the steady particle current induced
by the battery. Setting $j_{i}(t)=J$ for all $i$ in (2.3), we obtain the
conduction equation

$\frac{1}{\Lambda}\left(\tilde{\mu}_{i+1}-\tilde{\mu}_{i}+\phi\delta_{i,N}\right)=-\frac{J}{\sigma(\rho_{i}^{\rm
m})}.$ (2.28)

The uniformity of $J$ yields

$J=\phi\left(\Lambda\sum_{i=1}^{N}\frac{1}{\sigma(\rho_{i}^{\rm
m})}\right)^{-1}.$ (2.29)

When $\kappa_{\Lambda}\gg 1$ and $N\to\infty$, we can analyze the continuum
limit of (2.28). As shown in Appendix B, we find that the solution of the
deterministic equation is unique and that the chemical potential and pressure
are continuous at the interface in the space $\hat{x}=x/L$ with $L\to\infty$.
We can then conclude that the chemical potential at the interface is $\mu_{\rm
c}$. See Appendix B for the derivation of this result.

In contrast, when $\kappa_{\Lambda}\ll 1$ and $N\to\infty$, there is a family
of stationary solutions of the conduction equation (2.28) with
$\kappa_{\Lambda}\to 0$. In the equilibrium system, they correspond to
metastable states characterized by (2.19), as shown in Fig. 3(b), and the
metastable states were candidates for the most probable profile. We thus
expect that a family of stationary solutions $\rho_{i;X}^{\phi}$ for $\phi>0$
correspond to metastable states among which the most probable profile is
selected. In the remainder of this section, we express the metastable states
explicitly as a preliminary for the analysis in the next section.

We concentrate on small $\phi$ and ignore the contribution of $O(\phi^{2})$.
The solution $\rho_{i;X}^{\phi}$ should be given as a perturbation of the
equilibrium solution $\rho_{i;X}^{\phi=0}$ given in (2.20). We thus assume
that $|\rho_{i;X}^{\phi}-\rho^{\rm L}_{\rm c}|<|\rho_{i;X}^{\phi}-\rho^{\rm
G}_{\rm c}|$ for $1\leq i<NX$, and $|\rho_{i;X}^{\phi}-\rho^{\rm G}_{\rm
c}|<|\rho_{i;X}^{\phi}-\rho^{\rm L}_{\rm c}|$ for $NX<i\leq N$. Letting

$\sigma^{L}=\sigma(\rho^{\rm L}_{\rm c}),\qquad\sigma^{G}=\sigma(\rho^{\rm
G}_{\rm c}),$ (2.30)

the conduction equation (2.28) with $\kappa_{\Lambda}\to 0$ results in

$\displaystyle\mu(\rho_{i;X}^{\phi})=-\frac{1}{\sigma^{L}}JL\frac{i-1}{N}+\mu_{1:X}^{\phi}+O(\phi^{2}),\qquad(1\leq
i\leq NX),$ (2.31)
$\displaystyle\mu(\rho_{i;X}^{\phi})=-\frac{1}{\sigma^{G}}JL\frac{i-N}{N}+\mu_{N;X}^{\phi}+O(\phi^{2}),\qquad(NX<i\leq
N),$ (2.32)

where

$\mu_{1;X}^{\phi}=\mu_{N;X}^{\phi}+\phi.$ (2.33)

From the direct calculation, we also obtain the difference between the
chemical potentials of two adjacent sites over the interface as

$\mu(\rho_{X-1/(2N);X}^{\phi})-\mu(\rho_{X+1/(2N);X}^{\phi})=O(\phi/N).$
(2.34)

Finally, the condition

$\frac{1}{N}\sum_{i=1}^{N}\rho_{i;X}^{\phi}=\bar{\rho}$ (2.35)

uniquely determines $\mu_{1;X}^{\phi}$ and $\rho_{1;X}^{\phi}$ for a given
$X$. We can easily confirm that $\lim_{\phi\to 0}\rho_{i;X}^{\phi}$ is
equivalent to (2.20).

Figure 4:  Schematic profiles of $\mu_{X}^{\phi}(\hat{x})$ and
$\rho_{X}^{\phi}(\hat{x})$.

###  2.4 Singular continuum description

We introduce a continuum description using a real variable $\hat{x}=x/L$
defined on the interval $[0,1]$ by taking the limit $N\to\infty$. The
mathematical formalize of this limit is beyond the scope of the present paper.
However, to simplify the calculation in the subsequent sections, we naively
introduce the continuum description based on expected behaviors of
$\rho_{i}(t)$ for large $N$. For a discrete variable
${\bm{\rho}}=(\rho_{i}(t))_{i=1}^{N}$ at time $t$, we first define
$\rho_{N}(\hat{x},\hat{t})$ as the piece-wise linear function obtained by
connecting two consecutive points $(i/N,\rho_{i}(t))$ and
$((i+1)/N,\rho_{i+1}(t))$ for $0\leq i\leq N$ in the $(\hat{x},\rho)$ space,
where we set $\hat{t}=t/L^{2}$ for later convenience. For sufficiently large
$N$, we expect that there exists an almost continuous function
$\rho(\hat{x},\hat{t})$ such that
$|\rho(\hat{x},\hat{t})-\rho_{N}(\hat{x},\hat{t})|=O(1/N)$. The chemical
potential $\mu(\hat{x},\hat{t})$ in the continuum description is defined from
$\mu_{i}(t)$ by the same procedure and it is expected that
$\mu(\hat{x},\hat{t})=\mu(\rho(\hat{x},\hat{t}))$. The definition of
$j(\hat{x},\hat{t})$ is slightly different from $\rho(\hat{x},t)$ and
$\mu(\hat{x},t)$, because the current $j_{i}(t)$ is defined on the $i$-th bond
connecting the $i$ site and $i+1$ site as seen in (2.8) and (2.10).
Concretely, we define $j(\hat{x},\hat{t})$ by using the piece-wise linear
function obtained by connecting points $(i/N+1/2N,j_{i})$ in the $(\hat{x},j)$
space. Then, the continuum limit of (2.8) and (2.10) with $\kappa_{\Lambda}\to
0$ is expressed by

$\partial_{\hat{t}}\rho(\hat{x},\hat{t})+\partial_{\hat{x}}j(\hat{x},\hat{t})=0,$
(2.36)

and

$\displaystyle
j(\hat{x},\hat{t})=-\sigma(\rho(\hat{x},\hat{t}))\left[\frac{\partial\mu(\hat{x},\hat{t})}{\partial\hat{x}}-\phi\delta(\hat{x})\right]+\sqrt{\frac{2\sigma(\rho(\hat{x},\hat{t}))T_{\rm
eff}}{L}}\cdot\hat{\xi}(\hat{x},\hat{t}),$ (2.37)

where $\hat{\xi}(\hat{x},\hat{t})$ satisfies

$\left\langle\hat{\xi}(\hat{x},\hat{t})\xi(\hat{x}^{\prime},\hat{t}^{\prime})\right\rangle=\delta(\hat{x}-\hat{x}^{\prime})\delta(\hat{t}-\hat{t}^{\prime}).$
(2.38)

Apparently, (2.36) and (2.37), which we call the singular continuum
description, take the same form as a singular case of the standard continuum
description (2.2) and (2.3) with $\kappa\to 0$ by setting $\hat{x}=x/L$ and
$\hat{t}=t/L^{2}$. We here notice the order of limits in both the
descriptions. In the singular continuum description, we first take
$\kappa/\Lambda^{2}\to 0$, and then $N\to\infty$ with $L$ fixed. In the
standard continuum description, we take $\kappa/\Lambda^{2}\to\infty$ and
$N\to\infty$ with $L$ fixed. Then, as a singular case of the standard
continuum description, we consider the limit $\kappa\to 0$. The behavior of
the two descriptions are rather different, but we do not discuss the
difference anymore. In the argument below, we will consider only the singular
continuum description.

The stationary solution of the deterministic equation $\rho_{i;X}^{\phi}$ and
the corresponding chemical potential $\mu(\rho_{i;X}^{\phi})$ are expressed as
$\rho_{X}^{\phi}(\hat{x})$ and $\mu_{X}^{\phi}(\hat{x})$ in the singular
continuum description. Note that $\rho_{X}^{\phi}(\hat{x})$ is discontinuous
at $x=X$, while $\mu_{X}^{\phi}(\hat{x})$ is continuous at $x=X$ as is seen
from (2.34). The chemical potential at the interface $x=X$ is denoted by

$\mu_{X}^{\rm I}\equiv\mu_{X}^{\phi}(X).$ (2.39)

As shown in Fig. 4, $\rho_{X}^{\phi}(\hat{x})$ is a piece-wise continuous
function and that $\mu_{X}^{\phi}(\hat{x})$ is a piece-wise smooth function
with singular points at $\hat{x}=0$ and $\hat{x}=X$. Furthermore, we rewrite
(2.31) and (2.33) as

$\mu_{X}^{\phi}(\hat{x})=-\frac{1}{\sigma^{L}}JL\hat{x}+\mu_{X}^{\phi}(0)+O(\phi^{2})$
(2.40)

for $0\leq\hat{x}<X$, and

$\mu_{X}^{\phi}(\hat{x})=-\frac{1}{\sigma^{G}}JL(\hat{x}-1)+\mu_{X}^{\phi}(1)+O(\phi^{2})$
(2.41)

for $X<\hat{x}\leq 1$. We also have

$\mu_{X}^{\phi}(0)=\mu_{X}^{\phi}(1)+\phi.$ (2.42)

Here, from (2.29), we obtain

$JL=\phi\left(\frac{X}{\sigma^{\rm L}}+\frac{1-X}{\sigma^{\rm
G}}\right)^{-1}.$ (2.43)

Then, the density field $\rho_{X}^{\phi}(\hat{x})$ for $0\leq\hat{x}<X$ is
determined from (2.40) with
$\mu_{X}^{\phi}(\hat{x})=\mu(\rho_{X}^{\phi}(\hat{x}))$ and
$|\rho_{X}^{\phi}(\hat{x})-\rho^{\rm L}_{\rm
c}|<|\rho_{X}^{\phi}(\hat{x})-\rho^{\rm G}_{\rm c}|$. Similarly, the density
field $\rho_{X}^{\phi}(\hat{x})$ for $X<\hat{x}\leq 1$ is determined from
(2.41) with $\mu_{X}^{\phi}(\hat{x})=\mu(\rho_{X}^{\phi}(\hat{x}))$ and
$|\rho_{X}^{\phi}(\hat{x})-\rho^{\rm G}_{\rm
c}|<|\rho_{X}^{\phi}(\hat{x})-\rho^{\rm L}_{\rm c}|$. Note that
$\rho_{X}^{\phi}(\hat{x})$ and $\partial_{x}\mu_{X}^{\phi}(\hat{x})$ are
discontinuous at $\hat{x}=0$ and $\hat{x}=X$. For given $X$ and system
parameters $(\phi,\bar{\rho})$, $\mu_{X}^{\phi}(\hat{x})$ and
$\rho_{X}^{\phi}(\hat{x})$ are uniquely determined from the condition

$\int_{0}^{1}d\hat{x}~{}\rho_{X}^{\phi}(\hat{x})=\bar{\rho}.$ (2.44)

##  3 Variational function

To simplify the notation, we use $x$ and $t$ for $\hat{x}$ and $\hat{t}$ in
this and next sections. In the previous section, we determined the candidates
of steady density profile $\rho_{X}^{\phi}(x)$ characterized by $X$ in the
weak noise limit $T_{\rm eff}\to 0$. To determine the most probable density
profile among them for given system parameters $(\bar{\rho},\phi)$, in Sec.
3.1, we derive a variational function using the Zubarev-McLennan
representation of the steady state. The variational function includes a time
integral of the current at $x=0$. After confirming some basic issues and
assumptions in Sec. 3.2, we calculate the time integral of the current in Sec.
3.3. The result is presented in Sec. 3.4.

###  3.1 Stationary distribution

We consider the stationary distribution ${\cal P}_{\rm
ss}({\bm{\rho}};\bar{\rho},\phi)$ of density field ${\bm{\rho}}$. When
$\phi=0$, the stationary distribution is given by (2.14). However, the
stationary distribution for the system with $\phi>0$ is not generally
obtained. Nevertheless, in the linear response regime out of equilibrium,
there is a useful expression called the Zubarev-McLennan representation

${\cal P}_{{\rm ss}}({\bm{\rho}};\bar{\rho},\phi)={\cal P}_{{\rm
eq}}({\bm{\rho}};\bar{\rho})\exp\left[-\frac{\phi\left\langle
Q\right\rangle^{{\rm eq}}_{{\bm{\rho}}}+O(\phi^{2})}{T_{\rm eff}}\right]$
(3.1)

with

$Q=\int_{0}^{\infty}dt~{}j_{N}(t),$ (3.2)

where $j_{N}(t)$ is a fluctuating current at the $N$-th bond, which is defined
in (2.10). The $N$-th bond is the only bond on which the driving $\phi$ is
imposed. $\left\langle\ \right\rangle^{{\rm eq}}_{{\bm{\rho}}}$ denotes the
conditioned expectation for the equilibrium path ensemble provided that the
initial density profile is given by the specified ${\bm{\rho}}$ as the
argument of the stationary distribution Zubarev ; McLennan ; Crooks ; fh-
zubarev ; KN ; Maes-rep . See Appendix C for the derivation of (3.1).

For the general expression (3.1) with (3.2), we take the limit
$\kappa_{\Lambda}\to 0$ and consider the singular continuum description
introduced in the previous section. From (2.14) and (3.1), we have

${\cal P}_{{\rm ss}}({\bm{\rho}};\bar{\rho},\phi)=\exp\left[-\frac{1}{T_{\rm
eff}}({\cal F}({\bm{\rho}})+\phi\left\langle Q\right\rangle^{{\rm
eq}}_{{\bm{\rho}}}+O(\phi^{2})+{\rm
const})\right]\delta\left(\int_{0}^{1}dx\rho(x)-\bar{\rho}\right)$ (3.3)

with

$Q=\int_{0}^{\infty}dt~{}j(0,t).$ (3.4)

Note that $j_{N}(t)$ becomes $j(0,t)=j(1,t)$ in the singular continuum
description. The most probable density profile is identified as
$\rho_{X_{*}}^{\phi}(x)$ that minimizes $-\log{\cal P}_{\rm
ss}({\bm{\rho}}_{X}^{\phi};\bar{\rho},\phi)$. This is the variational
principle to determine $X_{*}$. From (3.3), we obtain the variational function
as

${\cal F}_{\rm ss}(X;\bar{\rho},\phi)={\cal
F}({\bm{\rho}}_{X}^{\phi})+\phi\left\langle
Q\right\rangle_{{\bm{\rho}}_{X}^{\phi}}^{\rm eq}+O(\phi^{2}).$ (3.5)

We then find $X_{*}$ as the special value of $X$ that minimizes the
variational function ${\cal F}_{\rm ss}(X;\bar{\rho},\phi)$ with
$(\bar{\rho},\phi)$ fixed. That is, $X_{*}$ is determined as

${\cal F}_{\rm ss}(X_{*};\bar{\rho},\phi)=\min_{X}{\cal F}_{\rm
ss}(X;\bar{\rho},\phi).$ (3.6)

We thus need to calculate $\left\langle
Q\right\rangle_{{\bm{\rho}}_{X}^{\phi}}^{\rm eq}$.

###  3.2 Preliminaries for the calculation $\left\langle
Q\right\rangle_{{\bm{\rho}}_{X}^{\phi}}^{\rm eq}$

Figure 5:  Schematic relationship between $\rho_{X}^{\phi}(x)$ and
$\rho_{X}^{\phi=0}(x)$.

To calculate $\left\langle Q\right\rangle_{{\bm{\rho}}_{X}^{\phi}}^{\rm eq}$,
we have to know the typical time evolution of $\rho(x,t)$ starting from
$\rho_{X}^{\phi}(x)$ at $t=0$ under the equilibrium condition $\phi=0$ in the
weak noise limit $T_{\rm eff}\to 0$, where $T_{\rm eff}\to 0$ is taken after
$t\to\infty$ is considered. Here, from (3.5), we find that only
$\phi$-independent terms of $\left\langle
Q\right\rangle_{{\bm{\rho}}_{X}^{\phi}}^{\rm eq}$ are necessary for the
calculation because $\phi$-dependent terms are absorbed into $O(\phi^{2})$. We
then notice expansions

$\mu_{X}^{\phi}(x)=\mu_{X}^{\phi=0}(x)+O(\phi),\qquad\rho_{X}^{\phi}(x)=\rho_{X}^{\phi=0}(x)+O(\phi).$
(3.7)

See Fig. 5 for illustration of $\rho_{X}^{\phi}(x)$ and
$\rho_{X}^{\phi=0}(x)$. See also the sentence involving (2.21) for
$\mu_{X}^{\phi=0}$. Using these relations, we find

$\left\langle Q\right\rangle_{{\bm{\rho}}_{X}^{\phi}}^{\rm eq}=\left\langle
Q\right\rangle_{{\bm{\rho}}_{X}^{\phi=0}}^{\rm eq}+O(\phi).$ (3.8)

That is, we study the typical time evolution of the density field starting
from $\rho_{X}^{\phi=0}(x)$. Since we study the case $\kappa_{\Lambda}\to 0$,
$\rho_{X}^{\phi=0}(x)$ represents a metastable profile that does not evolve in
time without noises. Nevertheless, weak noise slowly drives the metastable
profile $\rho_{X}^{\phi=0}(x)$ to the equilibrium state characterized by
(2.14). One may conjecture that the most probable time evolution from
$\rho_{X}^{\phi=0}(x)$ in the weak noise limit is described by the relaxation
process to the equilibrium profile $\rho_{X_{*}}^{\phi=0}(x)$ in Fig. 3(c)
from the metastable profile in Fig. 3(b). However, this is not correct because
of the space translational symmetry of equilibrium systems. Note that
metastable profiles for a given width $X$ of the liquid region form a one-
parameter family of profiles obtained by any space translation of
$\rho_{X}^{\phi=0}(x)$, and similarly equilibrium profiles also form a one-
parameter family. Thus, stochastic dynamics in this neutral direction are
equally probable, which are represented by Brownian motion of the profile with
the width of the liquid region fixed. In other words, most probable process is
not uniquely determined in the weak noise limit. We have to consider a
collection of trajectories that are equally probable and dominantly contribute
to $\left\langle Q\right\rangle_{{\bm{\rho}}_{X}^{\phi=0}}^{\rm eq}$. We call
this collection the highly probable path ensemble to distinguish it with the
most probable process.

Here, we concretely describe the highly probable path ensemble starting from
the initial condition $\rho_{X}^{\phi=0}(x)$. We assume that each trajectory
in the highly probable path ensemble satisfies the following two conditions.
First, the liquid or gas region is not separated into smaller pieces of liquid
or gas at any time $t$. That is, the number of interfaces in the system is
always two. Second, the time evolution occurs along the metastable states by
the influence of weak noise. This means that the interface positions
correspond to the slowest mode of the time evolution and the densities of the
liquid and gas regions are determined from the condition of the metastable
profile. We then consider the time evolution $\rho(x,t)$ as follows. Let
$D^{\rm L}(t)$ and $D^{\rm G}(t)$ denote the liquid region and gas region at
any time $t$. Recalling the relation (2.20) with (2.19) and (2.21), we express
$\rho(x,t)$ in the highly probable path ensemble as

$\rho(x,t)=\rho^{\rm L}(t)\chi(x\in D^{\rm L}(t))+\rho^{\rm G}(t)\chi(x\in
D^{\rm G}(t)),$ (3.9)

where $\rho^{\rm L}(t)$ and $\rho^{\rm G}(t)$ are determined from

$\mu(\rho^{\rm L}(t))=\mu(\rho^{\rm G}(t))$ (3.10)

and

$\rho^{\rm L}(t)|D^{\rm L}(t)|+\rho^{\rm G}(t)|D^{\rm G}(t)|=\bar{\rho}.$
(3.11)

It should be noted that we adopt the singular continuum description introduced
in the previous section. Using this form of the time evolution, we estimate
(3.2).

Because $\rho(x,t)$ in (3.9) is described by $D^{\rm L}(t)$ and $D^{\rm
G}(t)$, we explicitly express them in terms of the interface positions
$X_{\uparrow}(t)$ and $X_{\downarrow}(t)$ at time $t$, where they satisfy

$\rho(X_{\uparrow}(t)+\epsilon,t)-\rho(X_{\uparrow}(t)-\epsilon,t)>0,\quad\rho(X_{\downarrow}(t)+\epsilon,t)-\rho(X_{\downarrow}(t)-\epsilon,t)<0$
(3.12)

for small positive $\epsilon$. See Fig. 6(a). Then, the liquid region $D^{\rm
L}(t)$ and gas region $D^{\rm G}(t)$ are expressed as $D^{\rm
L}(t)=[X_{\uparrow}(t),X_{\downarrow}(t)]$ and $D^{\rm G}(t)=[0,1]\backslash
D^{\rm L}(t)$ if $X_{\uparrow}(t)<X_{\downarrow}(t)$, or $D^{\rm
G}(t)=[X_{\downarrow}(t),X_{\uparrow}(t)]$ and $D^{\rm L}(t)=[0,1]\backslash
D^{\rm G}(t)$ if $X_{\uparrow}(t)>X_{\downarrow}(t)$.

Now, we consider the time evolution of interface positions, $X_{\uparrow}(t)$
and $X_{\downarrow}(t)$, starting from $X_{\uparrow}(0)=0$ and
$X_{\downarrow}(0)=X$. However, because $0\leq X_{\uparrow}(t)\leq 1$ and
$0\leq X_{\downarrow}(t)\leq 1$, $X_{\uparrow}(t)$ and $X_{\downarrow}(t)$ are
not continuous functions of $t$. This would lead to a complicated calculation
of the accumulated current defined by (3.4). To describe the interface motion
using continuous functions, we introduce generalized coordinates
${\hat{X}}_{\uparrow}(t)\in\mathbb{R}$ and
${\hat{X}}_{\downarrow}(t)\in\mathbb{R}$ such that displacements of the
interface from the initial time $0$ to the time $t$ are given by
${\hat{X}}_{\uparrow}(t)-{\hat{X}}_{\uparrow}(0)$ and
${\hat{X}}_{\downarrow}(t)-{\hat{X}}_{\downarrow}(0)$. That is,
${\hat{X}}_{\uparrow}(t)$ and ${\hat{X}}_{\downarrow}(t)$ describe the
positions of the left and right interfaces of the liquid region in a
generalized coordinate space $\mathbb{R}$. The interface positions
$X_{\uparrow}(t)$ and $X_{\downarrow}(t)$ in the space $[0,1]$ are then
obtained as

$\displaystyle
X_{\uparrow}(t)={\hat{X}}_{\uparrow}(t)-\lfloor{\hat{X}}_{\uparrow}(t)\rfloor,\qquad
X_{\downarrow}(t)={\hat{X}}_{\downarrow}(t)-\lfloor{\hat{X}}_{\downarrow}(t)\rfloor,$
(3.13)

where $\lfloor\ \rfloor$ represents the floor function. The width of the
liquid and gas regions are then written as

$\displaystyle|D^{\rm
L}(t)|={\hat{X}}_{\downarrow}(t)-{\hat{X}}_{\uparrow}(t),\quad|D^{\rm
G}(t)|=1-|D^{\rm L}(t)|={\hat{X}}_{\uparrow}(t)-{\hat{X}}_{\downarrow}(t)+1$
(3.14)

irrespective of the sign of $X_{\downarrow}(t)-X_{\uparrow}(t)$. We also
introduce the center of the liquid region in the generalized coordinate space
as

$\hat{Y}(t)=\frac{{\hat{X}}_{\uparrow}(t)+{\hat{X}}_{\downarrow}(t)}{2}.$
(3.15)

Using the width $|D^{\rm L}(t)|$ and the center $\hat{Y}(t)$ of the liquid
region, we have

$\displaystyle{\hat{X}}_{\downarrow}(t)=\hat{Y}(t)+\frac{|D^{\rm
L}(t)|}{2},\qquad{\hat{X}}_{\uparrow}(t)=\hat{Y}(t)-\frac{|D^{\rm L}(t)|}{2}.$
(3.16)

We note that, in the weak noise limit, $|D^{\rm L}(t)|\to X^{\rm eq}$ as
$t\to\infty$, while $\hat{Y}(t)$ shows unbounded-free Brownian motion because
of the translation symmetry for the case $\phi=0$. This fact simplifies the
calculation of the accumulated current defined by (3.2).

At the end of this subsection, we discuss the difference between the most
probable process for the case $\kappa_{\Lambda}\gg 1$ and the highly probable
path ensemble for the case $\kappa_{\Lambda}\ll 1$. In the former case,
$\rho_{X}^{\phi=0}$ evolves to an equilibrium configuration in the
deterministic system, which is in contrast with the latter case where we
study. To obtain the accumulated current $Q$ for the former case, we have to
analyze the time-dependent solution of the deterministic equation, which is
out of the present paper.

###  3.3 Expression of $j(0,t)$

Figure 6:  (a) Density profile $\rho(x,t)$ whose interface positions are
$X_{\uparrow}(t)$ and $X_{\downarrow}(t)$. $\rho(x,t)=\rho^{\rm L}(t)$ in
$D^{\rm L}(t)=[X_{\uparrow}(t),X_{\downarrow}(t)]$, while $\rho(x,t)=\rho^{\rm
G}(t)$ in $D^{\rm G}(t)=[0,1]\backslash D^{\rm L}(t)$. (b) Example of
$\psi(x,t)$ defined in $-\infty<x<\infty$.

In this subsection, we calculate $j(x,t)$ in the singular continuum
description based on (2.36) and (2.37). In the argument below,
${\hat{X}}_{\uparrow}(t)$, ${\hat{X}}_{\downarrow}(t)$, $D^{\rm L}(t)$, and
$D^{\rm G}(t)$ are simply denoted by ${\hat{X}}_{\uparrow}$,
${\hat{X}}_{\downarrow}$, $D^{\rm L}$, and $D^{\rm G}$ if their
$t$-dependencies are clearly guessed.

For a given density profile $(\rho(x,t))_{x=0}^{1}$ at time $t$, we define

$\psi(x,t)\equiv\int_{0}^{x}dx^{\prime}(\rho(x^{\prime},t)-\bar{\rho})$ (3.17)

for $0\leq x\leq 1$. Because $\psi(0,t)=\psi(1,t)=0$, $\psi(x,t)$ can be
extended to a periodic function in $x$. That is, we define
$\psi(x,t)\equiv\psi(x-\lfloor x\rfloor,t)$ for any $-\infty<x<\infty$. See
Fig. 6(b) for the illustration. The time derivative of (3.17) leads to

$\partial_{t}\psi(x,t)=-j(x,t)+j(0,t)$ (3.18)

for any $x\in[0,1]$. We here integrate (3.18) over the liquid region and
divide by $\sigma^{\rm L}$. We repeat the same operation for the gas region.
Summing up the two results, we have a relation

$\frac{1}{\sigma^{\rm L}}\int_{D^{\rm
L}}dx\partial_{t}\psi+\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}}dx\partial_{t}\psi=-\frac{1}{\sigma^{\rm L}}\int_{D^{\rm
L}}dxj(x,t)-\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}}dxj(x,t)+\left(\frac{|D^{\rm L}|}{\sigma^{\rm L}}+\frac{|D^{\rm
G}|}{\sigma^{\rm G}}\right)~{}j(0,t).$ (3.19)

Recall the formula for $j(x,t)$ given in (2.37). Because $\int
dx~{}\xi(x,t)=0$ and $\mu(x,t)$ takes a constant value in the bulk regions,
the first and second terms on the right side of (3.19) turn out to be zero. We
thus obtain the expression for $j(0,t)$ as

$\displaystyle j(0,t)=\left[\frac{|D^{\rm L}|}{\sigma^{\rm L}}+\frac{|D^{\rm
G}|}{\sigma^{\rm G}}\right]^{-1}\left[\frac{1}{\sigma^{\rm L}}\int_{D^{\rm
L}}dx~{}\partial_{t}\psi+\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}}dx~{}\partial_{t}\psi\right].$ (3.20)

We remark that $j(0,t)$ is stochastic as specified by (2.37). In the right
side of (3.20), $D^{\rm L}(t)$, $D^{\rm G}(t)$, and the space integrals of
$\partial_{t}\psi(x,t)$ are affected by the Brownian motion of the interfaces.

The current $j(0,t)$ formulated in (3.20) contains the time derivative of
$\psi$ in the space integrals. We now transform (3.20) into a formula by
letting the time derivative be outside the integral. In the transformation
procedure, we need to pay attention to the time-dependent ranges $D^{\rm L}$
and $D^{\rm G}$ of the integrals. As the result, the current $j(0,t)$ is
expressed as

$j(0,t)=\Phi(t)-\Phi_{0}(t),$ (3.21)

where

$\displaystyle\Phi(t)\equiv\frac{d}{dt}\left\\{\left[\frac{|D^{\rm
L}|}{\sigma^{\rm L}}+\frac{|D^{\rm G}|}{\sigma^{\rm
G}}\right]^{-1}\left[\frac{1}{\sigma^{\rm L}}\int_{D^{\rm
L}}dx~{}\psi+\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}}dx~{}\psi\right]\right\\},$ (3.22)

and $\Phi_{0}(t)$ is determined from (3.20), (3.21), and (3.22). Concretely,
we perform the time-derivative of (3.22). In the time derivative of a function
of Brownian motion $X(t)$, we note the chain rule
$df(X(t))/dt=f^{\prime}(X(t))\circ dX/dt$, where the symbol $\circ$ represents
the Stratonovich product. We then have

$\Phi_{0}(t)=\left(\frac{1}{\sigma^{\rm L}}-\frac{1}{\sigma^{\rm
G}}\right)\left[\frac{|D^{\rm L}(t)|}{\sigma^{\rm L}}+\frac{|D^{\rm
G}(t)|}{\sigma^{\rm G}}\right]^{-1}|D^{\rm L}(t)|(\rho^{\rm
L}(t)-\bar{\rho})\circ\frac{d}{dt}\hat{Y}.$ (3.23)

See Appendix D for the derivation of $\Phi_{0}(t)$.

Due to the translational invariance for the case $\phi=0$, $\hat{Y}(t)$ shows
the free Brownian motion which is independent of $|D^{\rm L}(t)|$, $|D^{\rm
G}(t)|$, and $\rho^{\rm L}(t)$. Therefore, taking the path ensemble average
over noise realization, we have

$\displaystyle\left\langle\Phi_{0}(t)\right\rangle$
$\displaystyle=\left(\frac{1}{\sigma^{\rm L}}-\frac{1}{\sigma^{\rm
G}}\right)\left\langle{\left[\frac{|D^{\rm L}(t)|}{\sigma^{\rm
L}}+\frac{|D^{\rm G}(t)|}{\sigma^{\rm G}}\right]^{-1}|D^{\rm L}(t)|(\rho^{\rm
L}(t)-\bar{\rho})}\right\rangle\left\langle\frac{d}{dt}\hat{Y}\right\rangle$
$\displaystyle=0$ (3.24)

Combining this with (3.21), we obtain

$\left\langle j(0,t)\right\rangle^{\rm
eq}_{{\bm{\rho}}^{\phi=0}_{X}}=\left\langle\Phi(t)\right\rangle^{\rm
eq}_{{\bm{\rho}}^{\phi=0}_{X}}.$ (3.25)

###  3.4 Result of $\left\langle Q\right\rangle^{{\rm
eq}}_{{\bm{\rho}}_{X}^{\phi=0}}$

Let us define

$q(\tau)\equiv\int_{0}^{\tau}dt~{}\Phi(t).$ (3.26)

From (3.25), we have

$\left\langle Q\right\rangle^{{\rm
eq}}_{\rho_{X}^{\phi=0}}=\lim_{\tau\to\infty}\left\langle
q(\tau)\right\rangle^{{\rm eq}}_{\rho_{X}^{\phi=0}}.$ (3.27)

Substituting (3.22) into (3.26) with noting $X_{\uparrow}(0)=0$ and
$X_{\downarrow}(0)=X$, we obtain

$\displaystyle q(\tau)=$ $\displaystyle\left[\frac{|D^{\rm
L}(\tau)|}{\sigma^{\rm L}}+\frac{|D^{\rm G}(\tau)|}{\sigma^{\rm
G}}\right]^{-1}\left[\frac{1}{\sigma^{\rm L}}\int_{D^{\rm
L}(\tau)}dx~{}\psi+\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}(\tau)}dx~{}\psi\right]$ $\displaystyle\quad-\left[\frac{X}{\sigma^{\rm
L}}+\frac{1-X}{\sigma^{\rm G}}\right]^{-1}\left[\frac{1}{\sigma^{\rm
L}}\int_{0}^{X}dx~{}\psi+\frac{1}{\sigma^{\rm
G}}\int_{X}^{1}dx~{}\psi\right].$ (3.28)

From the piece-wise linear nature of $\psi(x,t)$, the integrals are calculated
by the trapezoidal rule as

$\displaystyle\int_{D^{\rm L}}dx~{}\psi(x,t)=\frac{1}{2}|D^{\rm
L}|\left(\psi({\hat{X}}_{\uparrow}(t),t)+\psi({\hat{X}}_{\downarrow}(t),t)\right),$
(3.29) $\displaystyle\int_{D^{\rm G}}dx~{}\psi(x,t)=\frac{1}{2}|D^{\rm
G}|\left(\psi({\hat{X}}_{\uparrow}(t),t)+\psi({\hat{X}}_{\downarrow}(t),t)\right),$
(3.30)

where $\psi(X_{\uparrow}(t)+1)=\psi(X_{\uparrow}(t))$ is applied. Using these
relations, we obtain

$\displaystyle\frac{1}{\sigma^{\rm L}}\int_{D^{\rm
L}}dx~{}\psi+\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}}dx~{}\psi=\frac{1}{2}\left(\frac{|D^{\rm L}|}{\sigma^{\rm L}}+\frac{|D^{\rm
G}|}{\sigma^{\rm
G}}\right)\left(\psi({\hat{X}}_{\uparrow}(t),t)+\psi({\hat{X}}_{\downarrow}(t),t)\right).$
(3.31)

Substituting this into (3.28), we have

$\displaystyle q(\tau)=$
$\displaystyle\frac{1}{2}\left(\psi({\hat{X}}_{\uparrow}(\tau),\tau)+\psi({\hat{X}}_{\downarrow}(\tau),\tau)\right)-\frac{1}{2}\left(\psi(0,0)+\psi(X,0)\right).$
(3.32)

Here,

$\displaystyle\psi(0,0)=0,\qquad\psi(X,0)=({\rho}^{\rm L}_{X}-\bar{\rho})X,$
(3.33)

and

$\displaystyle\psi({\hat{X}}_{\uparrow}(\tau),\tau)=\psi(X_{\uparrow}(\tau),\tau),\qquad\psi({\hat{X}}_{\downarrow}(\tau),\tau)=\psi(X_{\downarrow}(\tau),\tau).$
(3.34)

We thus obtain

$\lim_{\tau\to\infty}\left\langle q(\tau)\right\rangle^{{\rm
eq}}_{{\bm{\rho}}_{X}^{\phi=0}\to*}=\frac{1}{2}\lim_{\tau\to\infty}\left(\left\langle{\psi(X_{\uparrow}(\tau),\tau)}\right\rangle+\left\langle{\psi(X_{\downarrow}(\tau),\tau)}\right\rangle\right)-\frac{1}{2}({\rho}^{\rm
L}_{X}-\bar{\rho})X.$ (3.35)

Note that the positions of the interfaces, $X_{\uparrow}$ and
$X_{\downarrow}$, are uniformly distributed in the interval $x\in[0,1]$ as
$\tau\to\infty$. We thus conclude that

$\left\langle Q\right\rangle^{{\rm
eq}}_{{\bm{\rho}}_{X}^{\phi=0}}=C-\frac{1}{2}({\rho}^{\rm
L}_{X}-\bar{\rho})X,$ (3.36)

where $C$ is a constant independent of $X$. Using the relation

$\displaystyle\bar{\rho}=X{\rho}^{\rm L}_{X}+(1-X){\rho}^{\rm G}_{X},$ (3.37)

which comes from (2.44), we can also express (3.36) as

$\left\langle Q\right\rangle^{{{\rm
eq}}}_{{\bm{\rho}}_{X}^{\phi=0}}=C-\frac{1}{2}({\rho}^{\rm L}_{X}-{\rho}^{\rm
G}_{X})X(1-X).$ (3.38)

Substituting (3.38) into (3.5), we have

${\cal F}_{\rm
ss}(X;\bar{\rho},\phi)=\int_{0}^{1}dx~{}f(\rho_{X}^{\phi}(x))-\frac{\phi}{2}({\rho}^{\rm
L}_{X}-{\rho}^{\rm G}_{X})X(1-X)+\phi C.$ (3.39)

We note that $\rho_{X}^{\phi}(x)$ is uniquely determined for given
$(X,\phi,\bar{\rho})$ by $\mu(\rho_{X}^{\phi}(x))=\mu_{X}^{\phi}(x)$ with
(2.40) and (2.41) and that ${\rho}^{\rm L}_{X}$ and ${\rho}^{\rm G}_{X}$ are
functions of $(X,\bar{\rho})$. Thus, ${\cal F}_{\rm ss}(X;\bar{\rho},\phi)$ is
the variational function for determining the most probable interface position
$X_{*}$ for a given $(\bar{\rho},\phi)$. The most probable density profile is
expressed as $\rho_{X_{*}}^{\phi}(x)$.

##  4 Variational equation

In this section, using the variational function (3.39), we determine the most
probable value of $X$. This is regarded as an extension of the argument for
determining the equilibrium profile in the paragraph containing (2.24). Thus,
similarly to (2.24), we analyze the variational equation

$\left.\frac{d{\cal F}_{\rm ss}(X;\bar{\rho},\phi)}{dX}\right|_{X=X_{*}}=0$
(4.1)

with (3.39). The difference from (2.24) is the $O(\phi)$ terms in (3.39). As
shown in Fig. 4, one of the interfaces of $\rho_{X}^{\phi}(x)$ is located at
$x=0$ for any $X$, both the liquid and gas density profiles are sloped, and
the chemical potential profile $\mu_{X}^{\phi}(x)$ is piece-wise linear. Then,
as did for the equilibrium system, we first derive the chemical potential at
the interface $\mu_{X_{*}}^{\rm I}$ defined by (2.39) instead of directly
calculating $X_{*}$. Once we have $\mu_{X_{*}}^{\rm I}$, we obtain the most
probable profiles $\mu_{X_{*}}^{\phi}(x)$ and $\rho_{X_{*}}^{\phi}(x)$, and
the value of $X_{*}$. The determination of $\mu_{X_{*}}^{\rm I}$ is also
physically important because if $\mu_{X_{*}}^{\rm I}\neq\mu_{\rm c}$,
metastable states at equilibrium stably appear in the non-equilibrium phase
coexistence.

###  4.1 Preliminaries for the calculation

In the argument below, we ignore the contribution of $O(\phi^{2})$. We first
define

${\bar{\rho}}^{\rm
L}_{X}\equiv\frac{1}{X}\int_{0}^{X}dx\rho^{\phi}_{X}(x),\quad{\bar{\rho}}^{\rm
G}_{X}\equiv\frac{1}{1-X}\int_{X}^{1}dx\rho^{\phi}_{X}(x),$ (4.2)

and

${\bar{\mu}}^{\rm
L}_{X}\equiv\frac{1}{X}\int_{0}^{X}dx\mu^{\phi}_{X}(x),\quad{\bar{\mu}}^{\rm
G}_{X}\equiv\frac{1}{1-X}\int_{X}^{1}dx\mu^{\phi}_{X}(x).$ (4.3)

Because the density profile $\rho_{X}^{\phi}(x)$ and the chemical potential
profile $\mu_{X}^{\phi}(x)$ are linear in the respective regions, $[0,X]$ and
$[X,1]$, we obtain

$f\left({\bar{\rho}}^{\rm
L}_{X}\right)=\frac{1}{X}\int_{0}^{X}dxf(\rho^{\phi}_{X}(x)),\qquad
f\left({\bar{\rho}}^{\rm
G}_{X}\right)=\frac{1}{1-X}\int_{X}^{1}dxf(\rho^{\phi}_{X}(x)),$ (4.4)

and

${\bar{\mu}}^{\rm L}_{X}=\mu({\bar{\rho}}^{\rm L}_{X}),\qquad{\bar{\mu}}^{\rm
G}_{X}=\mu({\bar{\rho}}^{\rm G}_{X}).$ (4.5)

The pressures in the liquid and gas regions are characterized by $p^{\rm
L}_{X}\equiv p({\bar{\rho}}^{\rm L}_{X})$ and $p^{\rm G}_{X}\equiv
p({\bar{\rho}}^{\rm G}_{X})$, which are expressed by

$\displaystyle p^{\rm L}_{X}=\mu^{\rm L}_{X}{\bar{\rho}}^{\rm
L}_{X}-f({\bar{\rho}}^{\rm L}_{X}),\qquad p^{\rm G}_{X}=\mu^{\rm
G}_{X}{\bar{\rho}}^{\rm G}_{X}-f({\bar{\rho}}^{\rm G}_{X}).$ (4.6)

The first term on the right side of (3.39) is written as

$\int_{0}^{1}dx~{}f(\rho_{X}^{\phi}(x))=Xf\left({\bar{\rho}}^{\rm
L}_{X}\right)+(1-X)f\left({\bar{\rho}}^{\rm G}_{X}\right)+O(\phi^{2}),$ (4.7)

and the variational function (3.5) as

$\displaystyle{{\cal F}_{\rm ss}}(X;\bar{\rho},\phi)=Xf\left({\bar{\rho}}^{\rm
L}_{X}\right)+(1-X)f\left({\bar{\rho}}^{\rm
G}_{X}\right)-\frac{\phi}{2}\left({\rho}^{\rm L}_{X}-{\rho}^{\rm
G}_{X}\right)X(1-X)+\phi C.$ (4.8)

Since ${\bar{\rho}}^{\rm L}_{X}={\rho}^{\rm L}_{X}+O(\phi)$ and
${\bar{\rho}}^{\rm G}_{X}={\rho}^{\rm G}_{X}+O(\phi)$, we can rewrite (4.8) as

$\displaystyle{{\cal F}_{\rm ss}}(X;\bar{\rho},\phi)=Xf\left({\bar{\rho}}^{\rm
L}_{X}\right)+(1-X)f\left({\bar{\rho}}^{\rm
G}_{X}\right)-\frac{\phi}{2}\left({\bar{\rho}}^{\rm L}_{X}-{\bar{\rho}}^{\rm
G}_{X}\right)X(1-X)+\phi C.$ (4.9)

We emphasize that (4.9) is explicitly expressed as a function of $X$.

Noting that the chemical potential profile is piece-wise linear as shown in
Fig. 4, we estimate

$\displaystyle{\bar{\mu}}^{\rm L}_{X}=\frac{\mu_{X}^{\phi}(0)+\mu_{X}^{\rm
I}}{2},\qquad{\bar{\mu}}^{\rm G}_{X}=\frac{\mu_{X}^{\phi}(1)+\mu_{X}^{\rm
I}}{2}.$ (4.10)

From (4.10), $\phi=\mu_{X}^{\phi}(1)-\mu_{X}^{\phi}(0)$ is rewritten as

${\bar{\mu}}^{\rm L}_{X}-{\bar{\mu}}^{\rm G}_{X}=\frac{\phi}{2}.$ (4.11)

Furthermore, using (2.40) and (2.41), we can express ${\bar{\mu}}^{\rm L}_{X}$
and ${\bar{\mu}}^{\rm G}_{X}$ in terms of $\mu_{X}^{\rm I}$ as

$\displaystyle{\bar{\mu}}^{\rm L}_{X}=\mu_{X}^{\rm I}+\frac{JL}{\sigma^{\rm
L}}\frac{X}{2},\quad{\bar{\mu}}^{\rm G}_{X}=\mu_{X}^{\rm
I}-\frac{JL}{\sigma^{\rm G}}\frac{1-X}{2}.$ (4.12)

###  4.2 Steady state

We consider the variational equation

$\frac{d{\cal F}_{\rm ss}(X;\bar{\rho},\phi)}{dX}=0.$ (4.13)

Substituting (4.9) into (4.13), we have

$\displaystyle f({\bar{\rho}}^{\rm L}_{X})-f({\bar{\rho}}^{\rm
G}_{X})-\frac{\phi}{2}({\bar{\rho}}^{\rm L}_{X}-{\bar{\rho}}^{\rm
G}_{X})(1-2X)$ $\displaystyle+$ $\displaystyle{\bar{\mu}}^{\rm
L}_{X}\frac{d{\bar{\rho}}^{\rm L}_{X}}{dX}X+{\bar{\mu}}^{\rm
G}_{X}\frac{d{\bar{\rho}}^{\rm
G}_{X}}{dX}(1-X)-\frac{\phi}{2}\left(\frac{d{\bar{\rho}}^{\rm
L}_{X}}{dX}-\frac{d{\bar{\rho}}^{\rm G}_{X}}{dX}\right)X(1-X)$
$\displaystyle=$ $\displaystyle 0.$ (4.14)

Using (4.11), we rewrite the second line as

$\displaystyle[X{\bar{\mu}}^{\rm L}_{X}+(1-X){\bar{\mu}}^{\rm
G}_{X}]\left[X\frac{d{\bar{\rho}}^{\rm
L}_{X}}{dX}+(1-X)\frac{d{\bar{\rho}}^{\rm G}_{X}}{dX}\right].$ (4.15)

Here, taking the derivative of $\bar{\rho}=X{\bar{\rho}}^{\rm
L}_{X}+(1-X){\bar{\rho}}^{\rm G}_{X}$ in $X$, we obtain

$\displaystyle X\frac{d{\bar{\rho}}^{\rm
L}_{X}}{dX}+(1-X)\frac{d{\bar{\rho}}^{\rm G}_{X}}{dx}=-({\bar{\rho}}^{\rm
L}_{X}-{\bar{\rho}}^{\rm G}_{X}).$ (4.16)

We substitute this into (4.15) and combine it with (4.14). Then, the
variational equation (4.14) becomes

$\displaystyle f({\bar{\rho}}^{\rm L}_{X})-f({\bar{\rho}}^{\rm
G}_{X})-({\bar{\rho}}^{\rm L}_{X}-{\bar{\rho}}^{\rm
G}_{X})[(1-X){\bar{\mu}}^{\rm L}_{X}+X{\bar{\mu}}^{\rm G}_{X}]=0.$ (4.17)

The solution of (4.17) provides the most probable value $X_{*}$ of the
interface position. Therefore, we express

$\displaystyle f({\bar{\rho}}^{\rm L}_{X_{*}})-f({\bar{\rho}}^{\rm
G}_{X_{*}})-({\bar{\rho}}^{\rm L}_{X_{*}}-{\bar{\rho}}^{\rm
G}_{X_{*}})[(1-X_{*}){\bar{\mu}}^{\rm L}_{X_{*}}+X{\bar{\mu}}^{\rm
G}_{X_{*}}]=0.$ (4.18)

Now, we rewrite (4.18) as a different form using the chemical potential at the
interface position $X_{*}$. We subtract the equilibrium version of (4.18),
i.e., $f(\rho^{\rm L}_{\rm c})-f(\rho^{\rm G}_{\rm c})-(\rho^{\rm L}_{\rm
c}-\rho^{\rm G}_{\rm c})\mu_{\rm c}=0$, from (4.18). Noting
$f({\bar{\rho}}^{\rm L}_{X_{*}})-f(\rho^{\rm L}_{\rm c})=({\bar{\rho}}^{\rm
L}_{X_{*}}-\rho^{\rm L}_{\rm c})\mu_{\rm c}+O(\phi^{2})$, we obtain

$\displaystyle({\bar{\rho}}^{\rm L}_{X_{*}}-{\bar{\rho}}^{\rm
G}_{X_{*}})\left\\{\mu_{\rm c}-[(1-X_{*}){\bar{\mu}}^{\rm
L}_{X_{*}}+X_{*}{\bar{\mu}}^{\rm G}_{X_{*}}]\right\\}=0.$ (4.19)

Substituting (4.12) into this equation and using ${\bar{\rho}}^{\rm
L}_{X_{*}}\neq{\bar{\rho}}^{\rm G}_{X_{*}}$, we have

$\mu_{X_{*}}^{\rm I}=\mu_{\rm
c}-\frac{JLX_{*}(1-X_{*})}{2}\left(\frac{1}{\sigma^{\rm
L}}-\frac{1}{\sigma^{\rm G}}\right).$ (4.20)

Because $X_{*}=X^{\rm eq}+O(\phi)$, we can rewrite it as

$\mu_{X_{*}}^{\rm I}=\mu_{\rm c}-\frac{JLX^{\rm eq}(1-X^{\rm
eq})}{2}\left(\frac{1}{\sigma^{\rm L}}-\frac{1}{\sigma^{\rm G}}\right).$
(4.21)

Furthermore, combining (2.43) with (4.21), we finally obtain

$\mu_{X_{*}}^{\rm I}=\mu_{\rm c}+\frac{\phi}{2}\frac{(\sigma^{\rm
L}-\sigma^{\rm G})X^{\rm eq}(1-X^{\rm eq})}{\sigma^{\rm G}X^{\rm
eq}+\sigma^{\rm L}(1-X^{\rm eq})}.$ (4.22)

Recalling that $X^{\rm eq}$ is uniquely determined by $\bar{\rho}$ from
(2.18), we conclude that $\mu_{X_{*}}^{\rm I}$ is expressed in terms of the
system parameters. Thus, the chemical potential at the interface deviates from
$\mu_{\rm c}$ linearly with the voltage $\phi$. This means that metastable
states at equilibrium stably appear around the interface. The relation (4.22)
is the most important achievement of our theory.

Next, from (4.11) and (4.19), we obtain

$\displaystyle{\bar{\mu}}^{\rm L}_{X_{*}}=\mu_{\rm c}+X^{\rm
eq}\frac{\phi}{2},\quad{\bar{\mu}}^{\rm G}_{X_{*}}=\mu_{\rm c}-(1-X^{\rm
eq})\frac{\phi}{2},$ (4.23)

which also gives $\mu_{X_{*}}^{\phi}(0)$ and $\mu_{X_{*}}^{\phi}(1)$ using
(4.10). The result yields $\mu_{X_{*}}^{\phi}(x)$. Using
$\mu_{X_{*}}^{\phi}(x)=\mu(\rho_{X_{*}}^{\phi}(x))$, we have
$\rho_{X_{*}}^{\phi}(x)$. Finally, from $X_{*}{\bar{\rho}}^{\rm
L}_{X_{*}}+(1-X_{*}){\bar{\rho}}^{\rm G}_{X_{*}}=\bar{\rho}$, we can express
$X_{*}-X^{\rm eq}$ in terms of system parameters. In this manner, all
thermodynamic quantities are determined by the variational principle. As one
example, we discuss the pressure in the steady state. Using (4.6), we can
express (4.18) as

$\displaystyle p^{\rm L}_{X_{*}}-p^{\rm G}_{X_{*}}=\frac{\phi\bar{\rho}}{2}.$
(4.24)

Recalling that the local pressure is given by $p(x)=p(\rho(x))$, we define the
pressures at the left and right sides of the interface as

$\displaystyle p_{-}\equiv\lim_{x\to X_{*}^{-}}p(\rho(x)),\qquad
p_{+}\equiv\lim_{x\to X_{*}^{+}}p(\rho(x)).$ (4.25)

Then, using

$\rho^{\rm L/G}_{X_{*}}(\mu^{\rm L/G}_{X_{*}}-\mu^{\rm I}_{X_{*}})=p^{\rm
L/G}_{X_{*}}-p_{-/+},$ (4.26)

we can derive

$\displaystyle p_{-}=p_{\rm c}+(\mu_{X_{*}}^{\rm I}-\mu_{\rm c})\rho^{\rm
L}_{\rm c},\qquad p_{+}=p_{\rm c}+(\mu_{X_{*}}^{\rm I}-\mu_{\rm c})\rho^{\rm
G}_{\rm c},$ (4.27)

where $p_{\rm c}$ is the equilibrium coexistence pressure. This result
indicates that the pressure is not continuous at the interface. It should be
noted that this discontinuity occurs only in the singular continuum
description for $\kappa_{\Lambda}\ll 1$ but never occur in the continuum
description for $\kappa_{\Lambda}\gg 1$ as shown in Appendix B.

##  5 Global thermodynamics

In this section, we introduce an extended framework of thermodynamics based on
the variational function (4.9). We set $M\equiv\bar{\rho}L$, where $MA$
represents the total number of particles in the tube. In equilibrium
thermodynamics, the free energy function $F_{\rm eq}(L,M)$ is determined as

$F_{\rm eq}(L,M)={\cal F}_{\rm eq}(X_{*};\bar{\rho})$ (5.1)

using the variational function ${\cal F}_{\rm eq}(X;\bar{\rho})$ given by
(2.22). We then have the fundamental relation of thermodynamics

$dF_{\rm eq}=-pdL+\mu dM,$ (5.2)

where $p$ and $\mu$ are the equilibrium values of pressure and chemical
potential for a given $(L,M)$. We now attempt to extend the relations (5.1)
and (5.2) to non-equilibrium systems.

First, following the relation (5.1), we define the thermodynamic function
$F_{\rm ss}(L,M,\phi)$ as

$F_{\rm ss}(L,M,\phi)={\cal F}_{\rm ss}(X_{*};\bar{\rho},\phi).$ (5.3)

Using the expression of the variational function (4.9), we explicitly express
(5.3) as

$\displaystyle F_{{\rm ss}}(L,M,\phi)$
$\displaystyle=LX_{*}f({\bar{\rho}}^{\rm
L}_{X_{*}})+L(1-X_{*})f({\bar{\rho}}^{\rm
G}_{X_{*}})-\frac{\phi}{2}L({\bar{\rho}}^{\rm L}_{X_{*}}-{\bar{\rho}}^{\rm
G}_{X_{*}})X_{*}(1-X_{*}).$ (5.4)

The question here is to derive an extended form of (5.2) for $F_{{\rm
ss}}(L,M,\phi)$.

Then, using (4.6) and (4.11), we rewrite (5.4) as

$\displaystyle F_{{\rm ss}}(L,M,\phi)$ $\displaystyle=-LX_{*}p^{\rm
L}_{X_{*}}-L(1-X_{*})p^{\rm G}_{X_{*}}+(X_{*}{\bar{\mu}}^{\rm
L}_{X_{*}}+(1-X){\bar{\mu}}^{\rm G}_{X_{*}})M.$ (5.5)

Defining the average pressure and chemical potential as

$\displaystyle\bar{p}=X_{*}p^{\rm L}_{X_{*}}+(1-X_{*})p^{\rm
G}_{X_{*}},\qquad\bar{\mu}=X_{*}{\bar{\mu}}^{\rm
L}_{X_{*}}+(1-X_{*}){\bar{\mu}}^{\rm G}_{X_{*}},$ (5.6)

we further rewrite $F_{\rm ss}$ as a suggestive form

$\displaystyle F_{{\rm ss}}(L,M,\phi)$ $\displaystyle=-\bar{p}L+\bar{\mu}M,$
(5.7)

where $\bar{p}$ and $\bar{\mu}$ are functions of $(L,M,\phi)$.

Now, taking the derivative of $F_{\rm ss}$ in $L$, we have

$\frac{\partial F_{\rm ss}(L,M,\phi)}{\partial
L}=-\bar{p}-\frac{\partial\bar{p}(L,M,\phi)}{\partial
L}L+\frac{\partial\bar{\mu}(L,M,\phi)}{\partial L}M.$ (5.8)

Here, using (4.5), (4.6), (4.11), and (4.24), we can confirm

$\frac{\partial\bar{p}(L,M,\phi)}{\partial
L}=\rho\frac{\partial\bar{\mu}(L,M,\phi)}{\partial L}.$ (5.9)

Substituting this result into (5.8), we obtain

$\frac{\partial F_{\rm ss}(L,M,\phi)}{\partial L}=-\bar{p}.$ (5.10)

By repeating the similar calculation, we also have

$\frac{\partial F_{\rm ss}(L,M,\phi)}{\partial M}=\bar{\mu}.$ (5.11)

Finally, we define a new variable $\Psi$ by

$\Psi\equiv-\frac{\partial F_{\rm ss}(L,M,\phi)}{\partial\phi}.$ (5.12)

From (5.4), we have

$\Psi(L,M,\phi)=\frac{1}{2}L({\bar{\rho}}^{\rm L}_{X_{*}}-{\bar{\rho}}^{\rm
G}_{X_{*}})X_{*}(1-X_{*}).$ (5.13)

These three relations (5.10), (5.11), and (5.12) are summarized as an extended
form of the fundamental relation of thermodynamics

$\displaystyle dF_{\rm ss}=-\bar{p}dL+\bar{\mu}dM-\Psi d\phi.$ (5.14)

Therefore, $F_{{\rm ss}}(L,M,\phi)$ is regarded as a natural extension of the
equilibrium free energy $F_{\rm eq}(L,M)$ defined in (5.1).

As seen from this analysis, the extension of the variational function from
${\cal F}_{\rm eq}(X;\bar{\rho})$ to ${\cal F}_{\rm ss}(X;\bar{\rho},\phi)$ is
closely related to the extension of the thermodynamic function from $F_{\rm
eq}(L,M)$ to $F_{\rm ss}(L,M,\phi)$. Without analyzing specific stochastic
models, one can construct such an extended thermodynamic framework relying on
the consistency, uniqueness, and predictability. This phenomenological
argument, which is called global thermodynamics, was developed for heat
conduction systems exhibiting phase coexistence Global-PRL ; Global-JSP ;
Global-PRR . Furthermore, global thermodynamics was applied to the order-
disorder transition in heat conduction, and the prediction by global
thermodynamics was confirmed by numerical simulations KNS . Similarly, in the
present setup, starting from (5.3) and (5.14) without the explicit form of
(5.13), we can determine the variational function (4.9). See Global-PRR  for
the argument. Therefore, all the prediction made by global thermodynamics for
the present setup are the same as the theoretical result for the stochastic
model we study. In our research history, the variational function (4.9) was
first derived using global thermodynamics, and after that it was re-derived by
analyzing the stochastic model.

##  6 Numerical simulation

In this section, we perform numerical simulations of the discrete model and
compare numerical results with the theoretical predictions presented in the
previous sections. More explicitly, the time evolution of
$(\rho_{i})_{i=1}^{N}$ is defined by (2.8) accompanied with the current
$(j_{i})_{i=1}^{N}$ defined by (2.10). To obtain $(j_{i})_{i=1}^{N}$ using
(2.10), we determine $(\tilde{\mu}_{i})_{i=1}^{N}$ by (2.9) with

$\displaystyle\mu(\rho_{i})=(\rho_{i}-0.5)(\rho_{i}-1.5)(\rho_{i}-2.5)$ (6.1)

from (2.6). We adopt a simple form of the conductivity $\sigma(\rho)=\rho$,
where we have introduced a dimensionless time in this expression. For this
specific model, we have $\rho^{\rm L}_{\rm c}=2.5$, $\rho^{\rm G}_{\rm
c}=0.5$, and $\mu_{\rm c}=0$ from (2.16), and we thus obtain $\sigma^{\rm
L}=2.5$ and $\sigma^{\rm G}=0.5$ from (2.30).

Recalling that the independent parameters to be specified for numerical
determination of the steady state are $(\kappa_{\Lambda},T_{\rm
eff},\phi,\bar{\rho},N)$ as discussed around (2.13), we study the
$\kappa_{\Lambda}$ dependence of the steady state with fixing the other
parameter values as

$(T_{\rm eff},\phi,\bar{\rho},N)=(0.002,0.05,1.5,64).$ (6.2)

Since $A$ and $\Lambda$ are contained in the renormalized quantities
$\kappa_{\Lambda}$ and $T_{\rm eff}$, we do not need to specify the values of
$A$ and $\Lambda$, while the total volume and the total number of particles
are given by $LA$ and $\bar{\rho}LA$ with $L=\Lambda N$. To numerically solve
(2.8) and (2.10), we adopt the Heun method with a time step
$dt=0.01\Lambda^{2}$, where it should be noted that the time step $dt$ is
always coupled with $\Lambda^{2}$ for the time-discretized form of (2.8) and
(2.10).

In Fig. 7, we show the steady state for $\kappa_{\Lambda}=0.5$ and
$\kappa_{\Lambda}=1.5$, where the density profile $\rho_{i}$ and the chemical
potential profile $\tilde{\mu}_{i}$ are plotted for $i/N$. We remark that the
system reaches the steady state without dependence on initial conditions,
$(\rho_{i})_{i=1}^{N}$ at $t=0$. Note that $\tilde{\mu}_{i}=\mu_{\rm c}=0$ for
all $i$ for the equilibrium system with $\phi=0$, which is shown as the dotted
line in Fig. 7(b). It is observed that $\tilde{\mu}$ near the interface is
close to $\mu_{\rm c}=0$ for the case $\kappa_{\Lambda}=1.5$, while it clearly
deviates from $\mu_{\rm c}=0$ for the case $\kappa_{\Lambda}=0.5$. That is,
the metastable gas stably appears at the left-side of the interface for the
case $\kappa_{\Lambda}=0.5$.

Here, we determine the chemical potential at the interface, which is denoted
by $\mu^{I}$, more quantitatively from the numerical data
$(\rho_{i})_{i=1}^{N}$ and $(\tilde{\mu}_{i})_{i=1}^{N}$. In principle, we
first determine the interface position $X^{I}$ from the data
$(\rho_{i})_{i=1}^{N}$, and then read the value of the chemical potential at
the interface position from the data $(\tilde{\mu}_{i})_{i=1}^{N}$. In
practice, we use the linear interpolation of the data sets to systematically
estimate $\mu^{I}$ for several parameters. That is, we define a piece-wise
linear function $\rho(x)$ for $0\leq x\leq 1$ by connecting two consecutive
points $(i/N,\rho_{i})$ and $((i+1)/N,\rho_{i+1})$ for $0\leq i\leq N$ in the
graph of $(i/N,\rho_{i})_{i=1}^{N}$. We then define the interface position
$X^{I}$ as $\rho(X^{I})=1.5$, where $(\rho^{\rm G}_{\rm c}+\rho^{\rm L}_{\rm
c})/2=1.5$. Similarly, we define $\tilde{\mu}(x)$ from
$(\tilde{\mu}_{i})_{i=1}^{N}$. Using this $X^{I}$ and $\mu(x)$, we obtain
$\mu^{I}=\tilde{\mu}(X^{I})$. More explicitly, $\mu^{I}$ is determined as
follows. First, we find $i_{*}$ satisfying $\rho_{i_{*}-1}>1.5$ and
$\rho_{i_{*}}<1.5$. From the construction of $\rho(x)$, we obtain

$X^{I}=\frac{1.5-\rho_{i_{*}}}{\rho_{i_{*}-1}-\rho_{i_{*}}}\frac{i_{*}-1}{N}+\frac{\rho_{i_{*}-1}-1.5}{\rho_{i_{*}-1}-\rho_{i_{*}}}\frac{i_{*}}{N}.$
(6.3)

We then have

$\mu^{I}=\frac{1.5-\rho_{i_{*}}}{\rho_{i_{*}-1}-\rho_{i_{*}}}\tilde{\mu}_{i_{*}-1}+\frac{\rho_{i_{*}-1}-1.5}{\rho_{i_{*}-1}-\rho_{i_{*}}}\tilde{\mu}_{i_{*}}.$
(6.4)

Using this formula, we have $\mu^{I}=7.7\times 10^{-3}$ for the data of
$\kappa_{\Lambda}=0.5$, and $\mu^{I}=2.1\times 10^{-4}$ for the data of
$\kappa_{\Lambda}=1.5$. In Fig. 8, $\mu^{I}$ obtained by (6.4) are plotted for
several values of $\kappa_{\Lambda}$.

Now, we compare the numerical results with the theoretical predictions. We
developed the theory of the steady state in the weak noise limit $T_{\rm
eff}\ll 1$ and the macroscopic limit $N\gg 1$, with particularly focusing on
the two regimes $\kappa_{\Lambda}\gg 1$ and $\kappa_{\Lambda}\ll 1$. When
$\kappa_{\Lambda}\gg 1$, the chemical potential at the interface is $\mu_{\rm
c}=0$, as shown in Appendix B. When $\kappa_{\Lambda}\ll 1$, we have the
formula (4.22), where $\mu_{\rm c}=0$, $\sigma^{\rm L}=2.5$, and $\sigma^{\rm
G}=1.5$ were already determined for the model in the first paragraph of this
section. $X^{\rm eq}$ in the right-side of (4.22) is determined as $X^{\rm
eq}=1/2$ using (2.18) with $\bar{\rho}=1.5$, $\rho^{\rm L}_{\rm c}=2.5$,
$\rho^{\rm G}_{\rm c}=0.5$. By substituting these values into (4.22), we
obtain

$\displaystyle\mu_{X_{*}}^{\rm I}=\frac{\phi}{6}+O(\phi^{2}).$ (6.5)

The dotted lines in Fig. 8 represent the theoretical predictions $\mu^{\rm
I}/\phi=1/6$ for $\kappa_{\Lambda}\ll 1$ and $\mu^{\rm I}/\phi=0$ for
$\kappa_{\Lambda}\gg 1$. These are consistent with the numerical result in
Fig. 8.

It is quite interesting to elucidate the $\kappa_{\Lambda}$ dependence of
$\mu^{I}$ quantitatively. In particular, one may conjecture a phase transition
at some value of $\kappa_{\Lambda}$ in the limit $N\to\infty$. To investigate
the validity of this naive conjecture, we have to numerically study the
asymptotic behavior for $N\to\infty$, $T_{\rm eff}\to 0$ and $\phi\to 0$ in
more detail. From the theoretical viewpoint, we need to develop a calculation
method for thermodynamic properties of the system with finite
$\kappa_{\Lambda}$.

Figure 7: (a) Density profile $\rho_{i}$ for $i/N$. (b) Chemical potential
profile $\tilde{\mu}_{i}$ for $i/N$. $\kappa_{\Lambda}=0.5$ and
$\kappa_{\Lambda}=1.5$. The dotted line in (b) represents the equilibrium
profile of chemical potential $\tilde{\mu}_{i}=\mu_{\rm c}=0$.  Figure 8:
$\kappa_{\Lambda}$ dependence of $\mu^{I}$. Square symbols show the numerical
results for $N=64$. The two dotted lines represent the theoretical predictions
$\mu^{I}=\phi/6$ for $\kappa_{\Lambda}\ll 1$ and $\mu^{I}=0$ for
$\kappa_{\Lambda}\gg 1$.

##  7 Concluding remarks

We have derived the variational function determining the steady state for a
boundary-driven diffusive system with $\kappa_{\Lambda}\ll 1$. The result is
consistent with global thermodynamics, which is an extended framework of
thermodynamics. Before ending the paper, we present a few remarks.

The first remark is on the boundary conditions. It is natural to study a
system with different boundary conditions leading to the same most probable
profile. As more familiar boundary conditions, one considers the case that
chemical potentials at boundaries are fixed. However, as far as we attempted,
we could not evaluate the Zubarev-McLennan representation for this case. To
study the boundary condition dependence of the system is a future problem.

Second, in general, fluctuating hydrodynamics is regarded as a mesoscopic
model obtained by coarse-graining microscopic dynamics. Thus, it is a
significant problem to find relationship between microscopic dynamics and the
discrete fluctuating dynamics. As the first step of such studies, parameter
values of the model should be determined from the observation of microscopic
dynamics. In particular, it seems highly challenging to identify the value of
$\kappa_{\Lambda}$ from microscopic models.

Third, as a generalization of the present model, one may consider a discrete
fluctuating hydrodynamics describing liquid-gas phase coexistence in heat
conduction systems. One can numerically study the model by changing
$\kappa_{\Lambda}$. It is interesting to observe the deviation of the
interface temperature from the equilibrium coexistence temperature.
Furthermore, following the theoretical method presented in this paper, we may
develop a theory for calculating the deviation. We conjecture that the
deviation formula is equivalent to that predicted by global thermodynamics.

The most important future work is an experimental observation of non-
equilibrium phase coexistence in which metastable states are stable as the
influence of a non-equilibrium current. As shown in this paper, the phenomenon
is expected to occur in systems described by a discrete fluctuating
hydrodynamics. However, it is not obvious whether experimental systems are
described by a discrete fluctuating hydrodynamics. It would be interesting to
clarify an experimental condition for realizing a discrete fluctuating
hydrodynamics.

## Acknowledgments

The authors thank F. Kagawa, K. Saito, and Y. Yamamura for useful discussions
about experiments on non-equilibrium phase coexistence; M. Kobayashi, S.
Yukawa, and A. Yoshida for discussions on numerical simulations of non-
equilibrium phase coexistence; and K. Hiura, M. Itami, H. Nakano and Y.
Nakayama for discussions on fluctuating hydrodynamics. We also thank B.
Derrida, C. Maes, H. Spohn, and H. Tasaki for suggesting us to study boundary-
driven diffusive systems. This work was supported by JSPS KAKENHI Grant Number
JP22H01144 and JP23K22415.

Data Availability: Data sharing not applicable to this article as no datasets
were generated or analysed during the current study.

Conflict of Interest: The authors have no financial or proprietary interests
in any material discussed in this article.

##  Appendix A Derivation of (2.14)

In this section, we derive the the stationary distribution of ${\bm{\rho}}$
for the equilibrium system with $\phi=0$. In general, to analyze statistical
properties of ${\bm{\rho}}(t)$ obeying (2.8) with (2.10), it is convenient to
use a new variable ${\bm{\psi}}=(\psi_{i})_{i=1}^{N}$ for ${\bm{\rho}}$. The
introduction of the new variable is another purpose of this section. See
Appendix C for the analysis using ${\bm{\psi}}$.

We first define $\psi_{i}(t)$ by

$\frac{d\psi_{i}(t)}{dt}=-j_{i}(t)$ (A.1)

with $\psi_{i}(0)=\Lambda\sum_{j=1}^{i}\rho_{j}(0)$ at $t=0$. Substituting
(A.1) into (2.8) and integrating it in time, we have

$\rho_{i}(t)=\frac{\psi_{i}(t)-\psi_{i-1}(t)}{\Lambda}$ (A.2)

for any $t$. Substituting (A.2) into ${\cal F}({\bm{\rho}})$ given by (2.7),
we can define ${\cal F}({\bm{\psi}})$ from ${\cal F}({\bm{\rho}})$. Taking the
derivative of ${\cal F}({\bm{\psi}})$ in $\psi_{i}$, we obtain

$\displaystyle\frac{\partial{\cal F}}{\partial\psi_{i}}$
$\displaystyle=\frac{\partial{\cal F}}{\partial\rho_{i}}-\frac{\partial{\cal
F}}{\partial\rho_{i+1}}$ $\displaystyle=\tilde{\mu}_{i}-\tilde{\mu}_{i+1}.$
(A.3)

Using (A.1) and (A.3), we rewrite (2.10) as

$\frac{d\psi_{i}}{dt}=\frac{\sigma(\rho_{i}^{\rm
m})}{\Lambda}\left(-\frac{\partial{\cal
F}}{\partial\psi_{i}}-\phi\delta_{i,N}\right)+\sqrt{\frac{2\sigma(\rho_{i}^{\rm
m})T_{\rm eff}}{\Lambda}}\cdot\xi_{i}$ (A.4)

with

$\rho_{i}^{\rm m}=\frac{\psi_{i+1}-\psi_{i-1}}{2\Lambda}.$ (A.5)

Because $\sigma(\rho_{i}^{\rm m})$ is independent of $\psi_{i}$, the
multiplication of $\sigma(\rho_{i}^{\rm m})$ and $\xi_{i}$ is uniquely
determined independently of the multiplication rule. From (A.4), we obtain the
stationary distribution of ${\bm{\psi}}$ for the equilibrium system with
$\phi=0$ as

${\cal P}_{{\rm
eq}}({\bm{\psi}})=\frac{1}{Z^{\prime}}\exp\left[-\frac{1}{T_{\rm eff}}{\cal
F}({\bm{\psi}})\right],$ (A.6)

where $Z^{\prime}$ is the normalization constant. This gives the stationary
distribution of ${\bm{\rho}}$ as (2.14).

##  Appendix B Analysis of the Continuum model

In this section, we analyze the continuum model (2.2) and (2.3) which
corresponds to the discrete model (2.8) and (2.10) in the limit
$\kappa_{\Lambda}\to\infty$ and $N\to\infty$ with $L$ and $\kappa$ fixed.
Explicitly, we derive the phase coexistence condition (2.16) for the
equilibrium system with $\phi=0$ and the chemical potential at the interface
for the case $\phi>0$.

###  B.1 Equilibrium phase coexistence

Stationary solutions of (2.2) and (2.3) with $\phi=0$ and $T=0$ satisfy

$\partial_{x}\frac{\delta{\cal F}}{\delta\rho(x)}=0.$ (B.1)

Using (2.5), (B.1) is explicitly written as

$f^{\prime}(\rho)-\kappa\partial_{x}^{2}\rho=\mu_{\rm c},$ (B.2)

where $\mu_{\rm c}$ is a constant in $x$. Furthermore, multiplying
$(\partial_{x}\rho)$ with (B.2) and integrating in $x$, we obtain

$f(\rho)-\frac{\kappa}{2}(\partial_{x}\rho)^{2}-\mu_{\rm c}\rho=p_{\rm c},$
(B.3)

where $p_{\rm c}$ is also a constant in $x$.

We first consider necessary conditions under which there is a stationary and
spatially inhomogeneous solution, which we call the phase coexistence
solution, because it connects two stationary and spatially homogeneous
solutions that represent a liquid phase and a gas phase, respectively. To make
the argument clear, we take the limit $L\to\infty$ and we assume
$\rho(0)>\rho(\infty)$ for the phase coexistence solution. Because the phase
coexistence solution approaches to the stationary and spatially homogeneous
solutions $\rho(0)$ and $\rho(\infty)$ as $x\to 0$ and $x\to\infty$, we have

$\partial_{x}\rho(0)=0,\qquad\partial_{x}^{2}\rho(0)=0,\qquad\partial_{x}\rho(\infty)=0,\qquad\partial_{x}^{2}\rho(\infty)=0.$
(B.4)

Therefore, (B.2) and (B.3) lead to necessary conditions as

$\mu(\rho(0))=\mu(\rho(\infty))=\mu_{\rm c},$ (B.5)

and

$f(\rho(0))-\mu_{\rm c}\rho(0)=f(\rho(\infty))-\mu_{\rm c}\rho(\infty)=p_{\rm
c},$ (B.6)

which is further rewritten as

$\mu_{\rm c}=\frac{f(\rho(0))-f(\rho(\infty))}{\rho(0)-\rho(\infty)}.$ (B.7)

For the function $f(\rho)$ with two local minima, the conditions (B.5) and
(B.7) represents the common tangent line at the special values $\rho=\rho^{\rm
L}_{\rm c}$ and $\rho=\rho^{\rm G}_{\rm c}$. We set $\rho^{\rm L}_{\rm
c}>\rho^{\rm G}_{\rm c}$ without loss of generality. We thus identify
$\rho(0)=\rho^{\rm L}_{\rm c}$ and $\rho(\infty)=\rho^{\rm G}_{\rm c}$, and
the values of the constants $\mu_{\rm c}$ and $p_{\rm c}$ are also determined.
Note that (B.5) and (B.6) are regarded as the conditions giving $\rho^{\rm
L}_{\rm c}$, $\rho^{\rm G}_{\rm c}$, $\mu_{\rm c}$ and $p_{\rm c}$ by the form

$\mu(\rho^{\rm L}_{\rm c})=\mu(\rho^{\rm G}_{\rm c})=\mu_{\rm c},\qquad
p(\rho^{\rm L}_{\rm c})=p(\rho^{\rm G}_{\rm c})=p_{\rm c}.$ (B.8)

Now, suppose that $\bar{\rho}$ satisfies $\rho^{\rm G}_{\rm
c}<\bar{\rho}<\rho^{\rm L}_{\rm c}$, where $\rho^{\rm G}_{\rm c}$ and
$\rho^{\rm L}_{\rm c}$ are determined by (B.8). By setting $\rho(0)=\rho^{\rm
L}_{\rm c}$ and $\rho(\infty)=\rho^{\rm G}_{\rm c}$, we solve (B.2) with the
determined value of $\mu_{\rm c}$. We here notice that (B.2) is interpreted as
Newton’s equation describing the motion of a point particle under a potential
field, where $\rho$ and $x$ correspond to position and time, respectively.
$\kappa$ is interpreted as the mass, and the potential function $V(\rho)$ is
given by $V(\rho)\equiv\mu_{\rm c}\rho-f(\rho)$. (B.3) represents the energy
conservation for the equation of motion. From (B.8), we find that $\rho^{\rm
L}_{\rm c}$ and $\rho^{\rm G}_{\rm c}$ are local maximal points with the same
potential value. Therefore, the phase coexistence solution $\rho(x)$ with
$\rho(0)=\rho^{\rm L}_{\rm c}$ and $\rho(\infty)=\rho^{\rm G}_{\rm c}$ is
given by the heteroclinic orbit connecting the two maximal points with the
same potential value. This result corresponds to the statement involving
(2.16) in the main text.

###  B.2 Non-equilibrium phase coexistence

Stationary solutions of (2.2) and (2.3) with $\phi>0$ and $T=0$ satisfy

$\partial_{x}\frac{\delta{\cal F}}{\delta\rho(x)}=-\frac{J}{\sigma(\rho(x))},$
(B.9)

which corresponds to (2.28) in the continuum limit $\kappa_{\Lambda}\to\infty$
and $N\to\infty$, where $J$ is a constant given by

$J=\left(\int^{L}_{0}dx~{}\frac{1}{\sigma(\rho(x))}\right)^{-1},$ (B.10)

which corresponds to (2.29) in the continuum limit $\kappa_{\Lambda}\to\infty$
and $N\to\infty$. By analyzing (B.9), we determine the value of the chemical
potential at the interface.

First, to uniquely identify the interface position, we introduce a scaled
coordinate $\hat{x}=x/L$. Taking the limit $L\to\infty$, we find that the
interface width in the scaled coordinate space becomes zero. Thus, the
interface position $X$ in the $\hat{x}$ space is given by the discontinuous
point of $\rho(\hat{x})$ in the limit $L\to\infty$. We then define the
chemical potential at the interface by $\mu(\hat{x}=X)$. To determine the
value of $\mu(\hat{x}=X)$, we consider the generalized chemical potential
$\tilde{\mu}(x)$ given by

$\tilde{\mu}(x)=\mu(\rho(x))-\kappa\partial_{x}^{2}\rho,$ (B.11)

which corresponds to (2.12) in the limit $\kappa_{\Lambda}\to\infty$ and
$N\to\infty$. By integrating (B.9) in the range $[x_{1},x_{2}]$, we have

$\tilde{\mu}(x_{2})-\tilde{\mu}(x_{1})=-J\int_{x_{1}}^{x_{2}}dx~{}\frac{1}{\sigma(\rho(x))}$
(B.12)

for any $x_{1}$ and $x_{2}$. Even though $\sigma(\rho(\hat{x}))$ is
discontinuous at $\hat{x}=X$, the integration in the right-side of (B.12)
gives a continuous function in $x_{2}$ and $x_{1}$. Thus,
$\tilde{\mu}(\hat{x})$ is a continuous function in the limit $L\to\infty$.

Here, the key idea for the determination of $\mu(\hat{x}=X)$ is the
introduction of the generalized pressure $\tilde{p}$ satisfying

$\rho\partial_{x}\tilde{\mu}=\partial_{x}\tilde{p}.$ (B.13)

We can explicitly derive $\tilde{p}$ from (B.11) and (B.13) as

$\tilde{p}=p(\rho)-\kappa\rho\partial_{x}^{2}\rho+\frac{\kappa}{2}(\partial_{x}\rho)^{2},$
(B.14)

which was first obtained by van der Waals vanPress . By integrating (B.13) in
the range $[x_{1},x_{2}]$, we have

$\tilde{p}(x_{2})-\tilde{p}(x_{1})=-J\int_{x_{1}}^{x_{2}}dx~{}\frac{\rho(x)}{\sigma(\rho(x))},$
(B.15)

for any $x_{1}$ and $x_{2}$. We find from (B.15) that $\tilde{p}(\hat{x})$ is
a continuous function in the limit $L\to\infty$. Note that
$\tilde{p}(x)=p_{\rm c}$ for the equilibrium system, where $p_{\rm c}$ is the
constant given in (B.3).

Now, using the continuity of $\tilde{\mu}(\hat{x})$ and $\tilde{p}(\hat{x})$
at $\hat{x}=X$, we can determine the values of $\rho(\hat{x}=X-\epsilon)$ and
$\rho(\hat{x}=X+\epsilon)$ for small $\epsilon>0$ in the limit $L\to\infty$.
Let $\rho_{-}$ and $\rho_{+}$ be $\rho(\hat{x}=X-\epsilon)$ and
$\rho(\hat{x}=X-\epsilon)$ for $\epsilon\to 0+$ after taking the limit
$L\to\infty$. The continuity of $\tilde{\mu}(\hat{x})$ and
$\tilde{p}(\hat{x})$ at $\hat{x}=X$ leads to

$\mu(\rho_{-})=\mu(\rho_{+}),\qquad p(\rho_{-})=p(\rho_{+}).$ (B.16)

Recalling (B.8), we obtain

$\rho_{-}=\rho^{\rm L}_{\rm c},\qquad\rho_{+}=\rho^{\rm G}_{\rm c}.$ (B.17)

Therefore, the chemical potential at the interface is equal to $\mu_{\rm c}$.
The result is mentioned in the third paragraph of Sec. 2.3 in the main text.

Finally, we present a remark on the singular continuum description for the
case $\kappa_{\Lambda}\to 0$ in Sec. 2.4. The chemical potential
$\mu(\hat{x})$ is continuous at the interface position $\hat{x}=X$ where
$\rho(\hat{x})$ is discontinuous. In this case, however, $p(\hat{x})$ is
discontinuous at $\hat{x}=X$, as shown by (4.27). That is, (B.13) does not
hold at the interface. This is the most essential difference between the two
cases $\kappa_{\Lambda}\gg 1$ and $\kappa_{\Lambda}\ll 1$.

##  Appendix C Zubarev-McLennan representation

In this section, we derive the Zubarev-McLennan representation (3.1) in Sec.
3.1.

We study stochastic processes of ${\bm{\psi}}=(\psi_{i})_{1\leq i\leq N}$
defined by (A.1). The time evolution of ${\bm{\psi}}$ is described by (A.4).
Let $\hat{{\bm{\psi}}}=({\bm{\psi}}_{t})_{t=0}^{\tau}$ be a trajectory in the
time interval $[0,\tau]$. The path probability density ${\cal P}_{\rm
path}(\hat{{\bm{\psi}}})$ in the system with $\phi>0$ starting from a density
profile sampled from an equilibrium distribution ${\cal P}_{\rm
eq}({\bm{\psi}}_{0})$ is expressed as

$\displaystyle{\cal P}_{\rm path}(\hat{{\bm{\psi}}})$ $\displaystyle={\cal
P}_{\rm eq}({\bm{\psi}}_{0})$ $\displaystyle\times{\rm
const}\times\exp\left(-\frac{1}{4T_{\rm
eff}}\int_{0}^{\tau}dt\sum_{i=1}^{N}\frac{\Lambda}{\sigma(\rho_{i}^{\rm
m})}\left[\frac{d\psi_{i}}{dt}+\frac{\sigma(\rho_{i}^{\rm
m})}{\Lambda}\left(\frac{\partial{\cal
F}}{\partial\psi_{i}}+\phi\delta_{i,N}\right)\right]^{2}\right).$ (C.1)

For the time-reversed trajectory
$\hat{{\bm{\psi}}}^{\dagger}=({\bm{\psi}}_{\tau-t})_{t=0}^{\tau}$ of
$\hat{{\bm{\psi}}}$, we have

$\displaystyle{\cal P}_{\rm path}(\hat{{\bm{\psi}}}^{\dagger})$
$\displaystyle={\cal P}_{\rm eq}({\bm{\psi}}_{\tau})$ $\displaystyle\times{\rm
const}\times\exp\left(-\frac{1}{4T_{\rm
eff}}\int_{0}^{\tau}dt\sum_{i=1}^{N}\frac{\Lambda}{\sigma(\rho_{i}^{\rm
m})}\left[-\frac{d\psi_{i}}{dt}+\frac{\sigma(\rho_{i}^{\rm
m})}{\Lambda}\left(\frac{\partial{\cal
F}}{\partial\psi_{i}}+\phi\delta_{i,N}\right)\right]^{2}\right).$ (C.2)

The ratio of the two yields

$\displaystyle\frac{{\cal P}_{\rm path}(\hat{{\bm{\psi}}})}{{\cal P}_{\rm
path}(\hat{{\bm{\psi}}}^{\dagger})}=\exp\left[-\frac{1}{T_{\rm
eff}}\phi\int_{0}^{\tau}dt\frac{d\psi_{N}(t)}{dt}\right],$ (C.3)

where (A.6) has been substituted into ${\cal P}_{\rm eq}({\bm{\psi}}_{0})$ and
${\cal P}_{\rm eq}({\bm{\psi}}_{\tau})$.

To simplify the notation, we introduce the accumulated current

$Q_{\tau}(\hat{{\bm{\psi}}})\equiv\int_{0}^{\tau}dt~{}j_{N}(t).$ (C.4)

Using (A.1) and (C.4), we rewrite (C.3) as

$\displaystyle\frac{{\cal P}_{\rm path}(\hat{{\bm{\psi}}})}{{\cal P}_{\rm
path}(\hat{{\bm{\psi}}}^{\dagger})}=\exp\left[\frac{\phi Q_{\tau}}{T_{\rm
eff}}\right].$ (C.5)

The distribution of ${\bm{\psi}}$ at time $t$ is expressed as

$\displaystyle{\cal P}_{\tau}({\bm{\psi}})$ $\displaystyle=\int{\cal
D}\hat{{\bm{\psi}}}~{}{\cal P}_{\rm
path}(\hat{{\bm{\psi}}})\delta({\bm{\psi}}_{\tau}-{\bm{\psi}})$
$\displaystyle=\int{\cal D}\hat{{\bm{\psi}}}^{\dagger}\exp\left[\frac{\phi
Q_{\tau}(\hat{{\bm{\psi}}})}{T_{\rm eff}}\right]{\cal P}_{\rm
path}(\hat{{\bm{\psi}}}^{\dagger})\delta({\bm{\psi}}_{\tau}-{\bm{\psi}}),$
(C.6)

where we have used ${\cal D}\hat{{\bm{\psi}}}={\cal
D}\hat{{\bm{\psi}}}^{\dagger}$ and (C.5). When the path integration variable
is transformed, the right side of (C.6) is rewritten as

$\displaystyle\int{\cal D}\hat{{\bm{\psi}}}~{}\exp\left[\frac{\phi
Q_{\tau}(\hat{{\bm{\psi}}}^{\dagger})}{T_{\rm eff}}\right]{\cal P}_{\rm
path}(\hat{{\bm{\psi}}})\delta({\bm{\psi}}_{0}-{\bm{\psi}}),$ $\displaystyle=$
$\displaystyle\int{\cal D}\hat{{\bm{\psi}}}~{}\exp\left[-\frac{\phi
Q_{\tau}(\hat{{\bm{\psi}}})}{T_{\rm eff}}\right]{\cal P}_{\rm
path}(\hat{{\bm{\psi}}})\delta({\bm{\psi}}_{0}-{\bm{\psi}}),$ (C.7)

where we have used
$Q_{\tau}(\hat{{\bm{\psi}}}^{\dagger})=-Q_{\tau}(\hat{{\bm{\psi}}})$. We thus
have the relation

$\displaystyle{\cal P}_{\tau}({\bm{\psi}})={\cal P}_{\rm
eq}({\bm{\psi}})\left\langle\exp\left[-\frac{\phi Q_{\tau}}{T_{\rm
eff}}\right]\right\rangle_{{\bm{\psi}}}.$ (C.8)

Taking the limit $\tau\to\infty$, we have

$\displaystyle{\cal P}_{\rm ss}({\bm{\psi}})={\cal P}_{\rm
eq}({\bm{\psi}})\left\langle\exp\left[-\frac{\phi Q}{T_{\rm
eff}}\right]\right\rangle_{{\bm{\psi}}}$ (C.9)

with

$\displaystyle Q=\int_{0}^{\infty}dt~{}j_{N}(t).$ (C.10)

In the limit $T_{\rm eff}\to 0$, we estimate

$\displaystyle\left\langle\exp\left[-\frac{\phi Q}{T_{\rm
eff}}\right]\right\rangle_{{\bm{\psi}}}\simeq\exp\left[-\frac{\phi\left\langle
Q\right\rangle_{{\bm{\psi}}}}{T_{\rm eff}}\right].$ (C.11)

We then expand $\left\langle Q_{\infty}\right\rangle_{{\bm{\psi}}}$ in $\phi$,
we obtain

$\displaystyle{\cal P}_{\rm ss}({\bm{\psi}})={\cal P}_{\rm
eq}({\bm{\psi}})\exp\left[-\frac{\phi\left\langle
Q\right\rangle_{{\bm{\psi}}}^{\rm eq}+O(\phi^{2})}{T_{\rm eff}}\right].$
(C.12)

This is the Zubarev-McLennan representation of the steady state distribution.
By using (A.2), we obtain the stationary distribution of ${\bm{\rho}}$ as the
form (3.1).

##  Appendix D Derivation of (3.23)

In (3.21), we consider the decomposition of $j(0,t)$ into $\Phi(t)$ and
$\Phi_{0}(t)$. In this section, we calculate $\Phi_{0}(t)$.

###  D.1 Preliminaries for the calculation

We first note that (3.14) and (3.16) yield

$\displaystyle\frac{d{\hat{X}}_{\downarrow}}{dt}=\frac{d\hat{Y}}{dt}+\frac{1}{2}\frac{d|D^{\rm
L}|}{dt},\qquad\frac{d{\hat{X}}_{\uparrow}}{dt}=\frac{d\hat{Y}}{dt}-\frac{1}{2}\frac{d|D^{\rm
L}|}{dt},\qquad\frac{d|D^{\rm G}|}{dt}=-\frac{d|D^{\rm L}|}{dt}.$ (D.1)

From the chain rule of the derivative, we have

$\displaystyle\frac{d}{dt}\left[\frac{|D^{\rm L}|}{\sigma^{\rm
L}}+\frac{|D^{\rm G}|}{\sigma^{\rm G}}\right]^{-1}$
$\displaystyle=-\left[\frac{|D^{\rm L}|}{\sigma^{\rm L}}+\frac{|D^{\rm
G}|}{\sigma^{\rm G}}\right]^{-2}\circ\left[\frac{1}{\sigma^{\rm
L}}\frac{d|D^{\rm L}|}{dt}+\frac{1}{\sigma^{\rm G}}\frac{d|D^{\rm
G}|}{dt}\right]$ $\displaystyle=-\left(\frac{1}{\sigma^{\rm
L}}-\frac{1}{\sigma^{\rm G}}\right)\left[\frac{|D^{\rm L}|}{\sigma^{\rm
L}}+\frac{|D^{\rm G}|}{\sigma^{\rm G}}\right]^{-2}\circ\frac{d|D^{\rm
L}|}{dt},$ (D.2)

and

$\displaystyle\frac{d}{dt}\left[\frac{1}{\sigma^{\rm L}}\int_{D^{\rm
L}}dx~{}\psi+\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}}dx~{}\psi\right]=\frac{d}{dt}\left[\frac{1}{\sigma^{\rm
L}}\int_{{\hat{X}}_{\uparrow}(t)}^{{\hat{X}}_{\downarrow}(t)}dx~{}\psi+\frac{1}{\sigma^{\rm
G}}\int_{{\hat{X}}_{\downarrow}(t)}^{{\hat{X}}_{\uparrow}(t)+1}dx~{}\psi\right]$
$\displaystyle=\frac{1}{\sigma^{\rm
L}}\left[\psi({\hat{X}}_{\downarrow}(t),t)\circ\frac{d{\hat{X}}_{\downarrow}}{dt}-\psi({\hat{X}}_{\uparrow}(t),t)\circ\frac{d{\hat{X}}_{\uparrow}}{dt}\right]+\frac{1}{\sigma^{\rm
G}}\left[\psi({\hat{X}}_{\uparrow}(t),t)\circ\frac{d{\hat{X}}_{\uparrow}}{dt}-\psi({\hat{X}}_{\downarrow}(t),t)\circ\frac{d{\hat{X}}_{\downarrow}}{dt}\right]$
$\displaystyle\qquad+\frac{1}{\sigma^{\rm L}}\int_{D^{\rm
L}}dx~{}\partial_{t}\psi+\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}}dx~{}\partial_{t}\psi$ $\displaystyle=\left(\frac{1}{\sigma^{\rm
L}}-\frac{1}{\sigma^{\rm
G}}\right)\left[\left(\psi({\hat{X}}_{\downarrow}(t),t)-\psi({\hat{X}}_{\uparrow}(t),t)\right)\circ\frac{d\hat{Y}}{dt}+\frac{1}{2}\left(\psi({\hat{X}}_{\downarrow}(t),t)+\psi({\hat{X}}_{\uparrow}(t),t)\right)\circ\frac{d|D^{\rm
L}|}{dt}\right]$ $\displaystyle\qquad+\frac{1}{\sigma^{\rm L}}\int_{D^{\rm
L}}dx~{}\partial_{t}\psi+\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}}dx~{}\partial_{t}\psi.$ (D.3)

###  D.2 Relation between $\Phi(t)$ and $j(0,t)$

Noting the two derivatives (D.2) and (D.3), we define

$\displaystyle\Phi_{1}$ $\displaystyle\equiv-\left(\frac{1}{\sigma^{\rm
L}}-\frac{1}{\sigma^{\rm G}}\right)\left[\frac{|D^{\rm L}|}{\sigma^{\rm
L}}+\frac{|D^{\rm G}|}{\sigma^{\rm G}}\right]^{-2}\left[\frac{1}{\sigma^{\rm
L}}\int_{D^{\rm L}(t)}dx~{}\psi+\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}(t)}dx~{}\psi\right]\circ\frac{d|D^{\rm L}|}{dt},$ (D.4)
$\displaystyle\Phi_{2}$ $\displaystyle\equiv\left(\frac{1}{\sigma^{\rm
L}}-\frac{1}{\sigma^{\rm G}}\right)\left[\frac{|D^{\rm L}|}{\sigma^{\rm
L}}+\frac{|D^{\rm G}|}{\sigma^{\rm G}}\right]^{-1}$
$\displaystyle\times\left[\left(\psi({\hat{X}}_{\downarrow}(t),t)-\psi({\hat{X}}_{\uparrow}(t),t)\right)\circ\frac{d\hat{Y}}{dt}+\frac{1}{2}\left(\psi({\hat{X}}_{\downarrow}(t),t)+\psi({\hat{X}}_{\uparrow}(t),t)\right)\circ\frac{d|D^{\rm
L}|}{dt}.\right]$ (D.5)

Then, the definition of $\Phi(t)$ in (3.22) leads to

$\displaystyle\Phi(t)=\Phi_{1}(t)+\Phi_{2}(t)+j(0,t)$ (D.6)

by using (3.20). Here, from the piece-wise linear nature of $\psi(x,t)$, we
have

$\displaystyle\frac{1}{\sigma^{\rm L}}\int_{D^{\rm
L}(t)}dx~{}\psi+\frac{1}{\sigma^{\rm G}}\int_{D^{\rm
G}(t)}dx~{}\psi=\frac{1}{2}\left[\frac{|D^{\rm L}|}{\sigma^{\rm
L}}+\frac{|D^{\rm G}|}{\sigma^{\rm
G}}\right]\left(\psi({\hat{X}}_{\downarrow}(t),t)+\psi({\hat{X}}_{\uparrow}(t),t)\right).$
(D.7)

See also (3.31) for the same equation. Using this relation, we obtain

$\Phi_{1}+\Phi_{2}=\left(\frac{1}{\sigma^{\rm L}}-\frac{1}{\sigma^{\rm
G}}\right)\left[\frac{|D^{\rm L}(t)|}{\sigma^{\rm L}}+\frac{|D^{\rm
G}(t)|}{\sigma^{\rm
G}}\right]^{-1}\left(\psi({\hat{X}}_{\downarrow}(t),t)-\psi({\hat{X}}_{\uparrow}(t),t)\right)\circ\frac{d\hat{Y}}{dt}.$
(D.8)

Finally, we note that

$\displaystyle\psi({\hat{X}}_{\downarrow}(t),t)-\psi({\hat{X}}_{\uparrow}(t),t)$
$\displaystyle=\int_{D^{\rm L}}dx~{}\partial_{x}\psi,$ $\displaystyle=|D^{\rm
L}(t)|(\rho^{\rm L}(t)-\bar{\rho}).$ (D.9)

Substituting this into (D.8), we obtain (3.23) where
$\Phi_{0}\equiv\Phi_{1}+\Phi_{2}$.

## References

  * (1) J. R. Thome, Boiling in microchannels: a review of experiment and theory, Int. J. Heat and Fluid flow 25, 128-139 (2004).
  * (2) E. Ben-Jacob and P. Garik, The formation of patterns in nonequilibrium growth, Nature 343, 523-530 (1990).
  * (3) G. Ahlers, L. I. Berge, and D. S. Cannell, Thermal convection in the presence of a first-order phase change, Phys. Rev. Lett. 70, 2399 (1993).
  * (4) A. Yoshida, N. Nakagawa, and S.-i. Sasa, Heat-induced liquid hovering in liquid-gas coexistence under gravity, ArXiv:2310.058171.
  * (5) J.-Q. Zhong, D. Funfschilling, and G. Ahlers, Enhanced heat transport by turbulent two-phase Rayleigh-Benard convection, Phys. Rev. Lett. 102, 124501 (2009).
  * (6) S. Weiss and G. Ahlers, Nematic-isotropic phase transition in turbulent thermal convection, J. Fluid Mech. 737, 308-328 (2013).
  * (7) P. Urbana, D. Schmoranzerb, P. Hanzelkaa, K. R. Sreenivasanc, and L. Skrbekb, Anomalous heat transport and condensation in convection of cryogenic helium, Proc. Nat. Acad. Sci. 110, 8036-8039 (2013).
  * (8) M. E. Cates and J. Tailleur, Motility-Induced Phase Separation, Annual Review of Condensed Matter Physics 6, 219-244 (2015).
  * (9) N. Nakagawa and S.-i. Sasa, Liquid-gas transitions in steady heat conduction, Phys. Rev. Lett. 119, 260602 (2017).
  * (10) N. Nakagawa and S.-i. Sasa, Global thermodynamics for heat conduction states, J. Stat. Phys. 177, 825-888 (2019).
  * (11) N. Nakagawa and S.-i. Sasa, Unique extension of the maximum entropy principle to phase coexistence in heat conduction, Phys. Rev. Res. 4, 033155 (2022).
  * (12) M. Kobayashi, N. Nakagawa, and S.-i. Sasa, Control of metastable states by heat flux in the Hamiltonian Potts model, Phys. Rev. Lett. 130, 247102 (2023).
  * (13) S.-i. Sasa, N. Nakagawa, M. Itami, and Y. Nakayama, Stochastic order parameter dynamics for phase coexistence in heat conduction, Phys. Rev. E 103, 062129 (2021).
  * (14) R. Schmitz, Fluctuations in nonequilibrium fluids, Physics Reports 171, 1-58 (1988).
  * (15) G. L. Eyink, Dissipation and large thermodynamic fluctuations, J. Stat. Phys. 61, 533-572 (1990).
  * (16) H. Spohn, Large scale dynamics of interacting particles (Springer-Verlag, Berlin Heidelberg, 1991).
  * (17) L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Macroscopic fluctuation theory, Rev. Mod. Phys. 87, 593 (2015).
  * (18) Y. Kuramoto, Effects of diffusion on the fluctuations in open chemical systems, Prog. Theor. Phys. 52, 711 (1974).
  * (19) Y. Yanagisawa and S.-i. Sasa, Phase coexistence in a weakly stochastic reaction-diffusion system, arXiv:2403.19198
  * (20) Y. Kado and S.-i. Sasa, Microscopic singularity of an entropic force in stochastic order parameter dynamics, Phys. Rev. Lett. 132, 057101 (2024).
  * (21) D. N. Zubarev, Nonequilibrium Statistical Thermodynamics, (Consultants Bureau, New York, 1974).
  * (22) J. A. McLennan, Phys. Fluids 3, 493 (1960); Introduction to Non-equilibrium Statistical Mechanics (Prentice-Hall, 1988).
  * (23) G. E. Crooks, Path-ensemble averages in systems driven far from equilibrium, Phys. Rev. E 61, 2361 (2000).
  * (24) S.-i. Sasa, A perturbation theory for large deviation functionals in fluctuating hydrodynamics, J. Phys. A: Math. Theor. 41, 045006 (2008).
  * (25) T. S. Komatsu and N. Nakagawa, Expression for the stationary distribution in nonequilibrium steady states, Phys. Rev. Lett. 100, 030601 (2008).
  * (26) C. Maes and K. Netočný, Rigorous meaning of McLennan ensembles, J. Math. Phys. 51, 015219 (2010).
  * (27) J. S. Rowlinson, Translation of J. D. van der Waals’ “The thermodynamik theory of capillarity under the hypothesis of a continuous variation of density”, J. Stat. Phys. 20 197 (1979).