Bose–Einstein-like condensation due to diffusivity edge under periodic confinement

Systems in which detailed balance is broken at the microscopic scale are commonly referred to as active matter [1, 2]. This definition encompasses various processes, which often result in self-propulsion of the microscopic units. When coupled to other mechanisms, activity generally triggers novel physics as it possibly leads to nontrivial types of emergent self-organization [3–5].

One of the many fascinating properties of active systems is their ability to phase separate even in the absence of explicit attractive interactions. A good example of such a feature is the motility-induced phase separation, which emerges when persistent motion is coupled to local motility inhibition [6–10], and can lead to the formation of close-packed ordered structures [11–13]. Dilute systems with short-range velocity alignment also exhibit phase separation at the onset of macroscopic orientational order [14, 15]. Clustering is, moreover, known to arise when the interaction between active particles is induced by a self-generated scalar field (concentration, temperature, etc) [16–21], and from hydrodynamics [11, 22], resulting in effective long-range interactions.

In the absence of long-range orientational order [23], the long-time mean field description of the aforementioned systems is commonly achieved via a conservation law for the density field ρ, with generic drift and diffusion contributions. The effective mobility and diffusion coefficients that result from coarse-graining are then generally explicit functions of ρ. Recently, a new class of scalar active matter was introduced in which the consequences of the existence of a diffusivity edge at a critical concentration ρ_c (i.e. diffusion vanishes when ρ ⩾ ρ_c) was examined [24]. It was discovered that when confined in a harmonic potential, systems falling into this class undergo a transition formally akin to Bose–Einstein condensation (BEC), thus providing a new non-equilibrium mechanism for the emergence of clustering.

There are many examples in which structure formation in active matter results in a pattern formation that involves the selection of a characteristic length-scale, such that the cluster sizes are limited and do not scale with the system size [17–19, 22, 25]. Such microscopic confinement can be modelled at the mean field level by introducing an effective potential which provides multiple sites for condensation. Moreover, periodic potential landscapes are commonly used to manipulate driven colloidal systems. Such periodic potentials are expected to lead to cluster-lattices. Here, we study the phenomenology that arises from a diffusivity edge in such configurations.

We consider a sinusoidal egg-crate confinement in arbitrary dimension d, and identify two limiting regimes for the system. For deep potentials, the system behaves similarly to the case of a single harmonic trap case treated in reference [24]. However, we find that the existence of finite energy barriers between neighbouring condensates quantitatively modifies the transition. A generalized thermodynamic description shows that the overall phenomenology of BEC is always preserved. However, for shallow potentials we observe quantitative differences as compared to the classical BEC description. Most notably, we find that the exponent associated with the scaling of the condensate fraction with respect to an effective temperature is non-universal, and depends on how the diffusion scales with ρ_c − ρ.

The rest of the paper is organized as follows. We introduce the model in section 2 and characterize the phenomenology of the condensation transition in section 3. Section 4 is devoted to the development of the generalized thermodynamics associated with the phenomenology of the system, section 5 deals with the possibility of observing degenerate solutions and examines their stability, and section 6 concludes the paper.

We start by introducing the formalism that will be used throughout the paper. In the mean field approach considered here, the particle density field $\rho \left(\boldsymbol{r},t\right)$ obeys the following conservation law

$$\begin{equation}{\partial }_{t}\rho +\nabla \cdot \boldsymbol{J}=0,\quad \boldsymbol{J}=-M\left(\rho \right)\rho \nabla U-D\left(\rho \right)\nabla \rho ,\end{equation} \tag{ 1 }$$

where U(r) denotes the external confining potential. The dynamics conserves the total number of particles $N=\int \;{\mathrm{d}}^{d}\boldsymbol{r}\;\rho \left(\boldsymbol{r},t\right)$ in the accessible d-dimensional space at all times. The mobility M(ρ) and diffusion coefficient D(ρ) are in general density-dependent. Their ratio in the zero-density limit defines a tuning parameter

$$\begin{equation}{k}_{\text{B}}{T}_{\text{eff}}\equiv \frac{D\left(\rho \to 0\right)}{M\left(\rho \to 0\right)},\end{equation} \tag{ 2 }$$

which gives a measure of the fluctuations at the particle level, and can be assimilated to an effective temperature for the system. Because this study aims at describing systems that are non-equilibrium in essence, the fluctuation-dissipation theorem (FDT) can be broken for finite densities, namely,

$$\begin{equation}\frac{D\left(\rho \right)}{M\left(\rho \right)}\ne \frac{D\left(\rho \to 0\right)}{M\left(\rho \to 0\right)}.\end{equation} \tag{ 3 }$$

This feature can be interpreted as collective inhibition or activation caused by the interplay of activity and, for instance, interactions. In particular, for sufficiently large densities we assume the existence of a diffusivity edge in the system, defined as D(ρ)/M(ρ) = 0 for ρ ⩾ ρ_c. The non-local effects due to hydrodynamic interactions in the presence of broken FDT are neglected in our work [26].

The steady-state solutions of equation (1) are computed by setting the current to zero (J = 0), leading to

$$\begin{equation}\frac{\mathrm{d}U}{\mathrm{d}\rho }=-\frac{D\left(\rho \right)}{M\left(\rho \right)\rho },\end{equation} \tag{ 4 }$$

which can be formally used to obtain ρ(U). The normalization condition in the stationary state can be written as $N=\int \;{\mathrm{d}}^{d}\boldsymbol{r}\;\rho \left(U\left(\boldsymbol{r}\right)\right)=\int \;\mathrm{d}U\;g\left(U\right)\rho \left(U\right)$, where g(U) is the relevant density of states.

Since D(ρ)/M(ρ) ⩾ 0, we surmise that ρ is a decreasing function of U. We denote ρ₀ as the maximal value that ρ takes in the ground state U = 0. When ρ₀ < ρ_c, ρ(U) can be obtained simply by integrating and inverting equation (4), in which case ρ₀ is then determined from the density normalization. When the effective temperature decreases, ρ₀ increases until it reaches the maximally allowed value of ρ_c. The transition temperature T_c is defined as the value taken by T_eff when ρ₀ = ρ_c. For T_eff ⩽ T_c, ρ is thus not a smooth function at U = 0 (see figure 1) we obtain ρ(U → 0⁺) = ρ_c and the value ρ takes in the ground state is undefined, which reflects the formation of a condensate.

