Bosonic bright soliton in the mixture of repulsive Bose-Einstein condensate and polarized ultracold fermions under influence of the pressure evolution
Pavel A. Andreev K. V. Antipin Mariya Iv. Trukhanova
cond-mat.quant-gas







andreevpa@physics.msu.ru
Department of General Physics, Faculty of physics, Lomonosov Moscow State University, Moscow, Russian Federation, 119991.


kv.antipin@physics.msu.ru
Department of Quantum Statistics and Field Theory,
Faculty of physics, Lomonosov Moscow State University, Moscow, Russian Federation, 119991.


mar-tiv@yandex.ru
Department of General Physics, Faculty of physics, Lomonosov Moscow State University, Moscow, Russian Federation, 119991.
Russian Academy of Sciences, Nuclear Safety Institute (IBRAE), B. Tulskaya 52, Moscow, Russian Federation, 115191.




Repulsive Bose-Einstein condensate,
where the short-range interaction is included up to the third order by the interaction radius,
demonstrates existence of a bright soliton in a narrow interval of parameters.
This soliton is studied here for the boson-fermion mixture,
where spin-1/2 fermions are consider in the regime of full spin polarization.
Influence of fermions via the boson-fermion interaction
is considered up to the third order by the interaction radius.
Fermions themselves are considered by hydrodynamic model including the pressure evolution equation.
Interaction between fermions is considered.
The first order by the interaction radius gives zero contribution in the Euler equation and the pressure evolution equation,
but the third order by the interaction radius provides nonzero contributions in both equations.
Repulsive (attractive) boson-fermion interaction leads to the bright (dark) fermionic soliton.


03.75.Hh, 03.75.Kk, 67.85.Pq
                             










Bosonic bright soliton in the mixture of repulsive Bose-Einstein condensate
and polarized ultracold fermions under influence of the pressure evolution
    Mariya Iv. Trukhanova
    December 30, 2023
======================================================================================================================================================











§ INTRODUCTION


Solitons <cit.>, <cit.>, <cit.>,
vorticities <cit.>, <cit.>, <cit.>,
skyrmions <cit.> are fundamental nonlinear excitations in quantum gases.
The quantum droplets formation is also detected
in Bose-Einstein condensate of rare-earth atoms
due to the dipole-dipole interaction and quantum correlations <cit.>, <cit.>, <cit.>.
This work is focused on solitons in boson-fermion mixtures  <cit.>
studied in terms of hydrodynamic model.
The boson-fermion mixtures are experimentally obtained in different combinations, such as
^7Li-^6Li <cit.>, <cit.>,
^23Na-^6Li <cit.>,
and ^87Rb-^40K <cit.>.
If Bose-Einstein condensate (BEC) is considered in terms of the Gross-Pitaevskii equation <cit.>
it reveals the bright soliton at the attraction between bosons
or the dark soliton in the repulsive BECs.
Noninteracting degenerate fermions also demonstrate existence of soliton <cit.>,
where pair of a bright and dark solitons for trapped fermions is obtained.
However, the extended model with nonlocal short-range interaction considered up to the third order by the interaction radius
leads to formation of a bright soliton (earlier called the bright-like soliton) in repulsive BECs <cit.>.
The bright soliton in repulsive BECs is studied for the anisotropic short-range interaction <cit.>,
which is also considered up to the third order by the interaction radius.

Fermionic bright solitons at attractive boson-fermion interaction are found in Refs.
<cit.>, <cit.>.
Fermionic dark solitons are obtained in <cit.>.
The dark solitons in the Fermi superfluid within the Bogoliubov-de Gennes theory of the BEC to BCS crossover
using the random-phase approximation are considered in <cit.>,
where the decay of the soliton via the snake instability is studied.


BEC-BCS crossover in the two component fermion systems is actively studying branch of the field of ultracold fermions.
Sanner et. al. <cit.> considers two-component fermions with strong repulsive interaction between fermions with different spin projections.
It is demonstrated that pairing instability is faster than the ferromagnetic instability in this regime.
However, the average and weak repulsion shows stable coexisting of fermions with different spin projections.
Interaction of fermions of the same spin projection is not discussed there.
Metastable Stoner-like ferromagnetic phase supported by strong repulsion of fermions
with opposite spin projections in excited scattering states is studied in Ref. <cit.>.
Mass-imbalanced Fermi-Fermi mixture of ^161Dy and ^40K is created to study the strongly interacting regime  <cit.>.
The existence of the second sound in a system of one-dimensional fermions with repulsive interactions is reported in <cit.>,
where hydrodynamic equations for four conserved macroscopic characteristics of the fluid are presented
for the Luttinger liquid with linear excitation spectrum <cit.>, <cit.>, <cit.>.
This theoretical result corresponds to earlier experimental observation of the first and second sounds in ^6Li atoms <cit.>.

Single state ultracold fermions are also currently studied.
Single spin state degenerate fermions ^6Li confined
by the heavy bosonic atoms ^133Cs at attractive interspecies interaction is created <cit.>.
Rakshit et. al. <cit.> demonstrates that
higher order corrections to the standard mean-field energy can lead
to a formation of Bose-Fermi liquid droplets for attractive BEC and spin-polarized fermions.
A hydrodynamic approach including set of two nonlinear Schrodinger equations is used in <cit.>.
Weakly interacting single-component two-dimensional dipolar Fermi gas is considered in Ref. <cit.> to study the zero sound.


Mixtures of two BECs are studied either.
Bright-bright, dark-bright and dark-dark types of solitons are theoretically obtained in two-component BEC <cit.>, <cit.>.



A possibility that the boson-fermion interaction in the third order by the interaction radius approximation leads
to formation of new soliton is discussed in Ref. <cit.>,
where it is corresponds to a step beyond the mean-field approximation.
However, here we report an advanced version of hydrodynamic of fermions.
Therefore, result of Ref. <cit.> can be readdressed in terms of novel model,
where less assumptions are made for the force field in the Euler equation.
Moreover, the pressure evolution equation is included here.
The presence of the pressure evolution equation becomes especially important since the force field contains the pressure tensor.
Furthermore, the second interaction constants for each type of interatomic potential are not represented via the scattering length,
but they are considered as independent constants as it follows from their definitions.

Soliton dynamics beyond the mean-field approximation is studied in Refs.
<cit.>, <cit.>, <cit.>, <cit.>.
It is the dark solitons in BoseEinstein condensates <cit.>
and bright-dark solitons in two component quantum gases <cit.>.
Many-body quantum dynamics including the correlation dynamics is considered there.



Microscopic derivation of the quantum hydrodynamic model for boson-fermion mixture, where fermions are located in the single spin state is presented.
General structure of equations for balance of the particle number, momentum and the momentum flux is obtained.
The two-particle short-range boson-boson, boson-fermion, and fermion-fermion interactions are included in the model.
General form of the interaction in the momentum balance equations and the momentum flux evolution equations are found.
Weakly interacting limit of the interaction terms for bosons being in the BEC state and the degenerate fermions is derived up to the third order by the interaction radius.
In this regime the mixture is characterized by five interaction constants.
Small amplitude nonlinear evolution of collective excitations is considered
to find bright-bright and bright-dark solitons in the boson-fermion mixture.
These solitons exist purely due to the interaction constants existing in the third order by the interaction radius.
These solitons are solitons of concentrations of bosons and fermions.
The velocity field of bosons and fermions together with the diagonal elements of pressure reveal the soliton structure either.



