---
author:
- 'P. Karpov'
- 'F. Piazza'
date: 'September 30, 2019'
title: 'Supplemental Material for “Crystalline droplets with emergent topological color-charge in many-body systems with sign-changing interactions”'
---

Density-wave phase transition: mean-field theory {#Appendix Density-wave transition}
================================================

In this Supplementary Section we construct a mean-field theory for density-wave phase transition and find the temperature dependence of the order parameter. This can be done if the number of atoms inside a $D$-dimensional sphere with radius $\xi$ is much greater than one: $n \xi^D \gg 1$.

Let $n_A$ and $n_B$ be the average concentrations of particles per cite at the two checkerboard sublattices $A$ and $B$. We define the sublattice imbalance order parameter as $$\begin{aligned}
m(T)= \frac{n_A-n_B}{2n}
\label{eq:order_param_Z2}\end{aligned}$$ Recall that $n = N/L^D$ is average number of particles per site in the uniform phase. We start from the density-wave phase $n_A\approx 2n$, $n_B \approx 0$, so $m\approx 1$ and want to describe a transition to the uniform state $n_A\approx n_B\approx n$, $m\approx 0$.

Consider a particle in the effective mean field potential created by other particles. The particle can be treated as an effective two-level system, where two states correspond to sublattices $A$ and $B$. The ground-state energy of the chosen particle is reached when it sits at the sublattice $A$: $$\begin{aligned}
E_0 = -n_A U_A + n_B U_B\end{aligned}$$ The particle is in the excited state, when it sits at the sublattice $B$, then its energy is $$\begin{aligned}
E_1 = -n_B U_A + n_A U_B\end{aligned}$$ The excitation energy of the particle is $$\begin{aligned}
E_1-E_0 = (n_A-n_B) [U_A+U_B]=2n\, m(T) [U_A+U_B]\end{aligned}$$ Let $n_0$ be the occupation number for the considered particle of sublattice $A$ and $n_1$ for sublattice $B$, so only one of $n_0, n_1$ is non-zero and $n_0+n_1=1$. Now we are dealing with a simple statistical-mechanics problem of a classical two-level system with Hamiltonian $H=n_0 E_0 + n_1 E_1$, but with additional self-consistency conditions $$\begin{aligned}
& \langle n_0 \rangle = n_A/2n \nonumber\\
& \langle n_1 \rangle = n_B/2n\end{aligned}$$ where average is defined as $\langle O \rangle = (O_0 \cdot e^{-E_0/T} + O_1 \cdot e^{-E_1/T})/Z$, and $Z$ is the single-particle partition function $Z=e^{-E_0/T}+e^{-E_1/T}$. Subtracting these two equations we get the self-consistency condition in a compact form $$\begin{aligned}
m(T) = 2\langle n_0 \rangle -1
\label{eq:self-consistency}\end{aligned}$$ Substituting $\langle n_0 \rangle =  e^{-E_0/T}/Z$ to (\[eq:self-consistency\]) we get $$\begin{aligned}
m = \frac{2}{1+e^{-2n[U_A+U_B] m/T}} - 1
\label{eq:self-consistency-full}\end{aligned}$$ Figure \[fig:plot\_m\_T\_analytical\] shows the dependence $m(T)$, implicitly defined by eq. (\[eq:self-consistency-full\]) (this equation is used in order to construct Fig. 2a of the main text). Expanding its RHS to the first order in $m$ we find $$\begin{aligned}
T_{\mathrm{DW}} = \frac{N [U_A + U_B ]}{L^D}
\label{eq:T_CDW_mean_field}\end{aligned}$$ The third-order expansion gives also the mean-field critical exponent $m(T)\sim |T-T_c|^{1/2}$).

Comparing the mean-field value of $T_{\mathrm{DW}}$ (\[eq:T\_CDW\_mean\_field\]) with the estimate obtained from the free energy arguments (\[eq:T\_CDW-arbitrary-D-N\]) we see that they indeed differ only by a numerical factor $2\ln2\approx 1.4$, with mean-field result giving higher value of the critical temperature.