Zoom In Zoom Out Reset image size

In most of this work we will consider for simplicity a step profile for the ratio of diffusion over mobility

$$\begin{equation}\frac{D\left(\rho \right)}{M\left(\rho \right)}=\begin{cases}_{\text{B}}{T}_{\text{eff}}\quad \hfill & \rho {< }{\rho }_{\mathrm{c}}\hfill \\ 0\quad \hfill & \rho {\geqslant}{\rho }_{\mathrm{c}}\hfill \end{cases}.\end{equation} \tag{ 5 }$$

Therefore, denoting β ≡ 1/k_BT_eff for convenience, the density is given by the Boltzmann weights

$$\begin{equation}\rho \left(U\right)=\begin{cases}_{0}\;\mathrm{exp}\left(-\beta U\right)\quad \quad \hfill & {T}_{\text{eff}}{ >}{T}_{\mathrm{c}}\hfill \\ {N}_{\mathrm{c}}\delta \left(U\right)g{\left(U\right)}^{-1}+{\rho }_{\mathrm{c}}\;\mathrm{exp}\left(-\beta U\right)\quad \hfill & {T}_{\text{eff}}{\leqslant}{T}_{\mathrm{c}}\hfill \end{cases},\end{equation} \tag{ 6 }$$

where the contribution N_cδ(U)g(U)⁻¹ ensures the overall normalization in the condensed phase.

In this work we consider the sinusoidal potential in d dimensions, as sketched in figure 1 for d = 2, and defined by

$$\begin{equation}U\left(\boldsymbol{r}\right)=\sum _{k=1}^{d}\tilde{U}\left({r}_{k}\right),\quad \tilde{U}\left({r}_{k}\right)=\frac{{U}_{\mathrm{b}}}{2d}\left[1-\mathrm{cos}\left(\frac{\pi {r}_{k}}{{r}_{\mathrm{b}}}\right)\right],\end{equation} \tag{ 7 }$$

where the r_k's denotes the Cartesian coordinates of the position r. The system is thus divided into identical cells of volume ${\left(2{r}_{\mathrm{b}}\right)}^{d}$, each separated by an energy barrier U_b. Assuming an equal partitioning of the particles over the cells, and using the fact that the density can be factorized (as can be seen by combining equations (6) and (7)), we find

$$\begin{equation}N=\begin{cases}_{0}{L}^{d}\;\mathrm{exp}\left(-\frac{\beta {U}_{\mathrm{b}}}{2}\right){I}_{0}^{d}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)\quad \hfill & {T}_{\text{eff}}{ >}{T}_{\mathrm{c}}\hfill \\ {N}_{\mathrm{c}}+{\rho }_{\mathrm{c}}{L}^{d}\;\mathrm{exp}\left(-\frac{\beta {U}_{\mathrm{b}}}{2}\right){I}_{0}^{d}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)\quad \hfill & {T}_{\text{eff}}{\leqslant}{T}_{\mathrm{c}}\hfill \end{cases},\end{equation} \tag{ 8 }$$

where L denotes the linear system size and ${I}_{\nu }\left(x\right)={\int }_{0}^{\pi }\mathrm{d}s\;\mathrm{exp}\left[x\;\mathrm{cos}\left(s\right)\right]\mathrm{cos}\left(\nu s\right)/\pi$ is the modified Bessel function of the first kind of integer rank ν, which has the following asymptotic forms

$$\begin{equation}\;\;\;\;\;\;\;\;\;\;\;\;{I}_{\nu }\left(x\right)x\to \infty {\sim }\frac{{\mathrm{e}}^{x}}{\sqrt{2\pi x}}\left(1+\mathrm{O}\left({x}^{-1}\right)\right),\quad {I}_{\nu }\left(x\right)\underset{x\to 0}{\sim }\frac{1}{\nu !}{\left(\frac{x}{2}\right)}^{\nu }\left(1+\mathrm{O}\left({x}^{2}\right)\right).\end{equation} \tag{ 9 }$$

Hence, in the strong and weak confinement limits the ground state density below the diffusivity edge obeys

$$\begin{equation}{\rho }_{0}\underset{{k}_{\text{B}}{T}_{\text{eff}}\ll {U}_{\mathrm{b}}}{\sim }n{\left(\frac{2\pi {k}_{\text{B}}{T}_{\text{eff}}}{k}\right)}^{-\frac{d}{2}},\quad {\rho }_{0}\underset{{k}_{\text{B}}{T}_{\text{eff}}\gg {U}_{\mathrm{b}}}{\sim }\frac{N}{{L}^{d}}\left(1+\mathrm{O}\left(\beta {U}_{\mathrm{b}}\right)\right),\end{equation} \tag{ 10 }$$

where $n=N{\left(2{r}_{\mathrm{b}}/L\right)}^{d}$ denotes the number of particles in each cell, and $k={\pi }^{2}{U}_{\mathrm{b}}/\left(2d{r}_{\mathrm{b}}^{2}\right)$ measures the effective potential stiffness in the ground state. Because edge effects vanish when the barrier height U_b is much larger than the effective temperature, the expression given in (10) for k_BT_eff ≪ U_b is identical to the one derived in reference [24] for an infinite harmonic trap. On the other hand, when k_BT_eff ≫ U_b the system is dominated by fluctuations and the density reaches a uniform profile.

When ρ₀ becomes larger than the diffusivity edge, some of the particles form a condensate in the ground state. In this case, the normalization of the density profile is given by the second line of equation (8). In the high energy barrier limit, the condensate fraction can be approximated by

