Equation of state and self-bound droplet in Rabi-coupled Bose mixtures

Laser induced transitions between internal states of atoms have been playing a fundamental role to manipulate atomic clouds for many decades. In absence of interactions each atom behaves independently and their coherent quantum dynamics is described by the Rabi model. Since the experimental observation of Bose condensation in dilute gases, static and dynamical properties of multicomponent quantum gases have been extensively investigated. Moreover, at very low temperatures quantum fluctuations crucially affect the equation of state of many-body systems. Here we study the effects of quantum fluctuations on a Rabi-coupled two-component Bose gas of interacting alkali atoms. The divergent zero-point energy of gapless and gapped elementary excitations of the uniform system is properly regularized obtaining a meaningful analytical expression for the beyond-mean-field equation of state. In the case of attractive inter-particle interaction we show that the quantum pressure arising from Gaussian fluctuations can prevent the collapse of the mixture with the creation of a self-bound droplet. We characterize the droplet phase and discover an energetic instability above a critical Rabi frequency provoking the evaporation of the droplet. Finally, we suggest an experiment to observe such quantum droplets using Rabi-coupled internal states of K39 atoms.

In atomic physics, laser beams can stimulate transitions among different hyperfine states. Inspired by remarkable experiments with both fermionic 1,2 and bosonic clouds 3,4 , in the recent years an extensive theoretical research was devoted to understand static and dynamical properties of quantum mixtures with artificial coupling beween their internal states. Concerning fermionic mixtures, for example, in the attempt to search for itinerant ferromagnetism driven by Rabi coupling, it was shown that a critical coupling frequency marks the transition of a two-state Fermi gas to a ferromagnetic phase. A detailed investigation in three spatial dimension was performed by Conduit 5 and, very recently, for a two-dimensional Fermi gas 6 . On the other side, for bosonic atoms at temperatures below the transition to the superfluid phase, coupling of hyperfine states offers the possibility to address fascinating phenomena such as the internal Josephson effect 7-9 emulating a space dependent double well potential, analogues of the Hawking radiation 10,11 , non-abelian gauge potentials 12 like magnetic monopoles 13,14 , Rashba spin-orbit coupling [15][16][17][18] , or they can be used for applications to quantum metrology [19][20][21] and for the quantum simulation of spin models with short or long-range interactions [22][23][24][25] .
In this article we study the effects of a Rabi coupling on a two-component Bose mixture deriving the corresponding beyond-mean-field equation of state. To achieve this result we perform a non-trivial regularization of Gaussian fluctuations, which have a divergent zero-point energy due to both gapless and gapped elementary excitations. In particular, we obtain a meaningful analytical formula for the ground-state energy of the Bose mixture as a function of Rabi coupling and scattering lengths. Setting the Rabi frequency to zero in our formula one recovers Larsen's equation of state 26 . In the case of attractive inter-particle interaction we investigate the conditions for the formation of a self-bound droplet finding that its density profile and collective oscillations crucially depend on the interplay between Rabi coupling and interaction strengths. A similar equation of state, albeit in absence of internal coupling, has been recently used by Petrov 27,28 . He shows that, in the case of negative inter-component scattering length, quantum fluctuations can arrest the collapse of the mixture inducing the

