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

Umberto Giuriato, Giorgio Krstulovic, Davide Proment

► To cite this version:

Umberto Giuriato, Giorgio Krstulovic, Davide Proment. Clustering and phase transitions in a 2D superfluid with immiscible active impurities. Journal of Physics A: Mathematical and Theoretical, 2019, 52 (30), pp.305501. 10.1088/1751-8121/ab2607. hal-02413891

HAL Id: hal-02413891

https://hal.science/hal-02413891v1

Submitted on 28 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**PAPER • OPEN ACCESS**

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

To cite this article: Umberto Giuriato et al 2019 *J. Phys. A: Math. Theor.* **52** 305501

View the [article online](#) for updates and enhancements.

You may also like

- Single hole doped strongly correlated ladder with a static impurity
  S Gayen

- Electronic and Shallow Impurity States in Semiconductor Heterostructures Under an Applied Electric Field
  Hai-Yang Zhou, Shi-Wei Gu and Yao-Ming Shi

- (Invited) Large Mobility Modulation Due to Discrete Impurities in Nanowires
  Nobuyuki Sano
# Clustering and phase transitions in a 2D superfluid with immiscible active impurities

Umberto Giuriato<fnref n="1" />,<fnref n="3" /> ID, Giorgio Krstulovic<fnref n="1" /> ID and Davide Proment<fnref n="2" /> ID

E-mail: umberto.giuriato@oca.eu

Received 18 January 2019, revised 8 May 2019

Accepted for publication 31 May 2019

Published 1 July 2019

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

Phase transitions of a finite-size two-dimensional superfluid of bosons in presence of active impurities are studied by using the projected Gross-Pitaevskii model. Impurities are described with classical degrees of freedom. A spontaneous clustering of impurities during the thermalization is observed. Depending on the interaction among impurities, such clusters can break due to thermal fluctuations at temperatures where the condensed fraction is still significant. The emergence of clusters is found to increase the condensation transition temperature. The condensation and the Berezinskii-Kosterlitz-Thouless transition temperatures, determined numerically, are found to strongly depend on the volume occupied by the impurities: a relative increase up to a 20% of their respective values is observed, whereas their ratio remains approximately constant.

Keywords: superfluids, Bose-Einstein condensates, clustering, phase transitions, impurities

S Supplementary material for this article is available online

(Some figures may appear in colour only in the online journal)

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Nice, France</footnote>

<footnote n="2">² School of Mathematics, University of East Anglia, Norwich Research Park, Norwich, NR4 7TJ, United Kingdom</footnote>

<footnote n="3">³ Author to whom any correspondence should be addressed.</footnote>
# 1. Introduction

When a fluid composed of bosons is cooled down or the number of particles is increased, the system experiences a phase transition giving rise to a macroscopic state known as Bose-Einstein condensate (BEC) [1]. Since its first experimental observation by Anderson et al [2], BECs have been realized in systems of very different nature such as cold atomic gases [1, 3], solid-state quasiparticles [4, 5] and light in optical micro-cavities [6]. Whereas in three spatial dimensions a condensate is stable with respect to thermal fluctuations, in two dimensions such fluctuations can destroy the long-range order of the system. This is a general result in statistical field theory known as the Mermin-Wagner-Hohenberg theorem [7, 8]: it states that a continuous symmetry cannot be spontaneously broken in dimensions lower than three, otherwise large-scale Goldstone modes would have an infinite infrared contribution to the two-point correlator. This theorem assumes the thermodynamic limit, that is the system size being infinite. However, for a finite system, condensation can be recovered, having a transition temperature $T_λ$ that vanishes as the inverse of the logarithm of the system size.

Although condensation is formally forbidden in an infinite two-dimensional system, a peculiar phase transition of a different nature has attracted the attention of physicists and mathematicians since its independent discovery in the early 70s by Berezinskii, Kosterlitz and Thouless (BKT) [9–11]. The BKT transition is an infinite-order topological phase transition and manifests itself in systems that belong to the same universality class. By approaching the BKT transition temperature, $T_{BKT}$, from below, the system switches from a gas of bounded vortex-antivortex pairs to a gas of free vortices, moving from a quasi-ordered phase to a disordered one. The BKT transition has been observed in BECs made of dilute gases [12–16], exciton-polaritons [17], liquid helium films [18] and studied theoretically and numerically [19–22]; for a review on the topic, see for instance [3].

The purpose of this article is to study how the statistical mechanical properties of a two dimensional superfluid of bosons are affected by the presence of impurities. Particles and impurities have been used in superfluids since the early experiments in $^4$He [23] mainly with detection purposes: electrons, ions and neutral impurities such as hydrogen particles and exciters have been exploited to visualize quantized vortices, to study their dynamics and the statistics of superfluid (quantum) turbulence [24–26]. More recently, the investigation of the interaction between one or more impurities and superfluids has been the main topic of experiments in cold atoms [27, 28], superfluids of light [29, 30] and polaritons in semiconductor microcavities [31]. On the theoretical side, the dynamics of impurities in a BEC has been also addressed [32], as well as the properties of $^3$He and H impurities on thin $^4$He films [33, 34]. In addition, Ricaard and Roberts studied how a collection of impurities affects the ground-state of a BEC by using a mean field model [35]. In this last work, four phases were identified, depending on the value of the interaction couplings. In particular, if the scattering lengths between impurity fields are positive, impurities behave as localized objects, they separate from the condensate and present a hard-sphere repulsion between each other.

We investigate here impurity clustering and phase transitions occurring in a minimal model that mimics such situation: the Gross-Pitaevskii (GP) equation coupled with active immiscible impurities having classical degrees of freedom. Such model was introduced in [36] and recently used in two-dimensional numerical simulations to study impurity-impurity and impurity-vortex interactions [37, 38]. Finite-temperature BECs can be studied by using the projected GP equation, that is obtained by introducing a cut-off $k_{max}$ in Fourier space: this regularizes the classical mean-field ultra-violet divergence. The projected GP model is an effective model to study the condensation transition in two and three dimensions [21, 39–41] and superfluid vortex dynamics at finite temperature [42–44].
## 2. Theoretical model and numerical results

### 2.1. Model for impurities in a superfluid

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

$$ H = \int \left( \frac{\hbar^2}{2m} |\nabla \psi|^2 + \frac{g}{2} |\mathcal{P}_G[|\psi|^2]|^2 + \sum_{i=1}^{N_I} V_I(|\mathbf{x} - \mathbf{q}_i|) \mathcal{P}_G[|\psi|^2] \right) \mathrm{d}\mathbf{x} \\ + \sum_{i=1}^{N_I} \frac{\mathbf{p}_i^2}{2M_I} + \sum_{i<j}^{N_I} V_{\mathrm{rep}}(|\mathbf{q}_i - \mathbf{q}_j|), \qquad (1) $$

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

$$ i\hbar \frac{\partial \psi}{\partial t} = \mathcal{P}_G \left[ -\frac{\hbar^2}{2m} \nabla^2 \psi + g \mathcal{P}_G [ |\psi|^2 ] \psi + \sum_{i=1}^{N_I} V_I ( |\mathbf{x} - \mathbf{q}_i | ) \psi \right] \qquad (2) $$

$$ M_I \ddot{\mathbf{q}}_i = - \int V_I(|\mathbf{x} - \mathbf{q}_i|) \mathcal{P}_G [\nabla |\psi|^2] \mathrm{d}\mathbf{x} - \sum_{j \neq i}^{N_I} \nabla V_{\mathrm{rep}}(q_{ij}), \qquad (3) $$

where $q_{ij} = |\mathbf{q}_i - \mathbf{q}_j|$ and we have replaced $\mathbf{p}_i = M_I \dot{\mathbf{q}}_i$. The previous set of equations exactly conserves the Hamiltonian, the number of bosons $N = \int |\psi|^2 \mathrm{d}\mathbf{x}$ and momentum $\mathbf{P} = \int \frac{i\hbar}{2} (\psi \nabla \bar{\psi} - \bar{\psi} \nabla \psi) \mathrm{d}\mathbf{x} + \sum_i \mathbf{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 $V_{\mathrm{rep}}(r) = \epsilon r_{\mathrm{min}}^{12}/r^{12}$, as in [37]. We fix the energy $\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 $V_I(r) = V_0(1 - \tanh[(r^2 - a_I^2)/4\Delta l^2])/2$, where $a_I$ 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 $|\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 $\xi = \sqrt{\hbar^2/2g\rho_{\infty}}$.

We integrate the system (2) and (3) by using a pseudo-spectral code with $N_{\mathrm{res}}$ uniform grid points per direction of a squared domain of size $L = 2\pi$. We set $k_{\mathrm{max}} = N_{\mathrm{res}}/3$ so that the truncated system exactly conserves all the invariants (provided that initially $\mathcal{P}_G \psi = \psi$ and $\mathcal{P}_G V_I = V_I$) [44], $c = \rho_{\infty} = 1$, $V_0 = 10$ and $\epsilon = 0.00674$. 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 $L/\xi$, $a_I/\xi$, $N_I$, the
<figure><img src="image_4.png" /><figcaption>Figure 1. Snapshots of the fluid density during the GP temporal evolution of a state with 40 impurities (dark holes).</figcaption></figure>

relative mass $M = M_1/\rho_\infty\pi a_1^2$ and $\xi k_{\max}$. We set $L/\xi = 128$. The value of $\xi k_{\max}$ controls the strength of the nonlinear interactions and it is kept fixed to $\xi k_{\max} = 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 $\xi k_{\max} \sim 2$; such value applies to a gas of $^{87}$Rb atoms.

We start by presenting a long temporal evolution of a system having 40 impurities of mass $M = 0.1$, radius $a_1 = 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 ($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 ($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 ($t = 1691\xi/c$). Eventually, the system reaches thermal equilibrium where only one big cluster is observed in a bath of thermalized waves, ($t = 10965\xi/c$). A movie of the evolution is available in the supplementary data (stacks.iop.org/JPhysA/52/305501/mmedia).

## 2.2. Grand-canonical thermal states

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

$$ \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}} \quad (4) $$

$$ \frac{\partial \mathbf{q}_i}{\partial t} = -\frac{\partial F}{\partial \mathbf{q}_i} + \sqrt{\frac{2}{\beta} \xi_i^q}, \quad M_1 \frac{\partial \dot{\mathbf{q}}_i}{\partial t} = -\frac{\partial F}{\partial \dot{\mathbf{q}}_i} + \sqrt{\frac{2M_1}{\beta} \xi_i^q} \quad (5) $$

where $F = H - \mu N$ is the free energy of the system, $\mu$ is the chemical potential controlling the number of bosons and $\beta$ is the inverse temperature; $\hat{\xi}_i^q$ and $\xi_i^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 $P[\hat{\psi}_{\mathbf{k}}, \mathbf{q}_i, \dot{\mathbf{q}}_i] \propto e^{-\beta F}$. In the micro-canonical ensemble, $P[\hat{\psi}_{\mathbf{k}}, \mathbf{q}_i, \dot{\mathbf{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 $T = 1/k_N\beta$, with $k_N = L^2/N$ and $N = \pi k_{max}^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 TL^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_0$ defined as

$$n_0 = \frac{\langle |\int_{V'} \psi(\mathbf{x}) \, d\mathbf{x}|^2 \rangle_T}{\langle \int_{V'} |\psi(\mathbf{x})|^2 d\mathbf{x} \rangle_T}, \quad (6)$$

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

## 2.3. Clustering of impurities

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 $r_{min}$ with a fixed number of impurities $N_I = 31$ having radius $a_I = 4\xi$. The results are also shown in figure 2. In figures 2(b) and (c) snapshots of in the steady state 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_0$ can be observed, as it is clear in figure 2(a). When the repulsion among impurities is strong enough (blue markers, $r_{min} \ge 11\xi$), clusters are broken already at temperatures lower than $T_\lambda^0$ (see figure 2(b)) and impurities have no appreciable effect on $n_0$. On the other hand, for $r_{min} \le 10.5\xi$ (red markers) impurities remain clustered and lead to an increase of the $n_0$ at medium-high temperatures (see figure 2(c)).

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 $E_{I \leftrightarrow I}$. 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 $E_{I \leftrightarrow I}$ is eventually set to zero. Impurities are then repelled because of $V_{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 $E_{I \leftrightarrow I}$, as a function of their distance $\Delta q$ (dotted black line). As a reference, the figure also displays the repulsive potential $V_{rep}$ with $r_{min} = 2a_I$ (dotted-dashed golden line). The sum of both potentials is displayed in the same figure for different values of $r_{min}$: for sufficiently small

<footnote n="4">⁴ This definition gives the same result as $n_0 = \left( \frac{|\langle\psi\rangle_T|^2}{N} \right) \left( \frac{|\langle\psi\rangle_{T=0}|^2}{N} \right)^{-1} = \frac{|\langle\psi\rangle_T|^2}{|\langle\psi\rangle_{T=0}|^2}$. For numerical convenience we use the latter in our computations.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 2. (a) Condensed fraction as a function of the temperature for different values of $r_{min}$. Temperatures are expressed in units of the condensation temperature with no impurities $T_{\lambda}^0$. (b) Snapshots of the density field in the steady state at different temperatures in the case of high repulsion ($r_{min} = 11\xi$) among impurities. (c) The same of (b) but for low repulsion ($r_{min} = 8\xi$). Scans are performed with $N_I = 31$ and $a_I = 4\xi$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3. (a) Impurity-impurity interaction $E_{I \leftrightarrow I}$ as a function of their distance $\Delta q/\xi$ (dotted black line) and repulsive potential (dotted-dashed golden line) for $r_{min} = 8\xi$. Different markers correspond to total energy $E_{I \leftrightarrow I} + V_{rep}$ for different values of $r_{min}$. (b) Relative impurity distance $\delta^2(T)$ as a function of the normalized temperature $T/T_{cl}$. Scans are performed with $N_I = 31$ and $a_I = 4\xi$. The ratios $T_{cl}/T_{\lambda}$ are 7.41, 2.51, 1.5, 1.02, 0.53, 0.002 for $r_{min}$ from 8ξ to 12ξ respectively. The markers refer to the same legend as in figure 2.</figcaption></figure>

values of $r_{min}$, a potential well $\Delta U$ centered at $r_{min}$ appears. We thus expect a pair of impurities to split in a finite time at the temperature $T_{cl} \sim \Delta U/k_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

$$ \delta^2(T) = \frac{d^2(T) - d^2(0)}{d^2(0)}, \quad \text{with} \quad d^2(T) = \sum_{j=1}^{N_I} \langle |\mathbf{q}_j - \mathbf{q}_{cm}|^2 \rangle \quad \text{and} \quad \mathbf{q}_{cm} = \sum_{j}^{N_I} \frac{\mathbf{q}_j}{N_I}. \qquad (7) $$

Figure 3(b) displays $\delta^2(T)$ as a function of $T/T_{cl}$ for different values of $r_{min}$. A transition around $T_{cl}$ is indeed observed, where discrepancies are likely due to oversimplifications made in the estimation of $T_{cl}$, namely by neglecting the many-body impurity interactions and by
<figure><img src="image_7.png" /><figcaption>Figure 4. (a) Condensed fraction as a function of temperature for different values of the filling fraction Φ. Temperatures are expressed in units of the condensation temperature with no impurities $T_{\lambda}^{0}$. (b) Condensed fraction as a function of the temperature normalized with $T_{\lambda}(\Phi)$. Scans with $r_{p} = 4\xi$ and $r_{\min} = 2r_{p} = 8\xi$.</figcaption></figure>

using the interaction potential obtained at $T = 0$. For weak repulsion, even if the condensate vanishes, impurities still feel the density-mediated attraction.

## 2.4. Condensation and BKT transition temperatures in presence of impurity clusters

Studying different values of $r_{\min}$ allowed us to show that an increasing of $n_{0}$ 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 $r_{\min} = 2a_{l}$. In order to quantitatively characterize the change in $n_{0}$, we study how the condensation transition changes when varying the filling fraction

$$ \Phi = 1 - \frac{|\langle\psi\rangle_{T=0}|^2}{\langle|\psi|^2\rangle_{T=0}}, \qquad (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 Φ, 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)$.

The condensation temperature $T_{\lambda}(\Phi)$ is measured for different values of Φ 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 Φ growing up to 20%. We have checked by varying the number of impurities, their size and the parameter $r_{\min}$ for values lower than $2a_{l}$, that $n_{0}$ only depends on Φ 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)).

Finally, we briefly address the role of impurities in 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 $g_{1}(r) = \langle\psi(0)\psi^{*}(\mathbf{r})\rangle$ at the BKT transition temperature $T_{\text{BKT}}$. At
<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) Relative increments of $T_λ$ and $T_{BKT}$ as a function of the filling fraction $\Phi$. The relative increment of the ratio $\frac{T_{BKT}}{T_λ}$ is also shown. (b) Spatial correlation function $g_1(r)$ for three different temperatures (lower, close and higher than $T_{BKT}$) with and without impurities.</figcaption></figure>

$T < T_{BKT}$, it presents a power-law decay $g_1(r) \sim r^{-\alpha}$, where $\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_1(r)$ at different temperatures (lower, close and higher than $T_{BKT}$) where these two behaviors are clearly distinguishable. The BKT transition temperature can be thus determined by finding where $g_1(r)$ abruptly changes its behavior [46]. With no impurities and the parameters used in this article, the BKT transition takes place at $T_{BKT} = 0.83T_λ$. 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_λ^0 < T_{BKT}$. We do not address such limit in this article. The presence of impurities in the system modifies the decay of $g_1(r)$ by shifting $T_{BKT}$ to higher temperatures. Figure 5(a) also displays the relative increase $\Delta T_{BKT} = (T_{BKT}(\Phi) - T_{BKT}^0)/T_{BKT}^0$ of the BKT transition temperature for different filling fractions $\Phi$. Although $T_λ$ and $T_{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_{BKT}$ if impurities are not clustered.

The increase of the transition temperatures $T_λ$ and $T_{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_{rep}$, as long as it is sufficiently short-range to allow the formation of a large-size cluster at all $T < T_λ$. 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.
## 3. Discussion

In this article we studied thermal states of two-dimensional superfluids with active impurities. We demonstrated how the phase transitions are affected by the emergence of impurity clusters, opening up the possibility to raise the transition temperatures in experiments by doping superfluids with specific types of impurities. Such result rises new questions that would be interesting to address in detail. In particular, it is remarkable that the presence of impurities does not disorder the system inducing a loss of coherence. Is there a maximum value of the critical temperature that can be reached using impurities, or it will continue to increase until the impurities occupy the full domain? Could the modification of the condensation curve be rephrased as a competition between the full perimeter and the full area of the impurities? Moreover, this system presents a rich behaviour that, up our knowledge, has not yet been addressed in details. For instance, during the thermalisation dynamics, impurities cluster similarly to a diffusion-limited aggregation process [47]. Also, a complete study of the BKT transition, considering the opposite limit $T_λ < T_{BKT}$, needs to be investigated and might devise new interesting physics.

## Acknowledgments

The authors were supported by the cost-share Royal Society International Exchanges Scheme (ref. IE150527) in conjunction with CNRS. GK and DP acknowledge the Federation Doeblin for supporting DP during his sojourn in Nice. DP was supported by the UK Engineering and Physical Sciences Research Council (EPSRC) research grant EP/P023770/1. GK and UG were also supported by the ANR JCJC GIANTE ANR-18-CE30-0020-01. Computations were carried out on the Mésocentre SIGAMM hosted at the Observatoire de la Côte d'Azur. The authors acknowledge the anonymous referees for their interesting comments.

## ORCID iDs

Umberto Giuriato https://orcid.org/0000-0002-0211-2458

Giorgio Krstulovic https://orcid.org/0000-0002-9934-6292

Davide Proment https://orcid.org/0000-0002-9472-0097

## References

[1] Griffin A, Snoke D W and Stringari S 1996 *Bose-Einstein Condensation* (Cambridge: Cambridge University Press)

[2] Anderson M H, Ensher J R, Matthews M R, Wieman C E and Cornell E A 1995 Observation of Bose-Einstein condensation in a dilute atomic vapor *Science* **269** 198-201

[3] Hadzibabic Z and Dalibard J 2011 Two-dimensional Bose fluids: an atomic physics perspective *Riv. Nuovo Cimento* **34** 389

[4] Kasprzak J et al 2006 Bose-Einstein condensation of exciton polaritons *Nature* **443** 409

[5] Demokritov S O, Demidov V E, Dzyapko O, Melkov G A, Serga A A, Hillebrands B and Slavin A N 2006 *Bose-Einstein condensation of quasi-equilibrium magnons at room temperature under pumping* *Nature* **443** 430

[6] Klaers J, Schmitt J, Vewinger F and Weitz M 2010 Bose-Einstein condensation of photons in an optical microcavity *Nature* **468** 545

[7] Mermin N D and Wagner H 1966 Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic heisenberg models *Phys. Rev. Lett.* **17** 1133-6
[8] Hohenberg P C 1967 Existence of long-range order in one and two dimensions *Phys. Rev.* **158** 383–6

[9] Kosterlitz J M and Thouless D J 1973 Ordering, metastability and phase transitions in two-dimensional systems *J. Phys. C: Solid State Phys.* **6** 1181

[10] Kosterlitz J M 1974 The critical properties of the two-dimensional xy model *J. Phys. C: Solid State Phys.* **7** 1046

[11] Berezinskii V L 1971 Destruction of long-range order in one-dimensional and two-dimensional systems having a continuous symmetry group I. classical systems *Sov. Phys.—JETP* **32** 493–500

[12] Stock S, Hadzibabic Z, Battelier B, Cheneau M and Dalibard J 2005 Observation of phase defects in quasi-two-dimensional Bose–Einstein condensates *Phys. Rev. Lett.* **95** 190403

[13] Hadzibabic Z, Krüger P, Cheneau M, Battelier B and Dalibard J 2006 Berezinskii–Kosterlitz–Thouless crossover in a trapped atomic gas *Nature* **441** 1118

[14] Krüger P, Hadzibabic Z and Dalibard J 2007 Critical point of an interacting two-dimensional atomic bose gas *Phys. Rev. Lett.* **99** 040402

[15] Schweikhard V, Tung S and Cornell E A 2007 Vortex proliferation in the Berezinskii–Kosterlitz–Thouless regime on a two-dimensional lattice of Bose–Einstein condensates *Phys. Rev. Lett.* **99** 030401

[16] Cladé P, Ryu C, Ramanathan A, Helmerson K and Phillips W D 2009 Observation of a 2d bose gas: from thermal to quasicondensate to superfluid *Phys. Rev. Lett.* **102** 170401

[17] Nitsche W H, Kim N Y, Roumpos G, Schneider C, Kamp M, Höfling S, Forchel A and Yamamoto Y 2014 Algebraic order and the Berezinskii–Kosterlitz–Thouless transition in an exciton-polariton gas *Phys. Rev. B* **90** 205430

[18] Bishop D J and Reppy J D 1978 Study of the superfluid transition in two-dimensional $^4$He films *Phys. Rev. Lett.* **40** 1727–30

[19] Prokof'ev N, Ruebenacker O and Svistunov B 2001 Critical point of a weakly interacting two-dimensional bose gas *Phys. Rev. Lett.* **87** 270402

[20] Simula T P, Lee M D and Hutchinson D A W 2005 Transition from the Bose–Einstein condensate to the Berezinskii–Kosterlitz–Thouless phase *Phil. Mag. Lett.* **85** 395–403

[21] Nazarenko S, Onorato M and Proment D 2014 Bose–Einstein condensation and Berezinskii–Kosterlitz–Thouless transition in the two-dimensional nonlinear schrödinger model *Phys. Rev. A* **90** 013624

[22] Shukla V, Brachet M and Pandit R 2013 Turbulence in the two-dimensional fourier-truncated Gross–Pitaevskii equation *New J. Phys.* **15** 113025

[23] Donnelly R J 1991 *Quantized Vortices in Helium* II vol 2 (Cambridge: Cambridge University Press)

[24] Bewley G P, Lathrop D P and Sreenivasan K R 2006 Superfluid helium: visualization of quantized vortices *Nature* **441** 588

[25] Zmeev D E, Pakpour F, Walmsley P M, Golov A I, Guo W, McKinsey D N, Ihas G G, McClintock P V E, Fisher S N and Vinen W F 2013 Excimers $he_2^*$ as tracers of quantum turbulence in $^4$He in the $t = 0$ limit *Phys. Rev. Lett.* **110** 175303

[26] La Mantia M and Skrbek L 2014 Quantum turbulence visualized by particle dynamics *Phys. Rev. B* **90** 014519

[27] Spethmann N, Kindermann F, John S, Weber C, Meschede D and Widera A 2012 Dynamics of single neutral impurity atoms immersed in an ultracold gas *Phys. Rev. Lett.* **109** 235301

[28] Catani J, Lamporesi G, Naik D, Gring M, Inguscio M, Minardi F, Kantian A and Giamarchi T 2012 Quantum dynamics of impurities in a one-dimensional bose gas *Phys. Rev. A* **85** 023623

[29] Michel C, Boughdad O, Albert M, Larré P-É and Bellec M 2018 Superfluid motion and drag-force cancellation in a fluid of light *Nat. Commun.* **9** 2108

[30] Carusotto I 2014 Superfluid light in bulk nonlinear media *Proc. R. Soc. A* **470** 20140320

[31] Amo A, Lefrère J, Pigeon S, Adrados C, Ciuti C, Carusotto I, Houdré R, Giacobino E and Bramati A 2009 Superfluidity of polaritons in semiconductor microcavities *Nat. Phys.* **5** 805

[32] Klein A, Bruderer M, Clark S R and Jaksch D 2007 Dynamics, dephasing and clustering of impurity atoms in Bose–Einstein condensates *New J. Phys.* **9** 411

[33] Krotscheck E, Saarela M and Epstein J L 1988 Impurity states in liquid–helium films *Phys. Rev. B* **38** 111–28

[34] Krotscheck E and Zillich R E 2008 Hydrogen and $^3$He atoms on $^4$He surfaces: bound states and scattering features *Phys. Rev. B* **77** 094507

[35] Rica S and Roberts D C 2009 Induced interaction and crystallization of self-localized impurity fields in a Bose–Einstein condensate *Phys. Rev. A* **80** 013609
[36] Winiecki T and Adams C S 2000 Motion of an object through a quantum fluid *Europhys. Lett.* **52** 257

[37] Shukla V, Brachet M and Pandit R 2016 Sticking transition in a minimal model for the collisions of active particles in quantum fluids *Phys. Rev. A* **94** 041602

[38] Shukla V, Pandit R and Brachet M 2018 Particles and fields in superfluids: insights from the two-dimensional Gross–Pitaevskii equation *Phys. Rev. A* **97** 013627

[39] Davis M J, Morgan S A and Burnett K 2001 Simulations of bose fields at finite temperature *Phys. Rev. Lett.* **87** 160402

[40] Connaughton C, Josserand C, Picozzi A, Pomeau Y and Rica S 2005 Condensation of classical nonlinear waves *Phys. Rev. Lett.* **95** 263901

[41] Krstulovic G and Brachet M 2011 Dispersive bottleneck delaying thermalization of turbulent Bose–Einstein condensates *Phys. Rev. Lett.* **106** 115303

[42] Berloff N G and Youd A J 2007 Dissipative dynamics of superfluid vortices at nonzero temperatures *Phys. Rev. Lett.* **99** 145301

[43] Krstulovic G and Brachet M 2011 Anomalous vortex-ring velocities induced by thermally excited kelvin waves and counterflow effects in superfluids *Phys. Rev. B* **83** 132506

[44] Krstulovic G and Brachet M 2011 Energy cascade with small-scale thermalization, counterflow metastability, and anomalous velocity of vortex rings in fourier-truncated Gross–Pitaevskii equation *Phys. Rev. E* **83** 066311

[45] Gardiner C 2009 *Stochastic Methods: a Handbook for the Natural and Social Sciences* (Springer Series in Synergetics) (Berlin: Springer)

[46] Foster C J, Blakie P B and Davis M J 2010 Vortex pairing in two-dimensional bose gases *Phys. Rev. A* **81** 023623

[47] Witten T A and Sander L M 1983 Diffusion-limited aggregation *Phys. Rev. B* **27** 5686–97