$$\begin{equation}\frac{{N}_{\mathrm{c}}}{N}\underset{{k}_{\text{B}}{T}_{\text{eff}}\ll {U}_{\mathrm{b}}}{\sim }1-{\left(\frac{{T}_{\text{eff}}}{{T}_{\mathrm{c}}^{0}}\right)}^{\frac{d}{2}},\end{equation} \tag{ 11 }$$

where the effective transition temperature reads ${T}_{\mathrm{c}}^{0}=\frac{k}{2\pi {k}_{\text{B}}}{\left(n/{\rho }_{\mathrm{c}}\right)}^{\frac{2}{d}}$ [24]. In this limit, N_c/N takes a similar form as in the case of a free ideal Bose gas [27]. Defining $\overline{\rho }\equiv N/{L}^{d}$ as the average density of particles, the condensate fraction in the shallow potential limit reads

$$\begin{equation}\frac{{N}_{\mathrm{c}}}{N}\underset{{k}_{\text{B}}{T}_{\text{eff}}\gg {U}_{\mathrm{b}}}{\sim }\left(1-\frac{{\rho }_{\mathrm{c}}}{\overline{\rho }}\right)\left(1-\frac{{T}_{\mathrm{c}}^{\infty }}{{T}_{\text{eff}}}\right),\end{equation} \tag{ 12 }$$

where ${T}_{\mathrm{c}}^{\infty }={U}_{\mathrm{b}}{\left[2{k}_{\text{B}}\left(1-\overline{\rho }/{\rho }_{\mathrm{c}}\right)\right]}^{-1}$. The fact that the transition temperature diverges when $\overline{\rho }$ approaches ρ_c is due to the finiteness of the barrier height U_b, which leads to flat density profiles at high effective temperatures. As shown in the phase diagram of figure 2(a), the phase behaviour for k_BT_eff/U_b ≫ 1 is only set by the ratio $\overline{\rho }/{\rho }_{\mathrm{c}}$.

Zoom In Zoom Out Reset image size

The scaling of the condensate fraction as a function of k_BT_eff/U_b is shown in figure 2(b). If $\overline{\rho }{ >}{\rho }_{\mathrm{c}}$, the transition is suppressed, and N_c/N goes from 1 at vanishing T_eff to finite values when k_BT_eff/U_b → ∞. When $\overline{\rho }{< }{\rho }_{\mathrm{c}}$, there exists a finite effective temperature T_c for which N_c/N reaches 0, and above which no condensation occurs. If $\overline{\rho }\ll {\rho }_{\mathrm{c}}$, the transition happens at small effective temperatures and is similar to BEC. On the other hand, for ${\rho }_{\mathrm{c}}\gtrsim \overline{\rho }$ the exponent associated to the scaling of the condensate fraction is equal to −1 in any dimension (see equation (12)). Finally, in the particular case of $\overline{\rho }={\rho }_{\mathrm{c}}$, we find that ${N}_{\mathrm{c}}/N\sim {\left({k}_{\text{B}}{T}_{\text{eff}}/{U}_{\mathrm{b}}\right)}^{-1}$ at large k_BT_eff/U_b, such that the transition temperature is located exactly at infinity.

Although deriving an analogue to equation (8) for arbitrary functions D(ρ)/M(ρ) is out of the scope of this work, A shows how some progress can be achieved in the limit k_BT_eff ≫ U_b. Indeed, assuming that the diffusivity edge is approached as

$$\begin{equation}D\left(\rho \right)/M\left(\rho \right)\underset{\rho \to {\rho }_{\mathrm{c}}}{\sim }{\left(1-\rho /{\rho }_{\mathrm{c}}\right)}^{z-1},\end{equation} \tag{ 13 }$$

where z ⩾ 1, the associated scaling of the condensate fraction reads

$$\begin{equation}\frac{{N}_{\mathrm{c}}}{N}\underset{{k}_{\text{B}}{T}_{\text{eff}}\gg {U}_{\mathrm{b}}}{\sim }\left(1-\frac{{\rho }_{\mathrm{c}}}{\overline{\rho }}\right)\left[1-{\left(\frac{{T}_{\mathrm{c}}^{\infty }}{{T}_{\text{eff}}}\right)}^{\frac{1}{z}}\right],\end{equation} \tag{ 14 }$$

with ${k}_{\text{B}}{T}_{\mathrm{c}}^{\infty }/{U}_{\mathrm{b}}\sim {\left(1-\overline{\rho }/{\rho }_{\mathrm{c}}\right)}^{-z}$. For shallow potentials the condensate fraction exponent therefore takes a nonuniversal value, which is set by how fast the diffusivity edge is reached. This result is in clear departure from equation (11) for k_BT_eff ≪ U_b, where the exponent $\frac{d}{2}$ remains independent of the shape of D(ρ)/M(ρ) [24].

We now turn to the construction of a generalized thermodynamic formalism for the system. The average potential energy $\langle U\rangle \equiv \int \;{\mathrm{d}}^{d}\boldsymbol{r}\;U\left(\boldsymbol{r}\right)\rho \left(U\left(\boldsymbol{r}\right)\right)$ reads

$$\begin{equation}\langle U\rangle =\begin{cases}\frac{N{U}_{\mathrm{b}}}{2}\left(1-\frac{{I}_{1}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}{{I}_{0}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}\right)\quad \hfill & {T}_{\text{eff}}{ >}{T}_{\mathrm{c}}\hfill \\ \frac{\left(N-{N}_{\mathrm{c}}\right){U}_{\mathrm{b}}}{2}\left(1-\frac{{I}_{1}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}{{I}_{0}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}\right)\quad \hfill & {T}_{\text{eff}}{\leqslant}{T}_{\mathrm{c}}\hfill \end{cases}.\end{equation} \tag{ 15 }$$

A heat capacity can then be defined from the mean energy via C ≡ d⟨U⟩/dT_eff. For the present system, we find the following expressions after some algebra

