#  Control of metastable states by heat flux in the Hamiltonian Potts model

Michikazu Kobayashi  School of Environmental Science and Engineering, Kochi
University of Technology, Miyanoguchi 185, Tosayamada, Kami, Kochi 782-8502,
Japan  Naoko Nakagawa  Department of Physics, Ibaraki University, Mito
310-8512, Japan  Shin-ichi Sasa  Department of Physics, Kyoto University,
Kyoto, 606-8502 Japan

###### Abstract

The local equilibrium thermodynamics is a basic assumption of macroscopic
descriptions of the out of equilibrium dynamics for Hamiltonian systems. We
numerically analyze the Hamiltonian Potts model in two dimensions to study the
violation of the assumption for phase coexistence in heat conduction. We
observe that the temperature of the interface between ordered and disordered
states deviates from the equilibrium transition temperature, indicating that
metastable states at equilibrium are stabilized by the influence of a heat
flux. We also find that the deviation is described by the formula proposed in
an extended framework of the thermodynamics.

Introduction.— The macroscopic dynamics of a Hamiltonian system driven by non-
equilibrium boundary conditions is expected to be described by hydrodynamic
equations with local equilibrium thermodynamics [1, 2, 3, 4, 5, 6]. However,
there are exceptional cases, such as the shear flow near a liquid–gas critical
point [7], where suppression of critical fluctuations by the shear flow
modifies the thermodynamic properties [8], and as a result, the local
equilibrium thermodynamics is violated. A natural question is whether the
violation of the local equilibrium thermodynamics occurs except at critical
points. The aim of this Letter is to provide a definitive example by studying
non-equilibrium dynamics near the first-order transition point.

A first-order phase transition is a different type of singularity than the
critical point [9, 10, 11, 12]. A characteristic feature of a first-order
phase transition is the existence of hysteresis [13]. For an order-disorder
transition [14, 15], the observed transition temperature when a material is
cooled from the disordered state is lower than the transition temperature when
the same material is heated from the ordered state. Such transition
temperatures, generally both in cooling and heating, deviate from the
equilibrium transition temperature $T_{\mathrm{c}}$. Thus, the supercooled
disordered or the superheated ordered states are often observed as a transient
dynamical process. These observations may suggest that the metastable states
become steady states when the system sets up at non-equilibrium conditions. To
explore this possibility, we consider steady states in heat conduction where
two heat baths with temperatures $T_{1}$ and $T_{2}$ are attached at the left
and right sides of the system.

We assume that $T_{1}$ and $T_{2}$ satisfy $T_{1}\leq T_{\mathrm{c}}\leq
T_{2}$ to observe phase coexistence, where ordered and disordered states
appear at the low and high temperature sides, respectively, with a unique
interface separating the two phases. Our main question in this Letter is
whether the temperature of the interface is equal to $T_{\mathrm{c}}$. If the
local equilibrium thermodynamics is assumed to hold at each point of the
system, the interface temperature should be equal to the equilibrium
transition temperature. However, the validity of this assumption is not
obvious because of the existence of metastable states. To our best knowledge,
there have been no experimental studies on this question, while a rich variety
of non-equilibrium phase-coexistence phenomena have been studied including
flow boiling heat transfer and pattern formation in crystal growth [17, 18,
19, 20, 21, 22].

We study the interface temperature by numerical simulations of a model that
exhibits phase coexistence in steady heat conduction. Since there is no
standard model for order-parameter dynamics with conducting energy, in this
Letter, we propose the Hamiltonian Potts model in two dimensions. We first
confirm the coexistence of ordered and disordered states in an isolated system
by numerically solving the Hamiltonian equation. Then, imposing a heat flux at
the boundary with the total energy fixed, we produce the phase coexistence in
steady heat conduction. A remarkable property in this system is that the
interface temperature deviates from the equilibrium transition temperature.
This indicates that metastable states are stabilized. That is, metastability
is controlled by the heat flux. Furthermore, we find that the deviation is
well-fitted by the formula proposed in an extended framework of thermodynamics
which we call global thermodynamics [24, 25, 26]. It provides a quantitative
prediction of the phase coexistence in the steady heat conduction, in contrast
to other extended frameworks of thermodynamics [27, 28, 29, 30, 31, 32, 33,
34].

Model and observables.— We study the Hamiltonian Potts model with $q$-fold
symmetry in two dimensions. Let $\phi^{a}({\boldsymbol{r}})$, $a=1,2,\cdots
q-1$, be a $q-1$ components field defined in a rectangular region
$D\equiv[0,L_{x}]\times[0,L_{y}]$ with $L_{x}>L_{y}$. We express
$(\phi^{1},\cdots,\phi^{q-1})$ as $\phi$. The conjugate momentum field of
$\phi^{a}({\boldsymbol{r}})$ is expressed by $\pi^{a}({\boldsymbol{r}})$,
where $\pi$ represents $(\pi^{1},\cdots,\pi^{q-1})$. We assume the total
Hamiltonian ${\cal H}$ as

${\cal
H}(\phi,\pi)=\int_{D}d^{2}{\boldsymbol{r}}\left\\{\frac{1}{2}\sum_{a=1}^{q-1}\left[(\pi^{a})^{2}+\left|{\boldsymbol{\nabla}}\phi^{a}\right|^{2}\right]+V(\phi)\right\\},$
(1)

where the potential $V(\phi)$ possesses $q$ symmetric minima in
$\mathbb{R}^{q-1}$. Let $\mu_{k}$, $1\leq k\leq q$, be coordinates of the
$q$-vertices for the regular ($q-1$)-simplex in $\mathbb{R}^{q-1}$ as
illustrated in Fig. 1 (a) for the cases $q=2$, $3$, and $4$, where the
centroid of the simplex is located at the origin. The explicit form of
$\mu_{k}$ is given in Supplemental Material [35]. As a potential with minima
at $\mu_{k}$, we set

$V(\phi)=\frac{1}{2}\prod_{k=1}^{q}\sum_{a=1}^{q-1}(\phi^{a}-\mu_{k}^{a})^{2}.$
(2)

This potential is regarded as a continuous extension of the standard $q$-state
Potts model. We thus expect that the equilibrium statistical mechanics for
${\cal H}$ describes the same phase transitions as those observed in the
standard Potts model. As for the standard $q$-state Potts model, the model
exhibits a first-order transition at a temperature $T_{\mathrm{c}}$ for $q\geq
5$ in equilibrium [23]. In this Letter, we numerically study the case $q=11$
for $L_{x}=384$ and $L_{y}=64$. The system is spatially discretized with a
grid spacing $\Delta x=1/8$ [35]. The Boltzmann constant is set to unity.

We define the local kinetic energy density per unit degree of freedom and the
local order parameter as

$\displaystyle\hat{T}({\boldsymbol{r}})=\frac{\sum_{a=1}^{q-1}(\pi^{a}({\boldsymbol{r}}))^{2}}{2(q-1)},\quad\hat{m}(\boldsymbol{r})\equiv\sum_{a=1}^{q-1}\phi^{a}({\boldsymbol{r}})\mu_{1}^{a},$
(3)

where the direction of symmetry breaking is fixed to the $\mu_{1}$ direction
in our numerical simulations by choosing specific initial conditions [35].
Their average over the $y$-direction with $x$ fixed is expressed as

$\displaystyle[\hat{A}]_{x}=\frac{1}{L_{y}}\int_{0}^{L_{y}}dy\>\hat{A}({\boldsymbol{r}}),$
(4)

where $\hat{A}$ is $\hat{T}$ or $\hat{m}$.

Figure 1:  (a) Examples of $(q-1)$-simplexes with $q$-vertices in
$(q-1)$-dimensional space: a line segment with $\mu_{1}$ and $\mu_{2}$ for
$q=2$ (left), a regular triangle with $\mu_{1}$, $\mu_{2}$ and $\mu_{3}$ for
$q=3$ (middle), and a regular tetrahedron with $\mu_{1}$, $\mu_{2}$,
$\mu_{3}$, and $\mu_{4}$ for $q=4$ (right). (b) Phase diagram for the $q=11$
Hamiltonian Potts model (1). The system shows the first–order phase transition
at $T=T_{\mathrm{c}}$ when $T$ is changed, while the ordered and disordered
states coexist in the range $E_{1}<E<E_{2}$. $T_{\mathrm{c}}$, $E_{1}$, and
$E_{2}$ are numerically estimated to be $T_{\mathrm{c}}\simeq 0.15$,
$E_{1}/(L_{x}L_{y})\simeq 0.66$, and $E_{2}/(L_{x}L_{y})\simeq 1.09$ [35].

Equilibrium phase coexistence.— We first examine the equilibrium phase diagram
by numerically investigating the isothermal dynamics [35]. The first-order
transition is observed at $T=T_{\mathrm{c}}\simeq 0.15$, where the energy
density changes discontinuously at $T=T_{\mathrm{c}}$. The system is occupied
by the ordered and disordered states for $T<T_{\mathrm{c}}$ and
$T>T_{\mathrm{c}}$ , respectively. The phase diagram is shown in Fig. 1 (b).

With the knowledge of the phase diagram, we concentrate on an isolated system
with the total energy $E$ fixed. The time evolution is given by the
Hamiltonian equation

$\displaystyle\partial_{t}\phi^{a}=\frac{\delta{\cal
H}}{\delta\pi^{a}},\quad\partial_{t}\pi^{a}=-\frac{\delta{\cal
H}}{\delta\phi^{a}}.$ (5)

Note that

$\frac{d{\cal
H}}{dt}=\int_{D}d^{2}{\boldsymbol{r}}~{}\sum_{a=1}^{q-1}{\boldsymbol{\nabla}}\left(\pi^{a}{\boldsymbol{\nabla}}\phi^{a}\right).$
(6)

For conserving energy, we assume the Neumann boundary condition in the
$x$-direction:

$\left.\partial_{x}\phi\right|_{x=0,L_{x}}=0$ (7)

for any $y$, and periodic boundary conditions in the $y$-direction.

Figure 2:  (a), (b) Snapshot of the order parameter density
$\hat{m}(\boldsymbol{r})$ in equilibrium with $J=0$ (panel (a)) and the steady
state with $J=-0.00002$ (panel (b)) for $E/(L_{x}L_{y})=0.88$. (c) One-
dimensional order parameter profiles $m^{\mathrm{mc}}(x)$ and
$m^{\mathrm{ss}}(x)$. The interface positions $X^{\mathrm{mc}}$ and
$X^{\mathrm{ss}}$ are obtained as $m^{\mathrm{mc}}(X^{\mathrm{mc}})=0.5$ and
$m^{\mathrm{ss}}(X^{\mathrm{ss}})=0.5$, respectively. The interface thickness
$\xi^{\mathrm{ss}}$ is estimated by using $m^{\rm
ss}(x)=a-b\tanh((x-x_{0})/\xi^{\mathrm{ss}})$. Inset: Dependence of
$X^{\mathrm{mc}}$ and $X^{\mathrm{ss}}$ on the energy density
$E/(L_{x}L_{y})$.

We start from initial conditions with a single interface and solve the
Hamiltonian equation (5) until an equilibrium state is realized [35]. Figure
2(a) shows a snapshot of the local order parameter $\hat{m}(\boldsymbol{r})$
in equilibrium. It shows an interface parallel to the $y$-axis. We thus
discuss the one-dimensional profile by using the average over $y$ defined by
(4).

Let $\langle\hat{A}\rangle_{E}$ be the long-term average of
$\hat{A}({\boldsymbol{r}})$ after equilibration, which corresponds to the
expected value with respect to the micro-canonical distribution with $E$. We
show the order parameter profile
$m^{\mathrm{mc}}(x)=\langle[\hat{m}]_{x}\rangle_{E}$ in Fig. 2(c), which
exhibits a typical interface structure. Here, the superscript “mc” indicates
“micro-canonical”. The interface position $X^{\mathrm{mc}}$ defined as
$m^{\mathrm{mc}}(X^{\mathrm{mc}})=0.5$ decreases approximately linearly with
$E$ such that $X^{\mathrm{mc}}=L_{x}$ for $E=E_{1}$ and $X^{\mathrm{mc}}=0$
for $E=E_{2}$. There is no interface for $E<E_{1}$ and $E>E_{2}$.

The one-dimensional temperature profile
$T^{\mathrm{mc}}(x)=\langle[\hat{T}]_{x}\rangle_{E}$ is homogeneous in $x$
even for the phase coexistence observed in $E_{1}<E<E_{2}$. This homogeneous
temperature is equal to $T_{\mathrm{c}}$ obtained in isothermal systems [35].
All these results show that the phase coexistence observed as an equilibrium
state for the case $E\in[E_{1},E_{2}]$ is an important feature of the isolated
system and that the behavior is equivalent to the discontinuous change
observed in the isothermal system, as displayed in the phase diagram in Fig. 1
(b).

Non-equilibrium phase coexistence.— We now consider phase coexistence in heat
conduction. Recall that the position of the interface is uniquely determined
in isolated systems for given $E\in[E_{1},E_{2}]$, whereas it is neutral and
not under control in isothermal systems at $T=T_{\mathrm{c}}$. Even for heat
conduction, whether the energy conserves or not influences the fluctuation
properties in the steady state. Thus, energy-conserving heat-conduction could
be preferable to standard heat conduction for a detailed study of the steady
states. From the thermodynamic equivalence for heat conducting states as well
as equilibrium states [26, 35], we expect that the obtained results for the
energy-conserving heat-conduction systems will explain the phase coexistence
observed in standard heat-conducting systems.

We construct a heat-conducting system where the energy flows in at $x=L_{x}$
and flows out at $x=0$, while keeping the energy ${\cal H}$ constant. That is,
we impose a constant heat flux $(JL_{y},0)$ at $x=0$ and $(JL_{y},0)$ at
$x=L_{x}$. Recalling (6), we set

$\displaystyle\left.\partial_{x}\phi^{a}\right|_{x=0,L_{x}}=-\frac{\left.JL_{y}\pi^{a}\right|_{x=0,L_{x}}}{\int_{0}^{L_{y}}dy\left.\sum_{b=1}^{q-1}{(\pi^{b})^{2}}\right|_{x=0,L_{x}}},$
(8)

which is a non-equilibrium generalization of (7).

Figure 3:  One-dimensional temperature profile $T^{\mathrm{ss}}(x)$ for
$J=-0.00002$ and $E/(L_{x}L_{y})=0.88$. The red solid line represents the
interface position $X^{\mathrm{ss}}$. The orange dashed line represents the
position $x_{\ast}$ for $T^{\mathrm{ss}}(x_{\ast})$ being the equilibrium
transition temperature $T_{\mathrm{c}}$. We note that the position $x_{\ast}$
is not at the interface region, i.e.,
$x_{\ast}<X^{\mathrm{ss}}-\xi^{\mathrm{ss}}\simeq 0.444L$.

