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

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 of their respective values is observed, whereas their ratio remains approximately constant.

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 two-and three-dimensional 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 twopoint 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 70's by Berezinskii, Kosterlitz and Thouless (BKT) [9][10][11].The BKT transition is an infiniteorder 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 dipoles (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][13][14][15][16], exciton-polaritons [17], liquid helium films [18] and studied theoretically and nu-merically [19][20][21]; for a review on the topic, see for instance [3].
The purpose of this Letter is to study how the statistical mechanical properties of a two dimensional quantum fluid of bosons are affected by the presence of impurities.Particles and impurities have been used in superfluids since their early experiments in 4 He [22] mainly with detection purposes: electrons, ions and neutral impurities such as hydrogen particles and excimers have been exploited to visualize quantized vortices, to study their dynamics and the statistics of superfluid (quantum) turbulence [23][24][25].More recently, the investigation of the interaction between one or more impurities and a quantum fluid has been the main topic of experiments in cold atoms [26][27][28], quantum fluids of light [29,30] and polaritons in semiconductor microcavities [31].On the theoretical side, by using a quantum master equation, Klein et al. [32] studied the dynamics of impurities in a BEC and Rica & Roberts studied how a collection of impurities affects the ground-state of a BEC by using a mean field model [33].
We investigate here impurity clustering and phase transitions occurring in a minimal model describing this system: the Gross-Pitaevskii (GP) equation coupled with active immiscible impurities having classical degrees of freedom.Such model was introduced in [34] and recently used in two-dimensional numerical simulations to study impurity-impurity and impurity-vortex interactions [35,36].It also mimics one of the states observed in [33], where for large coupling values impurity fields separate from the condensate, giving rise to an immiscible phase.Finite-temperature BECs can be studied by using the truncated (or 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 truncated GP model is an effective model to study the condensation transition in two and three di-mensions [21,[37][38][39] and superfluid vortex dynamics at finite temperature [40][41][42].
We generalize the truncated GP model to include the dynamics of active impurities.The model is then described by the Hamiltonian where ψ is the collective wave-function of bosons having mass m, and g = 4πa s 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 q i and 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 rep is a repulsive potential between impurities.The Galerkin projector P G truncates the system acting in Fourier space as P G ψk = θ(k max − k) ψk with θ(•) the Heaviside function, ψk the Fourier transform of ψ(x) and k the wave vector.The equations of motion are directly obtained by varying (1) and read where q ij = |q i − q j | and we have replaced The previous set of equations exactly conserves the Hamiltonian, the number of bosons N = |ψ| 2 dx and momentum P = i 2 ψ∇ψ − ψ∇ψ dx + i p i .The impurities in the system feel an attractive force mediated by the quantum fluid density field [33,35].However, unlike the case of impurities described by classical fields [33], in the model ( 1) no repulsion mediated by the fluid exists.In order to mimic such a repulsion and avoid an unnatural overlap of the impurities, we consider the Lennard-Jones-like repulsive potential V rep (r) = r 12 min /r 12 as in [35].For the impurities potential we use a smoothed hat-potential , where a I sets the characteristic radius of the impurity and ∆l is a smoothing parameter.Finally, let us notice that in absence of impurities and at zero temperature, eq.( 2) can be linearized about a uniform density state |ψ| 2 = ρ ∞ /m, defining the phonon (sound) velocity c = gρ ∞ /m 2 with dispersive effects taking place for length scales smaller than the healing length ξ = 2 /2gρ ∞ .We integrate the system (2-3) by using a pseudospectral code with N res uniform grid points per direc-tion of a squared domain of size L = 2π.We set k max = N res /3 so that the truncated system exactly conserves all the invariants (provided that initially P G ψ = ψ and P G V I = V I ) [42], c = ρ ∞ = 1, V 0 = 10 and = 0.00674.As the healing length changes with temperature, we parametrize the solutions of (2-3) using its value taken at zero temperature.In thermal states, the only relevant dimensionless parameters are L/ξ, a I /ξ, N I , the relative mass M = M I /ρ ∞ πa 2 I and ξk max .We set L/ξ = 128.The value of ξk max controls the strength of the nonlinear interactions and it is kept fixed to ξk max = 2π/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 [37], that is ξk max ∼ 2; such value applies to a gas of 87 Rb atoms.
Let us start by presenting a long temporal evolution of a system having 40 impurities of mass M = 0.1, characteristic radius a I = 4ξ, initially located at random positions (avoiding overlaps) and having zero velocity.The quantum fluid density field of the initial condition is displayed in Fig. 1.a.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 (Fig. 1.b); 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 even further (Fig. 1.c).Eventually, the system reaches thermal equilibrium where only one big cluster is observed in a bath of thermalized waves, (Fig. 1.d).A movie of the full temporal evolution is available in the Supplemental Material.
The evolution illustrated in Fig. 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 very costly and has the inconvenience of not directly providing access to the conjugate thermodynamical variables: temperature and chemical potential in this case (here we only consider zero momentum states).To overcome these issues, in [42] a stochastic relaxation was introduced in order to efficiently generate thermal states in the grandcanonical ensemble.We make use of this approach adapting it to the Hamiltonian (1).The stochastic dynamics is ruled by where F = H − µN is the free energy of the system, µ is the chemical potential controlling the number of bosons and β is the inverse temperature; ξk , ξ q i and ξ q i are independent Gaussian white noises of unit variance.It can be shown by using the Fokker-Planck equation associated to the system (4-5), that the stationary probability distribution is given by the Gibbs grand-canonical distribution P[ ψk , q i , qi ] ∝ e −βF .In the micro-canonical ensemble, P[ ψk , q i , qi ] is also the stationary solution of the Liouville equation that describes the evolution of the phase-space distribution of the Hamiltonian system (1) [42].We define the temperature as T = 1/k N β, with k N = L 2 /N and N = πk 2 max the total number of Fourier modes.With this definition T is an energy per unit of surface such that at low temperatures F ≈ 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 ρ = mN/L 2 = 1 by dynamically adjusting the chemical potential [42].
We use (4-5) to study the effect of impurities on the quantum fluid condensed fraction n 0 defined as where • T stands for spatial and noise realizations average at temperature T .We first perform a temperature scan without impurities.The condensed fraction is shown in Fig. 2.a (solid blue line).As expected, n 0 decreases with temperature and then vanishes at the condensation temperature T λ [43].We then perform tem-perature scans varying the repulsive potential parameter r min with a fixed number of impurities N I = 31 having radius a I = 4ξ.The results are also shown in Fig. 2.a.Depending on the strength of the repulsion potential two different behaviors of n 0 can be observed.When the repulsion among impurities is strong enough (blue markers, r min ≥ 11ξ), clusters are broken at temperatures lower than T λ and impurities have no appreciable effect on the condensation curve (see Fig. 2.a.1).On the other hand, for r min ≤ 10.5ξ (red markers) impurities remain clustered and lead to an increase of the condensed fraction at medium-high temperatures (see Fig. 2.a.2).Impurities are attracted due to the density mediated interaction energy E I↔I computed in [33,35] and repelled because of V rep , generating cluster structures as the one observed in Fig. 1.d.However, if thermal fluctuations are large enough, the bound among impurities can be broken.In Fig. 2.b we compute the interaction energy between two impurities at zero temperature E I↔I , as a function of their distance ∆q (dotted black line).As a reference, the figure also displays the repulsive potential V rep with r min = 8ξ (dotted-dashed orange line).The sum of both potentials is displayed in the same figure for different values of r min : for sufficiently small values of r min , a potential well ∆U centered at r min appears.We thus expect a pair of impurities to split in finite time at the temperature T cl ∼ ∆U/k N as in a standard escape problem from a potential well [44].In order to quantify this clustering transition, we measure the average square distance between impurities and their center of mass q cm = N I j q j /N I :  many-body impurity interactions and by using the interaction potential obtained at T = 0. Note that for weak repulsion, even if the condensate vanishes, impurities still feel the density-mediated attraction.
The condensed fraction only changes when impurities remain clustered at all temperatures and hence this effect can not be simply explained by the local increase of density in regions not occupied by impurities.In order to quantitatively characterise the change in n 0 , we study how the condensation transition changes when varying the filling fraction Φ = 1 − | ψ T =0 | 2 / |ψ| 2 T =0 , which corresponds to the fraction of the total volume occupied by the impurities and accounts for the smooth density depletion near impurity boundaries.We fix r min = 2a I to ensure the formation of the cluster at all temperatures.In Fig. 3 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.
The condensation temperature T Φ λ is measured for different values of Φ: the relative increase ∆T λ = (T Φ λ − T 0 λ )/T 0 λ is displayed in Fig. 4. Remarkably, ∆T λ scales linearly with Φ growing up to 20%.We have checked by varying the number of impurities and their size, that n 0 only depends on Φ and T (data not shown).Despite the change on T λ , the condensed fraction curves collapse as expected to a single one, once plotted versus T /T Φ λ (see inset of Fig. 3).
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 g 1 (r) = ψ(0)ψ * (r) at the BKT transition temperature T BKT .At T < T BKT , it presents a power-law decay g 1 (r) ∼ r −α , where α depends linearly on the temperature; at high temperatures, it exhibits the standard exponential decay of disordered systems.In the inset of Fig. 4 we show g 1 (r) at two temperatures 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 [45].With no impurities and the parameters used in this Letter, the BKT transition takes place at T BKT = 0.83T λ .Note that because of the Mermin-Wagner-Hohenberg theorem [7,8], T λ vanishes as 1/ log L in the thermodynamic limit, so in principle for a very large system we could have T λ < T BKT .We do not address such limit in this Letter.The presence of impurities in the system modifies the decay of g 1 (r) by shifting T Φ BKT to higher temperatures.Figure 4 displays the relative increase ∆T BKT = (T Φ BKT − T 0 BKT )/T 0 BKT of the BKT transition temperature for different filling fractions Φ.Although T λ and T BKT both grow up to 20% when the Φ 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, decreasing its effective degrees of freedom [46].At a given temperature, with less active modes, the energy is smaller and thus a higher temperature is necessary to induce a transition.
In this Letter we studied thermal states of twodimensional quantum fluids 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 quantum fluids with specific types of 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.

FIG. 1 .
FIG. 1. (Color online) Snapshots of the density fluid during the GP temporal evolution of a state with 40 impurities.They occupy the region corresponding to the dark circles, where the density vanishes.(a) Initial condition.(b) Impurities start clustering at t = 815ξ/c.(c) Intermediate stage at t = 1691ξ/c.(d) Fully thermalized state at t = 10695ξ/c.The parameters used are M = 0.1 and aI = 4ξ.
7) Fig.2.c.displays δ 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

FIG. 2 .
FIG. 2. (Color online) (a): Condensed fraction as a function of the temperature for different values of rmin.The snapshots are density fields at T = 0.7 for rmin = 11ξ (a.1) and rmin = 8ξ (a.2). (b): Impurity-impurity interaction EI↔I as a function of their distance ∆q/ξ (dotted black line) and repulsive potential (dotted-dashed golden line) for rmin = 8ξ.Different markers correspond to total energy EI↔I + Vrep for different values of rmin.(c): Relative impurity distance (7) as a function of the normalized temperature T /T cl .Scans are performed with NI = 31 and aI = 4ξ.The ratios T cl /T λ are 7.41, 2.51, 1.5, 1.02, 0.53, 0.002 for rmin from 8ξ to 12ξ respectively.

FIG. 3 .
FIG. 3. (Color online): Condensed fraction as a function of temperature for different values of the filling fraction Φ. (Inset): Condensed fraction as a function of the temperature normalized with T λ .Scans with rp = 4ξ and rmin = 2rp = 8ξ.

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.Computations were carried out on the Mésocentre SIGAMM hosted at the Observatoire de la Côte d'Azur.