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

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–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–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.

In the recent theory works some improvement over the standard energy functional were suggested. The problem with the LHY energy for a two-component Bose–Bose mixture in the region of droplet formation is that it contains a small imaginary part which is commonly neglected. In [11, 12] this problem is cured by accounting for a Bose–Bose pairing energy. Contribution beyond the LHY energy are introduced in [13].

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 [14–16], but also in optical lattices [17]. Energy of quantum fluctuations depends on the density of Bogoliubov modes [18, 19] which, in turn, is a function of the dimensionality of the phase space involved. Bose–Bose droplets in lower dimensions [20] and at dimensional crossover [21, 22] were investigated theoretically. Similarly, dipolar droplets at the crossover were studied [23]. Quantum fluctuations in low-dimensional dipolar systems or systems at a dimensional crossover exhibit many interesting features [24]. In a strongly interacting regime novel quantum droplets are predicted [25]. 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 [26]. 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 [27]. Binary collisions in the context of quantum liquid droplets were investigated theoretically in 1D geometry [28]. Collisions of droplets composed of two hyperfine states of ³⁹K were studied in the experiment of [29].

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 [30, 31]. 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 [32–34]. To this end external trap geometry and interactions have to be appropriately tuned to trigger the so-called roton instability [35] appearing when a roton minimum in the excitation spectrum touches zero energy. Coherent arrays of droplets linked by a superfluid 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 [36] 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 [37] are predicted theoretically. For the review of the recent achievements in the field of ultradilute quantum droplets see [38].

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 (BdG) 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 [39]. The formalism developed for zero-energy solitonic modes occurred to be very useful in the description of greying of a dark soliton [40], Anderson localization of a soliton or localization of two interacting solitons [41, 42].

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. In the second part of the introduction, section 1.2, we give a compact physical explanation of the main results of our paper. Next, in section 2, we briefly introduce the main formalism of description of two-component Bose–Bose droplets. In section 3, 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 description can be also applied to dipolar droplets. A two component droplet is considered in section 4.1 where we find general expressions for the two phase modes as well as for the translational mode, section 4.3. 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, section 4.2. In the last section 5, we summarize our results and conclude.

In this section we give a simple physical summary of our main results. We compromise between detailedness and compactness of our discussion. For a precise description, including some notation, the reader must check the main part of the paper.

The goal of this paper is to determine the excitation modes of the system which are related to broken U(1) symmetries—the so called zero modes. We consider the case of a two-component droplet whose ground state is described by two wavefunctions, ψ_i , i = 1, 2. Standard procedure to find the excitation spectrum is to add a small perturbation ${\delta }_{i}\left(\mathbf{r},t\right)$ to the stationary state:

$$\begin{equation}{\psi }_{i}\to {\text{e}}^{-\text{i}{\varepsilon }_{i}t/\hslash }\left({\psi }_{i}\left(\mathbf{r}\right)+{\delta }_{i}\left(\mathbf{r},t\right)\right).\end{equation} \tag{ 1 }$$

and to analyze a response of the system by monitoring a time evolution. The normalization condition reads $\int \mathrm{d}\mathbf{r}\vert {\psi }_{i}\left(\mathbf{r}\right){\vert }^{2}={N}_{i}^{\left(0\right)}$. Knowledge of the form of elementary excitations field δ_i (r, t) allows for determination of the excitation energy defined as the difference:

$$\begin{equation}{H}_{\text{ex}}\left({\delta }_{i}\right)={H}_{\mathrm{d}}\left[{\psi }_{i}+{\delta }_{i},{\psi }_{i}^{{\ast}}+{\delta }_{i}^{{\ast}}\right]-{H}_{\mathrm{d}}\left[{\psi }_{i},{\psi }_{i}^{{\ast}},\right],\end{equation} \tag{ 2 }$$

where the energy functional H_d is expanded up to the second order terms in δ_i . The form of H_d is given in section 2.

Standard approach to excitation spectrum of an ultracold system is to diagonalize the BdG equations. It gives Bogoliubov quasiparticles. But these are not all possible excitations. The two functions ψ₁ and ψ₂ describing a droplet break five continuous symmetries of the Hamiltonian (2). Three symmetries are related to translation of the whole droplet along any 3D vector. Evidently in the uniform space the energy of the droplet does not depend on its location. This is not a case of a trapped system. Droplet's solution break also two U(1) symmetries. Each wavefunction can be multiplied by an arbitrary phase factor ${\psi }_{i}\to {\text{e}}^{-\text{i}{\theta }_{i}}{\psi }_{i}$ without affecting energy of the system.

In section 4.3 we show that the translation of the droplet along z-axis (translation along any other principal axis can be treated in the same fashion) leads to the following modification of droplet's wavefunctions:

$$\begin{equation}{\delta }_{i}^{\text{tr}}\left(\mathbf{r},t\right)=z\left(t\right){\partial }_{z}{\psi }_{i}-i\frac{{p}_{z}}{\hslash }\frac{m}{{M}_{\text{tr}}}z{\psi }_{i},\end{equation} \tag{ 3 }$$

where z(t) = z₀ + t(p_z /M_tr) and ip_z m/(ℏM_tr) are amplitudes of the zero mode, ∂_z ψ_i , and the adjoint mode, zψ_i , while ${M}_{\text{tr}}=m\left({N}_{1}^{\left(0\right)}+{N}_{2}^{\left(0\right)}\right)$ is droplet mass. Amplitude z(t) is the position of the center of mass of the droplet and p_z is its momentum.

Energy related to the translational mode excitation ${H}_{\text{tr}}={H}_{\text{tr}}\left({\delta }_{i}^{\text{tr}}\right)$ is equal to the energy of a free motion of the droplet with velocity p_z /M_tr:

$$\begin{equation}{H}_{\text{tr}}=\frac{{p}_{z}^{2}}{2{M}_{\text{tr}}}.\end{equation} \tag{ 4 }$$

This is rather obvious result but it visualizes the main lines of the approach used.

In addition to the translational modes there are two zero-energy phase modes related to two U(1) symmetries broken by a particular choice of phases of the ground state wavefunctions. Excitations corresponding to the phase shifts are missing from the Bogoliubov spectrum. Corresponding zero energy eigenvectors are obtained by applying the symmetry operation to the ground state wavefunctions, ${\psi }_{i}\left(\mathbf{r},0\right)\to {\psi }_{i}\left(\mathbf{r},0\right){\text{e}}^{-\text{i}{\theta }_{i}}\simeq {\psi }_{i}\left(\mathbf{r},0\right)-i{\theta }_{i}{\psi }_{i}\left(\mathbf{r},0\right)$. In section 4.1 we show that perturbation field of the droplet's wavefunctions caused by shift of phases has the form:

$$\begin{equation}{\delta }_{i}^{\theta }\left(\mathbf{r},t\right)=\sum\limits _{j=1}^{2}\left(-i{\vartheta }_{0}^{\left(j\right)}\left(t\right){\beta }_{i}^{\left(j\right)}{\psi }_{i}+{\Delta}{\mathcal{N}}_{0}^{\left(j\right)}\left({\alpha }_{1}^{\left(j\right)}\frac{\partial {\psi }_{i}}{\partial {N}_{1}}+{\alpha }_{2}^{\left(j\right)}\frac{\partial {\psi }_{i}}{\partial {N}_{2}}\right)\right).\end{equation} \tag{ 5 }$$

