Clustering and phase transitions in a 2D superfluid with immiscible active impurities

We generalize the projected GP model to include the dynamics of active impurities. The model is then described by the Hamiltonian

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\grad}{{\boldsymbol{\nabla}}} \newcommand{\Mp}{M_{\rm I}} \newcommand{\Vp}{V_{\rm I}} \newcommand{\Np}{{N_{\rm I}}} \newcommand{\p}{{\bf p}} \newcommand{\x}{{\bf x}} \displaystyle \nonumber H&=\int\left( \frac{\hbar^2}{2m} |\grad \psi |^2 +\frac{g}{2}|\mathcal{P}_{\rm G}[|\psi|^2]|^2 +\sum_{i=1}^{\Np} \Vp(| \x -{\bf q}_i |)\mathcal{P}_{\rm G}[|\psi|^2]\right) \mathrm{d} \x \nonumber \\ & \quad +\sum_{i=1}^{\Np}\frac{{\bf p}_i^2}{2 \Mp}+ \sum_{i<j}^{\Np}V_{\rm rep}(|{\bf q}_i-{\bf q}_j|), \label{Eq:HGP} \nonumber \end{align} \tag{ 1 }$$

where $\newcommand{\p}{{\bf p}} \psi$ is the collective wave-function of bosons having mass m, and $\newcommand{\p}{{\bf p}} g=4 \pi a_\mathrm{s} \hbar^2 /m$ being $a_{\rm s}$ the s-wave scattering length of bosons interaction. $\newcommand{\Np}{{N_{\rm I}}} \Np$ is the total number of impurities of mass $\newcommand{\Mp}{M_{\rm I}} \Mp$, that are described using their classical position and momentum ${\bf q}_i$ and ${\bf p}_i$, respectively. The strong repulsive potential $\newcommand{\Vp}{V_{\rm I}} \Vp$ determines the shape of the impurities by creating a large depletion in the fluid density. $V_{\rm rep}$ is a repulsive potential between impurities. The Galerkin projector $\mathcal{P}_{\rm G}$ truncates the system acting in Fourier space as $\newcommand{\p}{{\bf p}} \newcommand{\kmax}{k_{\rm max}} \mathcal{P}_{\rm G} \hat{\psi}_{\bf k}=\theta(\kmax-k)\hat{\psi}_{\bf k}$ with $\theta(\cdot)$ the Heaviside function, $\newcommand{\p}{{\bf p}} \hat{\psi}_{\bf k}$ the Fourier transform of $\newcommand{\x}{{\bf x}} \newcommand{\p}{{\bf p}} \psi(\x)$ and ${\bf k}$ the wave vector. The equations of motion are directly obtained by varying (1):

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dert}[1]{\frac{{\rm d}{#1}}{{\rm d}t}} \newcommand{\gd}{{\boldsymbol{\nabla}}^{2}} \newcommand{\dertt}[1]{\frac{\partial{#1}}{\partial t}} \newcommand{\alps}{{\frac{\hbar^2}{2m}}} \newcommand{\alp}{{\frac{\hbar}{2 m}}} \newcommand{\bet}{{g}} \newcommand{\Vp}{V_{\rm I}} \newcommand{\Np}{{N_{\rm I}}} \newcommand{\p}{{\bf p}} \newcommand{\x}{{\bf x}} \displaystyle {\rm i}\hbar\dertt{\psi}=\mathcal{P}_{\rm G} [- \alps \gd \psi + \bet\mathcal{P}_{\rm G} [|\psi|^2]\psi+ \sum_{i=1}^{\Np} \Vp(| \x -{\bf q}_i |)\psi]\label{Eq:GPE} \nonumber \end{align} \tag{ 2 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Mp}{M_{\rm I}} \newcommand{\Vp}{V_{\rm I}} \newcommand{\Np}{{N_{\rm I}}} \newcommand{\p}{{\bf p}} \newcommand{\x}{{\bf x}} \displaystyle \Mp\ddot{\bf q_i}= - \int \Vp(| \x -{\bf q}_i |) \mathcal{P}_{\rm G}[\nabla|\psi|^2]\, \mathrm{d} \x-\sum_{j\neq i}^{\Np}\nabla V_{\rm rep}(q_{ij}),\label{Eq:Particles} \nonumber \end{align} \tag{ 3 }$$

where $q_{ij}=|{\bf q}_i-{\bf q}_j|$ and we have replaced $\newcommand{\Mp}{M_{\rm I}} {\bf p}_i=\Mp\dot{\bf q}_i$. The previous set of equations exactly conserves the Hamiltonian, the number of bosons $\newcommand{\x}{{\bf x}} \newcommand{\p}{{\bf p}} N=\int |\psi|^2\mathrm{d} \x$ and momentum $\newcommand{\x}{{\bf x}} \newcommand{\p}{{\bf p}} \newcommand{\psib}{\overline{\psi}} {\bf P}=\int \frac{{\rm i}\hbar}{2}\left( \psi {\boldsymbol{\nabla}}\psib - \psib {\boldsymbol{\nabla}}\psi\right)\mathrm{d} \x+\sum_i {\bf p}_i$.

The impurities in the system feel an attractive force mediated by the superfluid density field [35, 37]. However, unlike the case of impurities described by classical fields [35], in the model (1) no repulsion mediated by the fluid exists. In order to mimic such a hard-sphere repulsion, we consider the Lennard–Jones-like potential $\newcommand{\rmin}{r_{\rm min}} \newcommand{\e}{{\rm e}} \newcommand{\rmi}{{\rm i}} V_{\rm rep}(r)=\epsilon \rmin^{12}/r^{12}$, as in [37]. We fix the energy $\newcommand{\e}{{\rm e}} \epsilon$ in order to set the minimum of the total interaction energy between impurities at zero temperature at a distance $r_\mathrm{min}$. Note that the specific shape of $V_\mathrm{rep}$ is not important, as long as it reproduces a hard-sphere repulsion. For the impurities potential we use a smoothed hat-potential $\newcommand{\Vp}{V_{\rm I}} \newcommand{\Rp}{a_{\rm I}} \Vp(r)=V_0(1-\tanh\left[(r^2 -\Rp^2)/4\Delta l^2\right])/2$, where $\newcommand{\Rp}{a_{\rm I}} \Rp$ sets the characteristic radius of the impurity and $\Delta l$ is a smoothing parameter. Finally, let us notice that in absence of impurities and at zero temperature, equation (2) can be linearized about a uniform density state $\newcommand{\p}{{\bf p}} |\psi|^2=\rho_\infty/m$, defining the phonon (sound) velocity $c=\sqrt{g\rho_\infty/m^2}$ with dispersive effects taking place for length scales smaller than the healing length $\newcommand{\x}{{\bf x}} \xi=\sqrt{\hbar^2/2g\rho_\infty }$.

We integrate the system (2) and (3) by using a pseudo-spectral code with $N_{\rm res}$ uniform grid points per direction of a squared domain of size $\newcommand{\p}{{\bf p}} L=2\pi$. We set $\newcommand{\kmax}{k_{\rm max}} \kmax=N_{\rm res}/3$ so that the truncated system exactly conserves all the invariants (provided that initially $\newcommand{\p}{{\bf p}} \mathcal{P}_{\rm G} \psi=\psi$ and $\newcommand{\Vp}{V_{\rm I}} \mathcal{P}_{\rm G} \Vp=\Vp$) [44], $c=\rho_\infty=1$, $V_0=10$ and $\newcommand{\e}{{\rm e}} \epsilon =0.006\,74$. As the healing length changes with temperature, we parametrize the solutions of (2) and (3) using its value taken at zero temperature. In thermal states, the only relevant dimensionless parameters are $\newcommand{\x}{{\bf x}} L/\xi$, $\newcommand{\x}{{\bf x}} \newcommand{\Rp}{a_{\rm I}} \Rp/\xi$, $\newcommand{\Np}{{N_{\rm I}}} \Np$, the relative mass $\newcommand{\p}{{\bf p}} \newcommand{\Mp}{M_{\rm I}} \newcommand{\Rp}{a_{\rm I}} \mathcal{M}=\Mp/\rho_\infty\pi\Rp^2$ and $\newcommand{\x}{{\bf x}} \newcommand{\kmax}{k_{\rm max}} \xi\kmax$. We set $\newcommand{\x}{{\bf x}} L/\xi=128$. The value of $\newcommand{\x}{{\bf x}} \newcommand{\kmax}{k_{\rm max}} \xi\kmax$ controls the strength of the nonlinear interactions and it is kept fixed to $\newcommand{\x}{{\bf x}} \newcommand{\p}{{\bf p}} \newcommand{\kmax}{k_{\rm max}} \xi\kmax=2\pi/3$. For this value, most of the excitations are phonons when the condensate fraction is large. For instance, it is compatible with the one used in [39], that is $\newcommand{\x}{{\bf x}} \newcommand{\kmax}{k_{\rm max}} \xi\kmax\sim2$; such value applies to a gas of ⁸⁷Rb atoms.

We start by presenting a long temporal evolution of a system having 40 impurities of mass $\mathcal{M}=0.1$, radius $\newcommand{\x}{{\bf x}} \newcommand{\Rp}{a_{\rm I}} \Rp=4\xi$, initially located at random positions (avoiding overlaps) and having zero velocity. The density field of the initial condition is displayed in figure 1 ($\newcommand{\x}{{\bf x}} t= 0 \xi/c$). Impurities correspond to dark holes. During the time evolution, the short-range interaction among impurities mediated by the fluid let them collapse into small clusters ($\newcommand{\x}{{\bf x}} t= 815 \xi/c$); waves with random phases are generated, populating small length scales and starting the thermalization process. This thermal noise induces the clusters to move in a stochastic way, and to grow further ($\newcommand{\x}{{\bf x}} t= 1691 \xi/c$). Eventually, the system reaches thermal equilibrium where only one big cluster is observed in a bath of thermalized waves, ($\newcommand{\x}{{\bf x}} t= 10\,965 \xi/c$). A movie of the the evolution is available in the supplementary data (stacks.iop.org/JPhysA/52/305501/mmedia).

Zoom In Zoom Out Reset image size

The evolution illustrated in figure 1 is an example of thermalization occurring in the micro-canonical ensemble, as the thermal state is achieved keeping all the invariants conserved. Such dynamical process is numerically costly and does not directly provide access to the conjugate thermodynamical variables: temperature and chemical potential (here we only consider zero momentum states). To overcome these issues, in [44] a stochastic relaxation was introduced in order to efficiently generate thermal states in the grand-canonical ensemble. We make use of this approach adapting it to the Hamiltonian (1). The stochastic dynamics is ruled by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bet}{{g}} \newcommand{\p}{{\bf p}} \newcommand{\x}{{\bf x}} \displaystyle \hbar\frac{\partial\hat{\psi}_\mathbf{k}}{\partial t}=-\frac{\partial F}{\partial \hat{\psi}_\mathbf{k}^*} + \sqrt{\frac{2\hbar}{\beta}}\hat{\xi}_\mathbf{k}\label{eq:stocpsi} \nonumber \end{align} \tag{ 4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bet}{{g}} \newcommand{\Mp}{M_{\rm I}} \newcommand{\q}{{\bf q}} \newcommand{\p}{{\bf p}} \newcommand{\x}{{\bf x}} \displaystyle \qquad \qquad \qquad \qquad \qquad \frac{\partial {\bf q}_{i}}{\partial t}=-\frac{\partial F}{\partial {\bf q}_{i}} + \sqrt{\frac{2}{\beta}}\boldsymbol{\xi}_i^q,\quad \Mp\frac{\partial \dot{{\bf q}}_i}{\partial t}=-\frac{\partial F}{\partial\dot{{\bf q}}_i} + \sqrt{\frac{2\Mp}{\beta}}\boldsymbol{\xi}_i^{\dot{q}}\label{eq:stocq} \nonumber \end{align} \tag{ 5 }$$

where $F=H-\mu N$ is the free energy of the system, $\mu$ is the chemical potential controlling the number of bosons and $\newcommand{\bet}{{g}} \beta$ is the inverse temperature; $\newcommand{\x}{{\bf x}} \hat{\xi}_\mathbf{k},\boldsymbol{\xi}_i^q$ and $\newcommand{\x}{{\bf x}} \boldsymbol{\xi}_i^{\dot{q}}$ are independent Gaussian white noises of unit variance. It can be shown by using the Fokker–Planck equation associated to (4) and (5), that the stationary probability distribution is given by the Gibbs grand-canonical distribution $\newcommand{\p}{{\bf p}} \newcommand{\bet}{{g}} \mathbb{P}[\hat{\psi}_\mathbf{k},{\bf q}_i,{\dot {\bf{q}}}_i]\propto {\rm e}^{-\beta F}$. In the micro-canonical ensemble, $\newcommand{\p}{{\bf p}} \mathbb{P}[\hat{\psi}_\mathbf{k},{\bf q}_i,{\dot {\bf{q}}}_i]$ is also the stationary solution of the Liouville equation that describes the evolution of the phase-space distribution of the Hamiltonian system (1) [44]. It is evident from (5) that varying the impurity masses modifies only the variance of impurity momenta in the steady state. Namely, the configurations of impurities in the steady state and the statistical properties of the thermalized system are independent of impurity masses. We define the temperature as $\newcommand{\bet}{{g}} T=1/k_\mathcal{N}\beta$, with $k_\mathcal{N}=L^2/\mathcal{N}$ and $\newcommand{\p}{{\bf p}} \newcommand{\kmax}{k_{\rm max}} \mathcal{N}=\pi \kmax^2$ the total number of Fourier modes. With this definition T is an energy per unit of surface such that at low temperatures $F\approx T L^2$, because of equipartition. With these choices, intensive quantities remain constant when increasing the system size. In addition, we fix the total density mass $\bar{\rho}=mN/L^2=1$ by dynamically adjusting the chemical potential [44]. We use (4) and (5) to study the effect of impurities on the superfluid condensed fraction n₀ defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\p}{{\bf p}} \displaystyle n_0=\frac{\langle\left| \int_{\mathcal{V}'} \psi({\bf x})\,\mathrm{d}{\bf x}\right|^2 \rangle_T}{\langle\int_{\mathcal{V}'}| \psi({\bf x})|^2\,\mathrm{d}{\bf x}\rangle_{T}} \label{Eq:CondFrac}, \nonumber \end{align} \tag{ 6 }$$

where $\mathcal{V}'$ is the domain excluding the region occupied by the impurities and $\langle \cdot\rangle$ stands for average over realizations at temperature T⁴.

We first perform a temperature scan without impurities. The condensed fraction is shown in figure 2(a) (solid blue line). The transition temperature $T_\lambda$ is the lowest temperature where the condensed fraction can be considered negligible. We estimate it in a consistent way adopting the following numerical protocol: we take the points around which $n_0^2(T)$ is close to zero and we perform a linear interpolation of it. $T_\lambda$ is then determined by finding the point where the linear fit vanishes. From now on, we indicate with $T_\lambda^0$ the transition temperature in the system without impurities. Then, we perform temperature scans varying the repulsive potential parameter $\newcommand{\rmin}{r_{\rm min}} \newcommand{\rmi}{{\rm i}} \rmin$ with a fixed number of impurities $\newcommand{\Np}{{N_{\rm I}}} \Np=31$ having radius $\newcommand{\x}{{\bf x}} \newcommand{\Rp}{a_{\rm I}} \Rp=4\xi$. The results are also shown in figure 2. In figures 2(b) and (c) snapshots of in the steady date are displayed respectively in the case of high and low repulsion among impurities. For both cases we report three different temperatures. Depending on the strength of the repulsion potential two different behaviors of n₀ can be observed, as it is clear in figure 2(a). When the repulsion among impurities is strong enough (blue markers, $\newcommand{\x}{{\bf x}} \newcommand{\rmin}{r_{\rm min}} \newcommand{\rmi}{{\rm i}} \rmin\geqslant11\xi$), clusters are broken already at temperatures lower than $T_\lambda^0$ (see figure 2(b)) and impurities have no appreciable effect on n₀. On the other hand, for $\newcommand{\x}{{\bf x}} \newcommand{\rmin}{r_{\rm min}} \newcommand{\rmi}{{\rm i}} \rmin\leqslant10.5\xi$ (red markers) impurities remain clustered and lead to an increase of the n₀ at medium-high temperatures (see figure 2(c)).

Zoom In Zoom Out Reset image size

It has been shown that impurities experience a short-range attractive force, mediated by the superfluid density [35, 37]. This interaction is characterized by a potential energy, denoted here $\newcommand{\Epp}{E_{\rm I \leftrightarrow I}} \Epp$. We compute this energy numerically by measuring the full GP free energy of the ground state with two impurities placed at a distance $\Delta q$, without the contribution of the repulsion. The constant contribution to $\newcommand{\Epp}{E_{\rm I \leftrightarrow I}} \Epp$ is eventually set to zero. Impurities are then repelled because of $V_{\rm rep}$, generating cluster structures as the one observed in figure 1(d). However, if thermal fluctuations are large enough, the bound among impurities can be broken. In figure 3(a) we compute the interaction energy between two impurities at zero temperature $\newcommand{\Epp}{E_{\rm I \leftrightarrow I}} \Epp$, as a function of their distance $\Delta q$ (dotted black line). As a reference, the figure also displays the repulsive potential $V_{\rm rep}$ with $\newcommand{\Rp}{a_{\rm I}} \newcommand{\rmin}{r_{\rm min}} \newcommand{\rmi}{{\rm i}} \rmin=2\Rp$ (dotted–dashed golden line). The sum of both potentials is displayed in the same figure for different values of $\newcommand{\rmin}{r_{\rm min}} \newcommand{\rmi}{{\rm i}} \rmin$: for sufficiently small values of $\newcommand{\rmin}{r_{\rm min}} \newcommand{\rmi}{{\rm i}} \rmin$, a potential well $\Delta U$ centered at $\newcommand{\rmin}{r_{\rm min}} \newcommand{\rmi}{{\rm i}} \rmin$ appears. We thus expect a pair of impurities to split in a finite time at the temperature $T_{\rm cl} \sim \Delta U / k_\mathcal{N}$, as in a standard escape problem from a potential well [45]. In order to quantify this clustering transition, we measure the average square distance between impurities and their center of mass

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Np}{{N_{\rm I}}} \newcommand{\q}{{\bf q}} \displaystyle \delta^2(T)=\frac{{\rm d}^2(T)-{\rm d}^2(0)}{{\rm d}^2(0)},\quad {\rm with }\quad {\rm d}^2(T)=\sum_{j=1}^\Np\langle |{\bf q}_j-{\bf q}_\mathrm{cm}|^2\rangle \quad {\rm and }\quad \mathbf{q}_\mathrm{cm}=\sum_j^{N_I}\frac{\mathbf{q}_j}{N_\mathrm{I}}. \label{Eq:delta2} \nonumber \end{align} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