We study the case $J=-0.00002$ [36]. Figure 2(b) is a snapshot of the order-
parameter density field $\hat{m}(\boldsymbol{r})$ in the steady state, which
can be compared with the equilibrium case shown in Fig. 2(a). The interface in
the steady state is shifted slightly to the left of the equilibrium interface.
In Fig. 2(c), we display the order-parameter profile
$m^{\mathrm{ss}}(x)=\langle[\hat{m}]_{x}\rangle_{E,J}$, where
$\langle\cdot\rangle_{E,J}$ represents the long-term average in the steady
state with $E$ and $J$. For the interface position $X^{\mathrm{ss}}$ defined
by $m^{\mathrm{ss}}(X^{\mathrm{ss}})=0.5$, we find that
$|X^{\mathrm{ss}}-X^{\mathrm{mc}}|/L\simeq 0.05$. The inset plots
$X^{\mathrm{ss}}$ for an energy density $E/(L_{x}L_{y})$. The relative
difference $|X^{\mathrm{ss}}-X^{\mathrm{mc}}|/L$ becomes smaller as $E$
approaches $E_{1}$ or $E_{2}$.

In Fig. 3, we plot the temperature profile
$T^{\mathrm{ss}}(x)=\langle[\hat{T}]_{x}\rangle_{E,J}$. An important feature
is that the position $x_{\ast}$ satisfying
$T^{\mathrm{ss}}(x_{\ast})=T_{\mathrm{c}}$ is inside the ordered regions,
i.e., $x_{\ast}<X^{\mathrm{ss}}-\xi^{\mathrm{ss}}$, where $\xi^{\mathrm{ss}}$
is the interface thickness as indicated in Fig. 2(c). That is, super-heated
ordered states appear in a region where $T(x)>T_{\mathrm{c}}$ and $m(x)>0.5$.
This means that the local equilibrium thermodynamics is violated due to the
heat flux. In other words, the temperature at the interface $\theta\equiv
T^{\mathrm{ss}}(X^{\mathrm{ss}})$ deviates from $T_{\mathrm{c}}$.

To study the violation more quantitatively, in Fig. 4(a), we plot $\theta$ as
a function of $E$ with $J=-0.00002$ fixed. It is notable that the interface
temperature $\theta$ deviates from the equilibrium transition temperature
$T_{\mathrm{c}}$ for $E\in[E_{1},E_{2}]$. The deviation becomes a maximum
around the mid point of $[E_{1},E_{2}]$, where the interface position is
around $L_{x}/2$. The super-heated regions disappear for $E\rightarrow E_{1}$
or $E_{2}$, where the system becomes occupied by an ordered or disordered
state. To examine the parameter dependence, we compare the numerical results
for the violation with the formula

$\displaystyle\begin{aligned} \theta^{\rm Th}=T_{\rm
c}+{|J|}\left(\frac{1}{\kappa^{\rm o}}-\frac{1}{\kappa^{\rm
d}}\right)\frac{X(L_{x}-X)}{2L_{x}}\end{aligned}$ (9)

which was proposed in the global thermodynamics framework [25], where
$\kappa^{\rm o}$ and $\kappa^{\rm d}$ are the heat conductivities of ordered
and disordered states, respectively, and $X$ is the position of the interface.
The formula (9) can be derived for the order-disorder transition [35]. In Fig.
4(a), we simultaneously plot $\theta$ and $\theta^{\rm Th}$, where the latter
is obtained by substituting $X^{\mathrm{mc}}$ into $X$ in (9). We find that
the interface temperature $\theta$ directly measured in numerical simulations
is in good agreement with $\theta^{\rm Th}$ determined by the formula (9).

Finally, we investigate the finite size effects. We find that the violation of
local equilibrium becomes weaker for smaller $L_{x}$ or larger $\Delta x$.
Instead of the original system consisting of $3073\times 512$ grid points with
$\Delta x=1/8$, we study smaller systems consisting of $1537\times 256$ grid
points and $768\times 128$ grid points, keeping the value of $\Delta x$ fixed.
The results are shown in Fig. 4 (b). The interface temperature deviates from
the theoretical curve and eventually reaches the equilibrium transition
temperature. This indicates that long-wavelength fluctuations play an
important role in the violation of local equilibrium properties. A similar
trend is observed for rough systems with larger values of $\Delta x$ for the
same system size. The interface temperature becomes closer to $T_{\mathrm{c}}$
as $\Delta x$ increases from $1/8$ to $1/2$. This implies that regularity of
the short wavelength fluctuations of the continuum fields is also necessary
for the violation, which is in stark contrast to the properties in
equilibrium. We do not observe the violation of local equilibrium with smaller
or rougher systems. These observations suggest that hydrodynamic fluctuations
on some scales, which generate long-range correlations, play an important role
for the violation of the local equilibrium at the interface.

Figure 4:  (a) Interface temperature as a function of the energy density
$E/(L_{x}L_{y})$ for $J=-0.00002$. The green points are numerical values
$\theta\equiv T^{\mathrm{ss}}(X^{\mathrm{ss}})$ and the red points are the
theoretical predictions $\theta^{\rm Th}$ in (9), substituting the numerically
determined $X^{\mathrm{mc}}$ for $X$. (b) Finite size effects of the
numerically determined $\theta$ for $(L_{x},L_{y},\Delta x)$ depicted in the
box with $J=-0.00002$. See text.

Concluding remarks.— In this Letter, we studied non-equilibrium phase
coexistence under a heat flux using the two-dimensional Hamiltonian Potts
model with $q=11$. We found that a super-heated ordered region stably appears
in the heat-conducting state. This indicates that the metastable states are
controlled by the heat flux. Our numerical results also quantitatively support
the validity of global thermodynamics. We hope this work triggers further
experimental and theoretical research, such as engineering thermodynamic
metastable states. In conclusion, we present three remarks.

Although we expect that the violation of local equilibrium thermodynamics is
generically observed in other models and experiments, we need to carefully
choose the system conditions. For example, liquid–gas coexistence under heat
conduction was studied using molecular dynamics simulations, and no violation
was observed [37, 38]. We conjecture that the system sizes were too small to
allow long-wavelength fluctuations that would lead to the violation of local
equilibrium thermodynamics. Regarding this point, we remark that long-range
correlation of hydrodynamic fluctuations in a sheared system is observed only
when the particle number exceeds a cross-over value of $10^{7}$ [39].
Therefore, we expect that the violation of local equilibrium thermodynamics
would be observed in large systems such that the long-range correlations are
observed as well.

From a theoretical viewpoint, an important first study is to derive the
stationary distribution of the heat conduction system. In the linear response
regime, the Zubarev–McLennan representation is a generalization of the micro-
canonical distribution, where the correction term is expressed in terms of
entropy production [3, 4], as discussed for phase coexistence in heat
conduction [40]. Developing a method for estimating this correction term
explicitly for microscopic or mesoscopic models is an important goal.

The most critical challenge is the observation of the violation of local
equilibrium thermodynamics in laboratory experiments. As an example, for the
liquid–gas coexistence for water at 1 atm pressure, the interface temperature
was calculated to be 96 ∘C when the temperatures of the heat baths attached to
the left and right sides were 95 ∘C and 105 ∘C [25]. This means that super-
cooled gas stably appears near the interface due to the heat flux. We believe
that observation of this phenomenon is of fundamental importance.

The work of M.K. was supported in part by the Osaka City University Advanced
Mathematical Institute (MEXT Joint Usage/Research Center on Mathematics and
Theoretical Physics, Grant No. JPMXP0619217849). The present study was
supported by JSPS KAKENHI Grant Numbers JP19KK0066, JP20K03765,
JP19K03647,JP19H05795, JP20K20425, and JP22H01144.

## References

  * [1] L. D. Landau and E. M. Lifshitz, Fluid Mechanics, (Pergamon Press, Oxford 1959).
  * [2] S. R. de Groot and P. Mazur, Non-equilibrium Thermodynamics, (Dover, New York, 1984).
  * [3] D. N. Zubarev, Nonequilibrium Statistical Thermodynamics, (Consultants Bureau, New York, 1974).
  * [4] J. A. Mclennan, Phys. Fluids 3, 493 (1960); Introduction to Non-equilibrium Statistical Mechanics (Prentice-Hall, 1988).
  * [5] S.-i. Sasa, Derivation of hydrodynamics from the Hamiltonian description of particle systems, Phys. Rev. Lett. 112, 100602 (2014).
  * [6] T. Hayata, Y. Hidaka, T. Noumi, and M. Hongo, Relativistic hydrodynamics from quantum field theory on the basis of the generalized Gibbs ensemble method, Phys. Rev. D 92, 065008 (2015).
  * [7] A. Onuki and K. Kawasaki, Nonequilibrium steady state of critical fluids under shear flow: A renormalization group approach, Annals of Physics 121, 456-528 (1979).
  * [8] A. Onuki, Phase Transition Dynamics, (Cambridge University Press, Cambridge, 2004).
  * [9] K. Binder, Theory of first-order phase transitions, Rep. Prog. Phys. 50 783 (1987).
  * [10] F. S. Bates, J. H. Rosedale, G. H. Fredrickson, and C. J. Glinka, Fluctuation-induced first-order transition of an isotropic system to a periodic State, Phys. Rev. Lett. 61, 2229-2232 (1988).
  * [11] L. Jin, Y. Shi, F. I. Allen, L.-Q. Chen, and J. Wu, Probing the critical nucleus size in the metal-insulator phase transition of CO2, Phys. Rev. Lett. 129, 245701 (2022).
  * [12] Q. Yang, R. Bi, D. Banerjee, and H. Nasrabadi, Direct observation of the vapor-liquid phase transition and hysteresis in 2 nm Nanochannels, Langmuir 38, 9790-9798 (2022).
  * [13] G. S. Agarwal and S. R. Shenoy, Observability of hysteresis in first-order equilibrium and nonequilibrium phase transitions, Phys. Rev. A 23 2719-2723 (1981).
  * [14] G. L. Wang, Z. K. Qin, and D. L. Lin, First-order phase transition in order-disorder ferroelectrics, Phys. Rev. B 40, 680-684 (1989).
  * [15] T. Riste and L. Dobrzynski, Nematic-isotropic transition: Thermal hysteresis and magnetic field effects, Phys. Rev. Lett. 74, 2737 (1995).
  * [16] P. D. Olmsted and P. Goldbard, Theory of the nonequilibrium phase transition for nematic liquid crystals under shear flow, Phys. Rev. A 41, 4578-4581 (1990).
  * [17] J. R. Thome, Boiling in microchannels: a review of experiment and theory, Int. J. Heat and Fluid flow 25, 128-139 (2004).
  * [18] E. Ben-Jacob and P. Garik, The formation of patterns in nonequilibrium growth, Nature 343, 523-530 (1990).
  * [19] 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).
  * [20] J.-Q. Zhong, D. Funfschilling, and G. Ahlers, Enhanced heat transport by turbulent two-phase Rayleigh-Benard convection, Phys. Rev. Lett. 102, 124501 (2009).
  * [21] S. Weiss and G. Ahlers, Nematic-isotropic phase transition in turbulent thermal convection, J. Fluid Mech. 737, 308-328 (2013).
  * [22] 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).
  * [23] F. Y. Wu, The Potts model, Rev. Mod. Phys. 54 235 (1982).
  * [24] N. Nakagawa and S.-i. Sasa, Liquid-gas transitions in steady heat conduction, Phys. Rev. Lett. 119, 260602 (2017).
  * [25] N. Nakagawa and S.-i. Sasa, Global thermodynamics for heat conduction states, J. Stat. Phys. 177, 825-888 (2019).
  * [26] 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).
  * [27] J. Keizer, Thermodynamics at nonequilibrium steady states, J. Chem. Phys. 69, 2609 (1978).
  * [28] B. C. Eu, Irreversible thermodynamics of fluids, Annals of Physics 140, 341-371 (1982).
  * [29] D. Jou, J. Casas-Vázquez, and G. Lebon, Extended irreversible thermodynamics, Rep. Prog. Phys. 51, 1105-1179 (1988).
  * [30] Y. Oono and M. Paniconi, Steady state thermodynamics, Prog. Theor. Phys. Suppl. 130, 29 (1998).
  * [31] S.-i. Sasa and H. Tasaki, Steady state thermodynamics, J. Stat. Phys. 125, 125–224 (2006).
  * [32] E. Bertin, K. Martens, O. Dauchot, and M. Droz, Intensive thermodynamic parameters in nonequilibrium systems, Phys. Rev. E 75, 031120 (2007).
  * [33] P. Pradhan, R. Ramsperger, and U. Seifert, Approximate thermodynamic structure for driven lattice gases in contact, Phys. Rev. E 84, 041104 (2011).
  * [34] R. Dickman, Failure of steady-state thermodynamics in nonuniform driven lattice gases, Phys. Rev. E 90, 062123 (2014).
  * [35] See Supplemental Material at [URL will be inserted by publisher] for derivations and detailed explanations.
  * [36] This value was chosen so that the system is in the linear response regime in $J$. For example, linear temperature profiles are observed in the ordered and disordered phases [35].
  * [37] A. Røsjorde, D. W. Fossmo, D. Bedeaux, S. Kjelstrup, and B. Hafskjold, Nonequilibrium molecular dynamics simulations of steady-state heat and mass transport in condensation: I. Local equilibrium, J. Coll. Int. Sci. 232, 178-185 (2000).
  * [38] S. Yukawa, private communication
  * [39] H. Nakano and Y. Minami, Molecular dynamics study of shear-induced long-range correlations in simple fluids, Phys. Rev. Res. 4, 023147 (2022).
  * [40] 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).

Supplemental Material for
“Control of metastable states by heat flux in the Hamiltonian Potts model”
Michikazu Kobayashi1, Naoko Nakagawa2, and Shin-ichi Sasa3

1School of Environmental Science and Engineering, Kochi University of
Technology, Miyanoguchi 185, Tosayamada, Kami, Kochi 782-8502, Japan

2Department of Physics, Ibaraki University, Mito 310-8512, Japan

3Department of Physics, Kyoto University, Kyoto, 606-8502 Japan

This Supplemental Material consists of three sections. In Sec. I, we review
global thermodynamics for the system we study. The goal of this section is to
derive the formula (9) in the main text. In Sec. II, we present details of
numerical simulations of the Hamiltonian equation (5) to support the results
of the main text. In Sec. III, we explain the isothermal dynamics by which the
equilibrium phase diagram for the temperature is determined.