This paper is organized as follows.
In Sec. II major steps of derivation of hydrodynamic equations from the Schrodinger equation are demonstrated.
In Sec. III the quantum hydrodynamic equations for BEC and ultracold spin-polarized fermions.
Both species are described up to the third order by the interaction radius is presented.
In Sec IV method of the approximate nonlinear solution is presented.
In Sec. V numerical analysis of obtained Korteweg-de Vries equation is shown.
In Sec. VI the obtained results are summarized.







§ ON DERIVATION OF HYDRODYNAMIC EQUATIONS FROM MICROSCOPIC QUANTUM DYNAMICS



Boson-fermion mixture consisting of N particles,
which is superposition of number of bosons N_b and number of fermions N_f,
is described by the Schrodinger equation with the following Hamiltonian

    Ĥ=∑_i=1^N(p̂^2_i/2m_i+V_ext(r_i,t))
    +1/2∑_i,j≠ iU(r_i-r_j) ,

where m_i is the mass of i-th particle,
p̂_i=-ħ∇_i is the momentum of i-th particle.
The last term in the Hamiltonian (<ref>)
is the short representation of
boson-boson
(1/2)∑_i=1^N_b∑_j=1,j≠ i^N_bU_bb(r_i-r_j)
interparticle interaction,
fermion-fermion
(1/2)∑_i=N_f^N_b+N_f∑_j=N_f,j≠ i^N_b+N_fU_ff(r_i-r_j)
interparticle interaction,
and boson-fermion
(1/2)[∑_i=1^N_b∑_j=N_f^N_b+N_f
+∑_i=N_f^N_b+N_f∑_j=1^N_b]U_bf(r_i-r_j)
=∑_i=1^N_b∑_j=N_f^N_b+N_fU_bf(r_i-r_j)
interparticle interaction.
The Schrodinger equation
ħ∂_tΨ(R,t)=ĤΨ(R,t)
with Hamiltonian (<ref>)
describes the evolution of wave function of full boson-fermion mixture Ψ(R,t),
where full configurational space R={R_b,R_f}
is the combination of configurational space of bosons R_b and configurational space of fermions R_f.


However, there are formulations of the many-body problem for fermions with no interaction between fermions in the same spin state
see for instance equation 23 in Ref. <cit.>.
It is due to the argument that interactions are strongly inhibited by the Pauli exclusion principle.
The antisymmetry of the many-particle wave function is the manifistation of the Pauli exclusion principle.
However, we are not used properties of the wave function at this step.
Mentioning here some of our results we point out that
the first order on the interaction radius (an analog of s-wave) contribution is equal to zero due to the antisymmetry of the many-particle wave function.
However, the third order on the interaction radius (an analog of p-wave) terms are non zero.


At this stage bosons and fermions have arbitrary distributions on quantum states.
Transition to near equilibrium states with zero temperature is made at later stage of derivation,
where truncation of the chain of equations is made.


Concentration of bosons is defined as the quantum mechanical average of the operator of concentration
which is the superposition of the delta functions <cit.>, <cit.>:

    n_b(r,t)=∫
    dR∑_i=1^N_bδ(r-r_i)Ψ^*(R,t)Ψ(R,t),

where dR=dR_b dR_f,
dR_b=∏_i=1^N_bdr_i is the element of volume in 3N_b dimensional configurational space,
with N_b is the number of bosons, and
dR_f=∏_i=N_b+1^N_b+N_fdr_i is the element of volume in 3N_f dimensional configurational space,
with N_f is the number of fermions.
We need to integrate over coordinates of all particles
since the wave function describes both species.

Definition of the concentration of fermion has similar structure

    n_f(r,t)=∫
    dR∑_i=N_b+1^N_b+N_fδ(r-r_i)Ψ^*(R,t)Ψ(R,t),

but operator of concentration of fermions contains coordinates of different set of particles.


Considering time evolution of each concentration n_a via the evolution of wave function Ψ(R,t)
find the continuity equation <cit.>, <cit.>:

    ∂_tn_a+∇·j_a=0,

where subindex a stands for b or f,
and the current j_a is defined via the many-particle wave function of the system:

    j_a(r,t)
    =∫ dR∑_i∈ N_aδ(r-r_i)×


    ×1/2m_i(Ψ^*(R,t)p̂_iΨ(R,t)+c.c.),

where c.c. is the complex conjugation.


Next derive equation for the current evolution.
Consider the time derivative of the current (<ref>) using the Schrodinger equation and some straightforward calculations.
As the result find the current evolution equation (it can be also called the momentum evolution equation)

    ∂_tj_a^α+∂_βΠ_a^αβ
    =-1/m_an_a∂_αV_ext+1/m_aF^α_a,int,

where

    Π_a^αβ=∫ dR∑_i∈ N_aδ(r-r_i) 1/4m_i^2
    [Ψ^*(R,t)p̂_i^αp̂_i^βΨ(R,t)


    +p̂_i^α *Ψ^*(R,t)p̂_i^βΨ(R,t)+c.c.]

is the momentum flux (containing the pressure tensor),
and

    F^α_a,int=-∑_a'=b,f∫ (∂^αU_aa'(r-r'))
    n_2,aa'(r,r',t)dr',

with the following expression for the two-particle concentration

    n_2,aa'(r,r',t)


    =∫
    dR∑_i∈ N_a,j∈ N_a',j≠ iδ(r-r_i)δ(r'-r_j)Ψ^*(R,t)Ψ(R,t) .