Coefficients ${\alpha }_{i}^{\left(j\right)}$, ${\beta }_{i}^{\left(j\right)}$ and M_i are to be chosen to assure orthonormality of different modes involved. Phases ${\vartheta }_{0}^{\left(i\right)}$ are a linear combination of phase shifts of the two components, θ_i :

$$\begin{equation}{\vartheta }_{0}^{\left(i\right)}={\alpha }_{1}^{\left(i\right)}{\theta }_{1}+{\alpha }_{2}^{\left(i\right)}{\theta }_{2},\end{equation} \tag{ 6 }$$

and they evolve in time according to:

$$\begin{equation}{\vartheta }_{i}\left(t\right)={\vartheta }_{i}\left(0\right)+t\frac{{\Delta}{\mathcal{N}}_{0}^{\left(i\right)}}{\hslash {M}_{i}}.\end{equation} \tag{ 7 }$$

Equation (7) describes time dependence of the position ${\vartheta }_{0}^{\left(i\right)}\left(t\right)$ of a fictitious particle of a mass M_i moving with a constant velocity ${\Delta}{\mathcal{N}}_{0}^{\left(i\right)}/\left(\hslash {M}_{i}\right)$. This velocity is proportional to ${\Delta}{\mathcal{N}}_{0}^{\left(i\right)}$:

$$\begin{equation}{\Delta}{\mathcal{N}}_{0}^{\left(i\right)}={\beta }_{1}^{\left(i\right)}\delta {N}_{1}+{\beta }_{2}^{\left(i\right)}\delta {N}_{2},\end{equation} \tag{ 8 }$$

where δN_i is a difference between actual number of atoms occupying the ith component, N_i and the number of atoms corresponding to the stationary ground state solution: $\delta {N}_{i}={N}_{i}-{N}_{i}^{\left(0\right)}$. The pair of variables (θ_i , δN_i ) describes the phase shift and deviation of a particle number of the ith component of the droplet. They form a pair of conjugated variables—analogons of position and momentum.

Excitation energy ${H}_{0}={H}_{\text{ex}}\left({\delta }_{i}^{\theta }\right)$ corresponding to the perturbation of the ground state wavefunctions by ${\delta }_{i}^{\theta }\left(\mathbf{r},t\right)$ is the following:

$$\begin{equation}{H}_{0}={H}_{1}+{H}_{2},\end{equation} \tag{ 9 }$$

where:

$$\begin{equation}{H}_{i}\left(\delta {\vartheta }_{0}^{\left(i\right)},{\Delta}{\mathcal{N}}_{0}^{\left(i\right)}\right)=\frac{{\left({\Delta}{\mathcal{N}}_{0}^{\left(i\right)}\right)}^{2}}{2{M}_{i}}=\frac{{\left({\beta }_{1}^{\left(i\right)}\delta {N}_{1}+{\beta }_{2}^{\left(i\right)}\delta {N}_{2}\right)}^{2}}{2{M}_{i}}.\end{equation} \tag{ 10 }$$

The Hamiltonians H₁ and H₂ have the form of Hamiltonians of free particles. Excitations described by H₁ have energies much smaller than those described by H₂. Thus variables $\left({\vartheta }_{0}^{\left(1\right)},{\Delta}{\mathcal{N}}_{0}^{\left(1\right)}\right)$ describe the soft mode while $\left({\vartheta }_{0}^{\left(2\right)},{\Delta}{\mathcal{N}}_{0}^{\left(2\right)}\right)$, the hard mode excitations. Quantum description of the hard and soft modes can be obtained by imposing canonical commutation relations: $\left[{\hat{\vartheta }}_{0}^{\left(i\right)},{\Delta}{\hat{\mathcal{N}}}_{0}^{\left(j\right)}\right]=i{\delta }_{ij}$. Convenient representation of the operators is ${\hat{\vartheta }}_{0}^{\left(i\right)}={\vartheta }_{0}^{\left(i\right)}$ and ${\Delta}{\hat{\mathcal{N}}}_{0}^{\left(i\right)}=i\frac{\partial }{\partial {\vartheta }_{0}^{\left(i\right)}}$. Analogously one can quantize variables describing different components of the system: $\left[{\hat{\theta }}_{i},\delta \hat{{\mathcal{N}}_{j}}\right]=i{\delta }_{ij}$.

Identification of canonical variables and excitation Hamiltonians, in particular expressions for masses M_i are the main results of our paper.

A two component droplet in a Bose–Bose mixture is formed by N₁ atoms of the first kind and N₂ atoms of the second kind. The system must be described by two wavefunctions, ψ₁ and ψ₂ 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, 46], i.e. stationary solutions of equation (14).

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 : ${\psi }_{i}\to {\text{e}}^{-\text{i}{\varepsilon }_{i}t/\hslash }\left({\psi }_{i}\left(\mathbf{r}\right)+{\delta }_{i}\left(\mathbf{r},t\right)\right)$:

$$\begin{equation}{\delta }_{i}\left(\mathbf{r},t\right)=\sum\limits _{\nu }\left({b}_{\nu }^{\left(i\right)}{u}_{B\nu }^{\left(i\right)}\enspace {\text{e}}^{-\text{i}{\varepsilon }_{B\nu }t/\hslash }+{b}_{\nu }^{\left(i\right){\ast}}{v}_{B\nu }^{\left(i\right){\ast}}\enspace {\text{e}}^{\text{i}{\varepsilon }_{B\nu }t/\hslash }\right).\end{equation} \tag{ 53 }$$

We now introduce the following two-component vectors, ${\bar{\delta \psi }}_{i}\left(\mathbf{r},0\right)={\left({\delta }_{i}\left(\mathbf{r},0\right),{\delta }_{i}^{{\ast}}\left(\mathbf{r},0\right)\right)}^{\text{T}}\quad 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):

$$\begin{equation}\tilde {\delta \psi }\left(\mathbf{r},0\right)=\left(\begin{matrix}\hfill {\bar{\delta \psi }}_{1}\left(\mathbf{r},0\right)\hfill \\ \hfill {\bar{\delta \psi }}_{2}\left(\mathbf{r},0\right).\hfill \end{matrix}\right)\end{equation} \tag{ 54 }$$

The two upper components, i.e. the field ${\bar{\delta \psi }}_{1}={\left({\delta }_{1},{\delta }_{1}^{{\ast}}\right)}^{\text{T}}$ are perturbations of the first wavefunction, ψ₁, and its complex conjugate, ${\psi }_{1}^{{\ast}}$, while the two lower components give the perturbation field of ψ₂, and ${\psi }_{2}^{{\ast}}$, ${\bar{\delta \psi }}_{2}={\left({\delta }_{2},{\delta }_{2}^{{\ast}}\right)}^{\text{T}}$. In the standard Bogoliubov approach generalized to a two component mixture [18, 19] the eigenvalue problem of the BdG operator $\mathcal{L}$ has to be solved to get eigenvectors and eigenenergies. $\mathcal{L}$ is the '4 × 4' matrix in this case (we do not introduce a separate notation for 2 × 2 BdG operator from the previous section):