$$\begin{equation}C=\begin{cases}\frac{N{k}_{\text{B}}}{4d}{\left(\beta {U}_{\mathrm{b}}\right)}^{2}\left[1-\frac{{I}_{1}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}{{I}_{0}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}\left(\frac{2d}{\beta {U}_{\mathrm{b}}}+\frac{{I}_{1}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}{{I}_{0}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}\right)\right]\quad \hfill & {T}_{\text{eff}}{ >}{T}_{\mathrm{c}}\hfill \\ \frac{\left(N-{N}_{\mathrm{c}}\right){k}_{\text{B}}}{4d}{\left(\beta {U}_{\mathrm{b}}\right)}^{2}\quad \hfill & \hfill \\ {\times}\left\{d+1-\frac{{I}_{1}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}{{I}_{0}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}\left[2d\left(1+\frac{1}{\beta {U}_{\mathrm{b}}}\right)-\left(d-1\right)\frac{{I}_{1}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}{{I}_{0}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}\right]\right\}\quad \hfill & {T}_{\text{eff}}{\leqslant}{T}_{\mathrm{c}}\hfill \end{cases}.\end{equation} \tag{ 16 }$$

The change in the heat capacity at the transition, ${\Delta}C\equiv C\left(T={T}_{\mathrm{c}}^{-}\right)-C\left(T={T}_{\mathrm{c}}^{+}\right)$ is then given by

$$\begin{equation}{\Delta}C=\frac{N{k}_{\text{B}}}{4}{\left({\beta }_{\mathrm{c}}{U}_{\mathrm{b}}\right)}^{2}\left[1+\frac{{I}_{1}\left(\frac{{\beta }_{\mathrm{c}}{U}_{\mathrm{b}}}{2d}\right)}{{I}_{0}\left(\frac{{\beta }_{\mathrm{c}}{U}_{\mathrm{b}}}{2d}\right)}\left(\frac{{I}_{1}\left(\frac{{\beta }_{\mathrm{c}}{U}_{\mathrm{b}}}{2d}\right)}{{I}_{0}\left(\frac{{\beta }_{\mathrm{c}}{U}_{\mathrm{b}}}{2d}\right)}-2\right)\right],\end{equation} \tag{ 17 }$$

with ${\beta }_{\mathrm{c}}\equiv {\left({k}_{\text{B}}{T}_{\mathrm{c}}\right)}^{-1}$. Generally, ΔC is nonzero such that the heat capacity experiences a discontinuous jump at the transition (see figure 3(b)). For BEC in free space, this feature appears only for d ⩾ 5 [27], while it can be affected by confinement [28, 29]. Similar features are expected for the diffusivity edge problem, where the shape of D(ρ)/M(ρ) in the vicinity of ρ_c could additionally play a role. These questions will be addressed in a separate publication [30].

Zoom In Zoom Out Reset image size

Using the asymptotic expansions of the modified Bessel functions (9), the analytical expressions for ⟨U⟩ and C can be obtained in the strong and weak confinement limits. As shown in figures 3(a) and (b), for k_BT_eff ≪ U_b their behaviour corresponds to that of an ideal Bose gas [24], while for k_BT_eff ≫ U_b, we obtain

$$\begin{equation}\langle U\rangle \underset{{k}_{\text{B}}{T}_{\text{eff}}\gg {U}_{\mathrm{b}}}{\sim }\frac{N{U}_{\mathrm{b}}}{2},\quad C\underset{{k}_{\text{B}}{T}_{\text{eff}}\gg {U}_{\mathrm{b}}}{\sim }\frac{N{k}_{\text{B}}}{8d}{\left(\beta {U}_{\mathrm{b}}\right)}^{2}{T}_{\text{eff}}{ >}{T}_{\mathrm{c}},\end{equation} \tag{ 18 }$$

$$\begin{equation}\langle U\rangle \underset{{k}_{\text{B}}{T}_{\text{eff}}\gg {U}_{\mathrm{b}}}{\sim }\frac{{\rho }_{\mathrm{c}}}{\overline{\rho }}\frac{N{U}_{\mathrm{b}}}{2},\quad C\underset{{k}_{\text{B}}{T}_{\text{eff}}\gg {U}_{\mathrm{b}}}{\sim }\left(d+\frac{1}{2}\right)\frac{{\rho }_{\mathrm{c}}}{\overline{\rho }}\frac{N{k}_{\text{B}}}{4d}{\left(\beta {U}_{\mathrm{b}}\right)}^{2}\quad {T}_{\text{eff}}{\leqslant}{T}_{\mathrm{c}}.\end{equation} \tag{ 19 }$$

In the limit of a shallow potential and a high effective temperature, ⟨U⟩ becomes independent of T_eff and scales linearly with U_b. The heat capacity thus vanishes as ${\left(\beta {U}_{\mathrm{b}}\right)}^{2}$. Note that the functions below T_c are proportional to the ratio ${\rho }_{\mathrm{c}}/\overline{\rho }$, which highlights the fact that particles in the condensate do not contribute to the total energy.

A thermodynamic entropy can be defined for the system as dS ≡ d⟨U⟩/T_eff. After some algebra, we find the following expressions

$$\begin{equation}\;\;\;\;\;\;\;\;\;\;\;\;S=\begin{cases}\frac{N{k}_{\text{B}}}{2}\left\{2\left[1-\mathrm{ln}\left(\frac{{\rho }_{0}}{{\rho }_{\mathrm{c}}}\right)\right]+\beta {U}_{\mathrm{b}}\left[1-\frac{{I}_{1}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}{{I}_{0}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}\right]\right\}\quad \hfill & {T}_{\text{eff}}{ >}{T}_{\mathrm{c}}\hfill \\ \frac{\left(N-{N}_{\mathrm{c}}\right){k}_{B}}{2}\left[2+\beta {U}_{\mathrm{b}}\left(1-\frac{{I}_{1}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}{{I}_{0}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right)}\right)\right]\quad \hfill & {T}_{\text{eff}}{\leqslant}{T}_{\mathrm{c}}\hfill \end{cases}.\end{equation} \tag{ 20 }$$