We also note that the expression in the square brackets in (\[eq:T\_CDW\_mean\_field\]) is the integrated “staggered” strength of the interaction potential: $U_{stag} = -\sum (-1)^{i_x+i_y} U_{i_x,i_y} = U_A + U_B$. Therefore the mean-field expression for the DW-transition critical temperature can be expressed as $$\begin{aligned}
T_{\mathrm{DW}} = \frac{N U_{stag}}{L^D}
\label{eq:T_CDW_mean_field}\end{aligned}$$

![Analytical mean-field temperature dependence of the sublattice imbalance order parameter; $m(T)$ is implicitly defined by eq. (\[eq:self-consistency-full\]). []{data-label="fig:plot_m_T_analytical"}](figAppA_plot_m_T_mean_field_Z2.png){width="0.45\linewidth"}

We can estimate the limits of applicability of the mean-field theory using the Ginzburg criterion $t_G = |T_G-T_{\mathrm{DW}}|/T_{\mathrm{DW}} \sim a^2/(n \xi^2)$, where the Ginzburg temperature $T_G$ is the temperature above which fluctuation effects become important, $\xi$ is the range of the potential. The mean-field theory works better and better with the growing number of particles in $\xi \times \xi$ region, $n (\xi/a)^2$, ultimately becoming exact when the latter quantity goes to infinity.

