Zero-energy modes of two-component Bose-Bose droplets

Bose-Bose droplets are self-bound objects emerging from a mixture of two interacting Bose-Einstein condensates when their interactions are appropriately tuned. During droplet formation three continuous symmetries of the system's Hamiltonian are broken: translational symmetry and two U1 symmetries, allowing for arbitrary choice of phases of the mean-field wavefunctions describing the two components. Breaking of these symmetries must be accompanied by appearance of zero-energy excitations in the energy spectrum of the system recovering the broken symmetries. Normal modes corresponding to these excitations are the zero-energy modes. Here we find analytic expressions for these modes and introduce Hamitonians generating their time evolution -- dynamics of the droplet's centers of mass as well as dynamics of the phases of the two droplet's wavefunctions. When internal types of excitations (quasiparticles) are neglected then the very complex system of a quantum droplet is described using only few"global"degrees of freedom - the position of the center of mass of the droplet and two phases of two wave-functions, all these being quantum operators. This gives the possibility of describing in a relatively easy way processes of interaction of these quantum droplets, such as collisions.


I. INTRODUCTION
Quantum droplets, self-bound systems of ultracold dilute atomic clouds, were predicted to be formed from atomic mixtures of two bosonic components when repulsive intraspecies interactions are almost perfectly balanced by interspecies attraction [1]. Experimentally droplets were quite unexpectedly discovered in ultracold samples of atoms with large magnetic moments, first in a dysprosium system [2][3][4][5] and later in samples of erbium atoms [6]. Soon, as the theory suggested, droplets were observed in a two-component mixture of bosonic Potassium [7][8][9].
A common feature of all these ultradilute quantum systems is the stabilizing role of quantum fluctuations, essential to provide self-binding [1]. Their contribution to the system energy has the form of the celebrated Lee-Huang-Yang term (LHY) [10] which in a weakly interacting system comes into play only when the mean field energy almost vanishes. To be precise, the mean field energy should be slightly negative indicating overall effective attraction. The system is then on the collapse side of the stability diagram. The energy of the quantum fluctuations prevents this collapse and stable droplets can be formed. In systems with large magnetic dipole-dipole interactions the droplets can be created even in a single component arrangement, [2]. This is because both attractive and repulsive interactions exist simultaneously at ultralow energies even in a monoatomic case. Their relative strength can be tuned to reach conditions convenient for droplet formation.
The stabilizing role of quantum fluctuations is essential in case of short range as well as dipolar interactions. In the latter case the LHY energy was investigated not only in uniform systems, [11][12][13], but also in optical lattices [14]. Energy of quantum fluctuations depends on the density of Bogoliubov modes [15,16] which, in turn, is a function of the dimensionality of the phase space involved. Bose-Bose droplets in lower dimensions [17] and at dimensional crossover [18,19] were investigated theoretically. Similarly, dipolar droplets at the crossover were studied [20]. Quantum fluctuations in low dimensional dipolar systems or systems at a dimensional crossover exhibit many interesting features [21]. In a strongly interacting regime novel quantum droplets are predicted [22]. Transition from bright solitons to self-bound droplets was observed experimentally in an axially extended geometry [8].
Ultradilute liquid droplets provide an unique platform to study many quantum phenomena, even such exotic processes as disruption of a white dwarf star by a black hole, [23]. Unlike trapped systems droplets can be pushed against each other to make them collide. Dynamics of such collisions may be very rich because of possible merging or fragmentation, the analogues of fusion and fission reactions being so far the domain of nuclear physics [24]. Binary collisions in the context of quantum liquid droplets were investigated theoretically in 1D geometry [25]. Collisions of droplets composed of two hyperfine states of 39 K were studied in the experiment of [26].
The superfluid character of quantum droplets introduces an additional degree of freedom, namely a phase. To observe the quantum effects related to the phase of a droplet wavefunction some kind of interference phenomenon should take place. Coupling of two superfluid systems with different phases is the essence of the Josephson effect and coherent oscillations. This kind of process leads to a transfer of nucleons between the two colliding nuclei [27,28]. On the other hand, superfluid coupling of droplets is at the heart of recently observed supersolid systems, a bizarre marriage of a rigid solid and a fluid. These states of matter were created in a weakly coupled droplet system [29][30][31]. To this end external trap geometry and interactions have to be appropriately tuned to trigger the so-called roton instability [32] appearing when a roton minimum in the excitation spectrum touches zero energy. Coherent arrays of droplets linked by a superfuid created experimentally are a great experimental exemplification of a supersolid system. Other phenomena related to the phase of droplets are also of interest. Droplets of swirling superfluids [33] or metastable clusters of droplets located on a ring and experiencing rich dynamics depending on relative phases of the droplets and the radius of the ring [34] are predicted theoretically.
The mean field description of a droplet involves five free parameters which do not affect its energy. These are the two phases of the droplet's wavefunctions and one 3D vector of the droplet position. This vector is typically set to zero as it is convenient to locate the droplet at the center of a coordinate system. However when two or more droplets are considered their positions might become dynamical variables. Similarly, arbitrary phases of the two wavefunctions are ignored in the ground state. The dynamics of phases have to be accounted for in addition to the standard Bogoliubov excitations when low energy excitations are considered. These additional excitations related to fluctuations of parameters not affecting the system energy are called zero-energy modes.
Zero-energy modes need some special treatment. They do not appear automatically as a result of diagonalization of the Bogoliubov-de Gennes operator. This operator is non-hermitian and its eigenvectors do not span the entire Hilbert space. The system we are studying here resembles to a large extent solitons in a 1D nonlinear medium.
Their global phases and positions are the parameters which do not affect the energy. Soliton's zero-energy modes are studied in [35]. The formalism developed for zero-energy solitonic modes occurred to be very useful in the description of greying of a dark soliton [36], Anderson localization of a soliton or localization of two interacting solitons [37,38].
The results of these earlier works cannot be automatically adopted to the droplet's case. Though conceptually the problem is very similar, the droplets are described by two coupled wavefunctions, and not by a single one. This fact brings quite a lot of additional issues, some of them of technical nature, though altogether they make the problem of zero-energy modes of a two component droplet nontrivial and give rise to quite a lot of interesting aspects.
In this paper we define a formalism to describe the dynamics of the center of mass of a droplet as well as phases of its wavefunctions. The paper is organized as follows. First, in Sec.II we briefly introduce the main formalism of description of two-component Bose-Bose droplets. In Sec.III we present the main idea of zero-energy modes in a simplified case where hard mode excitations are frozen and a single wavefunction can be assigned to the droplet. This allows to introduce the main ideas of our approach. A two component droplet is considered in Sec.IV where we find general expressions for the two phase modes as well as for the translational mode. In addition to this general formalism we present explicit forms of Hamiltonians generating dynamics of phases and fluctuations of particle number of Bose-Bose droplets. In the last section we summarize our results and conclude.