Suggested model includes the pressure evolution equation for fermions.
It is not required for bosons since bosons are considered below in the Bose-Einstein state.
Therefore, the pressure of bosons equals to zero and the bosons are completely described by the concentration and the velocity field.
Therefore, it is required to derive equation for the momentum flux evolution,
since the momentum flux has clear relation to the wave function (<ref>)
similarly to the concentrations and the currents (<ref>).
The pressure evolution will be extracted from momentum flux evolution below.
Similarly to derivation of the current evolution equation
consider the time derivative of the momentum flux (<ref>) using the Schrodinger equation <cit.>:

    ∂_tΠ_f^αβ+∂_γM_f^αβγ
    =-1/m_fj_f^β∂_αV_ext-1/m_fj_f^α∂_βV_ext


    -1/m_f∫[∂^βU(r-r')]j_2^α(r,r',t)dr'


    -1/m_f∫[∂^αU(r-r')]j_2^β(r,r',t)dr',

where

    M_f^αβγ=∫ dR∑_i∈ N_fδ(r-r_i) 1/8m_i^3[Ψ^*(R,t)p̂_i^αp̂_i^βp̂_i^γΨ(R,t)


    +p̂_i^α *Ψ^*(R,t)p̂_i^βp̂_i^γΨ(R,t)
    +p̂_i^α *p̂_i^γ *Ψ^*(R,t)p̂_i^βΨ(R,t)


    +p̂_i^γ *Ψ^*(R,t)p̂_i^αp̂_i^βΨ(R,t)+c.c.],

and

    j_2(r,r',t)=∫
    dR∑_i∈ N_f, j∈ N, j≠ iδ(r-r_i)δ(r'-r_j)×


    ×1/2m_i(Ψ^*(R,t)p̂_iΨ(R,t)+c.c.) .

If quantum correlations are dropped function j_2^α(r,r',t) splits on product of the current j_f^α(r,t) and the concentration n_f(r',t).

Equations (<ref>), (<ref>), (<ref>)
are fundamental equations for collection of bosons and fermions.
These equations contain a number of new functions which should be expressed via the basic hydrodynamic functions.
The truncation is to be made for the bosons being in the BEC state and fermions at zero temperature, but collected in the single spin state.

Moreover, it is necessary to present hydrodynamic equations
(<ref>), (<ref>), (<ref>) in more traditional form.
To this end, introduce the velocity field v_a=j_a/n_a.
This definition allows to represent the continuity equation in the traditional form.
However, other equations require more detail description which can be found in details in Refs. <cit.>, <cit.>.
the method of introduction of the velocity field includes
the analysis of deviation of velocities of quantum particles introduced as gradients ħ∇_iS/m of the phase of wave function Ψ(R,t)=a(R,t)exp( S(R,t)) from the velocity field.
This deviations are also includes the thermal effects and other mechanisms
(Fermi surface caused by the Pauli blocking)
of distribution of particles on quantum states with velocities shifted from the average velocity v_a.
This method provides the structure of the momentum flux tensor:

    Π_a^αβ=n_av_a^αv_a^β +p_a^αβ+T_a^αβ.

The first and second terms on the right-hand side of equation (<ref>)
have classical meaning and include the pressure tensor p^αβ.
The last term in equation (<ref>) has quantum nature and can be presented in the following approximate form

    T_a^αβ=-ħ^2/4m_a^2[∂_α∂_βn_a
    -∂_αn_a·∂_βn_a/n_a].

It is related to the quantum Bohm potential.


Equation (<ref>) appears for noninteracting bosons in the BEC state.
Its linear part (the first term) is straightforward for interacting bosons or interacting fermions.
While the second term in equation (<ref>) has no proper justification for fermions even if interaction is neglected.
Hence, the second term in equation (<ref>) for fermions is an approximate equation of state.

Representation for tensor M_f^αβγ (<ref>) similar to (<ref>)
can be found either.
It is given in Ref. <cit.> as set of four rather large equations (see equations 25-28).


Equations are obtained for particles with arbitrary spin,
but further analysis is made for spin-0 bosons and spin-1/2 fermions being in the single spin projection state.


Equations (<ref>)-(<ref>) are derived for arbitrary potential U_ij=U(r_i-r_j).
It is necessary to specify that neutral particles interact via the short-range potential.
To stress the small radius of interaction represent the coordinates of interacting particles r_i and r_j
via the relative distance and coordinate of their center of mass.
Next, we can expand the delta functions and the wave function on the small interparticle distance
r_ij=r_i-r_j
since potential U_ij goes to zero at the large interparticle distances.
The straightforward calculations for weakly interacting particles including symmetry between fermions,
between bosons, and the absence of symmetry between bosons and fermions (for more details see <cit.>, <cit.>).
Terms in the zeroth order on the interparticle distance cancel each other.
In the first order, there are nonzero terms for boson-boson and boson-fermion interactions,
which corresponds to the Gross-Pitaevskii approximation. 
The integral over interparticle distance contains the potential of interaction and gives the interaction constant for each interaction.
The fermion-fermion interaction term in the first order equals to zero due to the antisymmetry of the wave function.
Formally, we have the first interaction constant for fermions, but it is multiplied by the function which equal to zero.
The second order terms go to zero due to integration over the interparticle distance (its angular dependence).
The third order of expansion gives nonzero results for all three interactions.
The second interaction constant appears for each interaction.


The derivation of fundamental hydrodynamic equations is made by the many-particle quantum hydrodynamic method
<cit.>, <cit.>, <cit.>.
Paying attention to development of hydrodynamic methods mention the generalized hydrodynamics actively developing in recent years
<cit.>, <cit.>.






§ HYDRODYNAMIC EQUATIONS FOR BOSON-FERMION MIXTURE



In this regime we have two continuity equations:

    ∂_tn_b+∇· (n_bv_b)=0,

and

    ∂_tn_f+∇· (n_fv_f)=0.



We also have two Euler (momentum balance) equations.
The Euler equation for bosons

    m_bn_b(∂_t +v_b·∇)v^α_b
    -ħ^2/2m_bn_b∂^α√(n_b)/√(n_b)


    +g_b n_b∂^αn_b +1/2g_2b∂^α n_b^2 =-n_b∂^αV_ext


    -g_bf n_b∂^αn_f-g_2,bf/2n_b∂^α n_f

contains the boson-boson interaction in the first and third orders by the interaction radius,
the terms proportional to g_b and g_2b constants, correspondingly.
It includes the boson-fermion interaction
which are proportional to g_bf in the first order and g_2bf in the third order.

The Euler equation for fermions is

    m_fn_f(∂_t +v_f·∇)v^α_f
    -ħ^2/2m_fn_f∂^α√(n_f)/√(n_f)+∂^βp_f^αβ


    =-g_bf n_f∂^αn_b-g_2,bf/2n_f∂^α n_b -g_2fm_f^2/2ħ^2I_0^αβγδ∂^β(n_fp^γδ_f),

where all fermions are in the quantum states with the same spin projection.
Euler equation (<ref>) has contribution of the boson-fermion interaction.
Similar to the Euler equation for bosons (<ref>),
they are proportional to g_bf in the first order by the interaction radius and g_2bf in the third order.
The fermion-fermion interaction gives the single term in Euler equation (<ref>).
It appears in the third order by the interaction radius being proportional to g_2f constant.



The third order by the interaction radius approximation has similarity to the p-wave interaction.
The p-wave fermion-fermion interaction is studied in the boson-fermion mixtures at study of solitons in mixtures <cit.>.
Traditional p-wave approximation assumes an equation of state for the pressure in terms of the concentration.
However, our model gives more accurate analysis of pressure via the pressure evolution equation.

The boson-boson interaction in the third order by the interaction radius is presented by the nonlocal interaction term showing similarity to the models presented in Refs. <cit.>, <cit.>.

The pressure evolution equation for fermions is also a part of developed and applied hydrodynamic model

    ∂_tp_f^αβ +v_f^γ∂_γp_f^αβ +p_f^αγ∂_γv_f^β +p_f^βγ∂_γv_f^α
    +p_f^αβ∂_γv_f^γ


    =-m_f/8ħ^2g_2f{I_0^αγδμ
    [3 n_f^2 v_f^βv_f^δ∂^γv_f^μ
    +2n_fp_f^μδ(∂^γv_f^β-∂^βv_f^γ)]


    +I_0^βγδμ
    [3 n_f^2 v_f^αv_f^δ∂^γv_f^μ
    +2n_fp_f^μδ(∂^γv_f^α-∂^αv_f^γ)]},

where
I_0^αβγδ=
δ^αβδ^γδ+δ^αγδ^βδ+δ^αδδ^βγ.

Equations (<ref>)-(<ref>) contain the following interaction constants
g_b=∫ U_bbdr,
g_bf=∫ U_bfdr,
g_2b=(1/24)∫ r^2U_bbdr,
g_2bf=∫ r^2U_bfdr,
and
g_2f=∫ r^2U_ffdr.

The pressure evolution equation (<ref>) appears from equation (<ref>)
after extraction of the thermal part or other mechanisms of distribution of particles in the momentum space like the Pauli blocking for degenerate fermions.

Equation (<ref>) has no trace of the external potential and boson-fermion interaction.

The fermion-fermion interaction gives nonzero contribution in the third order by the interaction radius.
It consists of the structure of two terms
which is repeated twice to give it the form symmetric on free indexes
since the pressure tensor p_f^αβ is a symmetric tensor.
Let us have a closer look on each of two terms.
One is highly nonlinear and includes product of three velocities n_f^2 v_f^βv_f^δ∂^γv_f^μ.
Another term is proportional to the pressure tensor.
Moreover, it is proportional to the classical hydrodynamic vorticity
ε^γβδΩ^δ
=∂^γv_f^β-∂^βv_f^γ,
where
Ω^α=ε^αβγ∂_βv_γ is the vorticity of classic uncharged fluid.

The left-hand side of equation (<ref>) contains
the divergence of a third rank tensor ∂_γQ^αβγ
which is the average of product of three thermal velocities
(the thermal part of tensor M_f^αβγ (<ref>)),
while p^αβ is the average of product of two thermal velocities.
It is assumed to be equal to zero.
It is the equation of state obtained as extension of equilibrium value of tensor Q^αβγ.

The boson part of the model is developed in Refs. <cit.>, <cit.>.
The fermion part of the hydrodynamic model is derived in Refs. <cit.>, <cit.>.
Here, same is in Ref. <cit.>, the pressure evolution equation is considered in the long-wavelength limit,
so high order derivatives are neglected.
The interspecies interaction is addressed in terms of many-particle quantum hydrodynamic method in Refs. <cit.>, <cit.>.







§ PERTURBATION METHOD


Following papers <cit.> and <cit.>
we use the reductive perturbation method <cit.>, <cit.>
to study the.
According to this method all hydrodynamic values may be
represented as:


    n_b=n_0b+ε n_1b+ε^2 n_2b+... ,


    n_f=n_0f+ε n_1f+ε^2 n_2f+... ,


    v_b^x=ε v_1b+ε^2 v_2b+... ,


    v_f^x=ε v_1f+ε^2 v_2f+... ,

and

    p_f^ii=p_0f^ii+ε p_1f^ii+ε^2 p_2f^ii+... ,

where ii stands for xx, yy and zz
since all diagonal elements of the pressure tensor are involved in dynamics of longitudinal perturbations.

It is assumed that there are nonzero constant equilibrium concentrations and pressure of fermions.
The velocity fields are equal to zero in equilibrium.

We also performed the following "scaling" of variables:

    ξ = ε ^1/2(x-Vt)

and

    τ = ε ^3/2Vt .

The latter expression introduces so-called "slow" time.







 §.§ The first order perturbations


Substitute scaling of hydrodynamic functions (<ref>)-(<ref>)
and space-time variables (<ref>), (<ref>)
in basic equations (<ref>)-(<ref>).
Separate contributions appearing in different orders on parameter ε.
Extract equations in the lowest order on parameter ε and find
the continuity equation for bosons

    n_0b∂_ξv_1b-V∂_ξn_1b=0,

the Euler equation for bosons

    m_bV∂_ξv_1b-g_b∂_ξn_1b-g_bf∂_ξn_1f=0,

the continuity equation for fermions

    n_0f∂_ξv_1f-V∂_ξn_1f=0,

the Euler equation for fermions

    m_fn_0fV∂_ξv_1f+∂_ξp_1f^xx=-g_bfn_0f∂_ξn_1b


    -g_2fm_f^2/2ħ^2∂_ξ[n_0f(3p_1f^xx+p_1f^yy+p_1f^zz)+(3p_0f^xx+p_0f^yy+p_0f^zz)n_1f],

and equations for evolution of the elements of pressure tensor

    V∂_ξp_1f^xx-3p_0f^xx∂_ξv_1f=0,

and

    V∂_ξp_1f^yy-p_0f^yy∂_ξv_1f=0,

zz element is the same as yy element.

Fermi surface in equilibrium regime is assumed to be a sphere.
Therefore, we have p_0f^xx=p_0f^yy=p_0f^zz≡ p_0f.

Equations (<ref>)-(<ref>) show
that the perturbation of pressure in the direction of wave propagation is three times larger
than the perturbation in the perpendicular directions.








Equations (<ref>)-(<ref>) are uniform linear differential equations.
We obtain explicit relation between these functions using boundary conditions at infinity,
where all perturbations tend to zero.
They have nonzero solution if parameter V have the following form

    V^2=1/2[n_0bg_b/m_b+3p_0f/m_fn_0f+g_2f8p_0fm_f/ħ^2


    ±√((n_0bg_b/m_b-3p_0f/m_fn_0f-g_2f8p_0fm_f/ħ^2)^2
    +4n_0bn_0fg_bf^2/m_bm_f)],

"+" corresponds to perturbations in system of bosons affected by fermions,
"-" corresponds to perturbations in system of fermions affected by bosons.

Consider sign "+" in front of the square root and drop contribution of fermions
then we get

    V_b^2=n_0bg_b/m_b

velocity for nonlinear perturbations in BEC considered in Refs. <cit.>, <cit.>.
This velocity corresponds to the long-wavelength limit of the Bogoliubov spectrum.
The velocity square V^2 is positive for the repulsive interaction between bosons.
Consider the influence of fermions on the boson solution.
If boson-fermion interaction is small we can expand the square root.
After expansion assuming that the partial velocity of bosons (<ref>) dominates
over the partial velocity of fermions V_f^2=3p_0f/m_fn_0f+8 g_2fp_0fm_f/ħ^2
we find the following expression

    V^2=n_0bg_b/m_b+
    n_0bn_0fg_bf^2/m_bm_f/
    (n_0bg_b/m_b -3p_0f/m_fn_0f -g_2f8p_0fm_f/ħ^2).

General behavior of velocity (<ref>) shows that the chosen solution should have positive second term.
The sign of boson-fermion interaction does not affect the velocity of nonlinear perturbations.
It is corresponds to the general solution (<ref>).

Solution (<ref>) with sign - corresponds to the acoustic wave in fully spin polarized fermions.
Separation on the bosonic and fermionic branches is partially conventional.
If the the partial velocity of fermions V_f dominates over the partial velocity of bosons V_b (<ref>)
we have sign - for bosonic branch (find solution (<ref>) for small boson-fermion interaction) and sign + for fermionic branch.
However, study of the mixture for intermediate boson-fermion interaction does not allow so straightforward separation on bosonic and fermionic branches.
So, we keep studying nonlinear solution corresponding to sign "+" in (<ref>) and we conventionally call it the bosonic branch.
The second branch conventionally called the fermionic branch will be studied elsewhere.


Fermions in partially polarized regime demonstrate two acoustic wave and the spin wave with ω(k=0)≠0 <cit.>.
Hamiltonian of the nonlinear Pauli equation in <cit.> contains
the interaction term corresponding to the total energy of a two-component Fermi gas presented in <cit.>.











 §.§ The second order perturbations


Next order on ε contribution in hydrodynamic equations leads to the following set of nonlinear differential equations:
the continuity equation for bosons

    V∂_τn_1b-V∂_ξn_2b+n_0b∂_ξv_2b+∂_ξ(n_1bv_1b)=0;

the Euler equation for bosons

    m_bn_0bV∂_τv_1b-m_bn_0bV∂_ξv_2b-m_bn_1bV∂_ξv_1b
    -ħ^2/4m_b∂_ξ^3n_1b


    +m_bn_0bv_1bV∂_ξv_1b
    =-g_bn_0b∂_ξn_2b-g_bn_1b∂_ξn_1b
    -g_2bn_0b∂_ξ^3n_1b


    -g_bfn_0b∂_ξn_2f -g_bfn_1b∂_ξn_1f
    -1/2g_2bfn_0b∂_ξ^3n_1f;

the continuity equation for fermions

    V∂_τn_1f-V∂_ξn_2f+n_0f∂_ξv_2f+∂_ξ(n_1fv_1f)=0;

the Euler equation for fermions

    m_fn_0fV∂_τv_1f-m_fn_0fV∂_ξv_2f-m_fVn_1f∂_ξv_1f


    +m_fn_0fv_1f∂_ξv_1f -ħ^2/4m_f∂_ξ^3n_1f+∂_ξp_2^xx


    =-g_2fm^2/2ħ^2∂_ξ[n_0f(3p_2f^xx+p_2f^yy+p_2f^zz)


    +(3p_0f^xx+p_0f^yy+p_0f^zz)n_2f +(3p_1f^xx+p_1f^yy+p_1f^zz)n_1f]


    -g_bfn_0f∂_ξn_2b-g_bfn_1f∂_ξn_1b
    -1/2g_2bfn_0f∂_ξ^3n_1b;

equations for evolution of the elements of pressure tensor

    V∂_τp_1f^xx-V∂_ξp_2f^xx+v_1f∂_ξp_1f^xx
    +3p_0f^xx∂_ξv_2f+3p_1f^xx∂_ξv_1f=0,


    V∂_τp_1f^yy-V∂_ξp_2f^yy+v_1f∂_ξp_1f^yy
    +p_0f^yy∂_ξv_2f+p_1f^yy∂_ξv_1f=0,

and

    V∂_τp_1f^zz-V∂_ξp_2f^zz+v_1f∂_ξp_1f^zz
    +p_0f^zz∂_ξv_2f+p_1f^zz∂_ξv_1f=0.

All functions of the first order can be represented
via the first order perturbations for concentrations of bosons n_1b and fermions n_1f.
Presenting the second order hydrodynamic perturbations via the derivative of second order concentration of fermions ∂_ξn_2f
and first order perturbations for concentrations of bosons n_1b and fermions n_1f
find equation for three variables,
where coefficient in front of ∂_ξn_2f equals to zero
if expression (<ref>) for V^2 is included.




After described manipulations obtain equation for concentrations n_2f, n_1f, n_1b: 

    g_bfn_0f· g_bfn_0b∂_ξn_2f
    +∂_ξn_2f[8g_2fp_0m_f^2/ħ^2+3p_0/n_0f-m_fV^2](m_bV^2-g_bn_0b)


    +(m_bV^2-g_bn_0b)[2m_fV^2(∂_τn_1f+n_1f∂_ξn_1f/n_0f)
    -ħ^2/4m_f∂_ξ^3n_1f
    +n_1f∂_ξn_1f(6p_0/n_0f^2+20g_2fp_0m_f^2/n_0fħ^2)


    +g_bfn_1f∂_ξn_1b+1/2g_2bfn_0f∂_ξ^3n_1b]
    +g_bfn_0f[2m_bV^2∂_τn_1b+2m_bV^2n_1b∂_ξn_1b/n_0b-ħ^2/4m_b∂_ξ^3n_1b


    +g_bn_1b∂_ξn_1b+g_2bn_0b∂_ξ^3n_1b
    +g_bfn_1b∂_ξn_1f+1/2g_2bfn_0b∂_ξ^3n_1f]=0.


The first and second terms in equation (<ref>) contain all contribution of the second order functions
(in this case they are expressed via the second order concentration of fermions n_2f).

The coefficient in front of n_2f goes to zero
if explicit form of velocity V (<ref>) is used.
Therefore, equation (<ref>) reduces to equation relatively two functions n_1b and n_1f.

The first term in equation (<ref>)
and combined inside square brackets group of seven last terms present the contribution of bosons.
Other terms present the contribution of fermions.

Equation (<ref>) appears as the Euler equation for fermions.
So the contribution of bosons equal to zero if the interspecies interaction constant g_bf goes to zero.

The lowest order on ε analysis gives relation between concentrations of bosons and fermions

    n_1f=m_b/n_0bg_bf(V^2-n_0bg_b/m_b)n_1b.

This relation allows to get an equation for single function (n_1b for instance) from (<ref>).
Next, then n_1b is found we obtain the structure of soliton solution for n_1f using relation (<ref>).


Moreover, we can use solution for V^2 given by equation (<ref>) to analyze relation (<ref>).

    n_1f/n_1b=1/2m_b/n_0bg_bf[3p_0/m_fn_0f+g_2f8p_0m_f/ħ^2-n_0bg_b/m_b


    +√((3p_0/m_fn_0f+g_2f8p_0m_f/ħ^2-n_0bg_b/m_b)^2
    +4n_0bn_0fg_bf^2/m_bm_f)].

The right-hand side of equation (<ref>) is the product of two functions:
the interaction constant g_bf
and combination of parameters located in brackets.
The structure of the parameters in brackets can be expressed as follows: Ξ+√(Ξ^2+Λ^2).
The sign of this structure does not depend on signs and values of parameters Ξ, Λ being always positive.
Therefore, perturbations for the bosons and fermions have same sign if boson-fermion interaction is repulsive g_bf>0
or they have opposite signs if boson-fermion interaction is attractive g_bf<0.

Expression (<ref>) can be rewritten in different equivalent form

    n_1f=n_0fg_bf/m_f1/(V^2-3p_0f/m_fn_0f-8g_2fp_0fm_f/ħ^2)n_1b.








 §.§ Korteweg-de Vries equation for perturbations of bosons


Korteweg-de Vries (KdV) equation for concentration of bosons has the following structure

    ã∂_τn_1b+b̃n_1b∂_ξn_1b+c̃∂_ξ^3n_1b=0,

where we find coefficients

    ã =2m_bV^2[1 +
    n_0fg_bf/m_fn_0bg_bf/m_b/ (V^2 -3p_0/m_fn_0f-8p_0fg_2fm_f/ħ^2)^2],



    b̃=1/n_0b[ 2m_bV^2 +g_bn_0b +g_bfn_0bn_0fg_bf/m_f/V^2 -3p_0/m_fn_0f-8p_0fg_2fm_f/ħ^2


    +g_bfn_0bn_0fg_bf/m_fn_0bg_bf/m_f/(V^2 -3p_0/m_fn_0f-8p_0fg_2fm_f/ħ^2)^2


    +m_bn_0f/n_0b(m_f/m_b)^-2(V^2 +3p_0/m_fn_0f +10p_0fg_2fm_f/ħ^2)
    ×


    ×(n_0bg_bf/m_b)^3/(V^2 -3p_0/m_fn_0f-8p_0fg_2fm_f/ħ^2)^3],

and

    c̃=g_2bn_0b -ħ^2/4m_b +g_2bfn_0bn_0fg_bf/m_f/V^2 -3p_0/m_fn_0f-8p_0fg_2fm_f/ħ^2


    -ħ^2/4m_fn_0fg_bf/m_fn_0bg_bf/m_f/(V^2 -3p_0/m_fn_0f-8p_0fg_2fm_f/ħ^2)^2.

Coefficient ã is always positive since V^2 is positive for solution to exist.
However, condition V^2>0 gives a restriction on parameters.
For instance if we drop contribution of fermions V^2=g_bn_0b/m_b.
Hence, the interaction between bosons should be repulsive g_b>0.

Equation (<ref>) can be reduced to single variable after introduction of new variable ζ=τ-Uξ.
Afterwards KdV equation can be integrated.
As the result of integration find nonlinear perturbation of boson concentration in the first order:

    n_1b=3Uã/b̃1/cosh^2(1/2√(Uã/c̃)ζ).

Since coefficient ã is positive, solution (<ref>) can exist if coefficient c̃ is positive.
Sign of coefficient b̃ defines the type of soliton:
bright soliton for b̃>0 or the dark soliton for b̃<0.
















Purely for bosons the coefficient c̃ consists of two terms c̃_b=g_2bn_0b-ħ^2/4m_b^2.
Hence, it is positive for nonzero interaction between bosons in third order by the interaction radius.
Moreover, it requires the repulsive boson-boson interaction g_2b>0.
This conclusion is in agreement with condition V^2>0 requiring g_b>0.
For bosons, coefficient ã simplifies to ã_b=2m_bV^2>0.
Coefficient b̃ also appears in simple form in this limit:
b̃_b=g_b+2m_bV^2/n_0b>0.
The transition to bosons is made by limit g_bf→ 0.
It shows existence of a specific soliton solution with positive amplitude in system of bosons studied up to the TOIR.
Terms proportional to g_bf give contribution of fermions
which is discussed numerically below.

This bright soliton solution purely for bosons is obtained in Ref. <cit.> and <cit.>.
Physical picture behind the bosonic bright soliton solution for repulsive bosons demonstrates
deep relation to the soliton solution experimentally obtained in Ref. <cit.>.





§ NUMERICAL ANALYSIS OF THE BRIGHT-LIKE SOLITON



To perform analysis of the soliton properties and area of its existence represent KdV equation (<ref>)
and velocity of perturbation (<ref>) in dimensionless form.

The dimensionless velocity is given via the mass of bosons and concentration of bosons
since we consider the soliton
which originally exists in the boson subsystem:

    W^2≡m_b^2V^2/ħ^2n_0b^2/3
    =1/2{g+3/5(6π^2)^2/3N_0^2/3/M_0^2(1+8/3F)


    +√((g-3/5(6π^2)^2/3N_0^2/3/M_0^2(1+8/3F))^2 +4N_0/M_0l^2)},

where
N_0=n_0f/n_0b,
M_0=m_f/m_b,
g≡ m_bg_bn_0b^1/3/ħ^2,
G≡ m_bg_2bn_0b/ħ^2,
l≡ m_bg_bfn_0b^1/3/ħ^2,
L≡ m_bg_2bfn_0b/ħ^2,
and
F≡ m_fg_2fn_0f/ħ^2.
The explicit form of the equilibrium pressure for the degenerate fermions
with the full spin polarization is used in form of the Fermi pressure:
p_0f=(6π^2)^2/3ħ^2n_0f^5/3/5m_f.

Dimensionless KdV can be written as follows

    a∂_τ̃N+bN∂_ξ̃N+c∂_ξ̃^3N=0,

where N=n_1b/n_0b, τ̃=√(n_0b)τ, ξ̃=√(n_0b)ξ,

    a=m_bã/ħ^2n_0b^2/3=2W^2[1+N_0/M_0l^2/D_W^2],


    b=m_bb̃/ħ^2n_0b^1/3=2W^2+g
    +N_0/M_0l^2/D_W
    +N_0/M_0^2l^3/D_W^2


    +2N_0/M_0^2l^3/D_W^3[W^2+3/5(6π^2)^2/3N_0^2/3/M_0^2(1+10/3F)],

and

    c=m_bc̃/ħ^2
    =-1/4+G +N_0/M_0Ll/D_W
    -1/4N_0/M_0^3l^2/D_W^2 ,

with

    D_W=W^2-3/5(6π^2)^2/3N_0^2/3/M_0^2[1+8/3F] .



Change of sign of g_bf does not change the velocity of perturbation V (<ref>) and coefficient ã (<ref>).

Change of sign of g_bf does not affect coefficient c̃ (<ref>) since it is also required change of sign of g_2bf.
While coefficient c̃ contains terms without g_bf and g_2bf, term with square of g_bf, and term with product of g_bf on g_2bf.

Change of sign of g_bf modifies coefficient b̃ (<ref>).
If we consider parameters chosen in Fig. (<ref>) find
that for g=1 change of sign of g_bf modifies coefficient b̂ in two times: b̂_+/b̂_-≈2.
However, if parameter g equals to 3 or more change of sign of g_bf gives few percent modification of b̂.


Influence of the boson-boson interaction via g on the form of soliton is demonstrated in Fig. (<ref>).
Existence of the soliton requires relatively large boson-boson interaction to get G>1/4.
However, Fig. (<ref>) shows that increase of boson-boson interaction
(change of g from 1 to 10 at fixed G=1) leads to decrease of amplitude and width of the soliton.
Further increase of g at fixed G=1 leads to decrease of width with no change of amplitude.
However, formal increase of g up to 100 together with increase of G up to G=10 increases the width with no modification of amplitude.

The Fig. (<ref>) is made for relatively small influence of the fermions.
Increase the contribution of fermions and find their contribution in properties of perturbations.


First focus on properties of the velocity.
Dependence of the dimensionless velocity (<ref>) on the boson-boson interaction g
for different boson-boson l and fermion-fermion F interactions is presented in
Figs. (<ref>), (<ref>), (<ref>).
Each figure is made for different rations of masses M_0 and concentrations N_0.

Velocity square W^2 (<ref>) shows almost linear dependence on the boson-boson interaction g.
However, the presence of fermions change this dependence from linear to the superposition of the linear and the square root functions.

Main change of the dependence happens at small boson-boson interaction.
Figs. (<ref>), (<ref>), (<ref>) are made for relatively strong boson-fermion interaction l∼1.
The fermion-fermion interaction is considered in an interval from average F=0.1 to strong F=1 values.
Mass (concentration) increase of each species decreases (increases) the velocity V (<ref>).
This tendency conserves for dimensionless velocity (<ref>),
where increase of the mass (concentration) of fermions relatively the boson mass (concentration) decreases (increases) the velocity W
(compare corresponding lines in Figs. (<ref>), (<ref>), (<ref>)).
Relatively small influence of fermions is demonstrated in Fig. (<ref>),
where mass of fermions is relatively large M_0=9
while variation of F is noticeable at small g.
The increase of l increases the velocity W as it is seen from analytical dependence (<ref>).
The increase of fermion-fermion repulsion gives small increase of the velocity W at small g at fixed l.
Role of fermion-fermion interaction increases if the mass and concentration ratios are getting closer to 1
as it is demonstrated at transition to Figs. (<ref>) and (<ref>).

On coefficients b and c.
Coefficient c is  the symmetric function of l for the fixed α_L.
But coefficient b shows nonsymmetric dependence on boson-fermion interaction l.
The third term in (<ref>) is positive since D_W>0.
However, the fourth and last terms can be negative for attractive boson-fermion interaction.

Hence, the boson-fermion repulsion increases the amplitude of soliton.
The increase can be nonmonotonic
since W^2(l^2) are located in the denominator of the amplitude.

There is a competition between different terms defining the amplitude A=3Ua/bn_0b for the attraction between bosons and fermions.
The fourth and last terms become negative in this regime
while the third term is positive.
Hence, the sign of the amplitude change depends on parameters of the system (see Fig. <ref>).

Area of the soliton existence is restricted by the condition that
the width of soliton D∼√(c)/√(a) is real.
It means that coefficient c should be positive, since coefficient a>0 is positive for all parameters.

There is simple dependence of the width of soliton on G.
It is linear via the second term in c.
Focus on G=1.

The third term in c contains dependence on L.
This is positive term for positive product lL.
Let us to point out that D_W is positive for all parameters.
Since parameters l and L are related introduce the following relation L=α_L l,
where α_L<1 is a parameter
which does not depend on l or other parameters
and represents independent variation of interaction constant L.

Coefficient α_L is an independent parameter.
Hence, if L is fixed at change of l.
It means that parameter α_L changes to compensate contribution of l in L.


If we consider dependence of c on l at fixed L the third term in c (<ref>) plays crucial role
(for instance at m = 4, n = 2, g = 10, F = 0.01, G = 1, L = 0.01).
The dependence numerically appears as almost parabolic dependence (see Fig. <ref>)
in spite more complex analytical dependence via W^2(l).
This parabola has branches going below from the maximum located at positive value of l.

However, l and L are moments of the same potential of boson-fermion interaction.
So we use representation L=α_L l introduced above.
It changes dependence of c on l.
In this case, parameter c is the function of l^2.






Value G=1 is chosen, so the boson-fermion interaction shifts coefficient c from value c_0=0.75.
For small positive l at fixed L=0.01 the shift of c is positive (see Fig. <ref>).
There is value of l=l_0(g),
where the shift becomes equal to zero.
Value of l_0 becomes larger at larger boson-boson interaction g.
At further increase of l above l_0 the shift becomes negative c<0.75.
However, parameter c shows small deviation from G-1/4
and has positive value.
Therefore, presence of fermions does not destroy the soliton solution.

Figs. <ref> and <ref> show that
deviations of c from value G-0.25 are small.
Therefore, small values of G can be chosen down to G_min=0.26.

Consider behavior of c at fixed α_L.
Monotonic increase of c as function of l^2 is found at relatively large α_L=0.01
and relatively large g=10.
Small g=2 at large α_L=0.01 and different g at smaller α_L=0.001 lead to decrease of c at small l
which replaces by the increase of c at larger l.
The area of decrease of c from c_0=0.75 becomes wider and c_min becomes smaller at smaller α_L and smaller g
as it is presented in Fig. <ref>.
All of it is obtained for small fermion-fermion interaction F=0.01.
Area of larger F is presented for the width of soliton D on Fig. <ref>.

Relatively large fermion concentration and large repulsive fermion-fermion interaction can significantly decrease D_W.
So, the contribution of the last negative term in c (<ref>) can increase faster in compare with the third term in coefficient c.



Consider two parts of Fig. <ref>.
If F=1 there is no visible modification of D(l) at different α_L.
If the fermion-fermion repulsion is smaller F=0.3 there is increase of function D(l) with increase of α_L at g=2.
However, the stronger boson-boson repulsion hides any contribution of α_L.




Crucial role in the soliton existence plays the second constant of the boson-boson interaction G>0.25.
First interaction constants for the boson-boson and boson-fermion interactions defines properties of the solution.
The constants of the boson-fermion and fermion-fermion interaction existing in the third order on the interaction radius have small influence
if the boson-boson repulsive interaction is strong g≥1.
Condition G>0.25 also corresponds to this criterium.




§ CONCLUSION


Boson-fermion mixtures have been studied in terms of hydrodynamic model.
Boson-boson, fermion-fermion and boson-fermion interactions have been considered up to the third order by the interaction radius.
A stress has been made on the models of fermions,
where the pressure tensor has been considered as an independent function.
Hence, no equation of state has been used for perturbations of pressure,
but additional hydrodynamic equation for the pressure evolution is derived from the microscopic quantum model.
Equation of state can be used for the equilibrium pressure.


Developed model has been used to study the bright soliton in repulsing Bose-Einstein condensate fraction.
It exists due to the repulsive boson-boson interaction giving positive interaction constant in the third order by the interaction radius.
Formation of soliton in fermion fraction has been found.
It has been obtained that type of soliton of fermion concentration depends
on the sign of boson-fermion interaction constant in the first order by the interaction radius.
Hence, the boson-fermion repulsion (attraction) leads to bright (dark) soliton in fermion fraction.
Influence of the fermions on the properties of soliton in boson fraction is analyzed.


The obtained model contains the first order on the interaction radius including the boson-boson interaction corresponding to the Gross-Pitaevskii equation
and boson-fermion interaction (existing in well-known and sited above works on boson-fermion mixtures)
which are three-dimensional zeroth moments of the interaction potential.
However, the consideration of the interaction terms in the third order by the interaction radius introduces three additional interaction constants
which are the second moments of the interaction potential for boson-boson, boson-fermion, and fermion-fermion interactions.

It is possible to make an estimation of new constants via well-known constants (for boson-boson and boson-fermion interactions) as it is presented in <cit.>.
However, all constants are independent and introduce additional information about interaction potential.
Hence, the experimental study of properties of the found here solitons allows to study the interaction potential in more details.

Moreover, the found solitons presents interest by themselves
since they are examples of new nonlinear phenomena in ultracold mixtures.








§ ACKNOWLEDGEMENTS

The work of P.A. and M.T. is supported by the Russian Foundation for Basic Research (grant no. 20-02-00476).




17

Yefsah Nat 13 T. Yefsah, A. T. Sommer, M. J. H. Ku, L. W. Cheuk,
W. Ji, W. S. Bakr, and M. W. Zwierlein, Nature 499, 426 (2013).

Becker Nat Phys 08 C. Becker, S. Stellmer, P. Soltan-Panahi, S. Drscher, M. Baumert, E.-M. Richter, J. Kronjger, K. Bongs, K. Sengstock,
Nature Phys. 4, 496 (2008).

Tylutki NJP 16 M. Tylutki, A. Recati, F. Dalfovo, and S. Stringari, New J. Phys. 18, 053014 (2016).

Anderson PRL 01 B. P. Anderson, P. C. Haljan, C. A. Regal, D. L. Feder, L. A. Collins, C. W. Clark, and E. A. Cornell,
Phys. Rev. Lett. 86, 2926 (2001).


Guenther PRA 17 N.-E. Guenther, P. Massignan, and A. L. Fetter, Phys. Rev. A 96, 063608 (2017).

Gautam PRA 14 S. Gautam, A. Roy, and S. Mukerjee, Phys. Rev. A 89, 013612 (2014).



Lee Sc Ad 18 W. Lee, A. H. Gheorghe, K. Tiurev, T. Ollikainen, M. Mttnen and D. S. Hall, Science Advances 4, 3820 (2018).



Wachtler 1601 F. Wachtler and L. Santos, Phys. Rev. A 93, 061603R (2016).




Bisset PRA 16 R. N. Bisset, R. M. Wilson, D. Baillie, P. B. Blakie,
Phys. Rev. A 94, 033619 (2016).




Wachtler 1605 F. Wachtler and L. Santos, Phys. Rev. A 94, 043618 (2016).




Karpiuk PRL 04 T. Karpiuk, M. Brewczyk, S. Ospelkaus-Schwarzer, K. Bongs, M. Gajda, and K. Rzazewski, Phys. Rev. Lett. 93, 100401 (2004).


Truscott Science 01 A. G. Truscott, K. E. Strecker, W. I. McAlexander,
G. B. Partridge, and R. G. Hulet, Science 291, 2570 (2001).

Schreck PRL 02 F. Schreck, L. Khaykovich, K. L. Corwin,
G. Ferrari, T. Bourdel, J. Cubizolles, and C. Salomon, Phys. Rev. Lett. 88, 160401 (2002).

Hadzibabic PRL 02 Z. Hadzibabic, C. A. Stan, K. Dieckmann, S. Gupta, M. W. Zwierlein,
A. Gorlitz, and W. Ketterle, Phys. Rev. Lett. 88, 160401 (2002).

Roati PRL 02 G. Roati, F. Riboli, G. Modugno, and M. Inguscio, Phys. Rev. Lett. 89, 150403 (2002).


Dalfovo RMP 99 F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999).


Karpiuk JPB 02 T. Karpiuk, M. Brewczyk and K. Rzazewski, J. Phys. B: At. Mol. Opt. Phys. 35, L315 (2002).

Andreev MPL B 12 P. A. Andreev, L. S. Kuzmenkov, Mod. Phys. Lett. B 26, 1250152 (2012).


Andreev LP 19 P. A. Andreev, Laser Phys. 29, 035502 (2019).

Adhikari PRA 07 S. K. Adhikari and L. Salasnich, Phys. Rev. A 76, 023612 (2007).

Parker PRA 12 N. G. Parker, D. A. Smith, Phys. Rev. A 85, 013604 (2012).


Cetoli PRA 13 A. Cetoli, J. Brand, R. G. Scott, F. Dalfovo, and L. P. Pitaevskii, Phys. Rev. A 88, 043639 (2013).



Sanner PRL 12 C. Sanner, E. J. Su,W. Huang, A. Keshet, J. Gillen, and W. Ketterle, Phys. Rev. Lett. 108, 240404 (2012).

Valtolina Nat Ph 17 G. Valtolina, F. Scazza, A. Amico, A. Burchianti, A. Recati, T. Enss,
M. Inguscio, M. Zaccanti, and G. Roati, Nature Phys. 13, 704 (2017).





Ravensbergen PRA 18 C. Ravensbergen, V. Corre, E. Soave, M. Kreyer, E. Kirilov, and R. Grimm, Phys. Rev. A 98, 063624 (2018).

Matveev PRL 17 K. A. Matveev, and A. V. Andreev, Phys. Rev. Lett. 119, 266801 (2017).




Matveev PRB 12 K. A. Matveev, A. V. Andreev, Phys. Rev. B 85, 041102 (2012).



Luttinger JMP 63 J. M. Luttinger, J. Math. Phys. 4, 1154 (1963).


Haldane JPC 81 F. D. M. Haldane, J. Phys. C: Solid State Phys. 14, 2585 (1981).



Sidorenkov Nat 13
L. A. Sidorenkov, M. K. Tey, R. Grimm, Y.-H. Hou, L. Pitaevskii, and S. Stringari, Nature 498, 78 (2013).



DeSalvo PRL 17 B. J. DeSalvo, K. Patel, J. Johansen, and C. Chin, Phys. Rev. Lett. 119, 233401 (2017).





Rakshit SPP 19 D. Rakshit, T. Karpiuk, M. Brewczyk, M. Gajda, Sci. Post. Phys. 6, 079 (2019).

Lu PRA 13 Zhen-Kai Lu, S. I. Matveenko, and G. V. Shlyapnikov, Phys. Rev. A 88, 033625 (2013).



Csire PRA 10 G. Csire, D. Schumayer, and B. Apagyi, Phys. Rev. A 82, 063608 (2010).

Hoefer PRA 11 M. A. Hoefer, J. J. Chang, C. Hamner, and P. Engels, Phys. Rev. A 84, 041605(R) (2011).


Zezyulin EPJ D 13 K. V. Zezyulin, P. A. Andreev, and L. S. Kuzmenkov, Eur. Phys. J. D 67, 140 (2013).



Katsimiga NJP 17 01 G. C. Katsimiga, G. M. Koutentakis, S. I. Mistakidis, P. G. Kevrekidis,
and P. Schmelcher, New J. Phys. 19, 073004 (2017).


Katsimiga NJP 17 02 G. C. Katsimiga, S. I. Mistakidis, G. M. Koutentakis, P. G. Kevrekidis
and P. Schmelcher, New J. Phys. 19, 123012 (2017).


Katsimiga PRA 18 G. C. Katsimiga, S. I. Mistakidis, G. M. Koutentakis, P. G. Kevrekidis,
and P. Schmelcher, Phys. Rev. A 98, 013632 (2018).


Mistakidis NJP 18 S. I. Mistakidis, G. C. Katsimiga, P. G. Kevrekidis, and P. Schmelcher, New J. Phys. 20, 043052 (2018).




Giorgini RMP 08 S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).



MaksimovTMP 2001 L. S. Kuz'menkov, S. G. Maksimov, and V. V. Fedoseev, Theor.
Math. Fiz. 126, 136 (2001) [Theoretical and Mathematical
Physics 126, 110 (2001)].


Andreev PRA08 P. A. Andreev, L. S. Kuz'menkov, Phys. Rev. A 78, 053624 (2008).

Andreev 2001 P. A. Andreev, arXiv:2001.02764.







Ruggiero arxiv 19 P. Ruggiero, P. Calabrese, B. Doyon, J. Dubail, arXiv:1910.00570


Bertini PRL 16 B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Phys. Rev. Lett. 117, 207201 (2016).


Adhikari JPB05 S. K. Adhikari, Journal of Physics B 38, 3607 (2005).


Rosanov N. N. Rosanov, A. G. Vladimirov, D. V. Skryabin, W. J. Firth, Phys. Lett. A. 293, 45 (2002).

Braaten E. Braaten, H.-W. Hammer, and Shawn Hermans, Phys. Rev. A. 63, 063609 (2001).

Andreev 1912 P. A. Andreev, arXiv:1912.00843.




Washimi PRL 66 H. Washimi and T. Taniuti, Phys. Rev. Lett 17, 996 (1966).

Kalita PlasmaPhys 98 B. C. Kalita, R. Das, Plasma Phys. 5, 3588 (1998).

Andreev LPL 18 P. A. Andreev, Laser Phys. Lett. 15, 105501 (2018).

Jo Sc 09 G.-B. Jo, Y.-R. Lee, J.-H. Choi, C. A. Christensen, T. H. Kim, J. H. Thywissen, D. E. Pritchard,
and W. Ketterle, Science 325, 1521 (2009).


Wang NJP 14 Z. Wang, M. Cherkasskii, B. A Kalinikos, L. D. Carr, M. Wu, New J. Physics 16, 053048 (2014).