##  I Global thermodynamics

First of all, we introduce an external field $H\geq 0$ for formulating
thermodynamics of ordered and disordered states. While the numerical
computation in the main text was performed without imposing any external
field, we study a spontaneous symmetry-breaking state by considering the limit
$H\rightarrow 0^{+}$.

To set up the system, we notice the order parameter field
$\hat{m}({\boldsymbol{r}})$ defined as (3) and introduce a uniform field $H$
acting on $\hat{m}({\boldsymbol{r}})$ by assuming the Hamiltonian of the
system as

${\cal H}(\phi,\pi)-H{\hat{M}}(\phi).$ (S.1)

The Hamiltonian $\cal H$ without the external field is defined in (1) and the
total order parameter $\hat{M}$ is defined as the spatial integration

$\hat{M}=\int d^{2}{\boldsymbol{r}}~{}\hat{m}({\boldsymbol{r}}).$ (S.2)

We assume that thermodynamic properties are not qualitatively changed by
applying the external field for sufficiently small $H$. For instance, the
equilibrium transition temperature may continuously change with $H$, which we
write $T_{\mathrm{c}}(H)$. Below $T_{\mathrm{c}}(H)$ the ordered phase appears
as studied in the main text. From the symmetry,
$T_{\mathrm{c}}(H)=T_{\mathrm{c}}(0)+O(H^{2})$.

In Sec. I.1, we briefly review equilibrium thermodynamics. In Sec. I.2, we
formulate global thermodynamics for the bulk of ordered and disordered states
in steady heat conduction. In Sec. I.3, we present the variational principle
for determining the steady phase coexistence state in heat conduction. Section
I.4 is devoted to the derivation of (9) used in the main text. In Sec. I.5, we
explain the thermodynamic relations in heat conduction and the equivalence of
steady states of the energy-conserved system with the system conducted by two
heat baths.

###  I.1 Equilibrium thermodynamics

We review equilibrium thermodynamics for the model we study in the main text.

In the thermodynamic limit, $\cal H$ and $\hat{M}$ form the internal energy
$U$ and macroscopic order parameter $M$. Letting $E$ be the total energy of
the system, (S.1) leads to $E=U-HM$. Changing the value of $H$ to $H+dH$
requires work $d^{\prime}W=-MdH$. Since $d^{\prime}Q=TdS$, the energy
conservation forms the first law of thermodynamics as

$\displaystyle dE=TdS-MdH.$ (S.3)

Thus, the entropy $S$ is given as a function of $(E,H)$ for a fixed volume of
the system $V=L_{x}L_{y}$ with the relation

$\displaystyle dS=\frac{dE}{T}+\frac{M}{T}dH.$ (S.4)

Since $S(E,H,V)$ satisfies the extensivity

$\displaystyle S(\lambda E,H,\lambda V)=\lambda S(E,H,V)$ (S.5)

for any $\lambda>0$, we can derive

$\displaystyle dS=\frac{dE}{T}+\frac{M}{T}dH-\frac{f}{T}dV$ (S.6)

with $f=(E-TS)/V$ which corresponds to the free-energy density. Since we study
the case $H=0$ in the main text, the internal energy $U$ is equal to the total
energy $E$.

###  I.2 Heat conduction without phase coexistence

Next, we study the linear response regime in heat conduction systems without
phase coexistence. Specifically, we consider the case that the heat flux is in
parallel to the $x$-axis for simplicity. We then assume the existence of local
thermodynamic quantities $a({\boldsymbol{r}})$, which are independent of $y$.
We express them as $a(x)$ simply. The arguments below follow Sec. 3 of ref.
[1].

Let $s(x)$, $T(x)$, $e(x)$, and $m(x)$ be local entropy density, local
temperature, local energy density, and local order parameter density. These
local quantities satisfy the fundamental relation,

$\displaystyle ds(x)=\frac{de(x)}{T(x)}+\frac{{m}(x)}{T(x)}d{H},$ (S.7)

for each $x$. The total amounts of extensive quantities in heat conduction are
obtained by the integration of corresponding local quantities such that

$\displaystyle S=L_{y}\int_{0}^{L_{x}}dx~{}s(x),\quad
E=L_{y}\int_{0}^{L_{x}}dx~{}e(x),\quad M=L_{y}\int_{0}^{L_{x}}dx~{}m(x).$
(S.8)

Hereafter, we call these “global thermodynamic quantities.” We then define a
global temperature and a global free energy density to characterize the heat
conduction states as the mean value of the respective local values as

$\displaystyle\tilde{T}=\frac{1}{L_{x}}\int_{0}^{L_{x}}dx~{}T(x),$ (S.9)
$\displaystyle\tilde{f}=\frac{1}{L_{x}}\int_{0}^{L_{x}}dx~{}f(x).$ (S.10)

In the linear response regime, the quantities defined by the spatial
integration such as (S.8), (S.9) and (S.10) are estimated by the trapezoidal
rule. Thus the global quantities are represented by corresponding local
thermodynamic quantities at the mid-point $x=L_{x}/2$. This indicates that the
non-equilibrium entropy in (S.8) is approximated by the equilibrium entropy
function; i.e.,

$\displaystyle S=S(E,H,V)+O(\varepsilon^{2}),$ (S.11)

where $\varepsilon$ is a dimensionless small parameter proportional to $J$.
Straightforwardly, the global thermodynamic quantities in heat conduction
satisfy the fundamental relation of thermodynamics

$\displaystyle
dS=\frac{dE}{\tilde{T}}+\frac{M}{\tilde{T}}dH-\frac{\tilde{f}}{\tilde{T}}dV,$
(S.12)

with an error of $O(\varepsilon^{2})$. While global quantities $E$, $M$, and
$V$ are naturally defined for the heat conduction states, $\tilde{T}$ and
$\tilde{f}$ are given as non-equilibrium extensions of the temperature and the
free energy density through the thermodynamic relation (S.12).

###  I.3 Phase coexistence in heat conduction

Let us concentrate on the energy-conserving heat-conduction system. Choosing
the value of $E$ in a certain range, say $[E_{1},E_{2}]$, we observe the phase
coexistence with a planner interface. In equilibrium, the temperature at the
phase coexistence is kept at the transition temperature $T_{\mathrm{c}}(H)$.
Now, we impose the energy flux $J$ in such a way that the total energy $E$ is
conserved, where $J<0$ is assumed without loss of generality. We then observe
the phase coexistence of the ordered state in $0\leq x\leq X$ and the
disordered state in $X\leq x\leq L_{x}$, where $X$ is to the interface
position. The temperature field $T(x)$ and the order parameter field $m(x)$
satisfy the local thermodynamic relation in each phase, while there is no
established law for determining the interface temperature $\theta$. Here, let
us recall that the equilibrium state in the phase coexistence for given
$(E,H,V)$ is determined by the maximum entropy principle, where the
equilibrium state maximizes the total entropy. We extend this variational
principle to that for heat conduction states.

We start with global thermodynamics in each region introduced in Sec. I.2. The
global entropies are expressed as

$\displaystyle S^{\rm o}=S(E^{\rm o},H,L_{y}X),\quad S^{\rm d}=S(E^{\rm
d},H,L_{y}(L_{x}-X)),$ (S.13)

where we put the superscript “o” or “d” to indicate the quantities in the
ordered or disordered phase. From the fundamental relation (S.12), we have

$\displaystyle\left(\frac{\partial S^{\rm o}}{\partial E^{\rm
o}}\right)_{H,X}=\frac{1}{\tilde{T}^{\rm o}},\quad\left(\frac{\partial S^{\rm
d}}{\partial E^{\rm d}}\right)_{H,L-X}=\frac{1}{\tilde{T}^{\rm d}},$ (S.14)
$\displaystyle\frac{1}{L_{y}}\left(\frac{\partial S^{\rm o}}{\partial
X}\right)_{E^{\rm o},H}=-\frac{\tilde{f}^{\rm o}}{\tilde{T}^{\rm
o}},\quad\frac{1}{L_{y}}\left(\frac{\partial S^{\rm
d}}{\partial(L_{x}-X)}\right)_{E^{\rm d},H}=-\frac{\tilde{f}^{\rm
d}}{\tilde{T}^{\rm d}}.$ (S.15)

Letting $T_{1}=T(0)$, $T_{2}=T(L)$, and $\theta=T(X)$, the trapezoidal rule
yields

$\displaystyle\tilde{T}^{\rm o}=\frac{T_{1}+\theta}{2},\quad\tilde{T}^{\rm
d}=\frac{T_{2}+\theta}{2}.$ (S.16)

Note that $\theta$ is the interface temperature. $T_{1}$ and $T_{2}$ are
temperatures at the left and right boundaries of the system, respectively.

For determining steady states in phase coexistence, thermodynamics provides
variational principles. For the equilibrium case $J=0$, the equilibrium state
with $E$, $H$, and $V$ fixed maximizes the total entropy $S=S^{\rm o}+S^{\rm
d}$ with respect to $E^{\rm o}$ and $X$. Here, $E$ is the total amount of
energy of the system and $E^{\rm d}=E-E^{\rm o}$. The variational equations
lead to $\tilde{T}^{\rm o}=\tilde{T}^{\rm d}$ and $\tilde{f}^{\rm
o}=\tilde{f}^{\rm d}$. For heat conduction cases $J<0$, as derived in [2], the
entropy as the variational function is concluded as

$\displaystyle S=S^{\rm o}+S^{\rm d}+\phi\Psi$ (S.17)

with

$\displaystyle\phi=\frac{T_{2}-T_{1}}{T_{\mathrm{c}}(H)},$ (S.18)
$\displaystyle\Psi=\frac{E^{\rm d}X-E^{\rm
o}(L_{x}-X)}{2T_{\mathrm{c}}(H)L_{x}}.$ (S.19)

The steady state with $E$, $H$, $L_{x}$, and $\phi$ fixed is determined as the
state maximizing (S.17) with respect to $E^{\rm o}$ and $X$. We then obtain
the condition for the steady state as

$\displaystyle\tilde{T}^{\rm d}-\tilde{T}^{\rm o}=\frac{\phi
T_{\mathrm{c}}(H)}{2},$ (S.20) $\displaystyle\frac{\tilde{f}^{\rm
o}}{\tilde{T}^{\rm o}}-\frac{\tilde{f}^{\rm d}}{\tilde{T}^{\rm
d}}=\phi\frac{E}{2T_{\mathrm{c}}(H)L_{x}L_{y}}.$ (S.21)

Note that these two relations are the extension of the equilibrium balance
relations $\tilde{T}^{\rm o}=\tilde{T}^{\rm d}$ and $\tilde{f}^{\rm
o}=\tilde{f}^{\rm d}$. Due to the non-equilibrium contributions in the right-
hand sides of (S.20) and (S.21), we find that the steady states contain
metastable states in a region adjacent to the interface. See (9) and its
derivation in Sec. I.4.

###  I.4 Derivation of (9)

For simplicity of the derivation below, we introduce $\alpha\equiv-f/T$. Since
$\alpha=s-e/T$, the equilibrium fundamental relation (S.7) yields

$\displaystyle d\alpha=\frac{e}{T^{2}}dT+\frac{m}{T}dH.$ (S.22)

We also notice $E=E^{\rm o}+E^{\rm d}$, $E^{\rm o}=e_{\rm c}^{\rm
o}XL_{y}+O(\varepsilon)$, and $E^{\rm d}=e_{\rm c}^{\rm
d}(L_{x}-X)L_{y}+O(\varepsilon)$ with $e_{\rm c}^{\rm o/d}\equiv e^{\rm
o/d}(T_{\rm c}(H),H)$. Using $\alpha$ and these relations, (S.21) is written
as

$\displaystyle\tilde{\alpha}^{\rm d}-\tilde{\alpha}^{\rm
o}=\frac{\phi}{2T_{\mathrm{c}}(H)}\left(e_{\rm c}^{\rm
o}\frac{X}{L_{x}}+e_{\rm c}^{\rm
d}\frac{L_{x}-X}{L_{x}}\right)+O(\varepsilon^{2}).$ (S.23)

Since $\tilde{\alpha}^{\rm o/d}=\alpha^{\rm o/d}(\tilde{T}^{\rm o/d},H)$ with
the equilibrium functions $\alpha^{\rm o/d}(T,H)$, we have

$\displaystyle\tilde{\alpha}^{\rm d}-\tilde{\alpha}^{\rm o}=\alpha^{\rm
d}(T_{\mathrm{c}}(H),H)+\frac{e_{\rm c}^{\rm
d}}{T_{\mathrm{c}}(H)^{2}}(\tilde{T}^{\rm d}-T_{\mathrm{c}}(H))-\alpha^{\rm
o}(T_{\mathrm{c}}(H),H)-\frac{e_{\rm c}^{\rm
o}}{T_{\mathrm{c}}(H)^{2}}(\tilde{T}^{\rm
o}-T_{\mathrm{c}}(H))+O(\varepsilon^{2}).$ (S.24)

Substituting this into (S.23) and using the equilibrium balance $\alpha^{\rm
o}(T_{\mathrm{c}}(H),H)=\alpha^{\rm d}(T_{\mathrm{c}}(H),H)$, we obtain

$\displaystyle e_{\rm c}^{\rm o}\left(\tilde{T}^{\rm
o}-T_{\mathrm{c}}(H)+\frac{\phi
T_{\mathrm{c}}(H)}{2}\frac{X}{L_{x}}\right)-e_{\rm c}^{\rm
d}\left(\tilde{T}^{\rm d}-T_{\mathrm{c}}(H)-\frac{\phi
T_{\mathrm{c}}(H)}{2}\frac{L_{x}-X}{L_{x}}\right)=0.$ (S.25)

The global temperature $\tilde{T}$ for the phase coexistence, which is defined
by (S.9), is expressed as

$\displaystyle\tilde{T}=\tilde{T}^{\rm o}\frac{X}{L_{x}}+\tilde{T}^{\rm
d}\frac{L_{x}-X}{L_{x}}.$ (S.26)

Combining this with (S.20), we have

$\displaystyle\tilde{T}^{\rm o}=\tilde{T}-\frac{\phi
T_{\mathrm{c}}(H)}{2}\frac{L_{x}-X}{L_{x}},\quad\tilde{T}^{\rm
d}=\tilde{T}+\frac{\phi T_{\mathrm{c}}(H)}{2}\frac{X}{L_{x}},$ (S.27)

where we have used $\tilde{T}=T_{\mathrm{c}}(H)+O(\varepsilon)$. Substituting
(S.27) into (S.25), we obtain