We note that the same result can be derived from a Gibbs definition of generalized entropy, which is consistent with equation (4) as the relationship between the energy and the probability measure. In case with a large energy barrier, we find that the entropy exhibits an ideal Bose gas behaviour [24]. On the other hand, for k_BT_eff ≫ U_b the energy states are distributed uniformly in space, and S/L^d ≃ −k_Bρ[ln(ρ/ρ_c) − 1], with distributions $\rho =\overline{\rho }\;\left({T}_{\text{eff}}{ >}{T}_{\mathrm{c}}\right)$ and ρ = ρ_c(T_eff ⩽ T_c).

A remarkable feature of BEC concerns the divergence of the isothermal compressibility at the transition. A thermodynamic pressure can be defined for the system from a generalized Helmholtz free energy $\mathcal{F}\equiv \langle U\rangle -{T}_{\text{eff}}S$. The typical volume $\mathcal{V}$ of the confined system can be obtained from dimensional analysis: $\mathcal{V}\equiv N/{\rho }_{0}=\left(N-{N}_{\mathrm{c}}\right)/{\rho }_{\mathrm{c}}$. Using equation (8), this reads⁴

$$\begin{equation}\mathcal{V}={L}^{d}\;\mathrm{exp}\left(-\frac{\beta {U}_{\mathrm{b}}}{2}\right){I}_{0}^{d}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right).\end{equation} \tag{ 21 }$$

When k_BT_eff ≪ U_b, we find $\mathcal{V}\sim {\lambda }^{d}={\left(2\pi {k}_{\text{B}}{T}_{\text{eff}}/k\right)}^{\frac{d}{2}}$, which corresponds to the typical volume occupied by an ideal Bose gas confined in a harmonic potential [31], while for k_BT_eff ≫ U_b, it is simply given by the size of the system L^d.

Defining P_therm as the conjugate variable to $\mathcal{V}$ provides the following equations of state

$$\begin{equation}{P}_{\text{therm}}\equiv -{\left.\left(\frac{\partial \mathcal{F}}{\partial \mathcal{V}}\right)\right\vert }_{{T}_{\text{eff}},N}=\begin{cases}\frac{N{k}_{\text{B}}{T}_{\text{eff}}}{\mathcal{V}}\quad \hfill & {T}_{\text{eff}}{ >}{T}_{\mathrm{c}}\hfill \\ {\rho }_{\mathrm{c}}{k}_{\text{B}}{T}_{\text{eff}}\quad \hfill & {T}_{\text{eff}}{\leqslant}{T}_{\mathrm{c}}\hfill \end{cases}.\end{equation} \tag{ 22 }$$

Equations (22) are identical to those derived in [24], and stress again the similarities with BEC. These expressions can moreover be derived from a mechanical definition of the pressure, which we denote as P_mech. Thanks to the periodicity of the potential U, the forces that it induces do not create pressure difference between adjacent unit cells, such that a bulk can be defined for this system. Following previous works [32], the bulk mechanical pressure is defined as the average force per unit surface area exerted by the particles on a potential W which confines the system in a finite volume L^d. From the symmetries of the problem, the calculation is moreover carried out in one dimension. Assuming that the edge of the system corresponds to a dip of U⁵, W(r) is monotonously growing with r, satisfies W(r) = 0 for r ⩽ L and W(r → ∞) → ∞ (it is therefore assumed that U = 0 outside the sample). This way, P_mech reads

$$\begin{equation}{P}_{\text{mech}}={\int }_{L}^{\infty }\mathrm{d}r\;\rho \left(W\left(r\right)\right){\partial }_{r}W\left(r\right)={\int }_{0}^{\infty }\;\mathrm{d}W\;\rho \left(W\right),\end{equation} \tag{ 23 }$$

which after replacing ρ by its expression as function of the potential (6), leads to P_mech = P_therm ≡ P. As at equilibrium, and in a limited number of nonequilibrium cases [32], this result is moreover independent of the details of W.

Let us consider a typical isotherm of P as shown in figure 3(c). Defining ${\mathcal{V}}_{\mathrm{c}}\equiv N/{\rho }_{\mathrm{c}}$ as the volume below which condensation occurs, we find that for $\mathcal{V}{ >}{\mathcal{V}}_{\mathrm{c}}$ the system possesses an ideal gas equation of state and P scales like ${\mathcal{V}}^{-1}$. For $\mathcal{V}{\leqslant}{\mathcal{V}}_{\mathrm{c}}$ the pressure becomes independent of $\mathcal{V}$ and the corresponding isotherm exhibits a plateau, such that the isothermal compressibility of the system, ${\kappa }_{{T}_{\text{eff}}}=-{\mathcal{V}}^{-1}{\left(\frac{\partial P}{\partial \mathcal{V}}\right)}_{{T}_{\text{eff}}}^{-1}$, diverges at the threshold.

We end this section by computing the generalized chemical potential μ, defined as the conjugate variable to N. From equations (15) and (20), we find

$$\begin{equation}\mu \equiv {\left.\left(\frac{\partial \mathcal{F}}{\partial N}\right)\right\vert }_{{T}_{\text{eff}},\mathcal{V}}=\begin{cases}_{\text{B}}{T}_{\text{eff}}\;\mathrm{ln}\left(\frac{{\rho }_{0}}{{\rho }_{\mathrm{c}}}\right)\quad \hfill & {T}_{\text{eff}}{ >}{T}_{\mathrm{c}}\hfill \\ 0\quad \hfill & {T}_{\text{eff}}{\leqslant}{T}_{\mathrm{c}}\hfill \end{cases}.\end{equation} \tag{ 24 }$$

We thus find that μ vanishes at the transition where ρ₀ = ρ_c and remains identically 0 in the condensate phase, as shown in figure 3(d). From equations (6) and (24), the density profile outside the ground state thus takes the general form

$$\begin{equation}\rho \left(U{ >}0\right)={\rho }_{\mathrm{c}}\;\mathrm{exp}\left[\beta \left(\mu -U\right)\right],\end{equation} \tag{ 25 }$$

above and below the transition.