$$\begin{equation}\mathcal{L}=\left(\begin{matrix}\hfill {\mathcal{L}}_{11}\hfill & \hfill {\mathcal{L}}_{12}\hfill \\ \hfill {\mathcal{L}}_{21}\hfill & \hfill {\mathcal{L}}_{22}\hfill \end{matrix}\right),\end{equation} \tag{ 55 }$$

where the diagonal blocks are defined as:

$$\begin{equation}{\mathcal{L}}_{ii}=\left(\begin{matrix}\hfill {H}_{\text{GP}}^{\left(i\right)}+{n}_{i}\frac{\partial {\mu }_{i}}{\partial {n}_{i}}\hfill & \hfill {\psi }_{i}^{2}\frac{\partial {\mu }_{i}}{\partial {n}_{i}}\hfill \\ \hfill -{\psi }_{i}^{{\ast}2}\frac{\partial {\mu }_{i}}{\partial {n}_{i}}\hfill & \hfill -{H}_{\text{GP}}^{\left(i\right)}-{n}_{i}\frac{\partial {\mu }_{i}}{\partial {n}_{i}}\hfill \end{matrix}\right)\end{equation} \tag{ 56 }$$

and the off-diagonal blocks are:

$$\begin{equation}{\mathcal{L}}_{ij}=\left({\partial }_{{n}_{j}}{\mu }_{i}\right)\left(\begin{matrix}\hfill {\psi }_{i}{\psi }_{j}^{{\ast}}\hfill & \hfill {\psi }_{i}{\psi }_{j}\hfill \\ \hfill -{\psi }_{i}^{{\ast}}{\psi }_{j}^{{\ast}}\hfill & \hfill -{\psi }_{i}^{{\ast}}{\psi }_{j}\hfill \end{matrix}\right).\end{equation} \tag{ 57 }$$

In the equations above the ${H}_{\text{GP}}^{\left(i\right)}$ are EGP Hamiltonians:

$$\begin{equation}{H}_{\text{GP}}^{\left(i\right)}=-\frac{{\hslash }^{2}}{2m}{\Delta}+{\mu }_{i}-{\varepsilon }_{i}.\end{equation} \tag{ 58 }$$

The two functions ψ_i are the droplet's solutions of the set of the two coupled EGP equations:

$$\begin{equation}{H}_{\text{GP}}^{\left(i\right)}{\psi }_{i}=0.\end{equation} \tag{ 59 }$$