II. EXTENDED GROSS-PITAEVSKII EQUATIONS
Although Bose-Bose droplets are stabilized by beyond mean field effects, the corresponding Lee-Huang-Yang (LHY) energy can be easily incorporated into the mean field formalism. The energy density of the two-component Bose-Bose systems with additional LHY energy has the form: where n i are densities of the components and m is the mass of the atom. We assume that masses of atoms of both components are identical, m 1 = m 2 = m. This assumption is not crucial but allows for an explicit analytic form of the LHY term and, as a consequence, for approximate but analytic solutions. The first term in Eq.(1) accounts for mean field interactions where g ij are interaction strengths related to the scattering length a ij in the standard way, g ij = (4π 2 /m)a ij . Intraspecies interactions are repulsive, a ii > 0, but interspecies attraction, a 12 = a 21 < 0, slightly dominates this repulsion so the "effective scattering length" related to the effective strength of the mean field interactions, δa = − a 12 + √ a 11 a 22 > 0, is positive. We define δa in such a way that it enters the energy functional with a 'minus' sign in front, −δa. This convention will directly visualize the character of interactions. The mean field interactions alone would lead to a collapse of the system. The last contribution to the energy in Eq.(1) is the LHY term accounting for quantum fluctuations. It can prevent the mean-field collapse and stabilize the system for small negative values of −δa < 0.
The stationary extended Gross-Pitaevskii equations (EGP) can be obtained by minimizing the full energy density functional with both interaction and kinetic energies included: We introduced above the wavefunctions of the droplet's components ψ i , which are related to atomic densities n i = |ψ i | 2 , and are normalized to a number of particles of a given kind, dr|ψ i | 2 = N i . These constraints on the number of particles are accounted for by introducing two Lagrange multipliers ε i -the chemical potentials of the two species. The stationary EGP equations resulting from the minimization procedure, δH d /δψ * i = 0, are: where by µ i we denote the mean field interaction energy of one atom (from the i-th component) with all other particles in the droplet: Obviously the two equations Eqs.(4) are coupled because µ 1 (n 1 , n 2 ) and µ 2 (n 1 , n 2 ) depend on both involved densities. For given interaction strengths the only parameters which determine stationary solutions are atom numbers of the two species. The corresponding time dependent EGP equations are as follows: They can be reduced to stationary equations, Eqs. (4) via the standard substitution ψ i (r, t) = exp[−iε i t/ ]ψ i (r). Stationary solutions of these equations describing stable self-bound droplet exist only for some particular values of N 1 /N 2 which must be chosen from a well defined range. The best known standard estimate of this value [1] gives quite a good approximation to the ratio of densities at the ground state of the droplet: This approximation works best if a 11 ≃ a 22 . By imposing the constraint Eq.(7) the largest contribution to the energy functional Eq.(2) is eliminated -the first term accounting for hard mode excitations vanishes. This assumption limits the formalism to low energy states of the droplet. A more subtle analysis of the stability condition of two-component droplets is presented in [39][40][41] where a correction due to the hard mode contribution is accounted for. Some small deviations from condition of [1] are possible but they are not important for our further analysis.
Thus both densities are related. To avoid excitations of the hard mode any infinitesimal changes of the densities also have to be related, i.e dn 2 = dn 1 /s. The energy density becomes a function of one independent variable only, say n 1 . Then, an infinitesimal change of energy density due to small changes of atomic densities is dE = µ 1 dn 1 + µ 2 dn 2 = µ 1 + 1 s µ 2 dn 1 . The interpretation of this equation is straightforward. When removing one atom from the first component one has to remove also 1/s atoms from the second component to avoid excitation of an energetically costly hard mode. Consequently the energy density changes by an appropriate weighted sum of µ 1 and µ 2 . We might therefore think that one atom of the first kind corresponds to 1 + 1 s 'effective atoms of a droplet'. Alternatively, choosing n 2 as an independent variable, one atom of the second kind corresponds to 1 + s 'effective atoms of a droplet'. Let us observe that the number of droplet atoms is N = 1 + 1 s N 1 = (1 + s)N 2 . It is not surprising that it equals to N = N 1 + N 2 .
Note that to suppress energy of the hard mode one should set ψ 1 = √ sψ 2 . Both functions become thus proportional. They differ only by a normalization factor. This fact allows to introduce a condensate wavefunction ψ proportional to the wavefunctions ψ 1 and ψ 2 : Neglecting hard mode contribution to the energy and introducing a common wavefunction ψ forced us to set the relative phase of the two wavefunctions ψ 1 and ψ 2 to zero. This is a consequence of the singe wavefunction approximation. When substituting the two wavefunctions by a single one there is no room for two different phases anymore.
Only the amplitude and the global phase of the droplet wavefunction ψ are dynamical variables. The corresponding density of the droplet is n = |ψ| 2 . It is normalized to the total number of particles dr n = N = N 1 + N 2 .
Using Eqs.(8), (9) we get the total energy density in the single mode approximation: where interaction energy density is: Minimization of this functional leads to the single component EGP equation for the droplet wavefunction: The first term of Eq. (12) gives the kinetic energy of a droplet atom. The interaction energy µ and the chemical potential ε of a single atom have the form of a weighted sum of the potential energies or chemical potentials of the two species respectively: The above expressions as well as the EGP equation Eq. (12), are symmetric with respect to exchange of the two components, i.e. the simultaneous interchanging of indices 1 ⇐⇒ 2 and substituting s → 1/s.