The derivation of the above results relies on the assumption that the particles are divided evenly among cells of volume ${\left(2{r}_{\mathrm{b}}\right)}^{d}$. However, a vanishing diffusion coefficient may seem 'pathological' at the mean field level considered here, as the absence of fluctuations could potentially lead to atypical behaviour such as spontaneous symmetry breaking among the cells. It should thus be emphasized that the range of validity of the results derived here concerns all systems described by equation (1) for which the density-dependent hydrodynamic coefficients result from the integration of various microscopic processes, e.g. interactions, while fluctuations, albeit possibly weak, remain present. In this section, we show that the steady state solution characterized so far is in fact degenerate with respect to how condensed particles are distributed within the system, and is always stable to perturbations.

We first note from equation (6) that the density profile in the steady state is uniquely determined by the potential landscape. We consider a system made of two identical wells for simplicity, and denote their respective density profiles by ρ₁(U) and ρ₂(U). This system can be in three distinct states corresponding to zero, one, or two condensates. As long as the potential barrier U_b between the wells remains finite, the continuity of the density profile at their boundary imposes ρ₁(U_b) = ρ₂(U_b). Without condensation, the previous condition implies that the two ground state densities are equal: ρ_0,1 = ρ_0,2. If condensation occurs on both sides, this equality will automatically be satisfied (as the densities will both be equal to ρ_c). In the intermediate case where only one of the two wells, say 1, exhibits condensation, we will have ρ_0,2 = ρ_c. This argument is generalizable in a straightforward way to an arbitrary number of minima. Therefore, the only steady state solutions admitted by equation (1) in the presence of periodic confinement are those for which the density profile outside the condensate is identical for all cells.

For each cell i, let us denote the number of particles in the condensed and in the excited states as n_c,i and n_e,i respectively. From the preceding discussion, we conclude that n_e,i's are all equal for every value of 1 ⩽ i ⩽ M, where $M={\left(L/2{r}_{\mathrm{b}}\right)}^{d}$ is the number of cells. Moreover, the previous argument also implies that in order to observe steady state condensation in one or more minima, the total number of particles N must be sufficiently large such that the ground state densities of all the other wells are at least equal to ρ_c. In this case, using equation (8) and the particle number conservation, we obtain

$$\begin{equation}{N}_{\mathrm{c}}=\sum _{i=1}^{M}{n}_{\mathrm{c},i}=N-M{n}_{\mathrm{e}}=N-{\rho }_{\mathrm{c}}{L}^{d}\;\mathrm{exp}\left(-\frac{\beta {U}_{\mathrm{b}}}{2}\right){I}_{0}^{d}\left(\frac{\beta {U}_{\mathrm{b}}}{2d}\right).\end{equation} \tag{ 26 }$$

All the ${M}^{{N}_{\mathrm{c}}}$ configurations ${\left\{{n}_{\mathrm{c},i}\right\}}_{i}$ that satisfy this identity thus share the same total number of condensed particles N_c, and consequently the same global transition temperature T_c. Moreover, as shown in section 4 the particles in the condensates do not contribute to the global thermodynamic properties of the system (see equations (15), (16), (20)–(22), and (24)). All the configurations ${\left\{{n}_{\mathrm{c},i}\right\}}_{i}$ are thus thermodynamically equivalent to the particular case addressed previously where the distribution of particles among the cells was assumed to be uniform. Hence, the steady state solution of equation (1) shows degeneracy in the condensed phase.

We now investigate the stability of these solutions. The generalized free energy of the system introduced in section 4 can equivalently be written as a functional of the density as

$$\begin{align}\hfill \mathcal{F}\left[\rho \right]=& \int \;{\mathrm{d}}^{d}\boldsymbol{r}\;\left\{\rho \left(U\left(\boldsymbol{r}\right)\right)U\left(\boldsymbol{r}\right)-{k}_{\text{B}}{T}_{\text{eff}}{\rho }_{\mathrm{c}}{\Theta}\left(\rho \left(U\left(\boldsymbol{r}\right)\right)-{\rho }_{\mathrm{c}}\right)\right.\hfill \\ \hfill & \left.+{k}_{\text{B}}{T}_{\text{eff}}\rho \left(U\left(\boldsymbol{r}\right)\right)\left[\mathrm{ln}\left(\rho \left(U\left(\boldsymbol{r}\right)\right)/{\rho }_{\mathrm{c}}\right)-1\right]{\Theta}\left({\rho }_{\mathrm{c}}-\rho \left(U\left(\boldsymbol{r}\right)\right)\right)\right\},\hfill \end{align} \tag{ 27 }$$

where Θ(x) denotes the Heaviside function that is 1 for x > 0 and 0 for x ⩽ 0. This definition is consistent with what is introduced in reference [33] in a more general context, where it is shown that $\mathcal{F}$ exhibits similar properties to its equilibrium counterpart. In particular, we note that $\mathcal{F}$ obeys an H-theorem ($\mathrm{d}\mathcal{F}/\;\mathrm{d}t{\leqslant}0$), and stable steady state solutions of equation (1) correspond to the minima of $\mathcal{F}$ associated with constant particle number N. In appendix B, we show that solving $\delta \left(\mathcal{F}-\mu N\right)/\delta \rho =0$ leads to equation (25). Moreover, in steady state we obtain

$$\begin{equation}\frac{{\delta }^{2}\left(\mathcal{F}-\mu N\right)}{\delta \rho \left(\boldsymbol{r},t\right)\delta \rho \left({\boldsymbol{r}}^{\prime },t\right)}=\frac{{k}_{\text{B}}{T}_{\text{eff}}}{\rho \left(U\left(\boldsymbol{r}\right)\right)}\;{\Theta}\left({\rho }_{\mathrm{c}}-\rho \left(U\left(\boldsymbol{r}\right)\right)\right)\;{\delta }^{d}\left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right).\end{equation} \tag{ 28 }$$

In absence of condensation, the density ρ is lower than ρ_c everywhere such that $\mathcal{F}-\mu N$ is always strictly convex in ρ, and the corresponding solution is stable. Interestingly, as a consequence of the diffusivity edge when ρ → ρ_c, we observe that equation (28) vanishes in the condensate. This feature means that condensation, as well as the distribution of condensed particles among the cells, do not affect the stability of the overall solution.