Figure 3(b) displays $\delta^2(T)$ as a function of $T/T_{\rm cl}$ for different values of $\newcommand{\rmin}{r_{\rm min}} \newcommand{\rmi}{{\rm i}} \rmin$. A transition around $T_{\rm cl}$ is indeed observed, where discrepancies are likely due to oversimplifications made in the estimation of $T_\mathrm{cl}$, namely by neglecting the many-body impurity interactions and by using the interaction potential obtained at T  =  0. For weak repulsion, even if the condensate vanishes, impurities still feel the density-mediated attraction.

Studying different values of $r_\mathrm{min}$ allowed us to show that an increasing of n₀ occurs only when the fluid depletion, due to the presence of impurities, is confined to a large connected region at all temperatures. Therefore such effect can not be simply explained by the local increase of density in regions not occupied by impurities. In the following we consider hard-sphere impurities by fixing $\newcommand{\Rp}{a_{\rm I}} r_\mathrm{min}=2\Rp$. In order to quantitatively characterize the change in n₀, we study how the condensation transition changes when varying the filling fraction

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\p}{{\bf p}} \displaystyle \Phi=1 - \frac{|\langle \psi \rangle_{T=0}|^2} {\langle|\psi|^2\rangle_{T=0}}, \label{Eq:fillingfraction} \nonumber \end{align} \tag{ 8 }$$