Results
Microscopic theory for Rabi-coupled mixtures. We consider a Bose gas with two relevant hyperfine states in a volume L 3 , at temperature T and with chemical potential µ. In addition to the usual intra-and interstate contact interactions, transitions between the two states are induced by an an external coherent Rabi coupling of frequency ω R . We adopt the path integral formalism, where each component is described by a complex bosonic field ψ i ( = i 1,2). Given the spinor ψ ψ Ψ = ( , ) T 1 2 [35][36][37] , the partition function of the system reads: where the Euclidean action Ψ Ψ S[ , ] is given by being a ij the scattering length for collisions between component i and component j (specifically a 11 , a 22 , and a 12 ). All relevant thermodynamical quantities can be derived from the grand potential  Ω = − β ln( ) 1 . We work in the superfluid phase, where a U(1) gauge symmetry of each bosonic component is spontaneously broken. The presence of the Rabi coupling in the Euclidean action in equation (2) implies that only the total number of atoms is conserved. We can then set ψ τ where v i are the uniform order parameters of the two-component Bose-Einstein condensate, and η τ r ( , ) i are the fluctuation fields above the condensate. The mean-field plus gaussian approximation is obtained by expanding equation (2) up to the second order in η τ (r, ) . The corresponding beyond-mean-field grand potential is then given by 37,38 is the mean-field grand potential, while µ Ω v v ( , , ) g 1 2 is the grand potential of Gaussian quantum and thermal fluctuations.
In our scheme, the Bose-Einstein order parameters v i satisfy the saddle-point equations , leading to coupled equations for the uniform and constant fields v 1 and v 2 : . One finds a symmetric configuration where the two internal states are equally populated, a polarized phase with non-zero population imbalance, and an unstable phase when the attractive inter-state interaction overcomes the intra-state repulsion < − g g and equal intra-component interaction. The corresponding mean-field grand potential µ Ω v ( , ) 0 is then given by By solving equation (5) in the case of symmetric ground-state, we get the crucial relation between the order parameter and the chemical potential: 12 . In this case, equation (6) reduces to µ µ ω It is important to stress that we work in a regime where Rabi coupling cannot produce polarization in the ground state. However, as shown in the following section, it still influences the stability of balanced configuration, i.e. the region between the symmetric and unstable phase in the diagram reported in Fig. 1 (top panel), also when Gaussian fluctuations are taken into account.

Gaussian Fluctuations. To compute
for the symmetric ground state and for equal interaction strengths, we consider the quadratic terms in η i and η ⁎ i of equation (2). In reciprocal Fourier space one finds Here ω { } n n are the bosonic Matsubara frequencies and  ω − q ( , ) n is the × 4 4 inverse of the fluctuations propagator, whose definition is reported in the Methods. At zero temperature, the Gaussian grand potential corresponds to the zero-point energy of bosonic excitations and it reads 37,43 is the spectrum of elementary excitations, which can be obtained by diagonalizing 35,37,44 . The diagonal blocks of  are two-by-two identity matrices 1 2 , while the off-diagonal ones are the Pauli matrix σ z . The eigenvalues are the two branches of the Bogoliubov spectrum: where we set ε = a a / 12 , with a the intra-component scattering length and a 12 the inter-component scattering lengt, , the zero-temperature Gaussian grand potential is ultraviolet divergent. We employ the convergence-factor regularization 37,43,45 which generates proper counterterms in the zero-point energy completely removing the divergence. These counterterms can be determined by expanding the two branches of the Bogoliubov spectrum at high momenta. The zero-temperature beyond-mean-field grand potential is then given by equation (7) plus the regularized zero-point energy, namely  (4). In the symmetric ground state the two components appear with the same particle density ν ν | | = | | 1 2 2 2 , whereas in the polarized phase densities are unequal. Dotted lines represent the asymptotic phase boundaries of the polarized region for large g g / 12 and µ ω / R  ratios respectively.
the symmetric solution is unstable in the thermodynamic limit. The grey region for µ ω <  (12) and its mean field approximation µ Ω ( ) 0 (red solid line) of equation (7) as a function of the chemical potential µ for = . g g 0 9 12 within the symmetric phase.
In Fig. 1 (bottom panel) we plot the grand potential µ Ω( ) of equation (12), including gaussian fluctuations, as a function of the chemical potential for = . g g 0 9 12 . We compare it with the mean field approximation µ Ω ( ) 0 of equation (7). The energy density of the system is  where the number density n is obtained via = − µ ∂Ω ∂ n L 1 3 . In the limit of small Rabi-coupling, which is also the most relevant experimentally 10 (see below), it is possible to get an analytical result for the energy density. By as energy unit (then ω ω = E R R B ) and defining the diluteness parameter = n na 3 , up to the linear term in ω R , from equation (12) we obtain the scaled energy density of the mixture components with constant densities: one recover the Larsen's zero-temperature equation of state 26 . From equation (14) one finds that for ε | |>1 the uniform configuration is not stable. If ε > 1, at the mean field level, one expects phase separation or population imbalance 39 . Instead, if ε < − 1 the term proportional to ε + n [(1 ) ] 5/2 becomes imaginary. This imaginary term does not induce dynamical instability, but only dissipation. As for other sources of losses (for instance three-body recombination), to study short-time dynamics this dissipative term can be neglected if n is not too large. The resulting real energy density displays a characteristic n 5/2 dependence which competes with the negative mean-field contribution, opening the door to the possibility of observing a droplet phase for finite systems. This stabilization mechanism based on quantum fluctuations has been proposed for the first time in two-component mixtures without Rabi coupling 27,28 and recently applied to dipolar condensates [32][33][34] . For ε < −1 the equilibrium density is obtained upon the minimization of the energy density in equation (14) with respect to n neglecting the imaginary term: The solution − n is a local maximum, while the equilibrium value is given by + n which is a local minimum of the energy per particle. Moreover to obtain a real solution, Rabi frequency is limited by: 4 . For larger ω R there is only the absolute minimum with zero energy at = n 0.
Droplet phase. For a finite system of N of particles we define a space-time dependent complex field φ t ( , ) 2 is the space-time dependent local number density, and clearly ∫ = N d n t r r ( , ) 3 . The dynamics of φ t r ( , ) is driven by the following real-time effective action where  ND is obtained from equation (14) neglecting the imaginary term proportional to ε + n [(1 ) ] 5/2 . In the inset of Fig. 2 we plot the density profile of the stationary solution obtained by numerically solving with imaginary time-evolution the Gross-Pitaevskii equation associated to equation (16) varying the number of particles for ω π = /2 1 R kHz. The solution indeed corresponds to a self-bound spherical droplet whose radial width increases by increasing the number of atoms. For a very large number of atoms, the plateau of the density profile approaches the thermodynamic density given by equation (15). Instead, for a small number of atoms the self-bound droplet does not exist.
One can model the droplet by using a Gaussian wavefunction where the particle number is = ∼ N Nn. By inserting equation (17) in the rescaled version of equation (16), one gets six Euler-Lagrange equations for the 46 is a variational energy functional which is function of the width of the droplet only (see Methods). The variational stability diagram of the droplet phase is illustrated in Fig. 2. Upon increasing the atom number droplets stabilize. For small particle numbers we find a metastable region where σ σ σ    U( , , ) 1 2 3 has a local minimum with positive energy, the global minimum corresponding to zero energy for a dispersed gas with zero density. Interestingly, tuning the Rabi coupling to large values, as shown with the red dashed line for = N 1200 particles in Fig. 2, we move into the unstable phase. Therefore, differently from dipolar gases 32 or bosonic mixtures with attractive inter-species interactions 27 , where transition to the instability is driven by interactions, here, a direct coupling between the two components serves as an additional tunable knob to cross from a stable into an unstable phase.
The low-energy collective excitations of the self-bound droplet are investigated by solving the eigenvalues problem for the Hessian matrix of effective potential energy in equation (31). From the form of the variational ansatz we naturally describe the monopole (breathing) mode of frequency ω M and the quadrupole mode of frequency ω Q .
The upper panel of Fig. 3 displays monopole and quadrupole frequencies as a function of the number N of atoms in the droplet, fixing Rabi coupling and scattering lengths. The lower panel of Fig. 3 reports the collective frequencies as a function of the Rabi coupling and two different values of N. Both frequencies go to zero at the Rabi coupling above which the droplet evaporates.
The experimental observation of a droplet phase with Rabi coupled internal states is within experimental reach. A promising candidate is a gas of 39