We have studied the consequences of having a diffusivity edge for a system of particles embedded in a sinusoidal potential in arbitrary dimensions. This configuration leads to the formation multiple coexisting condensates. We have identified two asymptotic regimes that exhibit qualitatively different properties. At low effective temperatures or high potential barriers, the behaviour of the system is analogous to that of an ideal Bose gas in free space, similarly to the case treated in reference [24] considering a single harmonic trap. For shallow potentials, qualitative features such as the presence of a transition, at which the heat capacity is discontinuous, as well as the divergence of the isothermal compressibility and the vanishing of the chemical potential in the condensed phase, persist when the mean density stays lower than the threshold ρ_c.

We have, however, uncovered quantitative differences is this case. For example, we have found that for $D\left(\rho \right)/M\left(\rho \right)\sim {\left({\rho }_{\mathrm{c}}-\rho \right)}^{z-1}$ near ρ_c the scaling of the condensate fraction with the effective temperature takes an exponent of −z⁻¹ (see equation (14)). The exponent z is also expected to affect the scaling of other functions. For instance, using the results derived in appendix A, it is possible to show that $C\sim {\left(\beta {U}_{\mathrm{b}}\right)}^{1+{z}^{-1}}$ as βU_b → 0. A systematic study of the effect of z and of the shape of the potential is relegated to future publications [30].

In section 5, we have finally shown that the steady state solution of equation (1) is degenerate in the condensed phase, as multiple configurations sharing the same macroscopic properties but showing different distributions of condensed particles exist and are stable. This last result has interesting consequences, as it allows the system to show a certain degree of memory since the steady state solution may depend to some extent on the initial particle configuration. This mean field picture does not however explicitly take into account density fluctuations, which may be large enough to restore the symmetry among cells. These enthralling questions will be addressed in a forthcoming publication as well [30].

This section is devoted to the derivation of equation (14) describing the weak confinement behaviour of the condensate fraction assuming a general functional form of D(ρ)/M(ρ) near ρ_c. Below T_c, the potential can be formally written as βU(ρ) ≡ u(ρ/ρ_c), where this rescaled form satisfies

$$\begin{equation}u\left(s\right)=-{\int }_{1}^{s}\;\mathrm{d}t\;\frac{y\left(t\right)}{t}\end{equation} \tag{ A.1 }$$

with y(ρ/ρ_c) ≡ βD(ρ)/M(ρ). Denoting ρ_b ≡ ρ(U_b), from equation (A.1) the limit u(ρ_b/ρ_c) = βU_b → 0 is attained for ρ_b → ρ_c. In the following, we assume that the diffusivity edge is reached at ρ = ρ_c following a power law with an exponent z − 1 ⩾ 0:

$$\begin{equation}y\left(s\right)\underset{s\to {1}^{-}}{\sim }{y}_{0}{\left(1-s\right)}^{z-1},\end{equation} \tag{ A.2 }$$

where y₀ is a constant and y(s) = 0 for all s ⩾ 1. The Taylor expansion of u(s) in s = 1 reads

$$\begin{equation}u\left(s\right)=\sum _{n=0}^{\infty }\frac{{u}^{\left(n+1\right)}\left(1\right)}{\left(n+1\right)!}{\left(s-1\right)}^{n+1},\end{equation} \tag{ A.3 }$$

where u⁽ⁿ⁺¹⁾ stands for the (n + 1)^th derivative of u, and is obtained from equation (A.1) as

$$\begin{equation}{u}^{\left(n+1\right)}\left(s\right)=-\sum _{p=0}^{n}\left(\genfrac{}{}{0pt}{}{n}{p}\right){y}^{\left(p\right)}\left(s\right){\left({s}^{-1}\right)}^{\left(n-p\right)},\end{equation} \tag{ A.4 }$$

$$\begin{equation}{u}^{\left(n+1\right)}\left(s\right)\underset{s\to {1}^{-}}{\sim }{y}_{0}\sum _{p=0}^{\mathrm{min}\left(z-1,n\right)}\left(\genfrac{}{}{0pt}{}{z-1}{p}\right)n!{\left(-1\right)}^{n+1}{\left(1-s\right)}^{z-1-p}.\end{equation} \tag{ A.5 }$$

It is clear from equation (A.5) that u⁽ⁿ⁺¹⁾(s) will cancel when s → 1⁻ for all n < z − 1, and that u^(z+i)(1) = y₀(z + i − 1)!(−1)^z+i for all i ⩾ 0. Inserting this expression in equation (A.3), we get

$$\begin{equation}u\left(s\right)={y}_{0}\sum _{n=0}^{\infty }\frac{{\left(1-s\right)}^{z+n}}{z+n},\end{equation} \tag{ A.6 }$$

which corresponds as expected to u(s) = −y₀ ln(s) for z = 1. We then get at leading order for s ≲ 1

$$\begin{equation}u\left(s\right)\underset{s\to {1}^{-}}{\sim }{y}_{0}\frac{{\left(1-s\right)}^{z}}{z}+\mathrm{O}\left({\left(1-s\right)}^{z+1}\right).\end{equation} \tag{ A.7 }$$

Therefore, inverting this expression and coming back to the initial variables we find outside the ground state

$$\begin{equation}\rho \left(U{ >}0\right)\underset{\beta {U}_{\mathrm{b}}\to 0}{\sim }{\rho }_{\mathrm{c}}\left[1-{\left(\frac{\beta Uz}{{y}_{0}}\right)}^{\frac{1}{z}}\right].\end{equation} \tag{ A.8 }$$

Using equation (A.8) and the definition of the potential $U\left(\boldsymbol{r}\right)=\frac{{U}_{\mathrm{b}}}{2}\left[1-{d}^{-1}{\sum }_{k}\mathrm{cos}\left(\pi {r}_{k}/{r}_{\mathrm{b}}\right)\right]$, the normalization in the condensation phase thus obeys