III. ZERO-ENERGY MODES -SINGLE WAVEFUNCTION APPROXIMATION
A stable droplet solution of the stationary equations Eq.(4), if it exists, breaks three symmetries of the Hamiltonian of the system. These are the two U1 symmetries corresponding to the freedom of choice of phases of each of the two mean fields involved: ψ 1 → e −iθ1 ψ 1 , and ψ 2 → e −iθ2 ψ 2 , and translation of both fields by an arbitrary vector r 0 , As a result there should be three zero-energy modes in the excitation spectrum of the system. These excitations are analogues of the so called Goldstone modes [42,43]. The first two are the scalar modes describing changes of the wavefunction due to a phase shift of the droplet's wavefunctions [35]. The third one has 3D vector character -it describes a shift of the whole droplet. The energy of a droplet does not depend on its position, neither on the phases of its wavefunctions. But once fixed they do not necessarily stay constant. There are excitations which couple dynamically to the zero-energy modes. These are adjoint modes, [35]. Excitation of these modes increases the energy of the system and triggers dynamics of the zero-energy modes amplitudes. These dynamically coupled modes are different from Bogoliubov modes, they are missing in the eigenspectrum of the Bogoliubov-de Gennes operators (BdG). Our goal in the following is to find explicit expressions for the zero-energy modes. Our approach is similar to the approch of [44] where zero modes of a Bose-Einstein condensate were studied, however we follow closely the method of Ref. [35] where phase and translational zero energy modes of dark soliton excitations are analyzed.
Following [1] we argue that at low energies the hard mode is frozen. A simplified description based on a single EGP equation Eq.(12) reduces the number of zero-energy modes. The relative phase of the two wavefunctions is equal to zero and only a common phase of the two components and a center of mass position of the droplet are dynamical variables. We start our discussion with this simplified situation. This way we can introduce basic concepts and methods which will be explored further when studying the problem of zero-energy modes in its full extent. We introduce the extended Gross-Pitaevskii Hamiltonian H GP by rewriting Eq. (12) in the following way: and then Elementary excitations of a droplet can be found by studying a linear response of the system to a small perturbation, δ(r, t), of the stationary solution: By inserting the above expression into Eq.(15) and keeping terms which are linear in the small perturbation δ and δ * we get the time dependent Bogoliubov-de Gennes (BdG) equations: To get the stationary BdG equations we postulate the ansatz: which leads to the eigenvalue problem of the BdG operator L : where w Bν = (u Bν , v Bν ) T are right eigenvectors of the BdG operator L which is a '2 × 2' matrix: where n = |ψ| 2 . The operator L is non-hermitian. Its eigenvectors do not span the entire functional space [16,35,45]. Moreover its left eigenvectors (adjoint) w ad Bν are not simply hermitian conjugates of its right eigenvectors. Due to symmetries of L the left eigenvectors are of the form w ad The norm of the Bogoliubov modes is a standard scalar product of a right vector and its left (adjoint) partner, At this point we shall introduce our notational convention. The vector components w = (u, v) can be arranged in a row, as in the example discussed, or in a column. We will denote the latter arrangement by the transposition symbol Our discussion of the standard Bogoliubov modes of energies ε Bν = 0 does not apply to the zero-energy modes, ε Bν = 0, the case we focus on here. From now on, without any loss of generality, we assume that the droplet wavefunction ψ is real. So is the Bogoliubov-de Gennes operator L. If the phase of the function ψ is shifted, the function is modified as ψ → e −iϑ0 ψ. For infinitesimally small ϑ 0 the phase-shifted state is: Obviously the complex conjugate wavefunction gets the opposite phase, iϑ 0 ψ. This observations allows to find the zero-energy phase mode: Indeed, it is the eigenvector of the BdG operator corresponding to the zero eigenvalue: It is a rather simple exercise to check that v 0 is orthogonal to all Bogoliubov excitations w ad Bν |v 0 = 0. At this point a problem appears. Although the adjoint vector v ad 0 = (ψ, ψ) satisfies the equation: v ad 0 L = 0, the norm of the vector v 0 is zero: v ad 0 |v 0 = 0.
To span the basis of the functional space (δ, δ * ) we have to find yet another vector u 0 which is dynamically coupled to the zero-energy vector v 0 . Moreover, this coupled vector u 0 should have a unit overlap with the zero-energy mode, u ad 0 |v 0 = 1, and similarly v ad 0 |u 0 = 1. The vectors: v 0 and u ad 0 as well as vectors: u 0 and v ad 0 are the two pairs of adjoint vectors.
Let us assume that initially a perturbation of the stationary state is proportional to the zero-energy phase-mode, i.e. δ(r, 0) = −iϑ 0 ψ(r), where −iϑ 0 is some complex amplitude: where for convenience we introduced the two component vector of perturbations δψ = (δ, δ * ) T . Because v 0 is real, Eq.(24) can be satisfied only if ϑ 0 is a real number.
Evidently the initial state δψ(r, 0) = −iϑ 0 v 0 does not change with time, it remains constant: i ∂ t (−iϑ 0 v 0 ) = L(−iϑ 0 v 0 ) = 0. However, it might vary if there exists a vector which contributes to the zero-energy mode when evolved in time by the BdG operator L. We assume that the initial state of the field is proportional to some vector u 0 with a real amplitude ∆N 0 , i.e. we have δψ(r, 0) = ∆N 0 u 0 . Evolution couples this vector to the zero-energy mode according to Eq.(17), i.e. it obeys the following equation: Eq. (25) can be satisfied if u 0 = (u 0 , v 0 ) T meets the condition: where M is a parameter to be determined from the normalization condition, u ad 0 |v 0 = 1. Because according to our convention the ground state wavefunction ψ is real, we shall look for solutions of Eq.(26) assuming that u 0 (r) and v 0 (r) are real too. Adding the two equations of the set Eq. (26) one to the other, we get u 0 = v 0 and using this condition the following equation can be obtained when the two equations are substracted: The solution of this equation can be found by observing that ψ obeys the EGP equation, H GP ψ = 0. This equation differentiated with respect to the number of particles gives: Eq. (28) was obtained under the very general assumption that the EGP Hamiltonian is a function of density n = ψ 2 only, i.e. H GP = H GP (n). Comparing Eq. (27) and Eq.(28) we get: We can check that the zero-energy eigenvector v 0 has unit overlap with the dynamically coupled vector u 0 : 1 = u ad 0 |v 0 = dr(2ψ∂ N ψ), and the value of M is: Eqs. (25) and (26) give a linear time dependence of the amplitude of the zero-energy mode: ϑ 0 (t) = ϑ 0 + ∆N0 M t, as well as the perturbation field δ(r, t) The zero-energy mode contributes to the excitation field δψ(r, t) on equal footing with the Bogoliubov modes. Expansion of δψ(r, t) into normal modes eg. (18) must account for the zero-energy excitation: where the vector: is the eigenvector of the BdG operator L corresponding to the eigenvalue −ε Bν . The term proportional to w s Bν in Eq.(32) is necessary to meet the assumed form of the Bogoliubov expansion, Eq. (18).
The amplitudes of the normal vectors in the expansion Eq.(32) can be found by projecting the vector of perturbations δψ onto the corresponding adjoint vectors, i.e. b ν = w ad Bν |δψ , −iϑ 0 = u ad 0 |δψ and ∆N 0 = v ad 0 |δψ . The energy of the perturbed ground state of the system is equal to: H s [ψ + δ, ψ * + δ * ] where H s is given by Eq. (10). Expanding this equation up to the second order in a small perturbation δ we get: The linear term vanishes because ψ satisfies the EGP equation. The second order term is the Bogoliubov Hamiltonian: It is a sum of two terms: H B = ν ε Bν |b ν | 2 +H θ , where the first term is the energy of standard Bogoliubov excitations while the second one corresponds to the energy of the excitation of the mode coupled to the phase mode v 0 . After some algebra this term can be brought to the form: The zero-energy phase mode Hamiltonian H θ does not depend on the phase of the wavefunction, it depends on the amplitude of the coupled mode, ∆N 0 . Deviation of particle number from the equilibrium value, N = drψ 2 , is defined as: Keeping only terms linear in δ and δ * we have: Thus the amplitude ∆N 0 of the coupled mode describes the difference of the actual number of particles in a droplet from its equilibrium value. The symbol we used to denote the amplitude of the coupled mode was not accidental.
The time dependence of the phase, Eq.(31) and the form of the Hamiltonian H θ imply that (ϑ 0 , ∆N 0 ) are canonically conjugate variables describing one degree of freedom -in a full analogy to position and kinetic momentum of a particle. In the frame of the mean field approach the amplitudes (b ν , b * ν ), and (ϑ 0 , ∆N 0 ) are classical variables. We combined them in pairs to stress that two real functions are needed to describe every excitation mode.
The simplest way of switching to the quantum description of the excitations leads via the standard procedure of canonical quantization in analogy to [44]. To this end we shall impose bosonic commutation relations on the mean field excitations: [δ(r, 0)], δ(r ′ , 0)] = 0 The above commutation relations will be automatically satisfied if classical amplitudes of Bogoliubov excitations are replaced by annihilation and creation operators, b ν →b ν and b * ν →b † ν satisfying standard bosonic commutation relations [b ν ,b † µ ] = δ ν,µ . Moreover, amplitudes of the zero-energy modes have to be replaced by operators, ϑ 0 →θ and ∆N 0 → δN whose commutator is: [θ, δN ] = i. One can use the following representation of the two operators: where ϑ 0 is the phase of a condensate wavefunction, as can be seen from Eqs. (16), (21) and (22).
The quantum HamiltonianĤ θ = δN 2 /2M has no stationary solution even if bound from below except when fluctuations of particle number are set to zero. Instead, it predicts the standard quantum diffusion of droplet phase in full analogy to the quantum diffusion of a free particle wavepacket. If initially the wavepacket of a given number of atoms is prepared then the squared-variance of the phase of the droplet wavefunction will grow linearly in time: The same kind of diffusion was predicted in [44] for the phase of a Bose-Einstein condensate and in [35] where diffusion of phase and position of a dark soliton were studied. However, in a disordered potential the position of the soliton may be localized [36].