Eigenstates of the $\mathcal{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 $\mathcal{L}$ is non-Hermitian. Its right eigenvectors ${\tilde {W}}_{B\nu }={\left({\bar{w}}_{B\nu }^{\left(1\right)},{\bar{w}}_{B\nu }^{\left(2\right)}\right)}^{\text{T}}$ have adjoint partners corresponding to the same eigenvalue. These are the left eigenvectors, ${\tilde {W}}_{B\nu }^{\text{ad}}=\left({\bar{w}}_{B\nu }^{\left(1\right)\text{ad}},{\bar{w}}_{B\nu }^{\left(2\right)\text{ad}}\right)$:

$$\begin{equation}{\tilde {W}}_{B\nu }^{\text{ad}}\left(\mathcal{L}-{\varepsilon }_{B\nu }\right)=\left(\mathcal{L}-{\varepsilon }_{B\nu }\right){\tilde {W}}_{B\nu }.\end{equation} \tag{ 60 }$$

Above we combined the two-component vectors ${\bar{w}}_{B\nu }^{\left(i\right)}={\left({u}_{B\nu }^{\left(i\right)},{v}_{B\nu }^{\left(i\right)}\right)}^{\text{T}}$ and similarly ${\bar{w}}_{B\nu }^{\left(i\right)\text{ad}}=\left({u}_{B\nu }^{\left(i\right)\text{ad}},{v}_{B\nu }^{\left(i\right)\text{ad}}\right)$ to construct the two four-component vectors ${\tilde {W}}_{B\nu }={\left({u}_{B\nu }^{\left(1\right)},{v}_{B\nu }^{\left(1\right)},{u}_{B\nu }^{\left(2\right)},{v}_{B\nu }^{\left(2\right)}\right)}^{\text{T}}$ and ${\tilde {W}}_{B\nu }^{\text{ad}}=\left({u}_{B\nu }^{\left(1\right)\text{ad}},{v}_{B\nu }^{\left(1\right)\text{ad}},{u}_{B\nu }^{\left(2\right)\text{ad}},{v}_{B\nu }^{\left(2\right)\text{ad}}\right)=\left({u}_{B\nu }^{\left(1\right){\ast}},-{v}_{B\nu }^{\left(1\right){\ast}},{u}_{B\nu }^{\left(2\right){\ast}},-{v}_{B\nu }^{\left(2\right){\ast}}\right)$. The form of the adjoint eigenvectors follows form symmetries of $\mathcal{L}$. The standard scalar product is therefore defined as: $\langle {\tilde {W}}_{B\nu }^{\text{ad}}\vert {\tilde {W}}_{B\nu }\rangle =\int \mathrm{d}\mathbf{r}\enspace {\tilde {W}}_{B\nu }^{\text{ad}}\left(\mathbf{r}\right){\tilde {W}}_{B\nu }\left(\mathbf{r}\right)={\sum }_{i=1,2}\langle {\bar{w}}_{B\nu }^{\left(i\right)\text{ad}}\vert {\bar{w}}_{B\nu }^{\left(i\right)}\rangle ={\sum }_{i=1,2}\int \mathrm{d}\mathbf{r}\left(\vert {u}_{B\nu }^{\left(i\right)}{\vert }^{2}-\vert {v}_{B\nu }^{\left(i\right)}{\vert }^{2}\right)$.

The ground state solutions ψ_i of equations (59) 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 non-Hermitian BdG operator, $\mathcal{L}$.

To simplify the analysis, we take from now on ψ_i as real functions. The same arguments which led to equation (32) allows to find the form of the two, i = 1, 2, zero-energy modes of the BdG operator, $\mathcal{L}{\tilde {W}}_{0}^{\left(i\right)}=0$, which correspond to a shift of phases of the droplet's wavefunctions:

$$\begin{equation}{\tilde {W}}_{0}^{\left(1\right)}=\left(\begin{matrix}\hfill {\bar{w}}_{0}^{\left(1\right)}\hfill \\ \hfill \bar{0}\hfill \end{matrix}\right),\end{equation} \tag{ 61 }$$

and

$$\begin{equation}{\tilde {W}}_{0}^{\left(2\right)}=\left(\begin{matrix}\hfill \bar{0}\hfill \\ \hfill {\bar{w}}_{0}^{\left(2\right)}\hfill \end{matrix}\right).\end{equation} \tag{ 62 }$$

In the equations above the two-components vectors are: ${\bar{w}}_{0}^{\left(1\right)}={\left({\psi }_{1},-{\psi }_{1}\right)}^{\text{T}}$, ${\bar{w}}_{0}^{\left(2\right)}={\left({\psi }_{2},-{\psi }_{2}\right)}^{\text{T}}$, and $\bar{0}={\left(0,0\right)}^{\text{T}}$. Any linear combination of vectors ${\tilde {W}}_{0}^{\left(1\right)}$ and ${\tilde {W}}_{0}^{\left(2\right)}$ is also a zero-energy eigenvector:

$$\begin{equation}{\tilde {V}}_{0}^{\left(i\right)}={\beta }_{1}^{\left(i\right)}{\tilde {W}}_{0}^{\left(1\right)}+{\beta }_{2}^{\left(i\right)}{\tilde {W}}_{0}^{\left(2\right)}.\end{equation} \tag{ 63 }$$

These are two independent vectors of zero-modes equation (63) corresponding to the two broken U(1) symmetries. They have zero norm, $\langle {\tilde {W}}_{0}^{\left(i\right)\text{ad}}\vert {\tilde {W}}_{0}^{\left(j\right)}\rangle =0$. Our goal is to find two vectors ${\tilde {U}}_{0}^{\left(i\right)}={\left({u}_{01}^{\left(i\right)},{v}_{01}^{\left(i\right)},{u}_{02}^{\left(i\right)},{v}_{02}^{\left(i\right)}\right)}^{\text{T}}$ having unit overlap with the corresponding zero-energy phase modes, ${\tilde {V}}_{0}^{\left(i\right)}$:

$$\begin{equation}\langle {\tilde {U}}_{0}^{\left(i\right)\text{ad}}\vert {\tilde {V}}_{0}^{\left(j\right)}\rangle ={\delta }_{i,j}.\end{equation} \tag{ 64 }$$

The vectors ${\tilde {U}}_{0}^{\left(i\right)}$ are not eigenvectors of $\mathcal{L}$. They can be found using a physical argument based on analysis of the dynamics generated by the BdG operator. If initially $\tilde {\delta \psi }\left(\mathbf{r},0\right)={\Delta}{\mathcal{N}}_{0}^{\left(1\right)}{\tilde {U}}_{0}^{\left(1\right)}+{\Delta}{\mathcal{N}}_{0}^{\left(2\right)}{\tilde {U}}^{\left(2\right)}$, then the amplitude of zero-energy phase modes will vary in time according to:

$$\begin{equation}i\hslash \frac{\mathrm{d}}{\mathrm{d}t}\left(-i{\vartheta }_{0}^{\left(1\right)}\left(t\right){\tilde {V}}_{0}^{\left(1\right)}-i{\vartheta }_{0}^{\left(2\right)}\left(t\right){\tilde {V}}_{0}^{\left(2\right)}\right)=\mathcal{L}\left({\Delta}{\mathcal{N}}_{0}^{\left(1\right)}{\tilde {U}}_{0}^{\left(1\right)}+{\Delta}{\mathcal{N}}_{0}^{\left(2\right)}{\tilde {U}}_{0}^{\left(2\right)}\right).\end{equation} \tag{ 65 }$$

It is a simple exercise to check by inspection of analogous dynamical equations that the amplitudes ${\Delta}{\mathcal{N}}_{0}^{\left(i\right)}$ of the coupled vectors remain constant. Evidently (65) is satisfied if the vectors ${\tilde {U}}_{0}^{\left(i\right)}$ are solutions of the following equations:

$$\begin{equation}\mathcal{L}{\tilde {U}}_{0}^{\left(i\right)}=\frac{1}{{M}_{i}}{\tilde {V}}_{0}^{\left(i\right)}.\end{equation} \tag{ 66 }$$

Then direct comparison of both sides of equation (65) gives:

$$\begin{equation}{\vartheta }_{0}^{\left(j\right)}\left(t\right)={\vartheta }_{0}^{\left(j\right)}+t\frac{{\Delta}{\mathcal{N}}_{0}^{\left(j\right)}}{{M}_{j}\hslash }.\end{equation} \tag{ 67 }$$

The parameters M_j play the role of masses in the dynamical equations describing the evolution of ${\vartheta }_{0}^{\left(j\right)}$. Their values have to be chosen to ensure the normalization condition (64).

Upon using ${\psi }_{i}={\psi }_{i}^{{\ast}}$ the set of equations (66) simplifies. By adding the first two equations and repeating the same procedure for the third and fourth equation we get:

$$\begin{equation}{u}_{01}^{\left(i\right)}={v}_{01}^{\left(i\right)},\end{equation} \tag{ 68 }$$

$$\begin{equation}{u}_{02}^{\left(i\right)}={v}_{02}^{\left(i\right)},\end{equation} \tag{ 69 }$$

i.e ${U}_{0}^{\left(i\right)}$ has only two independent components. The set of four coupled equations (66) reduces therefore to only two equations:

$$\begin{equation}\mathcal{K}\left(\begin{matrix}\hfill {u}_{01}^{\left(i\right)}\hfill \\ \hfill {u}_{02}^{\left(i\right)}\hfill \end{matrix}\right)=\frac{1}{{M}_{i}}\left(\begin{matrix}\hfill {\beta }_{1}^{\left(i\right)}{\psi }_{1}\hfill \\ \hfill {\beta }_{2}^{\left(i\right)}{\psi }_{2},\hfill \end{matrix}\right)\end{equation} \tag{ 70 }$$

where the matrix $\mathcal{K}$ has the form:

$$\begin{equation}\mathcal{K}=\left(\begin{matrix}\hfill {H}_{\text{GP}}^{\left(1\right)}+2{n}_{1}\frac{\partial {\mu }_{1}}{\partial {n}_{1}}\hfill & \hfill 2{\psi }_{1}{\psi }_{2}\frac{\partial {\mu }_{1}}{\partial {n}_{2}}\hfill \\ \hfill 2{\psi }_{1}{\psi }_{2}\frac{\partial {\mu }_{2}}{\partial {n}_{1}}\hfill & \hfill {H}_{\text{GP}}^{\left(2\right)}+2{n}_{2}\frac{\partial {\mu }_{2}}{\partial {n}_{2}}.\hfill \end{matrix}\right)\end{equation} \tag{ 71 }$$

To solve the above equations let us observe that by differentiating the EGP equations (59) with respect to ${N}_{i}^{\left(0\right)}$:

$$\begin{equation}\frac{\partial }{\partial {N}_{i}^{\left(0\right)}}{H}_{\text{GP}}^{\left(j\right)}{\psi }_{j}=0,\end{equation} \tag{ 72 }$$

we get the following two sets (i = 1, 2) of the two coupled equations:

$$\begin{equation}\mathcal{K}\left(\begin{matrix}\hfill \frac{\partial }{\partial {N}_{i}^{\left(0\right)}}{\psi }_{1}\hfill \\ \hfill \frac{\partial }{\partial {N}_{i}^{\left(0\right)}}{\psi }_{2}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill \frac{\partial {\varepsilon }_{1}}{\partial {N}_{i}^{\left(0\right)}}{\psi }_{1}\hfill \\ \hfill \frac{\partial {\varepsilon }_{2}}{\partial {N}_{i}^{\left(0\right)}}{\psi }_{2}\hfill \end{matrix}\right).\end{equation} \tag{ 73 }$$

The two sets of equation (73) have the form of (70). This way, simply by comparison, we can write down the two independent solutions we are looking for: ${\bar{u}}_{0}^{\left(i\right)}\equiv \left({u}_{01}^{\left(i\right)},{u}_{02}^{\left(i\right)}\right)=\left(\frac{\partial {\psi }_{1}}{\partial {N}_{i}^{\left(0\right)}},\frac{\partial {\psi }_{2}}{\partial {N}_{i}^{\left(0\right)}}\right)$. In general any linear combination:

$$\begin{equation}{\bar{u}}^{\left(i\right)}\equiv {\alpha }_{1}^{\left(i\right)}{\bar{u}}_{0}^{\left(1\right)}+{\alpha }_{2}^{\left(i\right)}{\bar{u}}_{0}^{\left(2\right)},\end{equation} \tag{ 74 }$$

is also a solution of equation (73):

$$\begin{equation}\mathcal{K}{\bar{u}}^{\left(i\right)}={\alpha }_{1}^{\left(i\right)}\left(\begin{matrix}\hfill \frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}^{\left(0\right)}}{\psi }_{1}\hfill \\ \hfill \frac{\partial {\varepsilon }_{2}}{\partial {N}_{1}^{\left(0\right)}}{\psi }_{2}\hfill \end{matrix}\right)+{\alpha }_{2}^{\left(i\right)}\left(\begin{matrix}\hfill \frac{\partial {\varepsilon }_{1}}{\partial {N}_{2}^{\left(0\right)}}{\psi }_{1}\hfill \\ \hfill \frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}^{\left(0\right)}}{\psi }_{2}\hfill \end{matrix}\right)=\frac{1}{{M}_{i}}\left(\begin{matrix}\hfill {\beta }_{1}^{\left(i\right)}{\psi }_{1}\hfill \\ \hfill {\beta }_{2}^{\left(i\right)}{\psi }_{2}\hfill \end{matrix}\right),\end{equation} \tag{ 75 }$$