$\displaystyle(e_{\rm c}^{\rm d}-e_{\rm c}^{\rm
o})\left(\tilde{T}-T_{\mathrm{c}}(H)-\frac{\phi
T_{\mathrm{c}}(H)}{2}\frac{L_{x}-2X}{L_{x}}\right)=0.$ (S.28)

Note that $e_{\rm c}^{\rm d}-e_{\rm c}^{\rm o}\neq 0$ in the first-order
transition. The relation (S.28) together with (S.27) leads to

$\displaystyle\tilde{T}^{\rm o}=T_{\mathrm{c}}(H)-\frac{\phi
T_{\mathrm{c}}(H)}{2}\frac{X}{L_{x}},\quad\tilde{T}^{\rm
d}=T_{\mathrm{c}}(H)+\frac{\phi T_{\mathrm{c}}(H)}{2}\frac{L_{x}-X}{L_{x}}.$
(S.29)

Applying (S.16), the first relation of (S.29) is rewritten as

$\displaystyle T_{1}=2T_{\mathrm{c}}(H)-\theta-\phi
T_{\mathrm{c}}(H)\frac{X}{L_{x}}.$ (S.30)

Finally, let us recall that the heat flux is uniform in $x$ from the
conservation of energy, which is expressed as

$\displaystyle J=-\kappa^{\rm o}\frac{\theta-T_{1}}{X}=-\kappa^{\rm
d}\frac{T_{2}-\theta}{L_{x}-X},$ (S.31)

where $\kappa^{\rm o/d}$ is the heat conductivity at the transition
temperature, $\kappa^{\rm o/d}=\kappa^{\rm o/d}(T_{\mathrm{c}}(H),H)$. (S.31)
is transformed to

$\displaystyle T_{1}=\theta+J\frac{X}{\kappa^{\rm o}},$ (S.32)
$\displaystyle\phi=-\frac{J}{T_{\rm c}}\left(\frac{X}{\kappa^{\rm
o}}+\frac{L_{x}-X}{\kappa^{\rm d}}\right).$ (S.33)

We now substitute these two relations into (S.30). We then obtain

$\displaystyle\theta=T_{\mathrm{c}}(H)-J\left(\frac{1}{\kappa^{\rm
o}}-\frac{1}{\kappa^{\rm d}}\right)\frac{X(L_{x}-X)}{2L_{x}}.$ (S.34)

To this point, we have assumed $J<0$. Repeating the same argument for the case
$J>0$, we find that $\theta$ is symmetric for the transformation $J\to-J$. We
thus conclude (9).

###  I.5 Thermodynamic relations in phase coexistence

Thermodynamic entropy corresponds to the maximum value of the variational
function (S.17). With $\phi$ and $\Psi$ given in (S.18) and (S.19), and
applying the variational equations (S.20) and (S.21), we have the fundamental
relation of global thermodynamics as [2]

$\displaystyle
dS=\frac{dE}{\tilde{T}}+\frac{M}{\tilde{T}}dH-\frac{\tilde{f}}{\tilde{T}}dV+\Psi
d\phi.$ (S.35)

Thus, the global entropy is $S=S(E,H,V,\phi)$, i.e., a different function from
the equilibrium entropy.

As a more familiar setup for heat conduction, one may consider a system
between two heat baths of $T_{1}$ and $T_{2}$. Since the total energy $E$ is
not conserved, natural thermodynamic variables would be $(\tilde{T},H,V,\phi)$
considering a straightforward extension from that for equilibrium systems with
a heat bath. Here the global temperature $\tilde{T}$ is defined by (S.9) and
connected to the global temperatures in the ordered and disordered regions as

$\displaystyle\tilde{T}=\frac{X}{L_{x}}\tilde{T}^{\rm
o}+\frac{L_{x}-X}{L_{x}}\tilde{T}^{\rm d}.$ (S.36)

We can show the equivalence of the steady state in heat conduction between
such systems and the energy-conserved systems investigated in the main text.
More precisely, (9) results from the minimization of a variational function
corresponding to the spatial integration of free energy
$F=L_{y}\int_{0}^{L_{x}}f(x)dx$. It is written as $F=F^{\rm o}+F^{\rm d}$ with
$F^{\rm o}=F(\tilde{T}^{\rm o},H,L_{y}X)$ and $F^{\rm d}=F(\tilde{T}^{\rm
d},H,L_{y}(L_{x}-X))$, and the minimum value of $F$ satisfies a thermodynamic
relation

$\displaystyle dF=-Sd\tilde{T}-MdH+\tilde{f}dV-\tilde{T}\Psi d\phi,$ (S.37)

which indicates the non-equilibrium free energy is written as
$F=F(\tilde{T},H,V,\phi)$. The two fundamental relations (S.35) and (S.37) are
connected by the Legendre transformation with

$\displaystyle F=E-\tilde{T}S+O(\varepsilon^{2}),$ (S.38)

which is led from the spatial integration of the local thermodynamic relation

$\displaystyle f(x)=e(x)-T(x)s(x)$ (S.39)

using $\tilde{T}$ in (S.36) and the entropy in (S.17).

##  II Details of numerical simulations

###  II.1 Explicit form for $\mu_{k}$ of (2)

We explicitly derive coordinates $\mu_{k}^{a}$ ($1\leq k\leq q$, $1\leq a\leq
q-1$) of the $q$-vertex for the regular $(q-1)$-simplex in $\mathbb{R}^{q-1}$.
We first notice that $\mu_{k}^{a}$ should satisfy the following three
conditions:
(i) the centroid is located at the origin:

$\displaystyle\sum_{k=1}^{q}\mu_{k}^{a}=0,$ (S.40)

(ii) the length between the centroid and each vertex is one:

$\displaystyle\sum_{a=1}^{q-1}(\mu_{k}^{a})^{2}=1,$ (S.41)

(iii) the length between two different vertices is constant:

$\displaystyle\sum_{a=1}^{q-1}(\mu_{k}^{a}-\mu_{k^{\prime}}^{a})^{2}=K(1-\delta_{kk^{\prime}})$
(S.42)

for any $1\leq k,k^{\prime}\leq q$, where $K$ is a constant. We then find that
(S.40), (S.41), and (S.42) are satisfied by

$\displaystyle\mu_{k}^{a}=\begin{cases}\displaystyle\frac{\sqrt{q}}{\sqrt{q-1}}\delta_{k}^{a}-\frac{\sqrt{q}+1}{(q-1)^{3/2}},&1\leq
k\leq q-1,\\\ \displaystyle\frac{1}{\sqrt{q-1}},&k=q,\\\ \end{cases}$ (S.43)

with $K=2q/(q-1)$.

###  II.2 Numerical simulations of the Hamiltonian equation (5) with the
boundary condition (7)

We numerically solve the Hamiltonian equation (5). By imposing the boundary
condition (7), the system conserves the total amount of energy $E$. Hereafter,
we consider only the case with odd $q$.

First, we discretize a box of $(L_{x},L_{y})=(384,64)$ with a lattice spacing
length $\Delta x=1/8$ in both $x$ and $y$ directions. Then the system consists
of $N_{x}\times N_{y}$ grid points, i.e., $(N_{x},N_{y})=(3073,512)$. To
calculate the spatial derivative $\nabla^{2}\phi^{a}$ included in ${\delta\cal
H}/\delta\phi^{a}$, we use the following spectral collocation method [3]. We
first construct $(q-1)/2$-complex fields

$\displaystyle\Phi^{a}=\phi^{2a-1}+i\phi^{2a}$ (S.44)

from $(q-1)$-real fields $\phi^{a}$, where $1\leq a\leq(q-1)/2$. We evaluate
the spatial derivative
$\nabla^{2}\Phi^{a}=\nabla^{2}\phi^{2a-1}+i\nabla^{2}\phi^{2a}$ by the Fourier
transformation

$\displaystyle\begin{aligned}
&\nabla^{2}\Phi^{a}(x,y,t)=-\frac{1}{(N_{x}-1)N_{y}}\sum_{n_{x}=0}^{N_{x}-1}{}^{\prime}\sum_{n_{y}=-N_{y}/2+1}^{N_{y}/2}(k_{x}^{2}+k_{y}^{2})\tilde{\Phi}^{a}(n_{x},n_{y},t)\cos(k_{x}x)e^{ik_{y}y},\\\
&k_{x}=\frac{\pi n_{x}}{L_{x}},\quad k_{y}=\frac{2\pi
n_{y}}{L_{y}},\end{aligned}$ (S.45)

where we define

$\displaystyle\sum_{n=n_{i}}^{n_{f}}{}^{\prime}f(n)\equiv\frac{f(n_{i})+f(n_{f})}{2}+\sum_{n=n_{i}+1}^{n_{f}-1}f(n).$
(S.46)

The Fourier basis functions in (S.45) are cosine functions in $x$-direction to
satisfy the boundary condition (7). The Fourier component
$\tilde{\Phi}^{a}(n_{x},n_{y},t)$ in (S.45) is obtained by the inverse
transformation

$\displaystyle\begin{aligned}
\tilde{\Phi}^{a}(n_{x},n_{y},t)=2\sum_{I=0}^{N_{x}-1}{}^{\prime}\sum_{J=0}^{N_{y}-1}\Phi^{a}(x,y,t)\cos(k_{x}I\Delta
x)e^{-ik_{y}J\Delta x},\end{aligned}$ (S.47)

where $\Delta x$ is the lattice spacing length. The variables $\Phi$ defined
on the $N_{x}\times N_{y}$ grid points evolve according to the Hamiltonian
equation (5). To obtain the time evolution of them, we adopt the 4th-order
Strang splitting (symplectic integrator) method [4] with the temporal interval
$\Delta t=1/4096$.

We denote the solution of the Hamiltonian equation (5) with the boundary
condition (7) as $(\phi,\pi)_{E}(\boldsymbol{r},t)$, where $E$ is the amount
of the energy (the value of $\mathcal{H}$), which is determined by the initial
condition and conserved in the time evolution. Aiming at observing the phase
coexistence of the ordered and disordered states, we set the initial condition
as a connected configuration of two equilibrium configurations prepared at
different temperatures:

$\displaystyle(\phi,\pi)_{E}(x,y,0)=\begin{cases}(\phi_{1},\pi_{1})_{T_{1}}(x,y,t_{\rm
f})&0\leq x<L_{x}/2,\\\ (\phi_{1},\pi_{1})_{T_{2}}(x-L_{x}/2,y,t_{\rm
f})&L_{x}/2\leq x<L_{x},\\\ (\phi_{1},\pi_{1})_{T_{2}}(0,y,t_{\rm
f})&x=L_{x},\end{cases}$ (S.48)

where $(\phi_{1},\pi_{1})_{T}$ is the solution of (S.71) with the isothermal
condition at $T$, which will be explained in Sec. III. Thus, the energy $E$ is
determined depending on the chosen initial condition with $T_{1}$ and $T_{2}$.
We fix $T_{1}=0.63T_{\mathrm{c}}$, and change $T_{2}$ from
$0.63T_{\mathrm{c}}$ to $2.0T_{\mathrm{c}}$. Accordingly, the energy $E$
changes from $E/(L_{x}L_{y})=0.5$ to $E/(L_{x}L_{y})=1.2$.

Figure S.1:  (a) Time evolution of the averaged order parameter density
$[\hat{m}]_{x}$ according to the Hamiltonian equation (5) with the boundary
condition (7). The initial condition is given by a combination of two
equilibrium snapshots with different temperatures $T_{1}=0.63T_{\mathrm{c}}$
(left region) and $T_{2}=1.1T_{\mathrm{c}}$ (right region). (b) Dependence of
the interface position $\hat{X}$ on time $t$. The dashed lines show the
fitting with $X_{0}(1-\alpha e^{-t/t_{\ast}})$ for $\hat{X}(t)$. The energy
density is set to $E/(L_{x}L_{y})=0.88$.  Figure S.2:  (Orange dots) micro-
canonical temperature $\bar{T}_{E}$ as a function of energy density
$u=E/(L_{x}L_{y})$. (Blue line) canonical energy density $\bar{u}_{T}$ as a
function of temperature $T$. $T_{\rm c}$ is the transition temperature
determined from the canonical simulations.

Figure S.1 (a) shows an example of the time evolution of the one-dimensional
order-parameter profile $[\hat{m}]_{x}$ as a function of $t$. We see that the
interface shifts from the initial position and relax to an equilibrium
position. We expect that the interface dynamics is the slowest in the
equilibration process. Then, we estimate the relaxation time $t_{\ast}$ of the
interface as follows. Let the interface position $\hat{X}$ be the value of $x$
satisfying $[\hat{m}]_{x}=0.5$. When there are several values of $x$
satisfying $[\hat{m}]_{x}=0.5$, we choose $\hat{X}$ as their centroid. In Fig.
S.1 (b), we show the $t$-dependence of $\hat{X}$. We define the relaxation
time $t_{\ast}$ by using a functional form $\hat{X}\simeq X_{0}(1-\alpha
e^{-t/t_{\ast}})$ in the range $0.5\times 10^{6}\leq t\leq 2.0\times 10^{6}$.
We obtain $t_{\ast}\simeq 5.3\times 10^{5}$. Based on the result, we define
the long-term average of $\hat{A}({\boldsymbol{r}},t)$ as

$\displaystyle\left\langle{\hat{A}({\boldsymbol{r}})}\right\rangle_{E}\equiv\frac{1}{t_{\rm
f}-t_{0}}\int_{t_{0}}^{t_{\rm f}}dt\>\hat{A}({\boldsymbol{r}},t),$ (S.49)

where we set $t_{0}=2\times 10^{6}$ and $t_{\rm f}=t_{0}+3.2\times 10^{4}$.
Note that $t_{0}\simeq 4t_{\ast}$, which would be sufficient for equilibration
as shown in Fig. S.1 (b). We thus conjecture that
$\langle\hat{A}({\boldsymbol{r}})\rangle_{E}$ is approximately equal to the
expected value of ${\hat{A}({\boldsymbol{r}})}$ with respect to the micro-
canonical ensemble

$\displaystyle\rho_{E}^{\mathrm{mc}}(\phi,\pi)=\frac{1}{\Sigma}\delta(E-\mathcal{H}(\phi,\pi)),$
(S.50)

where $\Sigma$ is the normalization constant.

In Fig. S.2, we show the numerical results of the micro-canonical temperature
$\bar{T}_{E}$ given by

$\displaystyle\bar{T}_{E}=\frac{1}{L_{x}L_{y}}\int_{D}d^{2}{\boldsymbol{r}}\>\left\langle\hat{T}({\boldsymbol{r}})\right\rangle_{E},\quad\hat{T}({\boldsymbol{r}})=\frac{\sum_{a=1}^{q-1}(\pi^{a}({\boldsymbol{r}}))^{2}}{2(q-1)}$
(S.51)