A. Zero-energy phase modes -general approach
A two component droplet in a Bose-Bose mixture is formed by N 1 atoms of the first kind and N 2 atoms of the second kind. The system must be described by two wavefunctions, ψ 1 and ψ 2 which we assume to be normalized to the number of particles in the corresponding component. The amount of atoms in every component must be carefully chosen to get the stationary stable droplet [1,41], i.e. stationary solutions of Eqs.(4).
In order to set the notation we have to repeat some elements of the previous discussion. The excitation spectrum of a two component droplet can be found by studying the response of the system to small perturbations of the two stationary wavefunctions ψ i : ψ i → e −iεit/ (ψ i (r) + δ i (r, t)): We now introduce the following two-component vectors, δψ i (r, 0) = (δ i (r, 0), δ * i (r, 0)) T i = 1, 2 which can be combined into a single four-component vector of perturbations (the four-component and the two-component vectors are denoted by the 'tilde' or 'line' above the symbol respectively): The two upper components, i.e. the field δψ 1 = (δ 1 , δ * 1 ) T are perturbations of the first wavefunction, ψ 1 , and its complex conjugate, ψ * 1 , while the two lower components give the perturbation field of ψ 2 , and ψ * 2 , δψ 2 = (δ 2 , δ * 2 ) T . In the standard Bogoliubov approach generalized to a two component mixture [15,16] the eigenvalue problem of the BdG operator L has to be solved to get eigenvectors and eigenenergies. L is the '4 × 4' matrix in this case: where the diagonal blocks are defined as: and the off-diagonal blocks are: In the equations above the H The two functions ψ i are the droplet's solutions of the set of the two coupled EGP equations: Eigenstates of the L operator of nonzero energies ε Bν = 0, are the Bogoliubov modes -the excitations build up on the droplet's wavefunctions. Again, similarly to the one-component case, the BdG operator L is non-hermitian. Its right eigenvectors W Bν = (w Bν ) T have adjoint partners corresponding to the same eigenvalue. These are the left eigenvectors, W ad Bν = (w Above we combined the two-component vectors w Bν ) to construct the two four-component vectors W Bν = (u (1) . The form of the adjoint eigenvectors follows form symmetries of L. The standard scalar product is therefore defined as: The ground state solutions ψ i of equations Eqs.(49) corresponding to the two-components of the self-bound droplet break three continuous symmetries of the many-body Hamiltonian: translation of the droplet as a whole and phase shifts of the two droplet's wavefunctions. The problem we discussed in the previous section reappears. In the excitation spectrum one has to account for three zero-energy modes recovering the broken symmetries of the many-body Hamiltonian. They cannot be obtained by diagonalization of the nonhermitian BdG operator, L.
To simplify the analysis, we take from now on ψ i as real functions. The same arguments which led to Eq.(22) allows to find the form of the two, i = 1, 2, zero-energy modes of the BdG operator, L W (i) 0 = 0, which correspond to a shift of phases of the droplet's wavefunctions: and In the equations above the two-components vectors are: w is also a zero-energy eigenvector: These are two independent vectors of zero-modes Eq.(53) corresponding to the two broken U(1) symmetries. They have zero norm, W The vectors U (i) 0 are not eigenvectors of L. They can be found using a physical argument based on analysis of the dynamics generated by the BdG operator. If initially δψ(r, 0) = ∆N , then the amplitude of zero-energy phase modes will vary in time according to: It is a simple exercise to check by inspection of analogous dynamical equations that the amplitudes ∆N Then direct comparison of both sides of Eq.(55) gives: The parameters M j play the role of masses in the dynamical equations describing the evolution of ϑ Upon using ψ i = ψ * i the set of equations Eqs.(56) simplifies. By adding the first two equations and repeating the same procedure for the third and fourth equation we get: i.e U (i) 0 has only two independent components. The set of four coupled equations Eq.(56) reduces therefore to only two equations: where the matrix K has the form: To solve the above equations let us observe that by differentiating the EGP equations, Eq.(49), with respect to N i : we get the following two sets (i = 1, 2) of the two coupled The two sets of Eqs.(63) have the form of Eqs.(60). This way, simply by comparison, we can write down the two independent solutions we are looking for: u In general any linear combination: is also a solution of Eq.(63): provided that: 2 ψ 2 ) T and the vectors U (i) 0 must have unit overlap. If in addition we account for the constraints following from the orthogonality conditions, Eq.(54), we get: Finally, there are eight conditions: four involving masses Eq.(66), and four orthonormality conditions Eq.(67) which have to be solved for ten unknowns: eight coefficients α j and two masses M i . The solutions for α-s are: 1 , α and similarly, solutions for β-s have the form: 1 , β Finally the 'inverse of masses' M i of the two phase modes are: In the equations above c (i) are arbitrary scaling factors. Their presence reflects the freedom resulting from a smaller number of equations than number of unknowns. The scaling factors play no role as long as dynamics of zeroenergy modes are considered because they do not enter the Hamiltonians generating dynamics -scaling of masses is compensated by the appropriate scaling of kinetic momenta. However, dimensional analysis suggests that it would be more elegant if the coefficients α-s and β-s were dimensionless. Therefore in the following we set the values of the arbitrary coefficients to c (i) = (∂ε 1 /∂N 1 ) −1/4 (∂ε 2 /∂N 2 ) −1/4 . The amplitudes of zero-energy modes, −iϑ Utilizing orthogonality conditions Eq.(54), dr 2β we see that the amplitudes of the coupled modes are physical observables, namely deviations of particle number from the equilibrium value: Amplitudes ∆N (i) 0 of modes U (i) are linear combinations of these deviations of particle numbers in the first and second component of the classical fields ψ i . Analogous relations can be found for amplitudes ϑ (i) 0 of the zero-energy modes. To show this we remind that a small shift of the phases of droplets' wavefunctions ψ i (r, 0)e −iθi ≃ ψ i (r, 0) (1 − iθ i ) leads to the following perturbation field δ i (r, 0) = −iθ i ψ i . Combining this result with the assumed form of the perturbation field Eq.(70): Now using the orthogonality conditions Eq.(67) we find The amplitudes ϑ where: Using Eqs.(68) and (69) the explicit form of the Hamiltonians Eq.(76) is found: Eqs. (77) and (78)  Finally, the classical Hamiltonians above can be easily quantized by imposing canonical commutation relations on canonical variables. Evidently the pairs (θ 1 , δN 1 ) and (θ 2 , δN 2 ) are canonically conjugate variables. One can see here a close analogy to canonical position and momentum variables (x, p x ). Following the canonical quantization procedure we substitute the phases and atom number fluctuation variables by operators, θ i →θ i and δN i → δN i and impose canonical commutation relations: A convenient choice of representation of canonical operators is δN i = −i ∂ ∂θi .
Consistently we can easily verify that the phases of the soft and hard modesφ s =θ Quantum dynamics will lead to a phase diffusion of the wavepacket being a superposition of states with various atom numbers. The same effect takes place in the simplified one-component case studied in the previous subsection. Note that the mass of the soft mode is much larger than the mass of the hard mode. Therefore the phase of the soft mode will diffuse slowly in comparison with the phase of the hard mode.

B. Phase modes -Bose-Bose droplets
The considerations presented above are quite general. They can be applied to any system with broken U1 and translational symmetries. The results obtained were based on the assumption that interaction energy, E(n 1 , n 2 ), is a function of densities of both species involved, n i = |ψ i | 2 . Here we want to address the experimentally important situation of a mixture of two interacting Bose-Einstein condensates forming a self-bound droplet. In the case of such a system we can give explicit formulae for the Hamiltonians of the zero-energy modes.
The stationary state of the droplet is described by two wavefunctions, ψ i , which are eigenstates of the coupled EGP equations, Eq.(4). Eigenvalues, ε i (ψ 1 , ψ 2 ), account for contributions from kinetic and interaction energies. The stationary densities, n i = |ψ i | 2 , are almost constant in the bulk of a droplet and fall to zero at its surface, [1]. The bulk contribution dominates over the one from the surface for sufficiently large droplets. With decreasing number of atoms the bulk diminishes and eventually disappears as the droplet becomes very small -only a surface remains.
We now focus on the large droplet case, where we neglect the impact of surface. We model the large droplet by uniform densities n 1 , n 2 of volume V . Then its total energy is equal to E(n 1 , n 2 )V where E(n 1 , n 2 ) is given by Eq. (2). Here N 1 = n 1 V and N 2 = n 2 V are the numbers of particles in both components. At equilibrium we have ∂E ∂V N1,N2 = 0 which results in V (N 1 , N 2 ). With such a model it is questionable to use the solution of Eq. (60) given by Eq. (63). There to get the masses M i the EGP equations were differentiated, Eq.(61) with respect to N i under the silent assumption that the ground state solutions ψ i are differentiable. This is a case when kinetic energy is included and the functions change smoothly with atom number as well as in space, eventually falling exponentially to zero at infinity. When the kinetic energy term is neglected in the EGP Hamiltonian, H (i) GP , the stationary solutions have the shape of a step-like function with a step height and position depending on atom number. Its derivative with respect to atom number must have a contribution which is a Dirac-delta distribution at the step position. Dealing with such functions is a mathematically subtle task. Therefore, instead of using Eq. (63) we analytically solve Eqs. (60), after neglecting the kinetic energy term in the operator K. We arrive at a set of linear algebraic equations which together with orthogonality conditions, Eq. (54), lead to the Hamiltonians In the model studied here the derivatives ∂µ j /∂n i can be calculated directly and substituted into the Hamiltonians H 1 and H 2 : The contribution of the LHY term to the Hamiltonian H 1 is essential. H 1 is proportional to the small parameter δa/ √ a 11 a 22 ≪ 1. It describes the zero-energy excitations of the soft mode. This is not the case of the H 2 Hamiltonian. Here terms of the order of δa/ √ a 11 a 22 are negligible as compared to other contributions. We have already mentioned that there are some indications that the masses entering the zero-energy soft and hard mode Hamiltonians are very different. Explicit formulae confirm this statement. The mass related to the soft mode is inversely proportional to the small parameter δa/ √ a 11 a 22 , 1/M 1 ∼ δa/ √ a 11 a 22 , and is much larger than the mass of the hard mode M 1 ≫ M 2 . The dynamics of the soft mode phase is very slow, much slower than the evolution of the hard mode phase.

C. Translational modes
In addition to the two U(1) symmetries related to the freedom of choice of phases of the wavefunctions there is yet another continuous symmetry which is broken by the mean field solution. This is the translational symmetry of the entire droplet. Droplet position in space can be arbitrary. It does not affect the energy. We consider 3D droplets, therefore the translational zero-energy mode has a vectorial character. In other words there are three modes which correspond to translations along the three principal axes of a Cartesian coordinate system. Because space is isotropic we consider here translation along the z-axis only. The modes corresponding to translations in the two remaining principal directions can be obtained analogously.
The translational zero-energy modes can be found by applying a translation operator (along z-axis) to the stationary wavefunctions. In particular for an infinitesimally small shift of the droplet's position we have: ψ tr i (r + z 0 e z ) = e z0∂z ψ i (r) = ψ i (r) + z 0 ∂ z ψ i (r), which indicates that the zero-energy mode has the form V tr = (∂ z ψ 1 , ∂ z ψ * 1 , ∂ z ψ 2 , ∂ z ψ * 2 ) T .
The vector U tr can be defined as the one which couples dynamically to the zero-energy vector:

V. SUMMARY
In this paper we studied quantum droplets -self bound dilute liquid systems of ultracold atoms. The description of droplets in the frame of the mean field theory breaks three continuous symmetries. The energy of a droplet does not depend on the phases of the wavefunctions of the involved components. It also does not depend on the position of the droplet in uniform space. The excitation spectrum of a droplet should account for these zero-energy excitations. They recover broken symmetries of the many body Hamiltonian. The standard Bogoliubov-de Gennes approach does not account for these zero-energy modes.
Here we found the explicit form of the zero-energy modes as well as classical and quantum expressions of the Hamiltonians generating their time evolution. For every variable quantifying the particular symmetry operationthe 'canonical positions', i.e. the magnitudes of the phase shifts θ i and spatial translation r 0 there exists a canonical conjugated momentum, ∆N i and P. They correspond to deviation of particle numbers from their equilibrium values in the case of phase modes and to the kinetic momentum of the droplet in the case of the translational mode. Every pair of canonically conjugate variables describes a single degree of freedom.
The formalism developed here shall be very useful in the description of dynamical processes involving several droplets in situations when their relative phases become important. The analysis presented here is limited to situations when the perturbations of the stationary state are small. As we see the dynamical equations derived predict linear growth in time of the initial phases and droplet's position. This does not indicate that the approach is valid only for very short times. The zero-energy modes depend explicitly on the instantaneous phases or position of the droplet. Therefore one should bear in mind that we consider small perturbations of the instantaneous phases and/or position of a droplet.