provided that:

$$\begin{equation}{\alpha }_{1}^{\left(i\right)}\frac{\partial {\varepsilon }_{j}}{\partial {N}_{1}^{\left(0\right)}}+{\alpha }_{2}^{\left(i\right)}\frac{\partial {\varepsilon }_{j}}{\partial {N}_{2}^{\left(0\right)}}=\frac{{\beta }_{j}^{\left(i\right)}}{{M}_{i}}.\end{equation} \tag{ 76 }$$

The solutions (74) can be easily 'upgraded' to four-component vectors solving equation (66), ${\tilde {U}}_{0}^{\left(i\right)}={\left({u}_{01}^{\left(i\right)},{u}_{01}^{\left(i\right)},{u}_{02}^{\left(i\right)},{u}_{02}^{\left(i\right)}\right)}^{\text{T}}={\alpha }_{1}^{\left(i\right)}{\left(\frac{\partial {\psi }_{1}}{\partial {N}_{1}^{\left(0\right)}},\frac{\partial {\psi }_{1}}{\partial {N}_{1}^{\left(0\right)}},\frac{\partial {\psi }_{2}}{\partial {N}_{1}^{\left(0\right)}},\frac{\partial {\psi }_{2}}{\partial {N}_{1}^{\left(0\right)}}\right)}^{\text{T}}+{\alpha }_{2}^{\left(i\right)}{\left(\frac{\partial {\psi }_{1}}{\partial {N}_{2}^{\left(0\right)}},\frac{\partial {\psi }_{1}}{\partial {N}_{2}^{\left(0\right)}},\frac{\partial {\psi }_{2}}{\partial {N}_{2}^{\left(0\right)}},\frac{\partial {\psi }_{2}}{\partial {N}_{2}^{\left(0\right)}}\right)}^{\text{T}}$.

Zero-energy vectors ${\tilde {V}}_{0}^{\left(i\right)}={\left({\beta }_{1}^{\left(i\right)}{\psi }_{1},-{\beta }_{1}^{\left(i\right)}{\psi }_{1},{\beta }_{2}^{\left(i\right)}{\psi }_{2},-{\beta }_{2}^{\left(i\right)}{\psi }_{2}\right)}^{\text{T}}$ and the vectors ${\tilde {U}}_{0}^{\left(i\right)}$ must have unit overlap. If in addition we account for the constraints following from the orthogonality conditions, (64), we get:

$$\begin{equation}{\alpha }_{1}^{\left(i\right)}{\beta }_{1}^{\left(j\right)}+{\alpha }_{2}^{\left(i\right)}{\beta }_{2}^{\left(j\right)}={\delta }_{i,j}.\end{equation} \tag{ 77 }$$

Finally, there are eight conditions: four involving masses (76), and four orthonormality conditions (77) which have to be solved for ten unknowns: eight coefficients ${\alpha }_{j}^{\left(i\right)}$, ${\beta }_{j}^{\left(i\right)}$ and two masses M_i . The solutions for α-s are:

$$\begin{align*}\hfill \left({\alpha }_{1}^{\left(1\right)},{\alpha }_{2}^{\left(1\right)}\right)& =\frac{1}{{c}^{\left(1\right)}}\frac{1}{\sqrt{2}}\left(\frac{1}{\sqrt{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}^{\left(0\right)}}}},\frac{1}{\sqrt{\frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}^{\left(0\right)}}}}\right),\hfill \\ \hfill \left({\alpha }_{1}^{\left(2\right)},{\alpha }_{2}^{\left(2\right)}\right)& =\frac{1}{{c}^{\left(2\right)}}\frac{1}{\sqrt{2}}\left(\frac{1}{\sqrt{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}^{\left(0\right)}}}},-\frac{1}{\sqrt{\frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}^{\left(0\right)}}}}\right),\hfill \end{align*}$$

and similarly, solutions for β-s have the form:

$$\begin{align}\hfill \left({\beta }_{1}^{\left(1\right)},{\beta }_{2}^{\left(1\right)}\right)& ={c}^{\left(1\right)}\frac{1}{\sqrt{2}}\left(\sqrt{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}^{\left(0\right)}}},\sqrt{\frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}^{\left(0\right)}}}\right),\hfill \\ \hfill \left({\beta }_{1}^{\left(2\right)},{\beta }_{2}^{\left(2\right)}\right)& ={c}^{\left(2\right)}\frac{1}{\sqrt{2}}\left(\sqrt{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}^{\left(0\right)}}},-\sqrt{\frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}^{\left(0\right)}}}\right).\hfill \end{align} \tag{ 78 }$$

Finally the 'inverse of masses' M_i of the two phase modes are:

$$\begin{align}\hfill \frac{1}{{M}_{1}}& =\frac{1}{{\left({c}^{\left(1\right)}\right)}^{2}}\left(1+\frac{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{2}^{\left(0\right)}}}{\sqrt{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}^{\left(0\right)}}\frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}^{\left(0\right)}}}}\right),\hfill \\ \hfill \frac{1}{{M}_{2}}& =\frac{1}{{\left({c}^{\left(2\right)}\right)}^{2}}\left(1-\frac{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{2}}}{\sqrt{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}^{\left(0\right)}}\frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}^{\left(0\right)}}}}\right).\hfill \end{align} \tag{ 79 }$$

In the equations above c⁽ⁱ⁾ 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 zero-energy 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}^{\left(i\right)}={\left(\partial {\varepsilon }_{1}/\partial {N}_{1}^{\left(0\right)}\right)}^{-1/4}{\left(\partial {\varepsilon }_{2}/\partial {N}_{2}^{\left(0\right)}\right)}^{-1/4}$.

The amplitudes of zero-energy modes, $-i{\vartheta }_{0}^{\left(i\right)}$, and amplitudes of the modes coupled to them, ${\Delta}{\mathcal{N}}_{0}^{\left(i\right)}$, are some parameters characterizing the perturbation field, $\tilde {\delta \psi }$. They are related to variations of phase and number of particles in both components. Indeed, a deviation of atom number from the equilibrium value calculated directly from the definition (see equation (47) for details) using the explicit form of the perturbation

$$\begin{equation}{\delta }_{i}\left(\mathbf{r},t\right)=\sum\limits _{j=1}^{2}\left(-i{\vartheta }_{0}^{\left(j\right)}\left(t\right){\beta }_{i}^{\left(j\right)}{\psi }_{i}+{\Delta}{\mathcal{N}}_{0}^{\left(j\right)}{u}_{0i}^{\left(j\right)}\left(\mathbf{r}\right)\right)\end{equation} \tag{ 80 }$$

reads

$$\begin{equation}\delta {N}_{i}={N}_{i}\left(t\right)-{N}_{i}^{\left(0\right)}=\int \mathrm{d}\mathbf{r}\enspace {\psi }_{i}\left({\delta }_{i}+{\delta }_{i}^{{\ast}}\right)=\sum\limits _{j=1,2}{\Delta}{\mathcal{N}}_{0}^{\left(j\right)}\int \mathrm{d}\mathbf{r}\left(2{\psi }_{i}{u}_{0i}^{\left(j\right)\text{ad}}\right).\end{equation} \tag{ 81 }$$

Utilizing orthogonality conditions (64), $\int \mathrm{d}\mathbf{r}\enspace \left(2{\beta }_{1}^{\left(j\right)}{u}_{01}^{\left(i\right)\text{ad}}\left(\mathbf{r}\right){\psi }_{1}\left(\mathbf{r}\right)+2{\beta }_{2}^{\left(j\right)}{u}_{01}^{\left(i\right)\text{ad}}{\psi }_{1}\left(\mathbf{r}\right)\right)={\delta }_{i,j}$, we see that the amplitudes of the coupled modes are physical observables, namely deviations of particle number from the equilibrium value:

$$\begin{equation}{\Delta}{\mathcal{N}}_{0}^{\left(i\right)}={\beta }_{1}^{\left(i\right)}\delta {N}_{1}+{\beta }_{2}^{\left(i\right)}\delta {N}_{2}.\end{equation} \tag{ 82 }$$

Amplitudes ${\Delta}{\mathcal{N}}_{0}^{\left(i\right)}$ of modes ${\tilde {U}}^{\left(i\right)}$ 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 ${\vartheta }_{0}^{\left(i\right)}$ of the zero-energy modes. To show this we remind that a small shift of the phases of droplets' wavefunctions ${\psi }_{i}\left(\mathbf{r},0\right){\text{e}}^{-\text{i}{\theta }_{i}}\simeq {\psi }_{i}\left(\mathbf{r},0\right)\left(1-i{\theta }_{i}\right)$ leads to the following perturbation field δ_i (r, 0) = −iθ_i ψ_i . Combining this result with the assumed form of the perturbation field (80): ${\delta }_{i}=-i{\vartheta }_{0}^{\left(1\right)}{\beta }_{i}^{\left(1\right)}{\psi }_{i}-i{\vartheta }_{0}^{\left(2\right)}{\beta }_{i}^{\left(2\right)}{\psi }_{i}$, we get

$$\begin{equation}{\theta }_{i}={\vartheta }_{0}^{\left(1\right)}{\beta }_{i}^{\left(1\right)}+{\vartheta }_{0}^{\left(2\right)}{\beta }_{i}^{\left(2\right)}.\end{equation} \tag{ 83 }$$

Now using the orthogonality conditions (77) we find

$$\begin{equation}{\vartheta }_{0}^{\left(i\right)}={\alpha }_{1}^{\left(i\right)}{\theta }_{1}+{\alpha }_{2}^{\left(i\right)}{\theta }_{2}.\end{equation} \tag{ 84 }$$

The amplitudes ${\vartheta }_{0}^{\left(i\right)}$ of the zero-energy phase modes ${\tilde {V}}^{\left(i\right)}$ are combinations of phases θ_i of the wavefunctions of the two components.

Excitations of modes coupled to zero-energy modes contribute to the total energy of the system. We remind that the Bogoliubov Hamiltonian, ${H}_{\text{B}}=\frac{1}{2}\langle {\tilde {\delta \psi }}^{\text{ad}}\vert \mathcal{L}\vert \tilde {\delta \psi }\rangle$, (compare equation (45)) is a sum of the Hamiltonian of the standard Bogoliubov excitations and the Hamiltonian of the zero-energy modes, H_B = ∑_ν ɛ_Bν |b_ν |² + H₀. The contribution of these modes takes the form of the kinetic energy of a free particle having momentum ${\Delta}{\mathcal{N}}_{0}^{\left(i\right)}$ and mass M_i . There are two degrees of freedom, therefore there are two contributions to the energy:

$$\begin{equation}{H}_{0}={H}_{1}+{H}_{2},\end{equation} \tag{ 85 }$$

where:

$$\begin{equation}{H}_{i}=\frac{{\left({\Delta}{\mathcal{N}}_{0}^{\left(i\right)}\right)}^{2}}{2{M}_{i}}=\frac{{\left({\beta }_{1}^{\left(i\right)}\delta {N}_{1}+{\beta }_{2}^{\left(i\right)}\delta {N}_{2}\right)}^{2}}{2{M}_{i}}.\end{equation} \tag{ 86 }$$

Using equations (78) and (79) the explicit form of the Hamiltonians (86) is found:

$$\begin{equation}{H}_{1}=\frac{1}{4}\left(1+\frac{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{2}}}{\sqrt{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}}\frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}}}}\right){\left(\sqrt{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}}}\delta {N}_{1}+\sqrt{\frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}}}\delta {N}_{2}\right)}^{2},\end{equation} \tag{ 87 }$$

$$\begin{equation}{H}_{2}=\frac{1}{4}\left(1-\frac{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{2}}}{\sqrt{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}}\frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}}}}\right){\left(\sqrt{\frac{\partial {\varepsilon }_{1}}{\partial {N}_{1}}}\delta {N}_{1}-\sqrt{\frac{\partial {\varepsilon }_{2}}{\partial {N}_{2}}}\delta {N}_{2}\right)}^{2}.\end{equation} \tag{ 88 }$$