$$\begin{equation}N={N}_{\mathrm{c}}+{\rho }_{\mathrm{c}}{L}^{d}\left[1-{\left(\frac{\beta {U}_{\mathrm{b}}z}{2{y}_{0}}\right)}^{\frac{1}{z}}\mathcal{G}\left(z\right)\right],\end{equation} \tag{ A.9 }$$

where the function $\mathcal{G}\left(z\right)={\int }_{0}^{1}\mathrm{d}{x}_{1}\cdots {\int }_{0}^{1}\mathrm{d}{x}_{d}\;{\left[1-{d}^{-1}{\sum }_{k}\mathrm{cos}\left(\pi {x}_{k}\right)\right]}^{1/z}$ is nonzero and analytic but has no simple expression in general. Defining ${k}_{\text{B}}{T}_{\mathrm{c}}=z{U}_{\mathrm{b}}/\left(2{y}_{0}\right){\mathcal{G}}^{z}\left(z\right){\left(1-\overline{\rho }/{\rho }_{\mathrm{c}}\right)}^{-z}$, equation (A.9) is finally recast as equation (14).

In order to investigate the stability of the steady state solution of equation (1), in this section we calculate $\mathcal{F}\left[\rho +\delta \rho \right]-\mu N\left[\rho +\delta \rho \right]$, where $\mathcal{F}$ is defined in equation (27). Expanding up to second order in δρ, we find

$$\begin{align}\hfill \mathcal{F}\left[\rho +\delta \rho \right]-\mu N\left[\rho +\delta \rho \right]=& \mathcal{F}\left[\rho \right]-\mu N\left[\rho \right]+\int \;{\mathrm{d}}^{d}\boldsymbol{r}\;\delta \rho \left(\boldsymbol{r},t\right)\left\{U-\mu +{k}_{\text{B}}{T}_{\text{eff}}\right.\hfill \\ \hfill & \left.{\times}\;\left[\mathrm{ln}\left(\frac{\rho }{{\rho }_{\mathrm{c}}}\right){\Theta}\left({\rho }_{\mathrm{c}}-\rho \right)-\rho \left(\mathrm{ln}\left(\frac{\rho }{{\rho }_{\mathrm{c}}}\right)+\frac{{\rho }_{\mathrm{c}}}{\rho }-1\right)\delta \left(\rho -{\rho }_{\mathrm{c}}\right)\right]\right\}\hfill \\ \hfill & +\;\frac{1}{2}\int \;{\mathrm{d}}^{d}\boldsymbol{r}\int \;{\mathrm{d}}^{d}{\boldsymbol{r}}^{\prime }\;\delta \rho \left(\boldsymbol{r},t\right)\delta \rho \left({\boldsymbol{r}}^{\prime },t\right)\;{\delta }^{d}\left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right){k}_{\text{B}}{T}_{\text{eff}}\hfill \\ \hfill & {\times}\;\left\{\frac{1}{\rho }\;{\Theta}\left({\rho }_{\mathrm{c}}-\rho \right)-2\;\mathrm{ln}\left(\frac{\rho }{{\rho }_{\mathrm{c}}}\right)\delta \left(\rho -{\rho }_{\mathrm{c}}\right)+{\delta }^{\prime }\left({\rho }_{\mathrm{c}}-\rho \right)\rho \left(\mathrm{ln}\left(\frac{\rho }{{\rho }_{\mathrm{c}}}\right)+\frac{{\rho }_{\mathrm{c}}}{\rho }-1\right)\right\},\hfill \end{align} \tag{ B.1 }$$

where we have used ${\Theta}\left({\rho }_{\mathrm{c}}-\rho -\delta \rho \right)\approx {\Theta}\left({\rho }_{\mathrm{c}}-\rho \right)-\delta \rho \;\delta \left({\rho }_{\mathrm{c}}-\rho \right)+\frac{1}{2}\delta {\rho }^{2}\;{\delta }^{\prime }\left({\rho }_{\mathrm{c}}-\rho \right)$, with δ' formally being the distributional derivative of δ (and removed space and time dependencies of ρ and U for clarity). In order to investigate the fate of the δ terms of equation (B.1), we expand ρ = ρ_c(1 ± ɛ) which leads to

$$\begin{equation}\mathrm{ln}\left(\frac{\rho }{{\rho }_{\mathrm{c}}}\right)\underset{\varepsilon \to 0}{\sim }{\pm}\varepsilon ,\quad \mathrm{ln}\left(\frac{\rho }{{\rho }_{\mathrm{c}}}\right)+\frac{{\rho }_{\mathrm{c}}}{\rho }-1\underset{\varepsilon \to 0}{\sim }\frac{{\varepsilon }^{2}}{2}.\end{equation} \tag{ B.2 }$$

Hence, from the relations xδ(x) = x²δ'(x) = 0, the last terms of the second, third and fourth lines of equation (B.1) vanish. From equation (B.1), we can then read the following results in the steady state

$$\begin{equation}\frac{\delta \left(\mathcal{F}-\mu N\right)}{\delta \rho \left(\boldsymbol{r},t\right)}=U\left(\boldsymbol{r}\right)-\mu +{k}_{\text{B}}{T}_{\text{eff}}\;\mathrm{ln}\left(\frac{\rho \left(U\left(\boldsymbol{r}\right)\right)}{{\rho }_{\mathrm{c}}}\right){\Theta}\left({\rho }_{\mathrm{c}}-\rho \left(U\left(\boldsymbol{r}\right)\right)\right),\end{equation} \tag{ B.3 }$$

$$\begin{equation}\frac{{\delta }^{2}\left(\mathcal{F}-\mu N\right)}{\delta \rho \left(\boldsymbol{r},t\right)\delta \rho \left({\boldsymbol{r}}^{\prime },t\right)}=\frac{{k}_{\text{B}}{T}_{\text{eff}}}{\rho \left(U\left(\boldsymbol{r}\right)\right)}\;{\Theta}\left({\rho }_{\mathrm{c}}-\rho \left(U\left(\boldsymbol{r}\right)\right)\right)\;{\delta }^{d}\left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right).\end{equation} \tag{ B.4 }$$