as a function of the energy density $u=E/(L_{x}L_{y})$. We observe a plateau
in a range $0.66\lesssim u\lesssim 1.09$, which indicates the realization of
the phase coexistence. We also find good agreement of the data
$(u,\bar{T}_{E})$ with another data $(\bar{u}_{T},T)$, where $\bar{u}_{T}$ is
the energy density of an isothermal system with the temperature $T$. See Sec.
III for numerical simulations of the isothermal system.

###  II.3 Numerical simulations of the Hamiltonian equation (5) with the
boundary condition (8)

Figure S.3:  (a) Time evolution of the averaged order parameter density
$[\hat{m}]_{x}$ in the non-equilibrium dynamics under the heat conduction. (d)
Dependence of the interface position $\hat{X}$ on the time $t$. The dashed
lines show the fitting with $X_{0}(1-\alpha e^{-t/t_{\ast}})$ for
$\hat{X}(t)$. The energy density is set to $E/(L_{x}L_{y})=0.88$.

To study the system under the constant heat flux $(JL_{y},0)$, we solve the
Hamiltonian equation (5) with the non-equilibrium boundary condition (8).
Since the evolution equation (5) is unchanged from the previous section, we
apply the same numerical methods, i.e., spectral collocation and $4$-th order
Strang splitting methods. The only difference from Sec. II.2 is the boundary
condition, which is changed from (7) to (8). To apply (8), letting
$B_{x=0,L_{x}}$ be the right-hand side in (8), (S.45) should be modified as

$\displaystyle\nabla^{2}\Phi^{a}(x,y,t)=\frac{B_{L}^{a}-B_{0}^{a}}{L_{x}}-\frac{1}{(N_{x}-1)N_{y}}\sum_{n_{x}=0}^{N_{x}-1}{}^{\prime}\sum_{n_{y}=-N_{y}/2+1}^{N_{y}/2}(k_{x}^{2}+k_{y}^{2})\tilde{\Phi}^{a}(n_{x},n_{y},t)\cos(k_{x}x)e^{ik_{y}y},$
(S.52)

where $B^{a}_{x=0,L_{x}}$ is explicitly expressed as

$\displaystyle
B^{a}_{x=0,L_{x}}\equiv-\frac{\left.JL_{y}\Pi^{a}\right|_{x=0,L_{x}}}{\int_{0}^{L_{y}}dy\left.\sum_{b=1}^{q-1}\Pi^{b\ast}\Pi^{b}\right|_{x=0,L_{x}}},\quad\Pi^{a}\equiv\pi^{2a-1}+i\pi^{2a}.$
(S.53)

Accordingly, the Fourier component $\tilde{\Phi}^{a}(n_{x},n_{y},t)$ in (S.45)
is replaced by

$\displaystyle\begin{aligned}
\tilde{\Phi}^{a}(n_{x},n_{y},t)=2\sum_{I=0}^{N_{x}-1}{}^{\prime}\sum_{J=0}^{N_{y}-1}\left\\{\Phi^{a}(x,y,t)-B_{0}^{a}-\frac{(B_{L}^{a}-B_{0}^{a})(I\Delta
x)^{2}}{2L}\right\\}\cos(k_{x}I\Delta x)e^{-ik_{y}J\Delta x}.\end{aligned}$
(S.54)

We denote the solution of the Hamiltonian equation (5) with the boundary
condition (8) as $(\phi,\pi)_{E,J}(\boldsymbol{r},t)$.

To observe the effect of the heat flux clearly, we set the initial condition
as the final state of the equilibrium simulations performed in the previous
section:
$(\phi,\pi)_{E,J}(\boldsymbol{r},0)=(\phi,\pi)_{E}(\boldsymbol{r},t_{\rm f})$.
Figure S.3 (a) shows an example of the dynamics for the one-dimensional order-
parameter profile $[\hat{m}]_{x}$. We see that the interface slowly shifts
from the initial equilibrium position corresponding to the final position in
Fig. S.1 (a). This clearly shows that the heat flux changes the steady state.
The relaxation dynamics seems similar to the equilibrium one. Then, we
estimate the equilibrating time $t_{\ast}$ of the interface in Fig. S.3 (b)
with using the interface position $\hat{X}$. By fitting $\hat{X}$ in the range
$0\leq t\leq t_{0}=2.0\times 10^{6}$, we obtain $t_{\ast}\simeq 3.2\times
10^{5}$. Thus, we expect that the system reaches the non-equilibrium steady
state before $t_{0}$ and regard the expected value
$\langle\hat{A}({\boldsymbol{r}})\rangle_{E,J}$ in the steady state as

$\displaystyle\left\langle{\hat{A}({\boldsymbol{r}})}\right\rangle_{E,J}\equiv\frac{1}{t_{\rm
f}-t_{0}}\int_{t_{0}}^{t_{\rm f}}dt\>\hat{A}({\boldsymbol{r}},t)$ (S.55)

with setting $t_{\rm f}=t_{0}+1.6\times 10^{4}$.

###  II.4 Measurement of heat conductivity

Figure S.4:  One dimensional temperature profile $T^{\rm ss}(x)$ with
$J=-0.00002$ starting from the initial conditions (S.56) with (a) $T=0.9T_{\rm
c}$ (ordered phase with $u=0.580$) and (b) $T=1.1T_{\rm c}$ (disordered phase
with $u=1.14$).

To numerically evaluate $\theta^{\rm Th}$ in (9), we have to know the values
of heat conductivities $\kappa^{\rm o}$ and $\kappa^{\rm d}$ in ordered and
disordered states. We calculate them by solving the Hamiltonian equation (5)
with the non-equilibrium boundary condition (8). Here, the initial condition
is given by the final state of (S.71), that is,

$\displaystyle(\phi,\pi_{0})_{E,J}(x,y,0)=\begin{cases}(\phi_{1},\pi_{1})_{T}(x,y,t=t_{\rm
f})&0\leq x<L_{x}/2,\\\ (\phi_{2},\pi_{2})_{T}(x-L_{x}/2,y,t=t_{\rm
f})&L_{x}/2\leq x<L_{x},\\\ (\phi_{2},\pi_{2})_{T}(0,y,t=t_{\rm
f})&x=L_{x}.\end{cases}$ (S.56)

Note that the initial condition (S.56) is an equilibrium configuration at
temperature $T$. That is, it is an ordered state when $T<T_{\mathrm{c}}$ and a
disordered state when $T>T_{\mathrm{c}}$. We determine the local temperature
profile $T^{\rm ss}(x)$ by $\langle\hat{[}\hat{T}]_{x}\rangle_{E,J}$ with
(S.55).

For the calculation of $\kappa^{\rm o}$, we choose $T=0.9T_{\mathrm{c}}$ for
the initial configuration (S.56). The corresponding energy density is
$u=0.580$, where the equilibrium system with energy fixed does not show phase
coexistence. See Fig. S.2. As the value of $u$ is sufficiently smaller than
the saturation value $E_{1}/(L_{x}L_{y})$, the local temperature $T^{\rm
ss}(x)$ is lower than $T_{\mathrm{c}}$ everywhere when applying the small
constant heat flux $J=-0.00002$. See Fig. S.4 (a). The local temperature
$T^{\rm ss}(x)$ linearly increases from $x=0$ to $x=L_{x}$, which supports
that the value $J=-0.00002$ used in this work is so small for the system to be
in the linear-response regime. We finally determine the heat conductivity by
$\kappa^{\rm o}=-JL_{x}/(T^{\mathrm{ss}}(L_{x})-T^{\mathrm{ss}}(0))$.

For the calculation of $\kappa^{\rm d}$, we choose $T=1.1T_{\mathrm{c}}$ for
the initial configuration (S.56), where $u=1.14$. As shown in Fig. S.4 (b),
$T^{\rm ss}(x)$ is higher than $T_{\mathrm{c}}$ everywhere in $x$ with
$J=-0.00002$, and therefore the non-equilibrium system is in the disordered
state. The linear temperature profile confirms that the system is in the
linear-response regime. We determine the heat conductivity as $\kappa^{\rm
d}=-JL_{x}/(T^{\mathrm{ss}}_{E,J}(L_{x})-T^{\mathrm{ss}}_{E,J}(0))$.

Finally, we obtain $\kappa^{\rm o}=8.76JL_{x}/T_{\rm c}$ and $\kappa^{\rm
d}=24.8JL_{x}/T_{\rm c}$. Using these values, we plot $\theta^{\rm Th}$ in
Fig. 4 (a).

##  III Determination of the phase diagram in Fig. 1 (b)

In this section, we figure out the calculation for determining the equilibrium
phase diagram in Fig. 1 (b). We concentrate on the system at constant
temperature $T$. The equilibrium state is described by the canonical ensemble

$\displaystyle\rho_{T}^{\mathrm{can}}(\phi,\pi)=\frac{1}{\mathcal{Z}}\exp(-{\cal
H}(\phi,\pi)/T),$ (S.57)

where $\mathcal{Z}$ is the normalization constant satisfying the condition

$\displaystyle\int{\cal D}\phi{\cal
D}\pi\>\rho_{T}^{\mathrm{can}}(\phi,\pi)=1.$ (S.58)

Since it is not easy to obtain the canonical ensemble in numerical simulations
of the model that exhibits a first-order transition, we explain the numerical
scheme we adopted.

###  III.1 Langevin equation and numerical methods

As an evolution equation consistent with the Hamiltonian equation (5), we
adopt the Langevin equation

$\displaystyle\begin{aligned} \partial_{t}\phi^{a}=\frac{\delta{\cal
H}}{\delta\pi^{a}},\quad\partial_{t}\pi^{a}=-\frac{\delta{\cal
H}}{\delta\phi^{a}}-\gamma\pi^{a}+\xi^{a},\end{aligned}$ (S.59)

where $\xi^{a}({\boldsymbol{r}},t)$ is the white Gaussian noise satisfying

$\displaystyle\left\langle\xi^{a}({\boldsymbol{r}},t)\xi^{b}({\boldsymbol{r}}^{\prime},t^{\prime})\right\rangle=2\gamma
T\delta^{ab}\delta({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime})\delta(t-t^{\prime}).$
(S.60)

The canonical ensemble (S.57) is given as the stationary distribution of
(S.59). We set $\gamma=1$ without loss of generality.

We solve the Langevin equation (S.59) in a periodic box with
$(L_{x},L_{y},\Delta x)=(192,64,1/8)$ consisting of $N_{x}\times N_{y}$ grid
points with $(N_{x},N_{y})=(1536,512)$ separated by the lattice spacing
$\Delta x$. To calculate the spatial derivative $\nabla^{2}\phi^{a}$ included
in ${\delta\cal H}/\delta\phi^{a}$, we again use the spectral collocation
method. There is a slight difference from (S.45) because of the difference of
the boundary condition. Concretely, we use

$\displaystyle\begin{aligned}
&\nabla^{2}\Phi^{a}(x,y,t)=-\frac{1}{N_{x}N_{y}}\sum_{n_{x}=-N_{x}/2+1}^{N_{x}/2}\sum_{n_{y}=-N_{y}/2+1}^{N_{y}/2}(k_{x}^{2}+k_{y}^{2})\tilde{\Phi}^{a}(n_{x},n_{y},t)e^{i(k_{x}x+k_{y}y)},\\\
&k_{x}=\frac{2\pi n_{x}}{L_{x}},\quad k_{y}=\frac{2\pi
n_{y}}{L_{y}}.\end{aligned}$ (S.61)

Note that $\sum$ denotes the standard summation which is different from
$\sum^{\prime}$ defined in (S.46). The Fourier component
$\tilde{\Phi}^{a}(n_{x},n_{y},t)$ in (S.61) is obtained by the inverse-Fourier
transformation

$\displaystyle\begin{aligned}
\tilde{\Phi}^{a}(n_{x},n_{y},t)=\sum_{I=0}^{N_{x}-1}\sum_{J=0}^{N_{y}-1}\Phi^{a}(x,y,t)e^{-i(k_{x}I+k_{y}J)\Delta
x}.\end{aligned}$ (S.62)

The $N_{x}\times N_{y}$ grid points evolve according to the Langevin equation
(S.59). To obtain the time evolution of these points, we use the 2nd-order
Strang splitting:

$\displaystyle\begin{aligned} &\phi^{a}(\boldsymbol{r},t+\Delta
t/2)=\phi^{a}(\boldsymbol{r},t)+\frac{\Delta
t}{2}\pi^{a}(\boldsymbol{r},t),\\\ &\pi^{a}(\boldsymbol{r},t+\Delta
t)=\pi^{a}(\boldsymbol{r},t)-\Delta t\left(\frac{\delta{\cal
H}}{\delta\phi^{a}}(\boldsymbol{r},t+\Delta
t/2)+\pi^{a}(\boldsymbol{r},t)\right)+\sqrt{2T(\Delta x)^{2}(\Delta
t)}\Theta^{a}(\boldsymbol{r},t),\\\ &\phi^{a}(\boldsymbol{r},t+\Delta
t)=\phi^{a}(\boldsymbol{r},t+\Delta t/2)+\frac{\Delta
t}{2}\pi^{a}(\boldsymbol{r},t+\Delta t),\end{aligned}$ (S.63)

where $\Delta t=1/4096$ and $\Theta^{a}(\boldsymbol{x},t)$ are the random
numbers generated by the standard normal distribution. Below, we denote the
solution of the Langevin equation (S.59) as
$(\phi,\pi)_{T}(\boldsymbol{r},t)$.

###  III.2 Metastable and equilibrium states

Compared with the Hamiltonian equation (5) investigated in Sec. II.2, the
Langevin equation (S.59) shows a faster relaxation because of no energy
conservation. Since a typical relaxation time of the Langevin equation is
independent of system sizes, it is four orders of magnitude shorter than that
for the Hamiltonian equation where the relaxation time is proportional to
$L_{x}^{2}$. However, near the first-order transition, the system may relax to
the metastable state rather than the equilibrium state. Once the system is
trapped in the metastable state, it takes significant time to get out of it.
Since the observed state depends on the initial condition, the system shows
hysteresis.

Throughout the following sections, we fix the transient time as
$t_{0}=1.0\times 10^{3}$, and observe states from $t_{0}$ to $t_{\rm
f}=t_{0}+1.6\times 10^{4}$. The time $t_{0}$ is much longer than the
relaxation time to the equilibrium or metastable states, however, much shorter
than the first passage time to get out of the metastable state. The
observation period $t_{\rm f}-t_{0}$ is so long that the finite time
fluctuations can be ignored. In such a time scale, we cannot distinguish the
metastable state from the equilibrium stats. Thus, we first examine the
hysteresis and then proceed to the introduction of a numerical scheme for
providing the canonical ensemble. Our aim is to determine the equilibrium
energy density of the system