Droplet-formation phase transition: exactly-solvable toy model and numerical tests {#Appendix Collapse transition}
==================================================================================

In this Supplementary Section we, first, consider a toy model for the droplet-formation phase transition and.

Consider $N$ particles occupying $L$ lattice sites in 1D. If two particles sit on the same lattice site they interact with energy $-U_0<0$ (on-site attraction) and particles at different sites don’t interact with each other. Let $n_i$ be the occupation number of site $i$. The Hamiltonian of the system is $$\begin{aligned}
H = -\frac{U_0}{2}\sum_{i=1}^L n_i (n_i-1) \approx -\frac{U_0}{2}\sum_{i=1}^L n_i^2\end{aligned}$$ Such model is equivalent to the fully-connected (mean-field) Potts model with $L$ states, which possesses a phase transition from disorder (uniform occupation of all sites) to order (all particles occupy the single site): for $L=2$ it is of the second order and for $L>2$ it is of the first order [@Wu:1982]. The exact solution[@Wu:1982] gives the critical temperature: $$\begin{aligned}
T_{\mathrm{dr}}  = \frac{L-2}{L-1} \frac{N U_0}{2\ln(L-1)}
\label{eq:Potts}\end{aligned}$$ In case $L \gg 1$ this gives us $$\begin{aligned}
T_{\mathrm{dr}} \approx \frac{N U_0}{2 \ln L}
\label{eq:Tc=NU0/2lnL}\end{aligned}$$ The simplest model described here works for the case when on-site interaction is much greater then the interaction at nearest-neighboring sites, which corresponds to the case of infinitely many degenerate modes in the optical cavity. Alternatively it can be viewed as a coarse-grained model for interaction with range $\xi$, where one site corresponds a coarse-grained cell of the size $\xi$, so we should substitute $L\rightarrow L/\xi$. Additional substitution $L/\xi \rightarrow (L/\xi)^D$ trivially generalizes eq. (\[eq:Tc=NU0/2lnL\]) to the $D$-dimensional case, thus we arrive to the formula (2) of the main text.

In order to test eq. (2) of the main text we perform Monte Carlo (MC) simulations of the 1D system with the same Gaussian potential with the range $\xi$ as we use in the main text. Since eq. (2) works only qualitatively, we check the linear dependence of $T_{\mathrm{dr}}$ on the number of particles $N$ (Fig. \[fig:plots\_Tdr\]a) and the dependence of the form $T_{\mathrm{dr}} = A / \ln(B/\xi)$ (Fig. \[fig:plots\_Tdr\]b).

Simulation of interaction of color-charged droplets {#Appendix Droplets interaction}
===================================================

Derivation of the light-mediated interaction of atoms and mapping to the classical lattice gas problem {#Appendix Derivation of the model}
======================================================================================================

![Geometry of the system. Cavity axis is parallel to $z$; both lasers and the atomic gas are in the $xy$ plane. []{data-label="fig:geometry"}](figAppGeometry.png){width="0.5\linewidth"}

In this Supplementary Section we describe the setting we propose for the experimental realization of the interaction we use in the main text. We consider an ensemble of atoms in a multimode cavity pumped transverse to the cavity axis by a pair of lasers with a frequency $\omega_L$ and mode functions $\mathbf{\Omega}(\mathbf{r}) = \mathbf{e}_z \Omega(\mathbf{r}) \equiv
\mathbf{e}_z (\Omega_1(\tilde{x}) + \Omega_2(\tilde{y})) \equiv  \mathbf{e}_z
(\cos(2\pi \tilde{x}/\lambda) + \cos(2\pi \tilde{y}/\lambda))$ (here $\tilde{x}=(x-y)/\sqrt{2}$, $\tilde{y}=(x+y)/\sqrt{2}$, Fig. \[fig:geometry\]). The cavity possesses (nearly) degenerate modes with frequencies $\omega_{c\alpha}$ (where $\alpha$ is the mode index). The atoms are treated as two-level systems with the level splitting $\omega_A$ and the transition dipole moment connecting the ground and excited state being parallel to $\mathbf{e}_x+\mathbf{e}_z$ (spatial isotropy can be broken by external magnetic field), which allows to couple the atomic dipole transitions to both pump and cavity modes. For more realistic experimental settings with many-level atoms such coupling can be achieved by employing vector polarizability of the atoms[@LeKien:2013].

For low density of atoms we can neglect their contact interaction and describe the light-mediated interaction in a multimode cavity by the following effective action [@Gopalakrishnan:2009]: $$\begin{aligned}
S = S_{at} + S_{int} + S_{em}\end{aligned}$$ $$\begin{aligned}
S_{at} =  \int d\mathbf{r} d\tau \left[ \Psi_g^*(\mathbf{r},\tau) \left( \partial_{\tau} -\frac{\hbar\nabla^2}{2M} -\frac{\mu}{\hbar} \right) \Psi_g(\mathbf{r},\tau) +
    \Psi_e^*(\mathbf{r},\tau) \left( \partial_{\tau} -\frac{\hbar\nabla^2}{2M} -\frac{\mu}{\hbar} + \Delta_a \right) \Psi_e(\mathbf{r},\tau)  \right]\end{aligned}$$ $$\begin{aligned}
S_{int} = \int d\mathbf{r} d\tau \left[ \sum_{\alpha} i g_{\alpha}(\mathbf{r}) \Psi^*_e(\mathbf{r},\tau) \Psi_g(\mathbf{r},\tau) a_{\alpha}(\tau) + i\Omega(\mathbf{r}) \Psi^*_e (\mathbf{r},\tau) \Psi_g(\mathbf{r},\tau) + h.c.\right]\end{aligned}$$ $$\begin{aligned}
S_{em} = \int d\tau \sum_{\alpha} a^*_{\alpha}(\tau) (\partial_{\tau}+\Delta_{c\alpha}) a_{\alpha}(\tau)\end{aligned}$$ Here $\Psi_g$ and $\Psi_e$ are the bosonic second-quantized fields corresponding to the ground and excited states of the atoms with mass $M$. We work in the reference frame rotating at the laser frequency $\omega_L$ and denote $\Delta_a = \omega_A - \omega_L$, $\Delta_{c\alpha} = \omega_{c\alpha} - \omega_L$. $g_{\alpha} (\mathbf{r}) = g_{0} \Xi_{\alpha}(\mathbf{r})$ where $g_0$ is the atom-cavity coupling strength and $\Xi_{\alpha}(\mathbf{r})$ is the normalized mode function of the cavity mode $\alpha$. Integrating out excited state $\Psi_e$ and cavity modes $a_{\alpha}$ we get the following effective action $$\begin{aligned}
S_{\text{eff}}[\Psi_g,\Psi_g^*]&=\int d\mathbf{r} d\tau \Psi_g^*(\mathbf{r},\tau) \left( \partial_{\tau} -\frac{\nabla^2}{2M} +(V(\mathbf{r})-\mu) \right ) \Psi_g(\mathbf{r},\tau)-\nonumber \\
&+ \int d\tau d\mathbf{r} d\mathbf{r}'
|\Psi_g(\mathbf{r},\tau)|^2 U(\mathbf{r},\mathbf{r}')  |\Psi_g(\mathbf{r}',\tau)|^2
\label{eq:eff-action-simplified}\end{aligned}$$ The external optical potential $V(\mathbf{r})$ and the cavity-mediated interaction $U(\mathbf{r}, \mathbf{r}')$ are given by $$\begin{aligned}
V(\mathbf{r}) &= -\frac{1}{\Delta_a} |\Omega(\mathbf{r})|^2 \label{eq:ExternalPotential}\\
U(\mathbf{r},\mathbf{r}') &=- \frac{g_0^2}{\Delta_a^2} \Omega^*(\mathbf{r}) \Omega(\mathbf{r}') \sum_{\alpha}  \frac{\Xi^*_{\alpha}(\mathbf{r}) \Xi_{\alpha}(\mathbf{r}')}{\Delta_{c\alpha}} \label{eq:InteractionPotential}\end{aligned}$$ We consider the case of red-detuned laser frequency with respect to both atom splitting and to cavity nearly degenerate modes $\Delta_a, \Delta_{c\alpha}>0$. In this case the atoms are attracted to maxima of $|\Omega(\mathbf{r})|^2$; the interaction is attractive at small distances $\mathbf{r}'\rightarrow \mathbf{r}$ i.e. $U(\mathbf{r},\mathbf{r})<0$, which creates a possibility for a peculiar droplet-formation phase transition.

[*Mapping to a classical lattice-gas problem.*]{} We consider the classical limit of the action (\[eq:eff-action-simplified\]), which bounds of applicability we find below. For a large number of particles $N\gg1$ we can use the canonical ensemble ($N=const$) instead of the grand canonical one ($\mu=const$), and describe the system of atoms by a classical Hamiltonian $H = \sum_i V(\mathbf{r}_i) + \sum_{i ,j} U(\mathbf{r}_i, \mathbf{r}_j)$.

We study the case when the external optical potential $V$ confines all the particles to the plane $z=0$. For the pump mode function chosen as $\Omega(\mathbf{x}) = \Omega_0 (\cos(2\pi \tilde{x}/\lambda)+ \cos(2\pi \tilde{y}/\lambda))$, the external optical potential (\[eq:ExternalPotential\]) becomes $V(\mathbf{x})= -V_0 (\cos(2\pi \tilde{x}/\lambda)+\cos (2\pi \tilde{y}/\lambda))^2$, with $V_0 = \Omega_0^2/\Delta_a^2$. The interaction potential (\[eq:InteractionPotential\]) takes the form $U(\mathbf{x},\mathbf{x}') \sim - (\cos\frac{2\pi \tilde{x}}{\lambda}+\cos\frac{2\pi \tilde{y}}{\lambda}) (\cos\frac{2\pi \tilde{x}'}{\lambda}+\cos\frac{2\pi \tilde{y}'}{\lambda})$, see Fig. 1. The interaction potential $U(\mathbf{x}, \mathbf{x}')$ is sign-alternating and attractive at short distances, as considered in the main text. For $V_0 \gg|U|$ it is convenient to define the lattice version of the interaction potential $U_{i, j}=U(0,i \,\mathbf{a}_1+j  \,\mathbf{a}_2)$ where $\mathbf{a}_1, \mathbf{a}_2$ are the basis vectors of the optical lattice, with $|\mathbf{a}_1|=|\mathbf{a}_2|=a=\lambda/\sqrt{2}$ (Fig. 1a).

Experimental realization and results for the 4-color case. {#Appendix Z4}
==========================================================

In this Supplementary Section we show how to modify the experimental scheme in order to have 4 types of color-charges and present the results of the simulations in this case.

We consider a pair of transverse to the cavity axis lasers with orthogonally-polarized pump modes $\mathbf{\Omega}(\mathbf{r}) =  (\mathbf{e}_y\Omega_1(x) + \mathbf{e}_x\Omega_2(y)) \equiv  (\mathbf{e}_y\cos(2\pi x/\lambda) + \mathbf{e}_x\cos(2\pi y/\lambda))$, with each polarization being coupled to one of the two degenerate excited states $\Psi_{e,1}$, $\Psi_{e,2}$, which scatter the light into two families of orthogonally-polarized cavity modes $a_{\alpha,1}$, $a_{\alpha,2}$ degenerate with respect to polarization. Analogously to the previous section we have

$$\begin{aligned}
S_{at} =  \int d\mathbf{r} d\tau \left[ \Psi_g^*(\mathbf{r},\tau) \left( \partial_{\tau} -\frac{\hbar\nabla^2}{2M} -\frac{\mu}{\hbar} \right) \Psi_g(\mathbf{r},\tau) +
    \sum_{i=1,2}  \Psi_{e,i}^*(\mathbf{r},\tau) \left( \partial_{\tau} -\frac{\hbar\nabla^2}{2M} -\frac{\mu}{\hbar} + \Delta_a \right) \Psi_{e,i}(\mathbf{r},\tau)  \right]\end{aligned}$$

$$\begin{aligned}
S_{int} = \int d\mathbf{r} d\tau \left[ i \sum_{\alpha,i} g_{\alpha}(\mathbf{r}) \Psi^*_{e,i}(\mathbf{r},\tau) \Psi_g(\mathbf{r},\tau) a_{\alpha,i}(\tau) + i\sum_{i=1,2}\Omega_i(\mathbf{r}) \Psi^*_{e,i} (\mathbf{r},\tau) \Psi_g(\mathbf{r},\tau) + h.c.\right]\end{aligned}$$

$$\begin{aligned}
S_{em} = \int d\tau \sum_{\alpha,i} a^*_{\alpha,i}(\tau) (\partial_{\tau}+\Delta_{c\alpha}) a_{\alpha,i}(\tau)\end{aligned}$$

After integrating out excited state $\Psi_e$ and cavity modes $a_{\alpha}$ we get the effective action (\[eq:eff-action-simplified\]) with the external and interaction potentials $$\begin{aligned}
V(\mathbf{r}) &= -\frac{1}{\Delta_a} \left(|\Omega_1(\mathbf{r})|^2 +|\Omega_2(\mathbf{r})|^2 \right) \label{eq:ExternalPotential}\\
U(\mathbf{r},\mathbf{r}') &= -\frac{g_0^2}{\Delta_a^2} \left(\Omega_1^*(\mathbf{r}) \Omega_1(\mathbf{r}') + \Omega_2^*(\mathbf{r}) \Omega_2(\mathbf{r}')\right)\sum_{\alpha}  \frac{\Xi^*_{\alpha}(\mathbf{r}) \Xi_{\alpha}(\mathbf{r}')}{\Delta_{c\alpha}} \label{eq:InteractionPotential}\end{aligned}$$ Periodicity of the external potential $V$ and the sign-changing nature of the interaction $U$ is shown in Fig. \[fig:potential\_orthogonal\_polarizations\]. DW-state here breaks the translational $\mathbb{Z}_2 \times \mathbb{Z}_2$ symmetry.

![Contour plot of the external potential $V(x,y)$ (blue contours) and signs of the interaction potential $U(0,0,x,y)$ ($+,-,0$) for the 4-color case. The contours encircle the local minima of the potential $V$. The interaction potential $U$ favors the formation of density-wave corresponding to one of the four sublattices $A$, $B$, $C$, $D$ (marked by red, green, black, and blue colors) and endows the droplets with $\mathbb{Z}_2 \times \mathbb{Z}_2$ color-charge.[]{data-label="fig:potential_orthogonal_polarizations"}](figAppZ4_potential_orthogonal_polarizations.png){width="5cm"}

Figures \[fig:Z4\_modeling\_plots\]-\[fig:quench\_Z4\] show the results of the numerical simulations for the $\mathbb{Z}_2 \times \mathbb{Z}_2$ case analogous to Figures 3-5 for $\mathbb{Z}_2$ case presented in the main text.

\

[99]{}

F.Y. Wu, *Rev. Mod. Phys.* **54**, 235 (1982).

F. Le Kien, P. Schneeweiss, and A. Rauschenbeutel, Dynamical polarizability of atoms in arbitrary light fields: general theory and application to cesium, *Eur. Phys. J. D* [**67**]{}, 92 (2013).

S. Gopalakrishnan, B. Lev, and P. Goldbart, Emergent crystallinity and frustration with BECs in multimode cavities, *Nature Physics* [**5**]{}, 845 (2009).
