Statistical mechanics of one-dimensional quantum droplets
T. Mithun, S. I. Mistakidis, P. Schmelcher, P. G. Kevrekidis
SStatistical mechanics of one-dimensional quantum droplets
T. Mithun, S. I. Mistakidis, P. Schmelcher,
2, 3 and P. G. Kevrekidis Department of Mathematics and Statistics, University of Massachusetts, Amherst MA 01003-4515, USA Center for Optical Quantum Technologies, Department of Physics,University of Hamburg, Luruper Chaussee 149, 22761 Hamburg Germany The Hamburg Centre for Ultrafast Imaging, University of Hamburg,Luruper Chaussee 149, 22761 Hamburg, Germany
We study the statistical mechanics of quantum droplets in a one-dimensional ring geometry. Therelevant model of a modified two-component Gross-Pitaevskii equation under symmetry consider-ations, i.e. same particle numbers and equal intra-component interaction strengths, reduces to asingle-component equation. To determine the classical partition function thereof, we leverage thesemi-analytical transfer integral operator (TIO) technique. The latter predicts a distribution of theobserved wave function amplitudes and yields two-point correlation functions providing insights intothe emergent dynamics involving quantum droplets. We compared the ensuing TIO results withthe equilibrium properties of suitably constructed Langevin dynamics and with direct numericalsimulations of the original system’s modulational instability dynamics where the generated dropletsare found to coalesce. The results indicate good agreement between the distinct methodologiesaside from intermediate temperatures in the special limit where the droplet widens, approachinga plane-wave. In this limit, the distribution acquires a pronounced bimodal character at the TIOlevel, not captured by the Langevin dynamics, yet observed within the modified Gross-Pitaevskiiframework.
I. INTRODUCTION
The celebrated Gross-Pitaveskii model has proved par-ticularly successful for studying and describing a va-riety of macroscopic many-body phenomena in zero-temperature Bose-Einstein condensates (BECs) [1–3].On the other hand, more recently, a new type of mat-ter wave has emerged first in two-dimensional (2D) andthree-dimensional (3D) geometries [4, 5], namely theso-called quantum droplets. The theoretical basis isa two-component (binary) BECs, with intra-componentself-repulsion, yet also inter-component attraction thatslightly exceeds the self-repulsion. Here, the famous Lee-Huang-Yang quantum correction [6] comes into play toaccount for the averaged effect of quantum fluctuationsbeyond the mean-field description. It competes with themean-field effects [1, 2] and thus prevents the BEC col-lapse predicted in the mean-field realm. Notice that thisstabilization mechanism depends crucially on the sys-tem’s dimensionality; namely beyond mean-field fluctua-tions are attractive in one-dimension (1D) while being re-pulsive in higher spatial dimensions. Importantly, thesepredictions led to a sequence of experimental realizationsof such quantum droplets [7–10] in higher dimensions re-vealing, for instance, the transition from bright solitonsto quantum droplets [8], and thus rendering these statesan interesting emerging topic of study in the context ofBECs.Along this vein of experimental developments, it is re-markable that collisions of such droplet patterns havebeen experimentally observed recently (leading to merg-ers for slow collisions and quasi-elastic passage for fastones) [9]. Additionally and while the above examplesfocused on K droplets, heteronuclear binary BECs of Rb and K showcased very stable quantum droplets on time scales of the order of a second [10]. In par-allel to this steady stream of experimental demonstra-tions, theoretical studies have spearheaded a number ofparallel directions. These include but are not limited tothe exploration of quantum droplets with intrinsic vortic-ity [11], the impact of discreteness in the form of semidis-crete droplets with or without topological charge [12],or the exploration of 3D stable generalizations of suchstates [13]. Many of these developments have been re-cently summarized in Ref. [14]. Interestingly, the 1D dy-namical features of droplet states are far less explored andare currently mainly restricted to inelastic collisional as-pects of these configurations especially for high momentaand flat-top droplets [15] or their spontaneous generationdue to the modulation instability (MI) [16].The above efforts have mainly focused on the dynami-cal aspects of the droplets, and principally so at the zerotemperature setting. However, naturally, studying the fi-nite temperature dynamics of BEC systems is a topic ofbroad theoretical and experimental appeal [17]. This isparticularly intriguing in 1D where the role of quantumfluctuations, being inherently related to droplet forma-tion, is more prominent [18, 19] and three-body lossesare suppressed [20] compared to higher dimensions whilethe experimental probe of these configurations is still elu-sive. One of the commonly used methods to include thetemperature induced fluctuations is the so-called trun-cated Wigner method [21, 22]. Recently, a considerableamount of attention has been drawn towards the so-calledpositive P-method too [23, 24]. In general, the final taskof these methods is solving a stochastic equation that in-corporates the thermal fluctuations to the correspondingHamiltonian dynamical system. In statistical mechanics,a principal task consists of the evaluation of the parti-tion function corresponding to the energy functional of a r X i v : . [ c ond - m a t . qu a n t - g a s ] F e b the system at hand. A central tool to this effect in 1Dsettings consists of the so-called transfer integral opera-tor (TIO) method [25]. Indeed, this method renders theabove problem equivalent to solving the single particleSchr¨odinger equation. Moreover, it has been shown thata stochastic method represented by a suitable Langevinequation is comparable with the solution of TIO [26, 27].Further, recently it has been shown that the steady stateof molecular dynamics and the classical field method alsoreproduces the TIO results [27, 28].To the best of our knowledge, such techniques moti-vated by statistical mechanics, while being successful atthe level of single-component GP model [27], have yetto be explored in the quantum droplet realm of multi-component systems. Indeed, this is a central focus of thepresent work. More concretely, we study the statisticalmechanics of a quantum droplet in a 1D ring geometry.One of the peculiarities of a quantum droplet is its in-compressibility leading to a maximal critical density forsufficiently large particle numbers.This relevant system, as mentioned above, can be rep-resented with a modified binary Gross-Pitaveskii (MGP)equation, where the form of a Lee-Huang-Yang term isdetermined by the dimension of the system [5, 15, 29].For the case of particle-balanced components and equalintra-component interaction strengths, which will be ofprimary interest herein, the binary MGP equation re-duces to a single-component one. We determine the clas-sical partition function corresponding to our model bymapping the functional integration of the partition func-tion to a single-particle Schr¨odinger equation via the TIOtechnique. We then numerically verify the predictions ofTIO via a suitably crafted Langevin dynamics [26, 27].Our results indicate that the equilibrium properties of theLangevin dynamics are well in line with the TIO solutionat low and high temperatures for µ → µ , the limit wherethe droplets tend to disappear; the agreement is less ad-equate at intermediate temperatures. We also comparethe TIO results with the long-time dynamical evolutionof the original droplet system in the regime where thelatter falls into MI [30] and subsequently relaxes. The re-cent studies on the statistical properties of MI motivatethis analysis [31, 32]. Indeed, in Ref. [16], it has beenshown that the MI resulting from the small perturbationof a plane-wave state leads to the formation of quantumdroplet structures that undergo inelastic collisions. Wefind that, as a result of the inelastic collisions these gen-erated droplet structures coalesce and their equilibriumproperties are well matched by the TIO analysis.The work is organized as follows. We introduce theMGP model in Sec. II and then determine the classicalpartition function by using the TIO method in Sec. III.Section IV is devoted to developing and exploring therelevant Langevin dynamics. We report the results ofthe MI dynamics within the MGP framework in Sec. Vand provide an outlook in Sec. VI. II. THE MODIFIED GROSS-PITAEVSKIIFRAMEWORK
The starting point of our analysis will be the dimen-sionless MGP equation for a 1D binary quantum dropletin the form [5, 16], i ∂ψ ∂t = − ∂ ψ ∂z + ( P + GP − ) | ψ | ψ − (1 − G ) | ψ | ψ − P π (cid:112) P| ψ | + P − | ψ | ψ − µ ψ ,i ∂ψ ∂t = − ∂ ψ ∂z + ( P − + GP ) | ψ | ψ − (1 − G ) | ψ | ψ − π P (cid:112) P − | ψ | + P| ψ | ψ − µ ψ , (1)where the involved parameters stand for P ≡ (cid:114) g g , G = δgg P (1 + P ) and δg = g + g. (2)Here g > g > g < µ and µ are the chemical potentials of the individual componentsof the system. Additionally, g = √ g g and for theexistence of a quantum droplet δg (cid:28) g should hold.Thus, the parameter P quantifies the intra-componentinteraction imbalance, while G measures the deviationfrom the balance point of mean-field repulsion and at-traction where δg = 0. In a corresponding experiment, δg can be tuned using the Feshbach resonance tech-nique [7, 8]. Moreover, we remark that in Eq. (1), theunits of length, time and wave functions are expressed interms of (cid:126) / ( mg ), (cid:126) / ( mg ), and √ mg/ (cid:126) , respectively.The above coupled set of MGP Eqs. (1) can be reducedto a single-component equation under the assumptionthat P = 1 ( g = g ≡ g ) and ψ = ψ [5, 15]. Thissingle-component system has the form i ∂ψ∂t = − ∂ ψ∂z + δgg | ψ | ψ − √ π | ψ | ψ − µψ, (3)with the normalization condition (cid:90) + ∞−∞ n dz = N, n = | ψ ( z ) | , (4)where N is the scaled number of atoms and µ representsthe chemical potential. For δg/g >
0, Eq. (3) gives rise toan exact localized flat-top (FT) solution [15]. The latterrepresents a quantum droplet (QD) which originates fromthe balance between the effective cubic self-repulsion andquadratic attraction characterized by the central density n and the chemical potential µ , where n = 89 π (cid:18) gδg (cid:19) , µ = − π gδg , (5)respectively. In what follows we shall focus on the δg/g > δg/g <
0, see Ref. [16] for a more elaborated discussion.
III. STATISTICAL MECHANICS
To study the statistical mechanics of the models un-der consideration, as described by Eq. (1) and (3), weestablish the corresponding classical partition functionvia the TIO method, as indicated in the Introduction inline with Refs. [25–27]. In order to proceed, we considerthe QD on a ring with length L = M ∆ z , where M rep-resents the number of grid points and ∆ z denotes the(quite fine) spatial discretization. Then, the form of theclassical partition function reads Z = (cid:90) D ( ψ, ψ ∗ ) e − βF [ ψ,ψ ∗ ] , (6)where F = (cid:90) dz ( H − µN )= (cid:90) (cid:104) | ∂ z ψ | + 12 δgg n −
12 2 / π n / − µn (cid:105) dz, (7)is the free energy, β is the inverse temperature, H repre-sents the Hamiltonian corresponding to Eq. (1) or (3) and N is the total number of particles. In accordance with theTIO methodology, the functional integration can be re-duced to an eigenvalue problem [25, 33]. For the Eqs. (1)and (3), we arrive at the corresponding eigenvalue equa-tions (see Appendix A for a detailed derivation) (cid:34) − β (cid:16) δ δψ + δ δψ (cid:17) + V d ( ψ , ψ ) (cid:35) φ n ( ψ , ψ )= E n φ n ( ψ , ψ ) , (8)and (cid:34) − β δ δψ + V d ( ψ ) (cid:35) φ n ( ψ ) = E n φ n ( ψ ) , (9)regarding the two-component and the one-componentproblems respectively. Evidently, Eqs. (8) and (9)correspond to one- and two-dimensional single-particleSchr¨odinger equations respectively. In these expressions, φ n represent the corresponding single-particle eigenfunc-tions, being functionals of the ψ field, while E n denotethe eigenvalues. Importantly, the effective TIO poten-tials appearing in Eqs. (8) and (9), possess the form V d ( ψ , ψ ) = P + GP − | ψ | + P − + GP | ψ | + ( G − | ψ | | ψ | − π (cid:18) P| ψ | + | ψ | P (cid:19) / − µ ( | ψ | + | ψ | ) (10) and V d ( ψ ) = 12 δgg | ψ | −
12 2 / π | ψ | − µ | ψ | . (11)The shape of these effective potentials in the differentparameter regions of the system is crucial in order to un-derstand where bound state solutions, and thus droplet-like configurations, are prone to appear. For this rea-son, we next show V d ( ψ = u ) in Fig. 1 since the one-component case will be our focus in the following. A rel-evant discussion about the two-component effective po-tential [Eq. (10)] is provided in Appendix B. Accord-ing to Eq. (11) the minimum of V d ( ψ = u ) is givenby u = (3 b + (cid:112) b + 32 aµ ) / (8 a ) for µ >
0, while u = (3 b ± (cid:112) b + 32 aµ ) / (8 a ) for µ ≤ a = δgg and b =
12 2 / π . As a result for µ < μ =- μ μ = μ =+ μ - - - V d ( u ) Figure 1. Effective anharmonic potential V d ( u ), described byEq. (11), of the single-component case for δg/g = 0 .
05. Shownare different values of the chemical potential µ (see legend)with µ given by Eq. (5). We solve the Eqs. (8) and (9) numerically to determinethe corresponding eigenvalues and eigenvectors utilizingexact diagonalization. Notice that this is the only as-pect that renders the TIO approach semi-analytical. Themain contribution to the partition function in Eq. (6) willbe from the lowest TIO eigenvalue. The thermodynamicproperties of an equilibrium system can be well under-stood by resorting to the underlying probability distri-bution of amplitudes (PDA), P ( | φ | = u ). The latter canbe expressed in terms of the eigenfunction correspondingto the lowest eigenvalue E of Eq. (9) [27]. In particular P ( | φ | = u ) acquires the form P ( | φ | = u ) = 2 u | φ ( u ) | . (12)Moreover, the steady state properties of a thermody-namic system can be further characterized by employingthe two-point spatial correlation function C ( z ) [34, 35].Here, z denotes the relative distance between two distinctspatial locations. This correlation function can be writ-ten with respect to the eigenvalues and eigenfunctions ofEq. (9) as C ( z ) = (cid:104) ψ ( . ) ψ ( . + z ) (cid:105) = (cid:88) n | (cid:90) duφ ∗ n ( u ) uφ ( u ) | e − β | z | ( E n − E ) . (13)It provides a measure of the coherence among the parti-cles [34, 36]. Concretely, it is bounded from above andbelow taking as a maximal and minimal value the con-sidered particle number and zero respectively. If C ( z ) ismaximal for every z then the system is termed fully co-herent and it is characterized by quasi long-range order.Otherwise losses of coherence occur. In the thermody-namic equilibrium, the long-range order is expected tovanish and accordingly C ( z ) tends to zero for increasing z [37]. We finally remark that, for the settings consideredherein, the contributions of the (higher) excited states ofthe effective potentials (described by Eqs. (12) and (13))are neglected, as being exponentially smaller than thedominant lowest-lying states. IV. LANGEVIN DYNAMICS
It has been demonstrated that, for any positive tem-perature, the exact results obtained from the TIO in aBEC system can be compared with that of a Langevinequation with a Gaussian additive white noise term [27].This type of stochastic dynamics contributes towards re-laxing the system configuration to the free-energy mini-mum, while at the same time accounting for the thermalfluctuations arising around this minimum [38]. In the fol-lowing, we consider the Langevin equations correspond-ing to Eqs. (1) and (3), in line also with the earlier workof [27] for regular one-component BECs ∂ψ ∂t = − (cid:104) − ∂ ψ ∂z + ( P + GP − ) | ψ | ψ − (1 − G ) | ψ | ψ − P π (cid:112) P| ψ | + P − | ψ | ψ − µ ψ (cid:105) + ξ ( z, t ) ,∂ψ ∂t = − (cid:104) − ∂ ψ ∂z + ( P − + GP ) | ψ | ψ − (1 − G ) | ψ | ψ − π P (cid:112) P − | ψ | + P| ψ | ψ − µ ψ (cid:105) + ξ ( z, t ) , (14) and ∂ψ∂t = − (cid:104) − ∂ ψ∂z + δgg | ψ | ψ − √ π | ψ | − µ (cid:105) ψ + ξ ( z, t ) , (15)respectively. The term ξ i represents the white noise hav-ing a Gaussian distribution with the correlation (cid:104) ξ ∗ i ( z, t ) ξ i ( z (cid:48) , t (cid:48) ) (cid:105) = 2 β δ ( z − z (cid:48) ) δ ( t − t (cid:48) ) . (16)The Langevin equations are indeed found to relax to theequilibrium state at large evolution times. This is causedby the fact that the time average of the spatio-temporalcorrelation relaxes to its equilibrium spatial correlationor in other words the time average of the noise approachesits equilibrium distribution [39]. We solve Eqs. (14) and(15) numerically by using the xmds package [40]. Figure 2. Probability distribution function obtained fromTIO (solid lines) and within the Langevin dynamics for dif-ferent values of the inverse temperature β . Namely, β = 0 . β = 0 . β = 2 (redsquares). The other parameters refer to µ = µ + 0 . δg/g = 0 . Our focus here will be in the single-componentLangevin case. For completeness, we will present selectedresults for the two-component case in Appendix B. Notealso that for the numerical simulations of Eq. (15), to bepresented below, a sample of 1000 trajectories is foundto be sufficient for the convergence of the relevant prob-ability distributions, with the domain size being L = 30.We first fix the chemical potential µ = µ + 0 . P ( | ψ | = u )results of the Langevin dynamics for (long) total evolu-tion time t F ≥ β = 0 . β = 0 . β = 2 . Figure 3. Probability distribution function P ( | ψ | = u ) atdifferent total evolution times t F (see legend) for β = 0 . P ( | ψ | = u ) is observed for t F > µ = µ + 0 . δg/g = 0 . respectively; the solid lines represent the correspondingTIO solutions. Inspecting P ( | ψ | = u ) [Fig. 2] we canreadily deduce that as the temperature increases (from β = 2 to β = 0 . β = 0 . β = 0 . t F that allows to reach thesteady state, we determine the P ( | ψ | = u ) for different t F at β = 0 .
5, see in particular Fig. 3. Interestingly,we observe that the dynamics at initial times developsa bimodal distribution, being more proximal to the ex-pected picture from the TIO analysis, but as time in-creases the weight of one of the two humps reduces andeventually the distribution approaches a unimodal shape.We will further address this discrepancy when consider-ing the dynamics of the full model, Eq. (3), (rather thanthe modified Langevin one) in the following section. To gain additional insight on the dependence of the steadystate distribution at intermediate temperatures, we nowvary the chemical potential µ by fixing the inverse tem-perature parameter β = 0 .
5. The obtained results forthe P ( | ψ | = u ) are provided in Fig. 4. It shows thatthe TIO solution develops a bimodal distribution onlywhen µ → µ . As µ deviates from this FT QD limit, theTIO solution leads to a unimodal distribution and theresults of the Langevin dynamics match well with thesesolutions. This once again suggests that presumably theLangevin dynamics is not able to capture the delicate FTconfiguration of the QD at the intermediate temperaturelimit, as we will argue further below. Nevertheless, in allother settings, the TIO semi-analytical prediction is veryadequately captured by the modified Langevin dynamics. Figure 4. Probability distribution P ( | ψ | = u ) obtained fromthe Langevin dynamics for µ = µ + 0 .
25 (blue crosses), µ = µ + 0 . µ = µ + 0 .
01 (red squares), µ = µ + 0 .
005 (black triangles). The solid lines representthe corresponding TIO solutions. The other parameters are δg/g = 0 .
05 and β = 0 . To gain further insight into the emergent steady stateof the QDs, we additionally determine the two-point cor-relation function C ( | z | ) [Eq. (13)]. The latter is pro-vided in Fig. 5 for sufficiently large evolution times, i.e., t F ≥ C ( | z | ) → z . This fact also sup-ports the appearance of a steady state [37] for the droplet.Note that coherence losses are more dramatic for large β ,compare lower and upper panels of Fig. 5. On the otherhand, at the intermediate inverse temperature ( β = 0 . µ → µ limitas seen in the middle panel of Fig. 5. Concretely, the TIOmethod shows a vanishing coherence for larger values of z , whilst the Langevin dynamics predicts that the coher-ence of the system is maintained. Figure 5. Correlation function C ( z ) obtained from the TIOmethod (solid lines) and from the Langevin dynamics (dots)for different temperatures. Namely, β = 0 . β =0 . β = 2 (bottom panel) at µ = µ +0 . t F ≥ V. MI DYNAMICS IN THE MODIFIEDGROSS-PITAEVSKII APPROACH
Recently, it has been shown that the dynamics of aclassical field method (CFM) adequately traces the fea-tures of the TIO solution for the equilibrium statisti-cal mechanics of the single-component BEC case [27].Indeed, a key point of that work is the comparison ofthree different methods: the Langevin approach detailedabove with a molecular dynamics (MD) one [28] (basedon Hamiltonian dynamics of Klein-Gordon type) and theCFM. Considering the MD methodology as less directlyrelated to the system at hand, we focus here on theCFM as a complement of the Langevin approach pre-sented above. The core idea of the CFM is to obtain asteady state where dynamics is governed by the Gross-Pitaevskii equation and then compare this steady statesolution with the TIO solution. In this method, the ini-tial state is a non-equlibirium one where only a few modesin momentum space are excited. Moreover, it is wellknown that in such systems the MI phenomenon natu-rally creates non-equilibrium conditions when subjectedto a small perturbation (seeding the relevant instabil-ity) [30, 41]. Since the Gross-Pitaevskii equation with aneffective attractive interaction is modulationally unsta- ble and the MI dynamics of the model as described byEq. (3) has been recently studied in detail [16] we shallsubsequently focus on the equilibrium properties gener-ated by direct numerical simulations associated with theunstable dynamics stemming from the MI.As an initial condition we consider a homogeneous par-ticle density distribution subjected to a weak amplitudeperturbation of the form δψ = A p e iθ r . The strength A p of the perturbation is of the order of 10 − and θ r repre-sents the uniform noise distribution taking values in theinterval [0 , π ] [16]. We fix n = 9 and δg/g = 0 .
05 thatcorresponds to µ ≈ µ + 0 . n of the FT shaped droplet. The correspond-ing PDA is depicted in Fig. 6(b) (blue dotted points).Interestingly, the PDA exhibits a bimodal distribution,where the peak at large amplitude ( u ≈ .
8) representsthe large amplitude droplets, while peaks at small ampli-tude ( u ≈ .
92) are related to small-amplitude phononicexcitations [29] that are widespread within the spatialextension of the cloud.To compare this equilibrium state with the TIO so-lution, we considered different values of β and obtainedthat, e.g., for β = 0 .
62 the probability distribution func-tions are in close correspondence with each other asshown in Fig. 6(b) with the green solid line. Interest-ingly, the accordance between the PDA and TIO solu-tion predictions holds also for the corresponding corre-lation function in momentum space [42] i.e. (cid:104)| ψ ( k ) | (cid:105) =(1 / π ) (cid:82) dzC ( z ) e − ikz . The latter is routinely accessiblein current ultracold atom experiments via time-of-flightimaging [43]. This observable is provided in Fig. 6(c) at t F = 4000 and β = 0 .
62. As can be seen, it features acentral peak structure around k = 0 signaling the steadystate of the QD while its shape is in perfect agreementbetween the two approaches. To further corroborate thatthe same argumentation holds also for other values of thechemical potential, we exemplarily showcase in Fig. 7 thecases of µ = µ + 0 . µ = µ + 0 .
25. It becomes ap-parent that even for these values of the chemical poten-tial, the final state resulting from the MI dynamics is alarge excitation on top of a fluctuating background bear-ing small amplitudes. Hence, the corresponding PDAis a bimodal distribution. However, the weight of thesmall amplitude excitations in the corresponding PDAdistribution decreases with increasing µ , in line with therespective TIO prediction shown in Fig. 4. | ψ | 0 1000 2000 3000 4000t-40-20 0 20 40 z (a) Figure 6. MI dynamics for the system parameter values δg/g = 0 .
05 and n = 9. (a) Spatiotemporal evolution of | ψ ( z, t ) | . (b) Probability distribution function P ( | ψ | = u )and (c) correlation function in momentum space (cid:104)| ψ ( k ) | (cid:105) ata final evolution time of t F = 4000. The results shown in(b) and (c) are averaged ones over a sample of twelve dif-ferent initial conditions. The green solid line represents theTIO solution at β = 0 .
62, while, the blue dotted points areobtained through the direct simulation of the MGP Eq. (3).Here, µ = µ + 0 . | ψ |-40-20 0 20 40 z (a) | ψ | 0 1000 2000 3000 4000t-40-20 0 20 40 z (b) Figure 7. Time-evolution of | ψ ( z, t ) | explicating the MI dy-namics, the subsequent generation of QDs and their coales-cence when considering δg/g = 0 .
05 while (a) µ = µ + 0 . µ = µ + 0 . VI. CONCLUSIONS AND FUTURECHALLENGES
In the present work, we have sought to addressthe equilibrium statistical mechanics of one-dimensionalquantum droplets. This system is represented by a mod-ified two-component Gross-Pitaevskii equation and un-der the symmetry considerations (equal wave functionsand equal intra-component interactions strengths) it re-duces to a single-component equation. To determine theclassical partition function, we first reduce the functionalintegration of the partition function to a single-particleSchr¨odinger problem by using the transfer integral oper-ator (TIO) technique. The probability distribution of theamplitudes obtained from TIO for the single-componentcase shows unimodality in the limit µ → µ at large andsmall temperatures; however, at intermediate tempera-tures, the relevant probability exhibits a bimodal dis-tribution. On the other hand, for the values of µ farfrom µ , the probability distribution is always found tobe unimodal. We briefly also consider the binary systemin the Appendix below, which, in turn, leads to a two-dimensional system. Our results therein suggest similarlya predominantly unimodal distribution.We compared the results of the TIO analysis withthe equilibrium properties of suitably crafted Langevindynamics. The equilibrium properties of the Langevindynamics reproduce well the corresponding of the TIOsolutions except in the limit µ → µ at an intermedi-ate temperature for the single-component case, wherea bimodal distribution is featured within the TIO ap-proach. Finally, we contrasted the probability distri-bution of the modulational instability dynamics of thesingle-component scenario at large times with the TIOfindings, obtaining good agreement in the cases underconsideration for a suitable choice of the temperature.Very recently, significant attention has been drawn tothe thermodynamics of quantum droplets [44]. Beyondthe realm of the present study, several topics remainopen for future studies. An extension of thermodynamicconsiderations in higher dimensional settings would be atopic of particular interest for a variety of reasons. On theone hand, semi-analytical tools such as the TIO do notstraightforwardly generalize to higher dimensions, hencedifferent theoretical approaches would need to be broughtto bear. Moreover, even numerically the nonlinearities ofthe system being logarithmic in 2D and featuring a dif-ferent form than the one considered herein ( ∝ | ψ | ψ ) in3D could yield different outcomes as regards the proba-bility distributions and the correlation functions. On theother hand, the genuinely multi-component case wherethe two components are not forced to be equal and itssystematic consideration and direct comparison e.g. witha variational approach [36] is also of interest for furtherstudy. ACKNOWLEDGMENTS
We thank S. Flach and K. Burnett for fruitful com-ments and discussions. S. I. M. gratefully acknowledgesfinancial support in the framework of the Lenz-IsingAward of the University of Hamburg. The present paperis based on work that was supported by the US NationalScience Foundation under DMS-1809074 (PGK).
Appendix A: Derivation of the transfer integralproblem
Here, we briefly discuss the derivation of the single-particle Schr¨odinger Eq. (9), describing the equilibriumstate of the symmetric mixture, from the underlying clas-sical partition function of Eq. (6) [25]. For a ring of length L , we can discretize the free energy functional introducedin Eq. (7) as F = ∆ z M (cid:88) i =1 f ( ψ i , ψ i +1 ) = ∆ z M (cid:88) i =1 (cid:104) | ψ i +1 − ψ i ∆ z | + 12 δgg n i −
12 2 / π n / i − µn i (cid:105) , (A1)and we obtain the partition function [25] Z = M (cid:89) i =1 (cid:90) ∞−∞ d ˜ ψ (cid:48) d ˜ ψ i e − β ∆ zf ( ψ i ,ψ i +1 ) δ ( ˜ ψ − ˜ ψ (cid:48) ) , (A2)where d ˜ ψ i = (cid:113) β π ∆ z dψ i and M is the number of seg-ments within the ring of width ∆ z . We now expand the δ function in terms of a normalized set of (complete)eigenfunctions and obtain Z = (cid:88) n M (cid:89) i =1 (cid:90) ∞−∞ d ˜ ψ (cid:48) d ˜ ψ i φ n ( ˜ ψ (cid:48) ) e − β ∆ zf ( ψ i ,ψ i +1 ) φ n ( ˜ ψ ) . (A3)We can calculate this partition function from the eigen-value problem (cid:90) ∞−∞ d ˜ ψ i e − β ∆ zf ( ψ i ,ψ i +1 ) φ n ( ˜ ψ i ) = e − β ∆ zE n φ n ( ˜ ψ i +1 ) , (A4)where φ n and E n are the eigenfunctions and eigenvaluesof the transfer matrix equation. We then reduce thisequation to the single-particle Schr¨odinger problem asfollows. Performing the Taylor series expansion of theeigenfunction φ n ( ˜ ψ i ) around ( ˜ ψ i +1 ) we arrive at (cid:90) ∞−∞ d ˜ ψ i e − β ∆ zf ( ψ i ,ψ i +1 ) φ n ( ˜ ψ i ) = e − β ∆ z (cid:104) δgg n −
12 25 / π n / − µn (cid:105) × (cid:16) z β ∂ ∂ψ i +1 (cid:17) φ n ( ˜ ψ i +1 )= e − β ∆ zH φ n ( ˜ ψ i +1 ) ≡ e − β ∆ zE n φ n ( ˜ ψ i +1 ) (A5)where (cid:16) ∆ z β ∂ ∂ψ i +1 (cid:17) ≈ e ∆ z β ∂ ∂ψ i +1 . The effective Hamil-tonian reads H = − β ∂ ∂ψ + 12 δgg n −
12 2 / π n / − µn. (A6) Here, the second term corresponds to the mean-field con-tribution and the third one is the leading-order beyondmean-field correction accounting for quantum fluctua-tions. Figure 8. Two-dimensional probability distribution function P [R e ( ψ ) = u , R e ( ψ ) = u ] obtained from the Langevindynamics (left panels) and the TIO approach (right panels)for β = 0 . β = 0 . β = 2(bottom panels). Other parameters used are µ = µ +0 . δg/g = 0 .
05, and P = 1. The total evolution time is t F = 240. Appendix B: Two-component Langevin equations
For the numerical simulation of Eq. (14) we consider1000 trajectories with the domain size L = 8 ×
8. Wefirst set µ ≈ µ ≡ µ + 0 . P = 1. Further,by considering two different noise distributions ( ξ and ξ ) for the two components, we break the symmetry be-tween them. Figure 8 depicts the results of the Langevindynamics (left panels) and the TIO solutions (right pan-els). Comparing the corresponding 2D probability distri-bution functions P [R e ( ψ ) = u , Re( ψ ) = u ] betweenthe Langevin dynamics and the TIO findings for variousinverse temperatures we deduce that they are in a goodagreement. Particularly, as the temperature increasesthe distribution widens similarly to the single-componentcase. But in sharp contrast to the single-component sce-nario, the distribution always remains unimodal for thecase examples considered herein (although, by necessitywe have considered fewer cases in the present setting andcannot exclude the possibility that under different pa-rameters bimodality may arise). Furthermore, in orderto estimate the impact of an increasing population asym-metry we change P = 1 to P = 1 .
25. The correspondingPDAs for β = 0 . Figure 9. Two-dimensional probability distribution function P [ Re ( ψ ) = u , Re ( ψ ) = u ] as predicted from the Langevindynamics (left panel) and the TIO method (right panel) for β = 0 . P = 1 .
25 as well as µ = µ +0 . t F = 240.[1] C. J. Pethick and H. Smith, Bose–Einstein condensationin dilute gases (Cambridge university press, 2008).[2] S. Stringari and L. Pitaevskii, Bose-Einstein Condensa-tion (Oxford University Press, Oxford, United Kingdom,2003).[3] P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-Gonz´alez, The defocusing nonlinear Schr¨odinger equa-tion: from dark solitons to vortices and vortex rings (So-ciety for Industrial and Applied Mathematics, Philadel-phia, PA, 2015).[4] D. S. Petrov, “Quantum mechanical stabilization of a col-lapsing bose-bose mixture,” Phys. Rev. Lett. , 155302 (2015).[5] D. S. Petrov and G. E. Astrakharchik, “Ultradilutelow-dimensional liquids,” Phys. Rev. Lett. , 100401(2016).[6] T. D. Lee, K. Huang, and C. N. Yang, “Eigenvaluesand eigenfunctions of a bose system of hard spheres andits low-temperature properties,” Phys. Rev. , 1135(1957).[7] C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas,P. Cheiney, and L. Tarruell, “Quantum liquid dropletsin a mixture of bose-einstein condensates,” Science ,301 (2018). [8] P. Cheiney, C. Cabrera, J. Sanz, B. Naylor, L. Tanzi, andL. Tarruell, “Bright soliton to quantum droplet transitionin a mixture of bose-einstein condensates,” Phys. Rev.Lett. , 135301 (2018).[9] G. Ferioli, G. Semeghini, L. Masi, G. Giusti, G. Mod-ugno, M. Inguscio, A. Gallem´ı, A. Recati, and M. Fat-tori, “Collisions of self-bound quantum droplets,” Phys.Rev. Lett. , 090401 (2019).[10] C. D’Errico, A. Burchianti, M. Prevedelli, L. Salasnich,F. Ancilotto, M. Modugno, F. Minardi, and C. Fort,“Observation of quantum droplets in a heteronuclearbosonic mixture,” Phys. Rev. Research , 033155 (2019).[11] Y. Li, Z. Chen, Z. Luo, C. Huang, H. Tan, W. Pang,and B. A. Malomed, “Two-dimensional vortex quantumdroplets,” Phys. Rev. A , 063602 (2018).[12] X. Zhang, X. Xu, Y. Zheng, Z. Chen, B. Liu, C. Huang,B. A. Malomed, and Y. Li, “Semidiscrete quantumdroplets and vortices,” Phys. Rev. Lett. , 133901(2019).[13] Y. V. Kartashov, B. A. Malomed, L. Tarruell, andL. Torner, “Three-dimensional droplets of swirling su-perfluids,” Phys. Rev. A , 013612 (2018).[14] Z. Luo, W. Pang, B. Liu, Y. Li, and B. A. Malomed, “Anew form of liquid matter: quantum droplets,” Frontiersof Physics , 32201 (2021).[15] G. E. Astrakharchik and B. A. Malomed, “Dynamics ofone-dimensional quantum droplets,” Phys. Rev. A ,013631 (2018).[16] T. Mithun, A. Maluckov, K. Kasamatsu, B. A. Mal-omed, and A. Khare, “Modulational instability, inter-component asymmetry, and formation of quantumdroplets in one-dimensional binary bose gases,” Symme-try , 174 (2020).[17] A. Griffin, T. Nikuni, and E. Zaremba, Bose-condensedgases at finite temperatures (Cambridge UniversityPress, 2009).[18] L. Parisi, G. Astrakharchik, and S. Giorgini, “Liq-uid state of one-dimensional bose mixtures: a quantummonte carlo study,” Phys. Rev. Lett. , 105302 (2019).[19] L. Parisi and S. Giorgini, “Quantum droplets in one-dimensional bose mixtures: a quantum monte-carlostudy,” arXiv:2003.05231 (2020).[20] G. E. Astrakharchik and S. Giorgini, “Correlation func-tions of a lieb–liniger bose gas,” J. Phys. B: At. Mol. andOpt. Phys. , S1 (2006).[21] P. B. Blakie, A. Bradley, M. Davis, R. Ballagh, andC. Gardiner, “Dynamics and statistical mechanics ofultra-cold bose gases using c-field techniques,” Advancesin Physics , 363 (2008).[22] A. Sreedharan, S. Choudhury, R. Mukherjee,A. Streltsov, and S. W¨uster, “Collisions of soli-tary waves in condensates beyond mean-field theory,”Phys. Rev. A , 043604 (2020).[23] P. Deuar and P. Drummond, “First-principles quantumdynamics in interacting bose gases: I. the positive p rep-resentation,” J. Phys. A: Mathematical and General ,1163 (2006).[24] P. Deuar and P. Drummond, “First-principles quantumdynamics in interacting bose gases ii: stochastic gauges,”J. Phys. A: Mathematical and General , 2723 (2006).[25] D. J. Scalapino, M. Sears, and R. A. Ferrell, “Statisticalmechanics of one-dimensional ginzburg-landau fields,”Phys. Rev. B , 3409 (1972). [26] A. Khare, S. Habib, and A. Saxena, “Exact thermody-namics of the double sinh-gordon theory in 1+ 1 dimen-sions,” Phys. Rev. Lett. , 3797 (1997).[27] A. Nunnenkamp, J. N. Milstein, and K. Burnett, “Clas-sical field techniques for condensates in one-dimensionalrings at finite temperatures,” Phys. Rev. A , 033604(2007).[28] G. Aarts, G. F. Bonini, and C. Wetterich, “On thermal-ization in classical scalar field theory,” Nuclear PhysicsB , 403 (2000).[29] M. Tylutki, G. E. Astrakharchik, B. A. Malomed, andD. S. Petrov, “Collective excitations of a one-dimensionalquantum droplet,” Phys. Rev. A , 051601 (2020).[30] J. H. Nguyen, D. Luo, and R. G. Hulet, “Formation ofmatter-wave soliton trains by modulational instability,”Science , 422 (2017).[31] A. Gelash, D. Agafontsev, V. Zakharov, G. El, S. Ran-doux, and P. Suret, “Bound state soliton gas dynam-ics underlying the spontaneous modulational instability,”Phys. Rev. Lett. , 234102 (2019).[32] A. E. Kraych, D. Agafontsev, S. Randoux, and P. Suret,“Statistical properties of the nonlinear stage of modu-lation instability in fiber optics,” Phys. Rev. Lett. ,093902 (2019).[33] S. E. Trullinger and R. M. DeLeonardis, “Analytic statis-tical mechanics for a two-component-kink system,” Phys.Rev. B , 5522 (1980).[34] M. Naraschewski and R. Glauber, “Spatial coherence anddensity correlations of trapped bose gases,” Phys. Rev.A , 4595 (1999).[35] G. C. Katsimiga, G. M. Koutentakis, S. I. Mistakidis,P. G. Kevrekidis, and P. Schmelcher, “Dark–bright soli-ton dynamics beyond the mean-field approximation,”New J. Phys. , 073004 (2017).[36] S. I. Mistakidis, G. C. Katsimiga, P. G. Kevrekidis, andP. Schmelcher, “Correlation effects in the quench-inducedphase separation dynamics of a two species ultracoldquantum gas,” New J. Phys. , 043052 (2018).[37] R. Nandkishore and D. A. Huse, “Many-body localiza-tion and thermalization in quantum statistical mechan-ics,” Annu. Rev. Condens. Matter Phys. , 15 (2015).[38] G. Parisi, Y. S. Wu, et al., “Perturbation theory withoutgauge fixing,” Sci. Sin , 483 (1981).[39] P. C. Hohenberg and B. I. Halperin, “Theory of dynamiccritical phenomena,” Rev. Mod. Phys. , 435 (1977).[40] G. R. Dennis, J. J. Hope, and M. T. Johnsson, “Xmds2:Fast, scalable simulation of coupled stochastic partial dif-ferential equations,” Computer Physics Communications , 201 (2013).[41] P. Everitt, M. Sooriyabandara, M. Guasoni, P. Wigley,C. Wei, G. McDonald, K. Hardman, P. Manju, J. Close,C. Kuhn, et al., “Observation of a modulational insta-bility in bose-einstein condensates,” Phys. Rev. A ,041601 (2017).[42] S. I. Mistakidis, T. Wulf, A. Negretti, and P. Schmelcher,“Resonant quantum dynamics of few ultracold bosons inperiodically driven finite lattices,” J. Phys. B: At. Mol.and Opt. Phys. , 244004 (2015).[43] I. Bloch, J. Dalibard, and W. Zwerger, “Many-bodyphysics with ultracold gases,” Rev. Mod. Phys. , 885(2008).[44] G. De Rosi, G. E. Astrakharchik, and P. Massignan,“Thermal instability, evaporation and thermodynamicsof one-dimensional liquids in weakly-interacting bose-1