$\displaystyle\bar{u}_{T}=\frac{1}{L_{x}L_{y}}\int{\cal D}\phi{\cal
D}\pi\>\rho_{T}^{\mathrm{can}}(\phi,\pi){\cal H}(\phi,\pi),$ (S.64)

and the order parameter of the system

$\displaystyle\bar{m}_{T}=\frac{1}{L_{x}L_{y}}\int_{D}d^{2}{\boldsymbol{r}}\>\int{\cal
D}\phi{\cal
D}\pi\>\rho_{T}^{\mathrm{can}}(\phi,\pi){\hat{m}}({\boldsymbol{r}}),$ (S.65)

which provides a phase diagram.

Below, we examine the equilibrium and metastable states over a wide range of
$T$ in $[T_{\rm min},T_{\rm max}]$, where $T_{\min}=0.0015$ and
$T_{\max}=0.3$. Precisely, we examine $201$ temperature points,

$\displaystyle T_{k}=T_{\rm min}+k\Delta T,\quad(k=0,1,\cdots,200),$ (S.66)

where $\Delta T\equiv(T_{\rm max}-T_{\rm min})/200$. Using the numerical data
at various temperatures,we specify the equilibrium transition temperature
$T_{\mathrm{c}}$ in Sec. III.2.3. The results of the equilibrium simulations
are shown in Figs. S.5 (a) and (b), where $\bar{u}_{T}$ and $\bar{m}_{T}$
exhibits the discontinuous jump at $T=T_{\mathrm{c}}\simeq 0.150$. This
indicates that the first-order transition occurs at $T=T_{\mathrm{c}}$. Figure
S.5 (a) is the same as that displayed in Fig. 1 (b).

####  III.2.1 Heating experiment

We start with the initial condition $(\phi,\pi)_{T_{\rm
min}}(\boldsymbol{r},0)=(\mu_{1},0)$ at $T=T_{\rm min}$. Note that the
direction of the symmetry breaking is specified by this initial condition in
the following heating experiment. We numerically solve the Langevin equation
(S.59) up to $t=t_{\rm f}$ at $T_{\rm min}=T_{k=0}$, and calculate the time
average of observables over the time interval $t_{\rm f}-t_{0}$. That is, for

$\displaystyle\bar{u}_{T_{k}}^{\rm
heat}\equiv\frac{1}{L_{x}L_{y}}\frac{1}{t_{\rm f}-t_{0}}\int_{t_{0}}^{t_{\rm
f}}dt\>{\cal H}(\phi(t),\pi(t)),$ (S.67) $\displaystyle\bar{m}_{T_{k}}^{\rm
heat}\equiv\frac{1}{L_{x}L_{y}}\int_{D}d^{2}{\boldsymbol{r}}\>\frac{1}{t_{\rm
f}-t_{0}}\int_{t_{0}}^{t_{\rm f}}dt\>{\hat{m}}({\boldsymbol{r}},t),$ (S.68)

we obtain $(\bar{u}_{T_{0}}^{\rm heat},\bar{m}_{T_{0}}^{\rm heat})$. We then
increment the temperature from $T_{0}$ to $T_{1}$ and reset the time to be
$t=0$. The initial condition at $T=T_{1}$ is equal to the final state at
$T_{0}$, i.e.,
$(\phi,\pi)_{T_{1}}(\boldsymbol{r},0)=(\phi,\pi)_{T_{0}}(\boldsymbol{r},t_{\rm
f})$. We then numerically solve the Langevin equation (S.59) up to $t=t_{\rm
f}$. We repeat the same procedure changing the temperature from $T_{k-1}$ to
$T_{k}$ and setting
$(\phi,\pi)_{T_{k}}(\boldsymbol{r},0)=(\phi,\pi)_{T_{k-1}}(\boldsymbol{r},t_{\rm
f})$, where $1\leq k\leq 200$. We plot the temperature dependence of
$\bar{u}_{T_{k}}^{\rm heat}$ and $\bar{m}_{T_{k}}^{\rm heat}$ as displayed in
Figs. S.5 (c) and (d). At $T\simeq 0.166$, the system suddenly enters the
disordered state with $\bar{m}_{T}^{\rm heat}\simeq 0$. At the same
temperature, $\bar{u}_{T}^{\rm heat}$ also shows a discontinuous jump.

####  III.2.2 Cooling experiment

We start with an initial state which is obtained as the final state of the
heating experiment. We then decrease the temperature from $T_{k}$ to
$T_{k-1}$, where $k$ changes from $k=200$ to $k=1$. From each $T_{k}$, the
time is reset as $t=0$, and the initial condition at $T=T_{k-1}$ is
$(\phi,\pi)_{T_{k-1}}(\boldsymbol{r},0)=(\phi,\pi)_{T_{k}}(\boldsymbol{r},t_{\rm
f})$ in the cooling experiment. Similarly to the heating experiment, we define

$\displaystyle\bar{u}_{T_{k}}^{\rm
cool}\equiv\frac{1}{L_{x}L_{y}}\frac{1}{t_{\rm f}-t_{0}}\int_{t_{0}}^{t_{\rm
f}}dt\>{\cal H}(\phi(t),\pi(t)),$ (S.69) $\displaystyle\bar{m}_{T_{k}}^{\rm
cool}\equiv\frac{1}{L_{x}L_{y}}\int_{D}d^{2}{\boldsymbol{r}}\>\frac{1}{t_{\rm
f}-t_{0}}\int_{t_{0}}^{t_{\rm f}}dt\>{\hat{m}}({\boldsymbol{r}},t).$ (S.70)

At a glance, (S.69) and (S.70) are the same as (S.67) and (S.68). The only but
significant difference is the initial condition for each $T_{k}$. In the
cooling experiment, we obtain Figs. S.5 (e) and (f). They are obviously
different from Figs. S.5 (c) and (d). Through the cooling experiment, the
system always shows the disordered state with $\bar{m}_{T}\simeq 0$ and the
symmetry breaking never occurs, which suggests that the supercooled metastable
disordered state can exist even close to the zero temperature. These results
indicate the existence of a strong hysteresis in this system.

Figure S.5:  Temperature dependence of the energy density $\bar{u}_{T}$
(panels (a), (c), and (e)) and the order parameter $\bar{m}_{T}$ (panels (b),
(d), and (f)) as the function of $T$ for the equilibrium with the Langevin
equation (S.71) (panels (a) and (b)), the heating experiment with the Langevin
equation (S.59) (panels (c) and (d)), and the cooling experiment with the
Langevin equation (S.59) (panels (e) and (f)).

####  III.2.3 Equilibrium values

Now, we determine which quantity $(\bar{u}_{T}^{\rm heat},\bar{m}_{T}^{\rm
heat})$ or $(\bar{u}_{T}^{\rm cool},\bar{m}_{T}^{\rm cool})$ corresponds to
the equilibrium value. Here, we use the equilibration method developed by
Ohzeki and Ichiki [5]. We consider a duplicated system comprising $(\phi_{\rm
A},\pi_{\rm A})$ and $(\phi_{\rm B},\pi_{\rm B})$, each of which is defined in
the periodic box with $(L_{x},L_{y})=(192,64)$. The duplicated system evolves
by

$\displaystyle\begin{aligned} &\partial_{t}\phi^{a}_{\rm A}=\frac{\delta{\cal
H}(\phi_{\rm A},\pi_{\rm A})}{\delta\pi^{a}_{\rm A}}-\frac{\delta{\cal
H}(\phi_{\rm B},\pi_{\rm B})}{\delta\pi^{a}_{\rm
B}},\quad\partial_{t}\pi^{a}_{\rm A}=-\frac{\delta{\cal H}(\phi_{\rm
A},\pi_{\rm A})}{\delta\phi^{a}_{\rm A}}-\frac{\delta{\cal H}(\phi_{\rm
B},\pi_{\rm B})}{\delta\phi^{a}_{\rm B}}-\pi^{a}_{\rm A}+\xi^{a}_{\rm A},\\\
&\partial_{t}\phi^{a}_{\rm B}=\frac{\delta{\cal H}(\phi_{\rm B},\pi_{\rm
B})}{\delta\pi^{a}_{\rm B}}+\frac{\delta{\cal H}(\phi_{\rm A},\pi_{\rm
A})}{\delta\pi^{a}_{\rm A}},\quad\partial_{t}\pi^{a}_{\rm
B}=-\frac{\delta{\cal H}(\phi_{\rm B},\pi_{\rm B})}{\delta\phi^{a}_{\rm
B}}+\frac{\delta{\cal H}(\phi_{\rm A},\pi_{\rm A})}{\delta\phi^{a}_{\rm
A}}-\pi^{a}_{\rm B}+\xi^{a}_{\rm B},\end{aligned}$ (S.71)

where $\xi_{i}^{a}(\boldsymbol{r},t)$ satisfies

$\displaystyle\left\langle\xi_{i}^{a}({\boldsymbol{r}},t)\xi_{j}^{b}({\boldsymbol{r}}^{\prime},t^{\prime})\right\rangle=2\gamma
T\delta^{ab}\delta_{ij}\delta({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime})\delta(t-t^{\prime}).$
(S.72)

Writing down the Fokker-Planck equation, we can confirm that the stationary
distribution of (S.71) coincides with the canonical distribution (S.57). It is
also proved that the relaxation time of (S.71) to the true equilibrium is
shorter than that of the original system (S.59) [5].

We solve (S.71) by discretizing the system with the lattice spacing $\Delta
x=1/8$ and using the same methods as those for the original Langevin equation
(S.59). The initial condition at the temperature $T$ is given as a couple of
the final states obtained in the heating and cooling experiments,

$\displaystyle(\phi_{\rm A},\pi_{\rm A})_{T}(\boldsymbol{r},0)$
$\displaystyle=(\phi,\pi)_{T}(\boldsymbol{r},t_{\rm f})^{\rm heat},$ (S.73)
$\displaystyle(\phi_{\rm B},\pi_{\rm B})_{T}(\boldsymbol{r},0)$
$\displaystyle=(\phi,\pi)_{T}(\boldsymbol{r},t_{\rm f})^{\rm cool},$ (S.74)

where $(\phi,\pi)_{T}(\boldsymbol{r},t_{\rm f})^{\rm heat}$ and
$(\phi,\pi)_{T}(\boldsymbol{r},t_{\rm f})^{\rm cool}$ are the final states at
$T$ calculated in Sec. III.2.1 and Sec. III.2.2, respectively. Let $t_{0}$ and
$t_{\rm f}$ be the same value used in numerical simulations of the original
Langevin equation (S.59). We observe that the subsystems A and B become to
exhibit the same statistical state until $t=t_{0}$. We thus expect that the
obtained states $(\phi_{\rm A},\pi_{\rm A})_{T}(\boldsymbol{r},t)$ and
$(\phi_{\rm B},\pi_{\rm B})_{T}(\boldsymbol{r},t)$ for $t\geq t_{0}$ obey the
respective canonical distribution $\rho_{T}^{\mathrm{can}}(\phi_{\rm
A},\pi_{\rm A})$ and $\rho_{T}^{\mathrm{can}}(\phi_{\rm B},\pi_{\rm B})$
defined in (S.57). The equilibrium values $\bar{u}_{T}$ and $\bar{m}_{T}$
defined in (S.64) and (S.65) are then given by time averages in the interval
$t_{\rm f}-t_{0}$. Explicitly, we calculate them as

$\displaystyle\bar{u}_{T}=\frac{1}{L_{x}L_{y}}\frac{1}{t_{\rm
f}-t_{0}}\int_{t_{0}}^{t_{\rm f}}dt\>\frac{{\cal H}(\phi_{\rm A}(t),\pi_{\rm
A}(t))+{\cal H}(\phi_{\rm B}(t),\pi_{\rm B}(t))}{2},$ (S.75)
$\displaystyle\bar{m}_{T}=\frac{1}{L_{x}L_{y}}\int_{D}d^{2}{\boldsymbol{r}}\>\frac{1}{t_{\rm
f}-t_{0}}\int_{t_{0}}^{t_{\rm f}}dt\>\sum_{a=1}^{q-1}\frac{\phi_{\rm
A}^{a}({\boldsymbol{r}})+\phi_{\rm B}^{a}({\boldsymbol{r}})}{2}\mu_{1}^{a}.$
(S.76)

Since the duplicated system has $(\phi_{\rm A},\pi_{\rm A})$ and $(\phi_{\rm
B},\pi_{\rm B})$, we take the arithmetic means of the two in (S.75) and
(S.76). Figures S.5 (a) and (b) show $\bar{u}_{T}$ and $\bar{m}_{T}$ obtained
by (S.75) and (S.76), respectively.

In Fig. S.2, we also show the data of Fig. S.5 (a) in addition to the micro-
canonical temperature $\bar{T}_{E}$ as a function of $u=E/(L_{x}L_{y})$. The
both data show quantitatively the same behavior including the states with the
interface at $T=\bar{T}_{E}=T_{\rm c}$, which shows the equivalence of
thermodynamic ensembles (S.50) and (S.57). This result also demonstrates the
validity of the numerical method generating the canonical ensemble (S.57).

## References

  * [1] N. Nakagawa and S.-i. Sasa, Global thermodynamics for heat conduction states, J. Stat. Phys. 177, 825-888 (2019).
  * [2] 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).
  * [3] B. Fornberg, A practical Guide to Pseudospectral Methods, Cambridge University Press (1995).
  * [4] R. I. McLachlan and G. R. W. Quispel, Acta Numerica 11, 341 (2002).
  * [5] M. Ohzeki and A. Ichiki, Langevin dynamics neglecting detailed balance condition, Phys. Rev. E 92, 012105 (2015).


# Universal crosstalk of twisted light in random media

David Bachmann  Physikalisches Institut, Albert-Ludwigs-Universität Freiburg,
Hermann-Herder-Str. 3, D-79104 Freiburg, Germany  Asher Klug  School of
Physics, University of the Witwatersrand, Private Bag 3, Johannesburg 2050,
South Africa  Mathieu Isoard  Physikalisches Institut, Albert-Ludwigs-
Universität Freiburg, Hermann-Herder-Str. 3, D-79104 Freiburg, Germany
Vyacheslav Shatokhin  Physikalisches Institut, Albert-Ludwigs-Universität
Freiburg, Hermann-Herder-Str. 3, D-79104 Freiburg, Germany  EUCOR Centre for
Quantum Science and Quantum Computing, Albert-Ludwigs-Universität Freiburg,
Hermann-Herder-Str.3, D-79104 Freiburg, Germany  Giacomo Sorelli  Laboratoire
Kastler Brossel, Sorbonne Université, ENS-Université PSL, Collège de France,
CNRS; 4 place Jussieu, F-75252 Paris, France  Andreas Buchleitner
Physikalisches Institut, Albert-Ludwigs-Universität Freiburg, Hermann-Herder-
Str. 3, D-79104 Freiburg, Germany  EUCOR Centre for Quantum Science and
Quantum Computing, Albert-Ludwigs-Universität Freiburg, Hermann-Herder-Str.3,
D-79104 Freiburg, Germany  Andrew Forbes  School of Physics, University of the
Witwatersrand, Private Bag 3, Johannesburg 2050, South Africa