Equations (87) and (88) are our main results. These equations define a formalism allowing for studies of dynamical processes related to small variations of phase and number of particles of quantum two-component droplets. Excitations described by H₁ have energies much smaller than those described by H₂. This can be seen by observing that $\partial {\varepsilon }_{1}/\partial {N}_{2}\propto {\partial }^{2}\mathcal{E}/\left(\partial {n}_{1}\partial {n}_{2}\right)\sim -\vert {a}_{12}\vert {< }0$, see equation (12), therefore 1/M₁ is much smaller than 1/M₂. In the case of Bose–Bose mixtures the Hamiltonian H₁ corresponds to the soft zero-energy mode while the Hamiltonian H₂ corresponds to zero energy excitations of the hard mode. A detailed discussion is given in section 4.2 devoted to Bose–Bose-mixtures.

Finally, the classical Hamiltonians above can be easily quantized by imposing canonical commutation relations on canonical variables. Evidently the pairs $\left({\theta }_{1},\delta {N}_{1}\right)$ and $\left({\theta }_{2},\delta {N}_{2}\right)$ 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, ${\theta }_{i}\to {\hat{\theta }}_{i}$ and $\delta {N}_{i}\to \delta {\hat{N}}_{i}$ and impose canonical commutation relations:

$$\begin{equation}\left[{\hat{\theta }}_{i},\delta {\hat{N}}_{j}\right]=i{\delta }_{i,j}.\end{equation} \tag{ 89 }$$

A convenient choice of representation of canonical operators is $\delta {\hat{N}}_{i}=-i\frac{\partial }{\partial {\theta }_{i}}$.

Consistently we can easily verify that the phases of the soft ${\hat{\vartheta }}_{0}^{\left(1\right)}$ and hard modes ${\hat{\vartheta }}_{0}^{\left(2\right)}$ are canonically conjugate to fluctuations of the soft and hard mode occupations respectively, ${\Delta}{\hat{\mathcal{N}}}_{0}^{\left(1\right)}$ and ${\Delta}{\hat{\mathcal{N}}}_{0}^{\left(2\right)}$:

$$\begin{equation}\left[{\hat{\vartheta }}_{i},{\Delta}{\mathcal{N}}_{0}^{\left(j\right)}\right]=i{\delta }_{ij}.\end{equation} \tag{ 90 }$$

Again, a possible choice of representing the fluctuation operators is ${\Delta}{\mathcal{N}}_{0}^{\left(1\right)}=-i\frac{\partial }{\partial {\vartheta }_{1}}$ and ${\Delta}{\hat{\mathcal{N}}}_{0}^{\left(2\right)}=-i\frac{\partial }{\partial {\vartheta }_{2}}$.

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.

The considerations presented above are quite general. They can be applied to any system with broken U(1) and translational symmetries. The results obtained were based on the assumption that interaction energy, $\mathcal{E}\left({n}_{1},{n}_{2}\right)$, is a function of densities of both species involved, n_i = |ψ_i |². 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 (14). Eigenvalues, ɛ_i (ψ₁, ψ₂), account for contributions from kinetic and interaction energies. The stationary densities, n_i = |ψ_i |², 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₁, n₂ of volume V. Then its total energy is equal to $\mathcal{E}\left({n}_{1},{n}_{2}\right)V$ where $\mathcal{E}\left({n}_{1},{n}_{2}\right)$ is given by equation (12). Here ${N}_{1}^{\left(0\right)}={n}_{1}V$ and ${N}_{2}^{\left(0\right)}={n}_{2}V$ are the numbers of particles in both components. At equilibrium we have ${\left(\frac{\partial E}{\partial V}\right)}_{{N}_{1}^{\left(0\right)},{N}_{2}^{\left(0\right)}}=0$ which results in $V\left({N}_{1}^{\left(0\right)},{N}_{2}^{\left(0\right)}\right)$. With such a model it is questionable to use the solution of equation (70) given by (73). There to get the masses M_i the EGP equations were differentiated (71) with respect to ${N}_{i}^{\left(0\right)}$ 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}_{\text{GP}}^{\left(i\right)}$, 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 equation (73) we analytically solve equation (70), after neglecting the kinetic energy term in the operator $\mathcal{K}$. We arrive at a set of linear algebraic equations which together with orthogonality conditions (64) lead to the Hamiltonians

$$\begin{equation}{H}_{1}=\frac{1}{4V}\left(1+\frac{\frac{\partial {\mu }_{1}}{\partial {n}_{2}}}{\sqrt{\frac{\partial {\mu }_{1}}{\partial {n}_{1}}\frac{\partial {\mu }_{2}}{\partial {n}_{2}}}}\right){\left(\sqrt{\frac{\partial {\mu }_{1}}{\partial {n}_{1}}}\delta {N}_{1}+\sqrt{\frac{\partial {\mu }_{2}}{\partial {n}_{2}}}\delta {N}_{2}\right)}^{2},\end{equation} \tag{ 91 }$$

$$\begin{equation}{H}_{2}=\frac{1}{4V}\left(1-\frac{\frac{\partial {\mu }_{1}}{\partial {n}_{2}}}{\sqrt{\frac{\partial {\mu }_{1}}{\partial {n}_{1}}\frac{\partial {\mu }_{2}}{\partial {n}_{2}}}}\right){\left(\sqrt{\frac{\partial {\mu }_{1}}{\partial {n}_{1}}}\delta {N}_{1}-\sqrt{\frac{\partial {\mu }_{2}}{\partial {n}_{2}}}\delta {N}_{2}\right)}^{2}.\end{equation} \tag{ 92 }$$

In the model studied here the derivatives ∂μ_j /∂n_i can be calculated directly and substituted into the Hamiltonians H₁ and H₂:

$$\begin{equation}{H}_{1}\simeq \frac{4\pi {\hslash }^{2}}{m}\frac{1}{16V}\frac{\delta a}{\sqrt{{a}_{11}{a}_{22}}}{\left(\sqrt{{a}_{11}}\delta {N}_{1}+\sqrt{{a}_{22}}\delta {N}_{2}\right)}^{2}\end{equation} \tag{ 93 }$$

$$\begin{equation}{H}_{2}\simeq \frac{4\pi {\hslash }^{2}}{m}\frac{1}{2V}{\left(\sqrt{{a}_{11}}\delta {N}_{1}-\sqrt{{a}_{22}}\delta {N}_{2}\right)}^{2}.\end{equation} \tag{ 94 }$$