which corresponds to the fraction of the total volume occupied by the impurities. In figure 4(a) the condensed fraction is shown for different values of $\Phi$, obtained by varying the number of impurities. It is evident that the larger is the number of clustered impurities, the higher results the condensation transition temperature. We explicit the dependence of the transition temperature on the filling fraction as $T_\lambda(\Phi)$.

Zoom In Zoom Out Reset image size

The condensation temperature $T_\lambda(\Phi)$ is measured for different values of $\Phi$ following the same procedure explained in the previous section. In figure 5(a) the relative increase $\Delta T_\lambda=(T_\lambda(\Phi)-T_\lambda^0)/T_\lambda^0$ is displayed. Remarkably, $\Delta T_\lambda$ scales linearly with $\Phi$ growing up to $20\%$. We have checked by varying the number of impurities, their size and the parameter $r_\mathrm{min}$ for values lower than $2a_\mathrm{I}$, that n₀ only depends on $\Phi$ and T (data not shown). Despite the change on $T_\lambda$, the condensed fraction curves collapse as expected to a single one, once plotted versus $T/T_\lambda(\Phi)$ (see figure 4(b)).

Zoom In Zoom Out Reset image size

Finally, we briefly address the role of impurities in the the BKT transition. A detailed study will be left for a further work. This phase transition manifests through a change in the behavior of the correlation function $\newcommand{\p}{{\bf p}} g_1(r)=\langle \psi(0)\psi^*({\bf r})\rangle$ at the BKT transition temperature $T_{\rm BKT}$. At $T<T_{\rm BKT}$, it presents a power-law decay $\newcommand{\alp}{{\frac{\hbar}{2 m}}} g_1(r)\sim r^{-\alpha}$, where $\newcommand{\alp}{{\frac{\hbar}{2 m}}} \alpha$ depends linearly on the temperature; at high temperatures, it exhibits the standard exponential decay of disordered systems. In figure 5(b) we show g₁(r) at different temperatures (lower, close and higher than $T_\mathrm{BKT}$) where these two behaviors are clearly distinguishable. The BKT transition temperature can be thus determined by finding where g₁(r) abruptly changes its behavior [46]. With no impurities and the parameters used in this article, the BKT transition takes place at $T_{\rm BKT}=0.83T_\lambda^0$. Note that because of the Mermin–Wagner–Hohenberg theorem [7, 8], the condensation critical temperature vanishes as $1/\log{L}$ in the thermodynamic limit, so in principle for a very large system we could have $T_\lambda^0<T_{\rm BKT}$. We do not address such limit in this article. The presence of impurities in the system modifies the decay of g₁(r) by shifting $T_{\rm BKT}$ to higher temperatures. Figure 5(a) also displays the relative increase $\Delta T_{\rm BKT}=(T_{\rm BKT}(\Phi)-T_{\rm BKT}^0)/T_{\rm BKT}^0$ of the BKT transition temperature for different filling fractions $\Phi$. Although $T_\lambda$ and $T_{\rm BKT}$ both grow up to $20\%$ when the $\Phi$ is increased, their ratio remains almost constant. Let us remark that there are no important effects on $T_\mathrm{BKT}$ if impurities are not clustered.

The increase of the transition temperatures $T_\lambda$ and $T_{\rm BKT}$ can be explained by a simple phenomenological argument. Large objects in the system, such as the clustered impurities, modify the fluid wave-function boundary conditions. In particular, they impose effective Dirichlet boundary condition leading to symmetries in Fourier space, decreasing the number of active modes. At a given temperature, with less active modes, the energy is smaller and thus a higher temperature is necessary to induce a transition. We stress that this result is general and does not depend on the choice of the repulsion potential $V_\mathrm{rep}$, as long as it is sufficiently short-range to allow the formation of a large-size cluster at all $T<T_\lambda$. Our results could apply as well to two-component BECs, with the components having different condensation temperatures and strong repulsion between them. Finally, the same effect on the condensation curve will occur in three dimensions, as it comes from a geometrical effect.