(August 29, 2024)

###### Abstract

Structured light offers wider bandwidths and higher security for
communication. However, propagation through complex random media, such as the
Earth’s atmosphere, typically induces intermodal crosstalk. We show
numerically and experimentally that coupling of photonic orbital angular
momentum (OAM) modes is governed by a _universal_ function of a single
parameter – the ratio between the random medium’s and the beam’s transverse
correlation lengths, even in the regime of pronounced intensity fluctuations.

_Introduction.—_ Optical communication strives to answer the growing demand of
high-bandwidth links Lavery  _et al._ (2018); Puttnam  _et al._ (2021); Zhou
_et al._ (2021). In particular, photonic spatial degrees of freedom such as
_orbital angular momentum_ (OAM) offer an unbounded Hilbert space to encode
information, and provide intrinsic support for quantum key or entanglement
distribution protocols Pors  _et al._ (2011a); Leach  _et al._ (2012); Fickler
_et al._ (2012); Otte  _et al._ (2020), via high-dimensional multiplexing
which enables greater channel capacities Wang  _et al._ (2012) and enhanced
security Gröblacher  _et al._ (2006). While communication based on such
_twisted photons_ was successfully demonstrated for table-top Mirhosseini  _et
al._ (2015), indoor Vallone  _et al._ (2014) or short outdoor Sit  _et al._
(2017) channels, transport of photonic OAM through _complex random media_ ,
e.g., the atmosphere Pors  _et al._ (2011b); Malik  _et al._ (2012); Ibrahim
_et al._ (2013), water Ren  _et al._ (2016); Bouchard  _et al._ (2018) or
multimode fiber Zhou  _et al._ (2021); Matthès  _et al._ (2021), remains
challenging. Typically, stochastic fluctuations of the underlying media’s
refractive index induce phase distortions as well as intensity fluctuations
upon propagation, leading to transmission losses and _intermodal crosstalk_
Paterson (2005); Cox  _et al._ (2021), potentially disrupting communication.
Of particular interest for applications is crosstalk between spatial modes
with opposite OAM.

In this Letter, we establish the universal dependence of the mentioned
crosstalk on the distortion strength in generic disordered media, and identify
the ratio of the transverse length scales of the light beam and the media as
the sole parameter governing this behavior.

_The model.—_ We consider transverse optical modes $\psi({\bf r})$ arising
within the paraxial approximation as solutions of the _parabolic wave
equation_ , which describes free diffraction in vacuum Andrews (2005). The
fundamental properties and potential applications of such modes have been
intensively studied Forbes  _et al._ (2021). Structured light intrinsically
features distinct phase and intensity profiles which may be characterized by
their transverse _correlation length_ $\xi$ Leonhard  _et al._ (2015), cf.
Fig. 1$(a)$. For example, OAM-carrying _Laguerre-Gaussian_ (LG) modes Allen
_et al._ (1992),
$u_{p,\ell}(\rho,\varphi,z=0)=R_{p,|\ell|}(\rho)\exp\left(i\ell\varphi\right)$,
are described by their radial and azimuthal indices $p\in\mathbb{N}_{0}$ and
$\ell\in\mathbb{Z}$, respectively, which allow for high-dimensional
information encoding.  Figure 1:  $(a)$ LG beam ($p=0,\ \ell=3$) with
indicated transverse correlation length $\xi$ propagating through a single
realization of a random medium represented by a phase screen with correlation
length $r_{0}$. $(b)$ Experimental implementation of the split-step method via
special light modulators (SLMs). $(c)$ Crosstalk amplitudes of a source LG
mode ($p=0,\ \ell=3$) after propagation through a random medium and vacuum-
propagated LG modes ($p=0,\ \ell^{\prime}$) for increasing distortion strength
$w=w_{0}/r_{0}$ (from top to bottom). Original and opposed OAM, i.e.,
$\pm\ell$, are indicated by the solid and dashed black lines, respectively.
In particular, $\ell$ governs the beam’s OAM which manifests itself in form of
intertwined helical phase fronts along the propagation axis $z$, as
illustrated in Fig. 1$(a)$. Consequently, with increasing $\ell$ the phase
structure of such beams becomes finer, leading to shorter correlation lengths.
Likewise, conservation of angular momentum ensures that $\ell$ is preserved
upon free diffraction.

However, photonic phase fronts are fragile under inherent refractive index
fluctuations of the atmosphere or other random media. This sensitivity is
modeled by the _stochastic_ parabolic equation where free diffraction is
combined with noise induced by the random medium Andrews (2005); Segev  _et
al._ (2013). To find its solution, we employ the _split-step method_ Schmidt
(2010); Chatterjee and Mohamed (2014), which relies on segmenting the entire
propagation path into discrete, medium-induced phase modulations, i.e., _phase
screens_ [see Fig. 1$(a)$], interconnected by free diffraction. We implement
this split-step propagation numerically through Fourier optics methods Born
(1999); Goodman (2005); Schmidt (2010). In particular, the phase screens
feature the medium’s transverse correlation length. For atmospheric
_Kolmogorov turbulence_ this characteristic length scale is known as the
_Fried parameter_ $r_{0}$ Andrews (2005) and the ratio $w:=w_{0}/r_{0}$ with
the _beam waist_ $w_{0}$ quantifies the _distortion strength_. Similarly, the
_Rytov variance_ $\sigma_{\text{R}}^{2}$ quantifies intensity fluctuations and
distinguishes between _weak scintillation_ $\sigma_{\text{R}}^{2}<1$ and
_strong scintillation_ $\sigma_{\text{R}}^{2}\geq 1$ Ishimaru (1978). For a
faithful description of light propagation through random media, we ensure that
every propagation step, i.e., a single phase screen and subsequent free
diffraction, complies with the weak scintillation condition
($\sigma_{\text{R}}^{2}<1$). Furthermore, we represent more general random
media by _Gaussian noise_ based on normally distributed block matrices with
block or _pixel size_ $r_{0}$, independently altering both amplitude and phase
of given transverse modes. Examples of a perturbed LG beam are shown in Fig.
2, where we distinguish three regimes of distortion: weak ($w<1$), moderate
($w\approx 1$) and strong ($w>1$).

To benchmark our theoretical description, we realize a turbulent link
experimentally by downscaling realistic channels to tabletop dimensions
Rodenburg  _et al._ (2014). The setup, shown in Fig. 1$(b)$, allows to
experimentally simulate propagation through weak to strong scintillation. The
salient elements are divided into three stages. In the generation stage, a He-
Ne laser beam (wavelength $\lambda=633$ nm) was expanded and collimated before
being directed onto a reflective PLUTO-VIS HoloEye _spatial light modulator_
(SLM1) which generated the desired source mode. This mode then entered the
turbulent section of the setup where it passed through a thick turbulent
medium implemented by a two stage split-step propagation, with turbulence
phases programmed to SLM2 and SLM3, each followed by 1 m of free diffraction.
This 2 m laboratory system maps to a $L=200$ m real-world channel with
effective Rytov variances of up to $\sigma_{R}^{2}=1.8$ matching the
numerically investigated conditions. While for this intermediate channel’s
length mode spreading effects due to diffraction are not negligible, the
modes’ distortions remain dominated by medium-induced phase and intensity
fluctuations. In the final stage, the modal decomposition was performed
optically Pinnell  _et al._ (2020) with the aid of a match filter programmed
to SLM4 and a Fourier transforming lens, with the on-axis intensity measured
with a camera (CCD detector).

Figure 2: Intensity-weighted transverse phase profiles of a LG mode ($p=0,\
\ell=5$) after propagation through Kolmogorov turbulence (top) or Gaussian
noise (bottom) for three different distortion strengths $w=1/2,1,2$,
corresponding to $\sigma_{\text{R}}^{2}=0.24,0.76,2.42$ (left to right).

_Crosstalk among distorted modes.—_ Communication relies on the correct
identification of received modes despite partial distortion due to propagation
through disorder. The probability of identifying a source mode’s original OAM
$\ell$ after propagation as $\ell^{\prime}$ is quantified via the _crosstalk
matrix_

$C_{\ell,\ell^{\prime}}(L):=|\langle\tilde{u}_{\ell}(L),u_{\ell^{\prime}}(L)\rangle|^{2},$
(1)

Figure 3: Numerically obtained average (5000 realizations) crosstalk matrices
in $(a)$ Kolmogorov turbulence and $(b)$ Gaussian noise. For each panel the
crosstalk matrices (top) correspond to pure OAM modes, LG modes and pure
amplitude modes (left to right) in moderate distortion ($w=1$). Rows
$C_{\ell=5,\ell^{\prime}}$ (bottom) depict numerical (orange curve, error
bands correspond to one standard deviation) as well as experimental results
(blue bars, average of 60 realizations). The crosstalk matrices have been
normalized such that
$\sum_{\ell,\ell^{\prime}=-10}^{10}C_{\ell,\ell^{\prime}}=1$.

where $\tilde{u}_{\ell}(L)$ and $u_{\ell^{\prime}}(L)$ are the propagated
complex modes 111The full mode functions read $u_{p=0,\ell}(\rho,\phi,z=L)$,
where the radial indices $p=0$ and the transverse coordinates $\rho,\phi$ were
suppressed for brevity. across a random medium and a vacuum channel each of
length $L$, respectively, and $\langle\cdot,\cdot\rangle$ denotes the standard
scalar product. Consequently, higher transmission fidelities correspond to
weaker crosstalk, i.e., to matrices $C_{\ell,\ell^{\prime}}(L)$ with
prevailing diagonals. Figure 1(c) sketches the crosstalk amplitudes for a
fixed source LG mode with initial OAM $\ell=3$ after propagation through a
channel of length $L$ for three distortion regimes quantified by $w$. As
observed previously Klug  _et al._ (2021), under weak distortion (top row) the
detection probability of the original OAM $\ell=3$ (solid black line) is
highest and crosstalk is limited to neighboring modes. In contrast, under
strong distortion [bottom row of Fig. 1(c)], we find a $\pm\ell$-symmetric
distribution of the crosstalk amplitudes with disjoint maxima close to
$\pm\ell$ (solid and dashed black lines). In this case, the mode’s phase front
is destroyed, while its amplitude – which is on average independent of the
sign of OAM – partially survives. Finally, the moderate distortion regime [see
middle row of Fig. 1(c)] is of particular interest: It represents the
transition between rather faithful OAM transmission (top row) and dominant
phase destruction (bottom row). Both, the beam’s amplitude and phase profiles
are altered, but neither is completely destroyed, resulting in the observed
asymmetric doublet.

To elucidate the origin of the disorder-induced crosstalk between OAM modes,
we decompose the transmitted modes $\tilde{u}_{p=0,\ell}(\rho,\phi)$ into the
LG basis $u_{p=0,\ell}(\rho,\phi)$, as well as in the helical
$\exp(i\ell\phi)$ and amplitude $R_{p=0,|\ell|}(\rho)$ components of LG modes
alone. We refer to the former as _pure OAM_ and to the latter as _pure
amplitude modes_. The upper row in Fig. 3 shows the resulting average (5000
disorder realizations) crosstalk matrices in moderate, i.e., $w=1$, Kolmogorov
turbulence $(a)$, and under Gaussian noise $(b)$, for pure OAM, LG and pure
amplitude modes, from left to right. Crosstalk among pure OAM modes [top left
plots in Fig. 3$(a)$ and $(b)$] is distributed diagonally, evidence of the
absence of coherent coupling between modes with opposite OAM, confirming
previous findings Klug  _et al._ (2021). The remaining phase information
expressed by the diagonal’s broadening is due to the disorder-induced coupling
to neighboring modes $\ell^{\prime}=\ell\pm 1,2..$. On the other hand, pure
amplitude modes [top right plots in Fig. 3$(a)$ and $(b)$] feature crosstalk
matrices with pronounced, symmetric diagonal and _anti-diagonal_ structures.
This is unsurprising because the amplitude of LG modes does not depend on the
OAM’s sign. Ultimately, when considering LG modes [top middle plots in Fig.
3$(a)$ and $(b)$] both effects are combined: Their crosstalk matrix features
the salient diagonal originating from direct OAM coupling, as well as a
suppressed anti-diagonal descending from $\pm\ell$-amplitude overlap, in
agreement with Sorelli  _et al._ (2019). The damping of anti-diagonal matrix
elements precisely results from the competition between pure OAM and amplitude
crosstalk which, respectively, suppress and enhance the overlap of LG modes
with opposite OAM.

Plots in the lower row of Fig. 3 compare experimental measurements with
numerical results in form of crosstalk matrix rows which describe the
detection probabilities of various $\ell^{\prime}$ when transmitting $\ell=5$.
Experimental results (blue bars) qualitatively agree with numerical results
(orange curves). Slightly lower experimentally measured crosstalk for larger
$\ell^{\prime}$ may originate from aperture effects damping overall
transmission of these wider modes Andrews (2005). For LG crosstalk [bottom
middle plots of Fig. 3$(a)$ and $(b)$], we observe an _asymmetric doublet_
whose second diminished peak is systematically shifted towards the center.
While crosstalk between opposite $\ell$ is still amplified, this behavior
highlights that we do not observe direct $\pm\ell$-coupling of LG modes but
rather a tradeoff between destruction of OAM-preserving phase and surviving
$\pm\ell$-symmetric intensity. Moreover, the LG crosstalk [middle plots in
Fig. 3$(a)$ and $(b)$] is, as expected, well approximated by the product of
pure OAM and pure intensity crosstalk. Note that the observed average doublets
do not contradict with the conservation of momentum after ensemble-averaging:
The received modes are projected onto $p=0$ LG basis modes and, hence,
coupling to OAM carrying modes with higher radial indices is ignored.