The contribution of the LHY term to the Hamiltonian H₁ is essential. H₁ is proportional to the small parameter $\delta a/\sqrt{{a}_{11}{a}_{22}}\ll 1$. It describes the zero-energy excitations of the soft mode. This is not the case of the H₂ Hamiltonian. Here terms of the order of $\delta a/\sqrt{{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 $\delta a/\sqrt{{a}_{11}{a}_{22}}$, $1/{M}_{1}\sim \delta a/\sqrt{{a}_{11}{a}_{22}}$, and is much larger than the mass of the hard mode M₁ ≫ M₂. The dynamics of the soft mode phase is very slow, much slower than the evolution of the hard mode phase.

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:

$$\begin{equation}{\psi }_{i}^{\text{tr}}\left(\mathbf{r}+{z}_{0}{\mathbf{e}}_{z}\right)={\text{e}}^{{z}_{0}{\partial }_{z}}{\psi }_{i}\left(\mathbf{r}\right)={\psi }_{i}\left(\mathbf{r}\right)+{z}_{0}{\partial }_{z}{\psi }_{i}\left(\mathbf{r}\right),\end{equation} \tag{ 95 }$$

which indicates that the zero-energy mode has the form ${\tilde {V}}_{\text{tr}}={\left({\partial }_{z}{\psi }_{1},{\partial }_{z}{\psi }_{1}^{{\ast}},{\partial }_{z}{\psi }_{2},{\partial }_{z}{\psi }_{2}^{{\ast}}\right)}^{\text{T}}$.

The vector ${\tilde {U}}_{\text{tr}}$ can be defined as the one which couples dynamically to the zero-energy vector:

$$\begin{equation}i\hslash {\partial }_{t}\left({z}_{0}{\tilde {V}}_{\text{tr}}\right)=\mathcal{L}\left(i\frac{{p}_{z}}{\hslash }{\tilde {U}}_{tr}\right),\end{equation} \tag{ 96 }$$

where ip_z /ℏ is the initial amplitude. Note that p_z must be real if equation (96) is fulfilled.

Again, in the following we choose the phases of the stationary solutions of the EGP equations to be equal to zero which assures that the droplet's wavefunctions ψ₁ and ψ₂ are real. Therefore only two components of the zero-energy vector are essential, ${\tilde {V}}_{\text{tr}}={\left({\partial }_{z}{\psi }_{1},{\partial }_{z}{\psi }_{1},{\partial }_{z}{\psi }_{2},{\partial }_{z}{\psi }_{2}\right)}^{\text{T}}$. The dynamically coupled vector, ${\tilde {U}}_{\text{tr}}$, satisfies the following equation:

$$\begin{equation}\mathcal{L}{\tilde {U}}_{\text{tr}}=\frac{{\hslash }^{2}}{{M}_{\text{tr}}}{\tilde {V}}_{\text{tr}}.\end{equation} \tag{ 97 }$$

Consistently, equation (97) indicates that the initial amplitude of the zero-energy mode grows linearly in time as the adjoint mode is populated:

$$\begin{equation}{z}_{0}\left(t\right)={z}_{0}+\frac{{p}_{z}}{{M}_{\text{tr}}}t.\end{equation} \tag{ 98 }$$

Equation (96) form a set of four coupled equations. By adding both sides of the first and both sides of the second pair of this set we get:

$$\begin{equation}\left({H}_{\text{GP}}^{\left(1\right)}+2{n}_{1}\frac{\partial {\mu }_{1}}{\partial {n}_{1}}\right)\left({u}_{1}^{\text{tr}}+{v}_{1}^{\text{tr}}\right)+2\frac{\partial {\mu }_{1}}{\partial {n}_{2}}{\psi }_{1}{\psi }_{2}\left({u}_{1}^{\text{tr}}+{v}_{1}^{\text{tr}}\right)=0,\end{equation} \tag{ 99 }$$

and similar equations with interchanged indices. The solutions are quite simple, ${v}_{i}^{\text{tr}}=-{u}_{i}^{\text{tr}}$. The number of essential components of the adjoint vector reduces then to ${\left({u}_{1}^{\text{tr}},{u}_{2}^{\text{tr}}\right)}^{\text{T}}$. In turn, by subtracting the two first equations of set equation (96) and repeating the same operation for the second pair, we get:

$$\begin{equation}{H}_{\text{GP}}^{\left(i\right)}{u}_{i}^{\text{tr}}=\frac{{\hslash }^{2}}{{M}_{tr}}{\partial }_{z}{\psi }_{i}.\end{equation} \tag{ 100 }$$

It can be checked by inspection that solutions of equation (100) are:

$$\begin{equation}{u}_{i}^{\text{tr}}=-\frac{m}{{M}_{\text{tr}}}z{\psi }_{i}.\end{equation} \tag{ 101 }$$

Utilizing the normalization condition (here all four components of the involved vectors are needed), $\langle {\tilde {U}}_{\text{tr}}^{\text{ad}}\vert {\tilde {V}}_{\text{tr}}\rangle =1$, the 'translational mass', M_tr, can be determined:

$$\begin{equation}{M}_{\text{tr}}=m\left({N}_{1}^{\left(0\right)}+{N}_{2}^{\left(0\right)}\right).\end{equation} \tag{ 102 }$$

To not a big surprise, the translational mass is equal to the mass of the entire droplet.

The translational modes are orthogonal to the Bogoliubov modes and to the phase modes. The last fact follows from different spatial symmetries of the modes. The mean field perturbation which is related to the translation of the droplet has the form:

$$\begin{equation}{\delta }_{i}\left(\mathbf{r},t\right)={z}_{0}\left(t\right){\partial }_{z}{\psi }_{i}-i\frac{{p}_{z}}{\hslash }\frac{m}{{M}_{\text{tr}}}z{\psi }_{i}.\end{equation} \tag{ 103 }$$

It is a simple exercise to check that the amplitude of the adjoint translational mode is nothing else but the mean momentum of the center of mass of the droplet:

$$\begin{equation}{P}_{z}=\sum\limits _{i}\int \mathrm{d}\mathbf{r}\left({\psi }_{i}+{\delta }_{i}^{{\ast}}\right)\left(-i\hslash {\partial }_{z}\right)\left({\psi }_{i}+{\delta }_{i}\right)={p}_{z}.\end{equation} \tag{ 104 }$$

The Hamiltonian generating the motion of the droplet as a whole can be obtained the same way we got the phase zero-energy modes Hamiltonians. It is quite obvious that it must correspond to the energy of a free point-like particle with momentum p_z and mass of the droplet M_tr:

$$\begin{equation}{H}_{\text{tr}}=\frac{{p}_{z}^{2}}{2{M}_{\text{tr}}}.\end{equation} \tag{ 105 }$$

Although the result could be easily anticipated on the basis of general arguments, nevertheless obtaining it within the frame of the quite sophisticated formalism of zero-energy modes is a nice physical illustration of the approach. The quantum description of the free motion of a droplet can be obtained by substituting classical variables z₀ and p_z by operators ${z}_{0}\to {\hat{z}}_{0}$, and ${p}_{z}\to {\hat{p}}_{z}$ and imposing the canonical commutation relation $\left[{\hat{z}}_{0},{\hat{p}}_{z}\right]=i\hslash$.

The result obtained, in particular the quantized version of the translational Hamiltonian (105) predicts quantum diffusion of the position of the center of mass of the droplet. But characteristic time for this process is by many orders of magnitudes larger than the life time of the droplet. Droplet is a macroscopic object and no quantum effects in the dynamics of its center of mass are expected to be observed.