Discussion
We derived the beyond-mean-field grand potential of a Rabi-coupled bosonic mixture within the formalism of functional integration, and performing regularization of divergent Gaussian fluctuations. In the small Rabi-coupling regime we also obtained an analytical expression for the internal energy of the system. In the case of attractive inter-particle scattering length we have shown how the Gaussian terms of the internal energy help to stabilize the system against the collapse and that, for a finite number of atoms, a self-bound droplet is produced. Rabi coupling works as an additional tool to tune the stability properties of the droplet, inducing an energetic instability for large inter-component couplings. The evaporation of the droplet is also signaled by both the breathing and quadruple modes which vanish at a critical Rabi coupling. Notably, our predictions provide a benchmark for experimental observations of Rabi-coupled self-bound droplets in current experiments.  (31). We observe three phases: a stable droplet-phase region (light green) of spherical self-bound droplets, a metastable droplet phase (yellow) where the energy of the droplet is positive and larger than a uniform background with vanishing density, and an unstable (white region) for small particle number N or high Rabi coupling ω where droplet evaporate. Here we consider ε | + | = . 1 05 which corresponds to ω .  31 8 c kHz. In the inset we plot the three dimensional density profile n r ( ) of droplets from the numerical solution of the Gross-Pitaevskii equation for different particle numbers at ω π = /2 1 R kHz, from the metastable region = N 977 and gaussian density limit = ⋅ N 5 10 3 to the Thomas-Fermi regime = ⋅ N 2 10 4 and = N 10 5 where system density is roughly constant up to a critical droplet radius. Moving along the vertical axis, increasing the Rabi coupling, droplets become metastable and finally unstable. Red dashed line refers to a system of = N 1200 particles (see Methods).