Subsequently, we systematically quantify the found anti-diagonal contribution
of LG crosstalk for a range of distortion strengths $w$ and azimuthal indices
$\ell$. First, we note that anti-diagonal crosstalk becomes weaker for larger
$|\ell|$ [compare fading of anti-diagonals in the top row of Fig. 3$(a)$ and
$(b)$ for LG modes]. This observation is consistent with well-known results
that entangled photonic OAM qubit states with opposite $\ell$ become more
robust in turbulence as $\ell$ increases Smith and Raymer (2006). Careful
study of LG modes Leonhard  _et al._ (2015) within the weak scintillation
regime attributed this enhanced robustness to the finer transverse spatial
structure of such photons characterized by their phase correlation length
$\xi(\ell)$. The latter gives the average distance between two points in the
transverse profile of LG beams with a phase difference of $\pi/2$ and was
found analytically Leonhard  _et al._ (2015) (for generalizations of
$\xi(\ell)$, see Yang  _et al._ (2021)).

In particular, we investigate the ratio between anti-diagonal and diagonal
crosstalk, i.e., $C_{\ell,-\ell}$ over $C_{\ell,\ell}$, which quantifies
$\pm\ell$-coupling of OAM. This ratio is plotted in Fig. 4 versus the medium’s
transverse correlation length $r_{0}$ normalized by the LG beam’s phase
correlation length $\xi(\ell)$ for Kolmogorov turbulence (squares) and
Gaussian noise (triangles). Such rescaling leads in both cases to a
_universal_ , i.e. $\ell$-independent, crosstalk curve which was identified as
a _Gaussian decay_. It describes the transition from symmetric doublets (top
inset in Fig. 4) in strong turbulence, where sign-sensitive phase information
is destroyed, to dominant direct $\ell$-coupling in weak turbulence (bottom
inset). The Gaussian form of the universal curve can be related to the shape
of pure OAM coupling which is well approximated by normal distributions around
the transmitted $\ell$ as a function of $\ell^{\prime}$ [cf. bottom left plots
in Fig. 3$(a)$ and $(b)$], consistent with analytic findings for Kolmogorov
turbulence Klug  _et al._ (2021). Furthermore, the numerical results are in
agreement with experimental data (circles in Fig. 4) in the moderate and
strong turbulence regime; the measurement for weak turbulence (third point
from the left) shows a systematic offset due to the decreased signal-to-noise
ratio of the diminished second peak in this case. The universality of
crosstalk for both Kolmogorov turbulence and generic Gaussian noise suggests
that it can also be found in other random media.

_Discussion.—_ The universal behavior of rescaled crosstalk amplitudes between
modes of opposite OAM was previously predicted in the regime of weak
scintillation, i.e., within the single phase screen model, for Kolmogorov
Leonhard  _et al._ (2015) and non-Kolmogorov turbulence Bachmann  _et al._
(2019) alike, as well as under deterministic perturbation of twisted photons
by angular apertures Sorelli  _et al._ (2020). Moreover, it was established
that in the weak scintillation regime bipartite entanglement of twisted
photonic qubit states with opposite OAM depends on the same parameters as
crosstalk – and, hence, exhibits universal behavior as well Leonhard  _et al._
(2015); Bachmann  _et al._ (2019); Sorelli  _et al._ (2020).

However, entanglement of twisted photons depends on intensity fluctuations
Roux  _et al._ (2015) and this additional sensitivity renders rescaling into a
universal entanglement evolution impossible under strong scintillation. Our
present theoretical as well as experimental results show that _crosstalk_
among structured light modes is nonetheless universal under strong
scintillation conditions. In other words, crosstalk is immune to intensity
fluctuations. Importantly, only the characteristic transverse correlation
length of the medium-induced errors needs be known in order to understand the
impact of modal scattering of any structured light mode, with clear impact in
imaging and communication through noisy channels. In the case of LG modes, our
results show that the ubiquitous superposition of $\pm\ell$, often used as
quantum qubit states and vectorial classical states for communication Ndagano
_et al._ (2017); Vallone  _et al._ (2014); Ibrahim  _et al._ (2013); Leonhard
_et al._ (2015); Smith and Raymer (2006); Bachmann  _et al._ (2019), is in
fact close to the least desirable combination in distorting media.

Figure 4:  Ratio of average (5000 realizations) crosstalk between two OAM-
opposed, i.e. $\ell$ and $-\ell$, LG modes ($w_{0}$ fixed) versus ratio of the
medium’s and beam’s transverse coherence lengths in numerically simulated
Kolmogorov turbulence (squares) and Gaussian noise (triangles) for a range of
$|l|$ (see color bar legend) with fitted Gaussian curves (black, parameters
see in footnote fit ). Circles with error bars represent experimentally
measured average (60 realizations) crosstalk in Kolmogorov turbulence. Insets
illustrate rows (at $\ell=5$) of corresponding crosstalk matrices (same axes
as in the bottom rows of Fig. 3). Error bands in the insets or bars in the
main plot correspond to one standard deviation.

_Conclusion.—_ In this Letter, we have studied the crosstalk of twisted
photons in Kolmogorov turbulence and Gaussian noise. In both cases, we have
numerically as well as experimentally confirmed stronger crosstalk among LG
modes of opposed OAM. Instead of originating from direct OAM cross coupling,
this behavior was identified as a tradeoff effect between matching intensity
patterns and destroyed phase information. Moreover, we have uncovered a
universal Gaussian crosstalk decay complying with experimental measurements by
setting the random media’s transverse correlation length into relation with
the LG beam’s phase structure.

These new insights on how crosstalk among twisted photons arises may lead to
novel venues for communication. In particular, we conclude that optimizing the
beam’s transverse length for given distortion strengths will diminish induced
crosstalk, independent of the precise nature of the underlying random medium
itself. As we have demonstrated here, this may be achieved by increasing
$\ell$, but we further anticipate adaptation of the beam waist or the LG
modes’ radial index $p$ to be fruitful in future considerations.

_Acknowledgements.—_ The authors acknowledge financial support of the
Studienstiftung des deutschen Volkes and the resources by the state of Baden-
Württemberg through bwHPC. This work was partially funded by French ANR under
COSMIC project (ANR-19-ASTR0020-01). V.S. and A.B. acknowledge partial funding
and support through the Strategiefonds der Albert-Ludwigs-Universität Freiburg
and the Georg H. Endress Stiftung.

## References

  * Lavery  _et al._ (2018) M. P. J. Lavery, M. M. Abadi, R. Bauer, G. Brambilla, L. Cheng, M. A. Cox, A. Dudley, A. D. Ellis, N. K. Fontaine, A. E. Kelly, C. Marquardt, S. Matlhane, B. Ndagano, F. Petruccione, R. Slavík, F. Romanato, C. Rosales-Guzmán, F. S. Roux, K. Roux, J. Wang, and A. Forbes, Nat. Photonics 12, 249 (2018).
  * Puttnam  _et al._ (2021) B. J. Puttnam, G. Rademacher, and R. S. Luís, Optica 8, 1186 (2021).
  * Zhou  _et al._ (2021) Y. Zhou, B. Braverman, A. Fyffe, R. Zhang, J. Zhao, A. E. Willner, Z. Shi, and R. W. Boyd, Nat. Commun. 12 (2021).
  * Pors  _et al._ (2011a) B.-J. Pors, F. Miatto, G. W. ’t Hooft, E. R. Eliel, and J. P. Woerdman, J. Opt. 13, 064008 (2011a).
  * Leach  _et al._ (2012) J. Leach, E. Bolduc, D. J. Gauthier, and R. W. Boyd, Phys. Rev. A 85 (2012).
  * Fickler  _et al._ (2012) R. Fickler, R. Lapkiewicz, W. N. Plick, M. Krenn, C. Schaeff, S. Ramelow, and A. Zeilinger, Science 338, 640 (2012).
  * Otte  _et al._ (2020) E. Otte, I. Nape, C. Rosales-Guzmán, C. Denz, A. Forbes, and B. Ndagano, J. Opt. Soc. Am. B 37, A309 (2020).
  * Wang  _et al._ (2012) J. Wang, J.-Y. Yang, I. M. Fazal, N. Ahmed, Y. Yan, H. Huang, Y. Ren, Y. Yue, S. Dolinar, M. Tur, and A. E. Willner, Nat. Photonics 6, 488 (2012).
  * Gröblacher  _et al._ (2006) S. Gröblacher, T. Jennewein, A. Vaziri, G. Weihs, and A. Zeilinger, New J. Phys. 8, 75 (2006).
  * Mirhosseini  _et al._ (2015) M. Mirhosseini, O. S. Magaña-Loaiza, M. N. O’Sullivan, B. Rodenburg, M. Malik, M. P. J. Lavery, M. J. Padgett, D. J. Gauthier, and R. W. Boyd, New J. Phys. 17, 033033 (2015).
  * Vallone  _et al._ (2014) G. Vallone, V. D’Ambrosio, A. Sponselli, S. Slussarenko, L. Marrucci, F. Sciarrino, and P. Villoresi, Phys. Rev. Lett. 113 (2014).
  * Sit  _et al._ (2017) A. Sit, F. Bouchard, R. Fickler, J. Gagnon-Bischoff, H. Larocque, K. Heshami, D. Elser, C. Peuntinger, K. Günthner, B. Heim, C. Marquardt, G. Leuchs, R. W. Boyd, and E. Karimi, Optica 4, 1006 (2017).
  * Pors  _et al._ (2011b) B.-J. Pors, C. H. Monken, E. R. Eliel, and J. P. Woerdman, Opt. Express 19, 6671 (2011b).
  * Malik  _et al._ (2012) M. Malik, M. O’Sullivan, B. Rodenburg, M. Mirhosseini, J. Leach, M. P. J. Lavery, M. J. Padgett, and R. W. Boyd, Opt. Express 20, 13195 (2012).
  * Ibrahim  _et al._ (2013) A. H. Ibrahim, F. S. Roux, M. McLaren, T. Konrad, and A. Forbes, Phys. Rev. A 88 (2013).
  * Ren  _et al._ (2016) Y. Ren, L. Li, Z. Wang, S. M. Kamali, E. Arbabi, A. Arbabi, Z. Zhao, G. Xie, Y. Cao, N. Ahmed, Y. Yan, C. Liu, A. J. Willner, S. Ashrafi, M. Tur, A. Faraon, and A. E. Willner, Sci. Rep. 6, 33306 (2016).
  * Bouchard  _et al._ (2018) F. Bouchard, A. Sit, F. Hufnagel, A. Abbas, Y. Zhang, K. Heshami, R. Fickler, C. Marquardt, G. Leuchs, R. w. Boyd, and E. Karimi, Opt. Express 26, 22563 (2018).
  * Matthès  _et al._ (2021) M. W. Matthès, Y. Bromberg, J. de Rosny, and S. M. Popoff, Phys. Rev. X 11, 021060 (2021).
  * Paterson (2005) C. Paterson, Phys. Rev. Lett. 94, 153901 (2005).
  * Cox  _et al._ (2021) M. A. Cox, N. Mphuthi, I. Nape, N. Mashaba, L. Cheng, and A. Forbes, IEEE J. Sel. Top. Quantum Electron. 27, 1 (2021).
  * Andrews (2005) L. Andrews,  _Laser beam propagation through random media_ (SPIE, Bellingham, 2005).
  * Forbes  _et al._ (2021) A. Forbes, M. de Oliveira, and M. R. Dennis, Nat. Photonics 15, 253 (2021).
  * Leonhard  _et al._ (2015) N. D. Leonhard, V. N. Shatokhin, and A. Buchleitner, Phys. Rev. A 91, 012345 (2015).
  * Allen  _et al._ (1992) L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, Phys. Rev. A 45, 8185 (1992).
  * Segev  _et al._ (2013) M. Segev, Y. Silberberg, and D. N. Christodoulides, Nat Photon 7, 197 (2013).
  * Schmidt (2010) J. Schmidt,  _Numerical simulation of optical wave propagation with examples in MATLAB_ (SPIE, Bellingham, 2010).
  * Chatterjee and Mohamed (2014) M. R. Chatterjee and F. H. A. Mohamed, Opt. Eng. 53, 126107 (2014).
  * Born (1999) M. Born,  _Principles of optics: Electromagnetic theory of propagation, interference and diffraction of light_ (Cambridge University Press, New York, 1999).
  * Goodman (2005) J. Goodman,  _Introduction to Fourier optics_ (Roberts & Co, Englewood, 2005).
  * Ishimaru (1978) A. Ishimaru,  _Wave Propagation and Scattering in Random Media_ (Academic Press, New York, 1978).
  * Rodenburg  _et al._ (2014) B. Rodenburg, M. Mirhosseini, M. Malik, O. S. Magaña-Loaiza, M. Yanakas, L. Maher, N. K. Steinhoff, G. A. Tyler, and R. W. Boyd, New J. Phys. 16, 033020 (2014).
  * Pinnell  _et al._ (2020) J. Pinnell, I. Nape, B. Sephton, M. A. Cox, V. Rodríguez-Fajardo, and A. Forbes, J. Opt. Soc. Am. A 37, C146 (2020).
  * Note (1) The full mode functions read $u_{p=0,\ell}(\rho,\phi,z=L)$, where the radial indices $p=0$ and the transverse coordinates $\rho,\phi$ were suppressed for brevity.
  * Klug  _et al._ (2021) A. Klug, I. Nape, and A. Forbes, New J. Phys. 23, 093012 (2021).
  * Sorelli  _et al._ (2019) G. Sorelli, N. Leonhard, V. N. Shatokhin, C. Reinlein, and A. Buchleitner, New J. Phys. 21, 023003 (2019).
  * Smith and Raymer (2006) B. J. Smith and M. G. Raymer, Phys. Rev. A 74 (2006).
  * Yang  _et al._ (2021) D. Yang, Z.-D. Hu, S. Wang, and Y. Zhu, Opt. Lett. 46, 5461 (2021).
  * Bachmann  _et al._ (2019) D. Bachmann, V. N. Shatokhin, and A. Buchleitner, J. Phys. A Math. Theor. 52, 405303 (2019).
  * Sorelli  _et al._ (2020) G. Sorelli, V. N. Shatokhin, and A. Buchleitner, J. Opt. 22, 024002 (2020).
  * Roux  _et al._ (2015) F. S. Roux, T. Wellens, and V. N. Shatokhin, Phys. Rev. A 92 (2015).
  * Ndagano  _et al._ (2017) B. Ndagano, B. Perez-Garcia, F. S. Roux, M. McLaren, C. Rosales-Guzman, Y. Zhang, O. Mouane, R. I. Hernandez-Aranda, T. Konrad, and A. Forbes, Nat. Phys. 13, 397 (2017).
  * (42) Fit of the form $f(x)=(1-a)\exp[-\frac{1}{2}\,(x/b)^{2}]+a$ with $a=0.083\pm 0.003,\ b=1.167\pm 0.005$ for Kolmogorov turbulence and $a=0.115\pm 0.004,\ b=1.552\pm 0.010$ for Gaussian noise .