Methods
Mean-field phase diagram. Our description of mean-field phase diagram starts from the mean-field free-energy density, where we identified = v n a a 2 The equilibrium configuration has to stationarize the energy density, namely In order to clarify the stability of these equilibrium points, one has to compute the determinant of the free-energy Hessian matrix (we assume an intra-species repulsion). Over the symmetric ground state, one finds and by imposing it to be positive, the following stability condition 39 In the symmetric ground-state, the density in terms of µ reads  µ ω = + + n g g /2 ( ) /( ) R ab , so we easily derive equation (7). In the polarized ground-state, since µ = n g / , the normalized imbalance equals The results of stability analysis of the stationary points of mean-field free energy are summarized in Fig. 1 (top panel). (8) is (see e.g. Refs 35,37 ), where: . The off-diagonal blocks of  ω q ( , ) n are given by

Convergence-factor regularization technique.
Among the available methods to regularize the zero-point energy of bosonic excitations 43 , we employ the convergence-factor technique. It consists in adding a factor ω + e i 0 n before performing the Matsubara summation contained in equation (8); this notation has to be intended as a limit procedure, i.e δ ω δ → + lim e i 0 n . As an example, we consider single-state bosons 37 with energy ε. The partition function is given by ^ † β ε µ ψ ψ = − − Tr[exp( ( ) )], and in the path-integral framework, one formally writes Since the path-integral relies on a time-axis discretization, this notation introduces an ambiguity 37 . We are not specifying on which time slice the field ψ τ ⁎ ( ) (corresponding to the operator ψ † ) acts. If ψ τ ( ) acts on the time slice τ i , then we can choose that ψ τ δ + ( ) acts on the τ + i 1 one, and δ → + 0 is needed to specify this prescription. In the Fourier space, it corresponds to the appearance of the convergence factor presented at the beginning of this subsection. If one instead chooses the opposite time-ordering, by taking the limit δ → − 0 , one can then verify that the corresponding partition functions differ only in a β ε µ − − e ( ) . factor. Then, from a partition function as the one in equation (28), with fields computed at the same time, we need to add a term such as ε µ − ( ) 1 2 in the grand potential to take this fact into account. This justifies the first two counterterm appearing in equation (13).
However, this equal-time prescription does not completely remove the ultraviolet divergences in the zero-point energy of bosonic excitations. The remaining ones are due to the presence of a gapped spectrum branch and to the zero-range approximation for the interaction potential. The grand potential of our system with Gaussian contribution is given by ) n is defined in equation (25). Gapped excitations induce a branch cut in the complex logarithm on the real q-axis at higher energy, corresponding to the single-particle continuum of states. This divergent contribution is removed as shown by Diener et al. 45 , giving rise to the third counteterm in equation (13). Finally, the fourth counterterm arises by making use of standard scattering theory at the second order 44   Indeed, one of the divergences encountered integrating the zero-point energy is due to the delta-shaped potential 37 . Its Fourier transform is constant for all momenta, while a reasonable interaction potential should fall at least as q 1/ 2 , giving back a finite contribution.
Variational and numerical analysis. The equation for σ  i is the classical equation of motion for a particle of coordinates σ σ σ σ =     ( , , ) T 1 2 3 moving in an effective potential given by the derivative of the potential energy per particle: where α = π 128 75 5 7/4 and γ = π 112 9 3 5/4 . The energy per particle of the ground state is simply σ where σ  m is the minimum of the effective potential energy. In absence of an external trapping, the system preserves its spherical symmetry, i.e. the critical point of the effective potential in equation (31) Figure 4 shows the energy per particle of the self-bound droplet: the numerical approach relying on imaginary-time evolution is in reasonable agreement with the variational one based on equation (17). Remarkably, above a critical Rabi frequency the internal energy of the droplet becomes positive, signaling that the droplet goes in a metastable configuration. Moreover, at a a slightly larger critical Rabi frequency the droplet evaporates.
Data availability. Data are available upon request. Requests should be addressed to either author.