The variational reduction for low-dimensional Fermi gases and Bose - Fermi mixtures: A brief review

TL;DR

This paper reviews recent theoretical advances in low-dimensional Fermi and Bose-Fermi gases, focusing on dimensional reduction techniques and their accuracy in modeling ground states and dynamics within a quasi-mean-field framework.

Contribution

It introduces a variational reduction method for low-dimensional Fermi and Bose-Fermi gases and compares reduced models with full 3D equations, demonstrating high accuracy.

Findings

Good agreement between reduced and full models for ground states

Effective 1D and 2D equations capture complex dynamics

Dimensional reduction simplifies analysis of degenerate gases

Abstract

We present a summary of some recent theoretical results for matter-wave patterns in Fermi and Bose-Fermi degenerate gases, obtained in the framework of the quasi-mean-field approximation. We perform a dimensional reduction from the three-dimensional (3D) equations of motion to 2D and 1D effective equations. In both cases, comparison of the low-dimensional reductions to the full model is performed, showing very good agreement for ground-state solutions. Some complex dynamical regimes are reported too for the corresponding 1D systems.

Tables

Table 1. Table 1: λ 1 subscript 𝜆 1 \lambda_{1} , λ 2 subscript 𝜆 2 \lambda_{2} , β 𝛽 \beta , and s F subscript 𝑠 F s_{\mathrm{F}} for three different regimes in the Fermi-gas dynamics [ 60 , 61 , 62 , 63 ] .

Regime	${\lambda }_1$	${\lambda }_2$	$\beta$	$s_F$
Polarized	1	1	1	0
BCS	2	4	1	1/2
Unitary	2	4	0.44	1/2

Equations

$$\mathcal{S}=\int dtd\mathbf{r}\mathcal{L},$$

$$\mathcal{L}=\frac{i\hslash }{2{\lambda }_1}\left({\Psi }^{\ast }\frac{\partial \Psi }{\partial t}-\Psi \frac{\partial {\Psi }^{\ast }}{\partial t}\right)-\frac{{\hslash }^2}{2{\lambda }_2{m}_{\mathrm{F}}}{\left|\nabla \Psi \right|}^2-U(\mathbf{r}){\left|\Psi \right|}^2-\frac{3}{5}\beta \frac{{\hslash }^2}{2m_{\mathrm{F}}}{C}_{\mathrm{F}}{\left|\Psi \right|}^{10/3}-\frac{g_{\mathrm{F}}}{2}{\left|\Psi \right|}^4,$$

$$\frac{i\hslash }{{\lambda }_1}\frac{\partial \Psi }{\partial t}=[-\frac{{\hslash }^2}{2{\lambda }_2{m}_{\mathrm{F}}}{\nabla }^2+U+g_{\mathrm{F}}{\left|\Psi \right|}^2+\frac{{\hslash }^2}{2m_{\mathrm{F}}}\beta C_{\mathrm{F}}{\left|\Psi \right|}^{4/3}]\Psi ,$$

$$U\left(\mathbf{r}\right)=\frac{1}{2}{m}_{\mathrm{F}}{\omega }_z^2{z}^2+U_{2D}\left({\mathbf{r}}_{\perp },t\right).$$

$$\Psi \left(\mathbf{r},t\right)=\frac{1}{{\pi }^{1/4}\sqrt{\xi ({\mathbf{r}}_{\perp },t)}}\exp \left(-\frac{z^2}{2(\xi ({\mathbf{r}}_{\perp },t){)}^2}\right)\varphi \left({\mathbf{r}}_{\perp },t\right).$$

$$\mathcal{S}{}_{2D}=\int dtdxdy{\mathcal{L}}_{2D},$$

$${\mathcal{L}}_{2D}$$

$$i\frac{\hslash }{{\lambda }_1}{\partial }_t\varphi$$

$$m_{\mathrm{F}}{\omega }_z^2{\xi }^4-\frac{2{\hslash }^2}{5m_{\mathrm{F}}}\beta C_{2\mathrm{D}}{\left|\varphi \right|}^{4/3}{\xi }^{4/3}-\frac{g_{\mathrm{F}}}{\sqrt{2\pi }}{\left|\varphi \right|}^2\xi -\frac{{\hslash }^2}{{\lambda }_2{m}_{\mathrm{F}}}=0.$$

$$U_{2D}=A\sum _{b=1}^2\sum _{a=1}^3{\sin }^2\left(\frac{2\pi }{\lambda }({\mathbf{k}}_{a,b}\cdot {\mathbf{r}}_{\perp })\right)+\frac{1}{2}{\omega }_x^2{x}^2+\frac{1}{2}{\omega }_y^2{y}^2,$$

$$U\left(\mathbf{r}\right)=U_{1\mathrm{D}}\left(z\right)+\frac{1}{2}m{\omega }_t^2{r}^2,$$

$$\Psi \left(\mathbf{r},t\right)=\frac{1}{{\pi }^{1/2}\sigma \left(z,t\right)}\exp \left(-\frac{r^2}{2\sigma {\left(z,t\right)}^2}\right)f\left(z,t\right),$$

$$i\frac{\hslash }{{\lambda }_1}{\partial }_{t}f=\left[-\frac{{\hslash }^2}{2{\lambda }_2{m}_{\mathrm{F}}}{\partial }_z^2+U_{1\mathrm{D}}+\frac{g_{\mathrm{F}}}{2\pi {\sigma }^2}{\left|f\right|}^2+\frac{{\hslash }^2}{2m_{\mathrm{F}}}\beta \frac{C_{1\mathrm{D}}}{{\sigma }^{4/3}}{\left|f\right|}^{4/3}+\frac{{\hslash }^2}{2m_{\mathrm{F}}{\lambda }_2{\sigma }^2}\frac{1}{2}{m}_{\mathrm{F}}{\omega }_t^2{\sigma }^2\right]f,$$

$$m_{\mathrm{F}}{\omega }_t^2{\sigma }^4-\frac{2}{5}\frac{{\hslash }^2}{m_{\mathrm{F}}}\beta C_{1\mathrm{D}}|f{|}^{4/3}{\sigma }^{2/3}-\frac{{\hslash }^2}{{\lambda }_2{m}_{\mathrm{F}}}-\frac{g_{\mathrm{F}}}{2\pi }|f{|}^2=0,$$

$$f_R\left(z,t_0\right)=\frac{1}{2}+\frac{1}{2}\sum _{j=1}^{N_s/2}\left[{\left(-1\right)}^{1+j}\left(\tanh \left(\frac{{z-z}_j}{{\Delta }_s}\right)+1\right)+{\left(-1\right)}^j\left(\tanh \left(\frac{z-z_{-j}}{{\Delta }_s}\right)+1\right)\right],$$

$$E_{ϵ}=\sum _{j=1}^{N_s}{v}_{j,ϵ}^2.$$

$$\mathcal{S}=\int dtd\mathbf{r}\left({\mathcal{L}}_{\mathrm{B}}+{\mathcal{L}}_{\mathrm{F}}+{\mathcal{L}}_{\mathrm{BF}}\right),$$

$${\mathcal{L}}_B=\frac{i\hslash }{2}\left({\Psi }_{\mathrm{B}}^{\ast }\frac{\partial {\Psi }_{\mathrm{B}}}{\partial t}-{\Psi }_{\mathrm{B}}\frac{\partial {\Psi }_{\mathrm{B}}^{\ast }}{\partial t}\right)-\frac{{\hslash }^2}{2m_{\mathrm{B}}}{\left|\nabla {\Psi }_{\mathrm{B}}\right|}^2-U_B(\mathbf{r}){\left|{\Psi }_{\mathrm{B}}\right|}^2-\frac{1}{2}{g}_{\mathrm{B}}{\left|{\Psi }_{\mathrm{B}}\right|}^4,$$

$${\mathcal{L}}_{\mathrm{F}}=\frac{i\hslash }{2{\lambda }_1}\left({\Psi }_{\mathrm{F}}^{\ast }\frac{\partial {\Psi }_{\mathrm{F}}}{\partial t}-{\Psi }_{\mathrm{F}}\frac{\partial {\Psi }_{\mathrm{F}}^{\ast }}{\partial t}\right)-\frac{{\hslash }^2}{2{\lambda }_2{m}_{\mathrm{F}}}{\left|\nabla {\Psi }_{\mathrm{F}}\right|}^2-U_{\mathrm{F}}(\mathbf{r}){\left|{\Psi }_{\mathrm{F}}\right|}^2-\frac{1}{2}{g}_{\mathrm{F}}{\left|{\Psi }_{\mathrm{F}}\right|}^4-\frac{3\beta C_{\mathrm{F}}{\hslash }^2}{10m_{\mathrm{F}}}{\left|{\Psi }_{\mathrm{F}}\right|}^{10/3},$$

$${\mathcal{L}}_{B\mathrm{F}}=-\frac{1}{2}{g}_{\mathrm{BF}}{\left|{\Psi }_B\right|}^2{\left|{\Psi }_{\mathrm{F}}\right|}^2.$$

$$i\hslash {\partial }_t{\Psi }_B=\left[-\frac{{\hslash }^2}{2m_B}{\nabla }^2+g_B{\left|{\Psi }_B\right|}^2+g_{\mathrm{BF}}{\left|{\Psi }_{\mathrm{F}}\right|}^2+U_B\right]{\Psi }_B,$$

$$\frac{i\hslash }{{\lambda }_1}{\partial }_t{\Psi }_{\mathrm{F}}=\left[-\frac{{\hslash }^2}{2{\lambda }_2{m}_{\mathrm{F}}}{\nabla }^2+g_{\mathrm{F}}{\left|{\Psi }_{\mathrm{F}}\right|}^2+g_{\mathrm{BF}}{\left|{\Psi }_B\right|}^2+U_{\mathrm{F}}+\frac{\beta C_{\mathrm{F}}{\hslash }^2}{2m_F}{\left|{\Psi }_{\mathrm{F}}\right|}^{4/3}\right]{\Psi }_{\mathrm{F}}.$$

$$C_s\left({\bar{n}}_B,{\bar{n}}_{\mathrm{F}}\right)=\frac{⟨{\bar{n}}_B{\bar{n}}_{\mathrm{F}}⟩}{\sqrt{⟨{\bar{n}}_B^2⟩⟨{\bar{n}}_{\mathrm{F}}^2⟩}},$$

$$U_{\mathrm{B}/\mathrm{F}}\left(\mathbf{r}\right)=\frac{1}{2}{m}_{\mathrm{B}/\mathrm{F}}{\omega }_{z,\mathrm{B}/\mathrm{F}}^2{z}^2+U_{2\mathrm{D},\mathrm{B}/\mathrm{F}}\left({\mathbf{r}}_{\perp }\right),$$

$${\Psi }_{\mathrm{B}/\mathrm{F}}\left(\mathbf{r},t\right)=\frac{1}{{\pi }^{1/4}\sqrt{{\xi }_{\mathrm{B}/\mathrm{F}}({\mathbf{r}}_{\perp },t)}}\exp \left(-\frac{z^2}{2{\left({\xi }_{\mathrm{B}/\mathrm{F}}({\mathbf{r}}_{\perp },t)\right)}^2}\right){\varphi }_{\mathrm{B}/\mathrm{F}}({\mathbf{r}}_{\perp },t)$$

$$\mathcal{S}=\int dtdxdy\left({\mathcal{L}}_{2D,\mathrm{B}}+{\mathcal{L}}_{2\mathrm{D},\mathrm{F}}+{\mathcal{L}}_{2\mathrm{D},\mathrm{BF}}\right),$$

$${\mathcal{L}}_{2\mathrm{D},\mathrm{B}}=i\frac{\hslash }{2}\left({\varphi }_{\mathrm{B}}^{\ast }{\partial }_t{\varphi }_{\mathrm{B}B}-{\varphi }_{\mathrm{B}}{\partial }_t{\varphi }_{\mathrm{B}}^{\ast }\right)-U_{2D,\mathrm{B}}{n}_{2\mathrm{D},\mathrm{B}}-e_{2\mathrm{D},\mathrm{B}},$$

$${\mathcal{L}}_{2\mathrm{D},\mathrm{F}}=i\frac{\hslash }{2{\lambda }_1}\left({\varphi }_{\mathrm{F}}^{\ast }{\partial }_t{\varphi }_{\mathrm{F}}-{\varphi }_{\mathrm{F}}{\partial }_t{\varphi }_{\mathrm{F}}^{\ast }\right)-U_{2\mathrm{D},\mathrm{F}}{n}_{2\mathrm{D},\mathrm{F}}-e_{2\mathrm{D},\mathrm{F}},$$

$${\mathcal{L}}_{2\mathrm{D},\mathrm{BF}}=-\frac{1}{{\pi }^{1/2}}\frac{g_{B\mathrm{F}}}{\sqrt{{\xi }_{\mathrm{B}}^2+{\xi }_{\mathrm{F}}^2}}{n}_{2\mathrm{D},\mathrm{B}}{n}_{2\mathrm{D},\mathrm{F}},$$

$$e_{2\mathrm{D},\mathrm{B}}=\frac{{\hslash }^2}{2m_{\mathrm{B}}}{\left|{\nabla }_{\perp }{\varphi }_{\mathrm{B}}\right|}^2+[\frac{g_{\mathrm{B}}}{\sqrt{8\pi }{\xi }_{\mathrm{B}}}{n}_{2\mathrm{D},\mathrm{B}}+\frac{{\hslash }^2}{4m_{\mathrm{B}}{\xi }_{\mathrm{B}}^2}+\frac{1}{4}{m}_{\mathrm{B}}{\omega }_{z,\mathrm{B}}^2{\xi }_{\mathrm{B}}^2]n_{2\mathrm{D},\mathrm{B}},$$

Full text

††journal: Condensed Matter

1 Introduction

Ultracold atomic gases have been widely explored from both experimental and theoretical point of view due to their ability to emulate many effects from condensed-matter physics and create novel states of quantum matter. Various results obtained in this area have been reviewed in many publications – see, in particular, Refs. [1, 2, 3, 4, 5, 6, 7, 8]. Important experimental tools, the application of which opens ways to the observation diverse novel phenomena in the quantum gases, are, inter alia, optical-lattice (OL) potentials, the use of the Feshbach resonance (FR) to control the strength of interactions between atoms, and the implementation of the effective spin-orbit coupling [9, 10, 11, 12, 13].

The effective spatial dimension of the setting in which quantum gases are created strongly affects the ensuing physics. The use of confining potentials makes it possible to reduce the dimension from 3D to 2D and 1D. In particular, the dimensional reduction of confined Bose gases can be approximated by means of the variational method [14, 15, 16, 17]. Recently, similar approaches for ultracold Fermi gases in confining potentials have been elaborated in Refs. [18, 19, 20, 21]. These reductions make it possible to study the complex dynamics and pattern formation in ultracold gases in 2D and 1D settings. In this context, the study of dark solitons in ultracold gases was reported in Bose-Einstein condensates (BECs) [22], and further developed later [23, 24]. For dark solitons in Fermi gases, several works have reported theoretical and experimental results [26, 27, 28, 29, 30, 31, 32]. The reduced 1D equation for Fermi gases was used for studies of interactions between dark solitons [21].

The earliest experimental studies of Bose-Fermi mixtures (BFMs) were performed with lithium isotopes [33, 34], as well as in 174Yb-6Li [35] and 87Rb-40K [36] settings. Much interest has been also drawn to heavy-atom mixtures, such as 87Sr-84Sr [37]. These isotopes, which are characterized by a large nuclear spin, have been proposed for the design of prototype quantum-information processing devices. The use of FRs in the mixtures plays a major role, as it allows one to control nonlinear interactions between the species. For the 87Rb-40K mixture, the FR has been observed in Ref. [38, 39], and a giant FR effect was reported in the 85Rb-6Li system [40]. Further, in the 6Li-133Cs mixture five FRs have been observed [41], and over 30 resonances are expected to be available in the 23Na-40K system [42]. Multiple heteronuclear FRs were reported in the triply quantum-degenerate mixture of bosonic 41K and two fermionic species, 40K and 6Li [43]. In a recently elaborated framework, the BFM is transformed into a strongly interacting isotopic mixture immersed into a Fermi sea, with the help of a wide $s$-wave resonance for the 41K-40K combination. Many theoretical works have addressed the dynamics of BFMs under various conditions [44, 45, 46, 47, 48, 49, 50]. To describe the ground state (GS) of the mixture, the quasi-mean-field theory may be a useful approach [51, 52, 53, 54, 55]. In this framework, the use of FRs was studied in 23Na-6Li, 87Rb-40K, 87Rb-6Li, 3He-4He, 173Yb-174Yb, and 87Sr-84Sr mixtures [56, 57]. Recently, effective 1D and 2D nonlinear Schrödinger equations have been derived for BFMs in cigar-shaped and disc-shaped configurations [58], using the variational approximation (VA) along the lines previously developed in Refs. [14, 21]. In addition, dark solitons in BFMs have been analyzed in Ref. [59]. Here, we address, in particular, dark solitons in the 7Li-6Li BFM, using the effective low-dimensional equations derived in Ref. [58].

The general aim of the present article is to present a brief review of the spatial reduction for Fermi gases and BFMs, based on the VA. In particular, we outline the procedure for implementing the 2D and 1D reduction, starting from the full 3D equations of motions. To test the accuracy of the approximations, we present a comparison of the results with full 3D numerical simulations. Using the corresponding effective equations, we address various dynamical settings, such as dark solitons and their interactions. In the case of BFMs, we consider the construction of GSs, varying the interaction strength. Finally, for the 1D situation, we address the formation of dark solitons in the mixture, and compare the corresponding 1D solution to results of the full numerical simulations, observing good agreement between them. The presentation is arranged as follows: the Fermi gases and BFMs are considered, severally,in Secs. 2 and 3, and the paper is concluded in Sec 4.

2 The Fermi Gas

We consider a dilute superfluid formed by $N$ fermionic atoms of mass $m_{\mathrm{F}}$ and spin $s_{\mathrm{F}}$, loaded into an optical trap at zero temeprature. We apply the local density approximation [4] to the description of this setting. The corresponding dynamical equations can be derived from the action functional

[TABLE]

where the Lagrangian density is

[TABLE]

$\Psi\left(\mathbf{r},t\right)$ being a complex order parameter, whose norm is equal to the number of particles. Here ${C_{\mathrm{F}}}={\left[{6{\pi^{2}}/\left({2{s_{\mathrm{F}}}+1}\right)}\right]^{2/3}}$ is a constant that depends on spin $s_{\mathrm{F}}$, ${g_{\mathrm{F}}}=4\pi{\hbar^{2}}({a_{\mathrm{F}}}/{m_{\mathrm{F}}})[2s_{\mathrm{F}}/(2s_{\mathrm{F}}+1)]$ with scattering length $a_{\mathrm{F}}$, which determines interactions of fermions belonging to different spin states (the interactions which are not forbidden by the Pauli principle) [3], and ${U(\mathbf{r})}$ is to an external potential applied to fermions.

Parameters $\lambda_{1}$, $\lambda_{2}$, $\beta$, and $s_{\mathrm{F}}$ in Eq. (2) correspond to three different regimes addressed in this article, which are listed in Table 1. It is relevant to mention that the spin polarization may affect some parameters, such as coefficient ${C_{\mathrm{F}}}$ [64].

Lagrangian density (2) gives rise to the following The Euler-Lagrange equation,

[TABLE]

which as an effective quasi-mean-field equation for the fermi gas under the consideration; note that it may be rewritten in the form of hydrodynamic equations [65, 66]. More details on the derivation of this equation are given in Appendix A. Below, we focus on the BCS (Bardeen-Cooper-Schrieffer) setting, referring to atoms of 6Li with mass $6$ a.u.

In numerical simulations we use the fourth-order Runge-Kutta method in time, and the centered second-order finite-difference method for handling the spatial discretization. In the next two subsections we reduce the full 3D equation to the corresponding 2D and 1D effective equations, using the VA proposed in Ref. [21].

2.1 The two-dimensional reduction

We derive effective 2D equations, applying the VA to the Fermi gas in the disk-shaped trap. For this purpose, we consider an external potential composed of two terms: the parabolic (harmonic-oscillator) one accounting for the confinement in the $z$ direction, transverse to the disk’s plane, and the in-plane potential, $U_{\mathrm{2D}}$:

[TABLE]

The initial ansatz assumes, as usual, the factorization of the 3D wave function into a product of functions of $z$ and $\mathbf{r}_{\bot}$, the former one being the Gaussian ground state of the harmonic-oscillator potential [14]:

[TABLE]

The Gaussian is subject to the unitary normalization, with transverse width $\xi$ considered as a variational parameter, while the 2D wave function, $\phi$, is normalized to the number of atoms. Therefore, the reduction from 3D to 2D implies that the system of equations should be derived for the pair of functions $\phi\left(\mathbf{r}_{\bot},t\right)$ and $\xi\left(\mathbf{r}_{\bot},t\right)$, using the reduced action functional, which is obtained by integrating the 3D action over the $z$-coordinate:

[TABLE]

where the respective Lagrangian density is

[TABLE]

$C_{\mathrm{2D}}\equiv{(3/5)^{1/2}}{(6/(2s_{\mathrm{F}}+1))^{2/3}}\pi$, the last two terms being produced by the reduction to 2D, the penultimate term corresponding to the spread in the confined dimension. Hence, the Euler-Lagrange equations, derived by varying the 2D action, which is generated by Lagrangian (7), with respect to $\phi$ and $\xi$ take the form of

[TABLE]

[TABLE]

Algebraic equation (9) for $\xi$ cannot be solved analytically, therefore we used the Newton’s method to solve it numerically. The necessity to find $\xi$ at each step of the integration is a numerical complication of a minimal cost compared to the 3D integration of the underlying equation (3). Note that a further simplifications can be achieved by assuming in Eq. (5) that the Gaussian width is a constant $\xi(\mathbf{r}_{\bot},t)=\xi_{0}$. In this case $\xi$, naturally, does not depend on $\phi$. Then, the solution of Eq. (9) with the density tending to zero can be calculated analytically and it is given by $\xi_{0}=\lambda_{2}^{-1/4}\sqrt{\hbar/m_{\mathrm{F}}\omega_{z}}$.

We consider a 2D potential consisting of the axisymmetric parabolic potential and the superposition of two triangular OLs:

[TABLE]

where $\{\mathbf{k}_{a,1}\}$ and $\{\mathbf{k}_{a,2}\}$ are triplets of unitary vectors of both triangular lattices, which are separated by a specific angle $\theta$. Here $A$ denotes the lattice’s amplitude, and $(\omega_{x},\omega_{y})$ are frequencies of the magnetic-optical trapping potential. In the absence of the OLs ($A=0$), we have verified the accuracy of the 2D reduction by comparing results generated by this approximation to those obtained by integrating the underlying 3D equation (3). The respective GS was found by means of the imaginary-time integration based on the fourth-order Runge-Kutta algorithm with $\Delta t=0.5$ $\mu$s. The spatial discretization for the simulations was performed with $\Delta x=0.25$ $\mu$m and $\Delta y=0.25$ $\mu$m. The comparison is displayed in panel (a) of Figure 1, where the radial-density profiles are plotted. We can observe excellent agreement between the reduced 2D and full 3D descriptions. This result suggests one to use the Eqs. (8) and (9) for studying 2D patterns. Panel (b) of Figure 1 shows a comparison of 3D full numerical simulations versus the VA, assuming a constant width $\xi_{0}$. One can observe that the latter approximation produces less accurate results, which is at least ten times worse than the VA with a density-dependent width.

Figure 2 shows the density as a function of coordinates $x$ and $y$ when the OLs are taken into account. We observe that this particular combination of the OLs (a superlattice) produces a pattern in the form of the superstructure, with the number of density peaks varying when the angle between the unitary vectors increases. Note that the multitude of different coexisting robust multi-peak patterns suggests that this setting has a potential for the use as a data-storage system.

2.2 The one-dimensional reduction

Next, we consider the system confined in one dimension, which implies a cigar-shaped configuration, elongated in the $z$ direction. In this case, the potential trap acting in transeverse plane is the harmonic oscillator in the transverse plane:

[TABLE]

where $r^{2}={x^{2}}+{y^{2}}$. It is assumed that the potential in the transverse direction is much stronger than the axial one. The simplest option is to adopt a Gaussian shape in the transverse plane, which represents the ground state of a the 2D harmonic oscillator , similar to what is adopted above in the case of the 2D reduction. As a result, the variable-separation assumption can be applied, defining the 3D wave function as [14, 15, 16]

[TABLE]

where $f$ is normalized to $N$, such that the 1D density is ${n_{1\mathrm{D}}}={\left|f\right|^{2}}$. Here $\sigma$ is the Gaussian width, which is a function of $z$ and time. After some algebra, similar to that performed above, one derives the Euler-Lagrange equations:

[TABLE]

[TABLE]

where $C_{1\mathrm{D}}=(3/5){(6\pi(2{s_{\mathrm{F}}}+1))^{2/3}}$. Similar to the 2D case, algebraic equation (14) is solved using the Newton’s method, and here too the quasi-BCS regime is addressed. We set $U_{1\mathrm{D}}=0$, $a_{\mathrm{F}}=-5$ nm, and $\omega_{t}=1000$ Hz, these parameters being in the range of experimental values [29]. Since Eqs. (13) and (14) produce results which agree well with the full 3D simulations [21], one can use the effective 1D equations to study more complex dynamical behavior, such as that of dark solitons [29, 32].

To generate a dark soliton, it is posible to consider the initial condition with zero imaginary part, $f_{I}(z,t=0)=0$, while the real part is given by $f_{R}(z,t=0)=f_{b}\tanh((z-z_{s})/\Delta_{s})$, where $f_{b}$ and $\Delta_{s}$ are the soliton’s amplitude and width, respectively. We have found that values of the square amplitude and width, $n_{b}=10$ particles/$\mu$m and $\Delta_{s}=0.8\leavevmode\nobreak\ \mu$m, respectively, can be chosen to minimize the background noise. If we consider a set of $N_{s}$ dark solitons, the initial condition for the imaginary part is again zero, $f_{I}=0$, while the real part can be cast in the form of

[TABLE]

where the positions of the solitons are $z_{j}$ and $z_{-j}$ on the positive and negative $z$ half-axes, receptively. Moreover, the widths of the solitons ($\Delta_{s}$) are considered the same, and that the number of initial solitons $N_{s}$ is even. This initial ansatz was normalized to secure the correct density of the wave function, $n_{b}=|f_{b}|^{2}$. Then, the system was simulated with the help of the standard fourth-order Runge-Kutta method with $\Delta t=0.095$ $\mu$s. The spatial discretization for the simulations was performed with $\Delta z=0.100$ $\mu$m. Figure 3 shows the shape of the initial conditions for the case of one and $N_{s}=18$ dark solitons.

In the case of two solitons, we have $z_{i}=z_{-i}=d/2$, where $d$ is the initial inter-core separation. Frame (a) of Fig. 4 shows the spatiotemporal diagram for two solitons at $d=4\mu$m. One clearly observes that both solitons separate in the course of the evolution. Frame (b) of Fig. 4 shows the speed taken after $90$ms of the evaluation as a function of different initial inter-core separations between the dark solitons. naturally, smaller initial separations generate higher speeds. In fact, at this fix time the speed follows the law $v_{s}\sim d^{\alpha}$, with $\alpha=-3.49$. Other features of the two-soliton interaction can be found in Ref. [32].

Figure 5 shows the spatiotemporal diagrams for the 1D density, $n_{1\mathrm{D}}$, for different numbers of dark solitons $N_{s}=(6,10,14,18)$. Similar to the case of two solitons, we observe that the solitons interact repulsively. To measure the strength of the interaction is provided by the distance between the central part and the positive-side border of the soliton gas, $\delta z_{e}=|z_{\mathrm{central}}-z_{\mathrm{bond}}|$.

Figure 6 shows $\delta z_{e}$ as a function of time for different values of $N_{s}$. We obverse that it increases monotonously with time, and its time derivative (speed) change as $N_{s}$ increases. Nevertheless, the speed tends to a limit with the increase of the number of solitons, so that there is no dramatic difference between $N_{s}=14$ and $N_{s}=18$. This happens because the interaction between the solitons has an effective range, as shown in frame (b) of Fig. 4, hence the solitons located near the edges interact weaker with the central ones.

To analyze the case of a large number of solitons, it is enough to take $N_{s}=18$. Frame (a) of Fig. 7 shows the speed at $t=90$ms of each soliton as a function of its initial position for different initial distances $d$. We observe that the central solitons have smaller speeds than their counterparts placed near the edges, so that the speed is given by $v_{s}\simeq\tanh(\gamma_{d}z_{d})$ with $\gamma_{d}=-0.01d+0.077$ in the range of Fig. 7. Frame (b) of Fig. 7 shows the speed of the soliton located near the positive edge at $t=90$ms. Similar to the two-soliton case, the speed decays with the increase of the initial distance, $v_{s}\sim d^{\alpha}$, with $\alpha=-3.385$.

Finally, we consider random initial positions of the solitons, with $N_{s}=18$. We define the initial positions as $z_{j,\epsilon}=z_{j,0}+\epsilon$, where $\epsilon$ is a random fluctuation, and $z_{j,0}$ are the soliton positions, with the mean distance between them $d=4\leavevmode\nobreak\ \mu$m, like in the symmetric case. Figure 8 shows the spatiotemporal diagrams of the 1D density $n_{1\mathrm{D}}$ for two different random realizations. In particular, we assume that $\epsilon$ takes random values in the ranges $[-\epsilon_{\max},\epsilon_{\max}]=[-0.4,0.4]\leavevmode\nobreak\ \mu$m and $[-0.8,0.8]\leavevmode\nobreak\ \mu$m, in panels (a) and (b), respectively. It is observed that the speed of the expansion is higher than in the absence of the randomness, because the interaction energy generates higher internal pressure in the gas of solitons. We analyze the influence of the random-fluctuation magnitude, $\epsilon$ on the dynamics. In particular, we calculate the sum of the squared velocities at the final moment of time,

[TABLE]

Panel (c) of Fig. 8 shows $E_{\epsilon}$ normalized to $E_{0}$ (the kinetic energy of the set of dark solitons with equidistant initial positions) as a function of $\epsilon_{\max}$. We can observe that $E_{\epsilon}$ strongly increases with the growth of $\epsilon_{\max}$, which naturally means that the gas of solitons expands faster when the fluctuations are stronger.

3 The Bose - Fermi Mixture

In this section we consider a dilute superfluid mixture formed by $N_{\mathrm{B}}$ bosonic atoms of mass $m_{\mathrm{B}}$, and $N_{\mathrm{F}}$ fermionic atoms of mass $m_{\mathrm{F}}$ and spin $s_{\mathrm{F}}$. The atoms interact through the pseudopotential, $\delta(\mathbf{r})$ [4]. We assume that bosons form a BEC, described by the Gross-Pitaevskii equation [4], while the local density approximation [4] applies to the description of the weakly interacting fermionic component. Accordingly, the dynamical equations can be derived from the functional,

[TABLE]

where $\mathcal{L}_{B}$ and $\mathcal{L}_{\mathrm{F}}$ are the Lagrangian densities of the Bose and Fermi components, while $\mathcal{L}_{\mathrm{BF}}$ accounts for the interaction between them [58]:

[TABLE]

[TABLE]

[TABLE]

Here {g_{B}}\equiv 4\pi{\hbar^{2}}{a_{B}}/$$m_{\mathrm{B}}, ${g_{\mathrm{F}}}\equiv 4\pi{\hbar^{2}}({a_{\mathrm{F}}}/{m_{\mathrm{F}}})[2S_{\mathrm{F}}/(2S_{\mathrm{F}}+1)]$, and {g_{\mathrm{BF}}}\equiv 4\pi{\hbar^{2}}$$a_{\mathrm{BF}}$$/$$m_{\mathrm{BF}} are three interaction parameters of the mixture, with $a_{\mathrm{B}}$, $a_{\mathrm{F}}$ and $a_{\mathrm{BF}}$ being the respective scattering lengths; ${m_{\mathrm{BF}}}\equiv{m_{B}}{m_{\mathrm{F}}}/({m_{B}}+{m_{\mathrm{F}}})$ is the reduced mass; and $U_{\mathrm{B}/\mathrm{F}}{(\mathbf{r})}$ are external potentials acting on bosons/fermions. Complex wave functions $\Psi_{\mathrm{B}/\mathrm{F}}\left(\mathbf{r},t\right)$ are normalized to the respective numbers of particles, $N_{\mathrm{B}/\mathrm{F}}$. The other parameters of the fermionic Lagrangian density are the same as in Sec. 2.

Varying action $\mathcal{S}$ with respect to $\Psi_{\mathrm{B}}^{\ast}$ and to $\Psi_{\mathrm{F}}^{\ast}$, we derive the following system of nonlinear Schrödinger equations for bosons and fermions:

[TABLE]

[TABLE]

We apply the formalism developed below to the 7Li-6Li mixture, with the same scattering parameter for both species, $a_{\mathrm{B}}=a_{\mathrm{F}}=5$nm. The use of isotopes of the same alkali element is suggested by the similarity of their electric polarizability, thus implying similar external potentials induced by an optical trap. Unless specified otherwise, in what follows below we consider configurations with fully polarized fermions. Note that the BCS and unitarity regimes involve more than one spin state of fermions, hence the magnetic trap will split the respective spin energy levels. For this reason, we assume the presence of the optical trap, which supports equal energy levels for all the spin states, making it possible to discriminate different regimes of the interaction in the BFM. In the BCS and unitarity regimes, we assume balanced populations of the two spin components.

Our analysis is first presented for the GS and dynamics of perturbations around it. In particular, for the GS we focus on determining the spatial correlation $C_{s}$ between the spatial particle densities in both species, defined as

[TABLE]

where ${{\bar{n}}_{\mathrm{B}/\mathrm{F}}}={n_{\mathrm{B}/\mathrm{F}}}-\left\langle{{n_{\mathrm{B}/\mathrm{F}}}}\right\rangle$, $\left\langle{}\right\rangle$ standing for the spatial average. For dynamical perturbations around the GS, a spatiotemporal correlation, which is defined by replacing the spatial average with the spatiotemporal average, is known as the Pearson coefficient $C_{s-t}$ [67]. We remark that when $C_{s}=1$ and $C_{s}=-1$ the mixture is fully synchronized and anti-synchronized, respectively; whereas, the mixture is not synchronized at $C_{s}=0$.

While numerical integration of this system in the 3D form is very heavy computationally, the effective dimension may be reduced to 1D or 2D when the system is tightly confined by a trapping potential. To this end, the VA is employed, making use, as above, of the factorization of the 3D wave function, which includes the Gaussian ansatz in the tightly confined transverse directions. As mentioned above too, the factorization has been widely used for Bose and Fermi systems separately, as it shown in Refs. [14] and [21], respectively. In the next two subsections we reduce the full 3D system to the corresponding 2D and 1D effective systems, using the VA proposed in Ref. [58].

3.1 The two-dimensional reduction

Similar to the case of the pure Fermi gas, we derive 2D equations for the disc-shaped configuration. Accordingly, the structure of the confinement potential is taken as

[TABLE]

where the second term corresponds to the strong harmonic-oscillator trap acting along the $z$ direction. The corresponding factorized ansatz is adopted as

[TABLE]

where $\phi_{\mathrm{B}/\mathrm{F}}$ is normalized to $N_{\mathrm{B}/\mathrm{F}}$, and ${\xi_{\mathrm{B}/\mathrm{F}}}\left({x,y,t}\right)$ are widths of the gas in the confined direction. Substituting the factorized ansatz (25) in action (17) and integrating over $z$, we arrive at the following expression for the effective 2D action:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

so that $n_{2\mathrm{D,B/F}}\equiv\left|{\phi_{\mathrm{B}/\mathrm{F}}}\left({x,y}\right)\right|^{2}$ are the 2D particle densities of the boson and fermion species, and ${e_{2\mathrm{D,B}}}$ and ${e_{2\mathrm{D},\mathrm{F}}}$ are their energy densities:

[TABLE]

[TABLE]

with $C_{2\mathrm{D},\mathrm{F}}\equiv{(3/5)^{1/2}}{(6/(2s_{\mathrm{F}}+1))^{2/3}}\pi$. The field equations for the 2D system are obtained by the variation of the action $S$ given by Eq. (26) with respect to variables $\phi_{\mathrm{B}}$ and $\phi_{\mathrm{F}}$:

[TABLE]

[TABLE]

Relations between $\xi_{\mathrm{B}/\mathrm{F}}$ and $\phi_{\mathrm{B}/\mathrm{F}}$ are produced by the Euler-Lagrange equations associated to $\xi_{\mathrm{B}/\mathrm{F}}$:

[TABLE]

[TABLE]

where ${\kappa_{I,\mathrm{F}}}\equiv{m_{\mathrm{F}}}\omega_{z,\mathrm{F}}^{2}+2{g_{\mathrm{BF}}}{n_{2\mathrm{D,B}}}/[{{\pi^{1/2}}{{({\xi_{\mathrm{B}}^{2}+\xi_{\mathrm{F}}^{2}})}^{3/2}}}]$. Thus, Eqs. (32)-(35) constitute a system of four 2D coupled equations produced by the reduction of the underlying 3D system (21). Note also that when $g_{BF}=0$, the system is decoupled and Eq. (32) corresponds to the dimensional reduction of the Gross-Pitaevskii equation. Equations (34) and (35) for $\xi_{\mathrm{B}/\mathrm{F}}$ can be solved numerically by dint of the Newton’s method. The basic external potential is taken as the harmonic-oscillator one: $U_{2\mathrm{D},\mathrm{B}/\mathrm{F}}=m_{\mathrm{B}/\mathrm{F}}\omega_{x,\mathrm{B}/\mathrm{F}}^{2}{x^{2}}/2+m_{\mathrm{B}/\mathrm{F}}\omega_{y,\mathrm{B}/\mathrm{F}}^{2}{y^{2}}/2$. The simulations were based on the fourth-order Runge-Kutta algorithm with $\Delta t=4.77$ $\mu$s. The spatial discretizations was performed with $\Delta x=1$ $\mu$m, $\Delta y=1$ $\mu$m and $\Delta z=0.05$ $\mu$m. The GS was found by means of the imaginary-time integration. We here focus on the case when the number of bosons is much greater than the number of fermions, viz., $N_{\mathrm{B}}=5\times 10^{4}$ and $N_{\mathrm{F}}=2.5\times 10^{3}$.

Frames (a) and (b) of Fig. 9 show the radial profile of both 2D bosonic and fermionic densities, $n_{2\mathrm{D},\mathrm{B}/\mathrm{F}}$, respectively. The panels for the bosonic and fermionic components are the left and right ones, respectively. Each density has been computed using the VA and the full 3D system. To obtain the 2D profile from the 3D simulations, Eqs. (21) and (22) were solved, and the 3D density was integrated along the $z$ axis, $\bar{n}_{\mathrm{2D,B/F}}=\int_{-\infty}^{+\infty}\left|\Psi_{\mathrm{2D,B/F}}(\mathbf{r})\right|dz$. We infer that the repulsive mixture concentrates the bosons at the center, while the attractive mixture concentrates both species at the center. Panels (c) and (d) of Figure 9 show the radial dependence of the width for both bosonic and fermionic component, respectively. We observe that only the width of the fermionic density profile varies significantly with the change of the scattering length of the inter-species interaction, which is a consequence of a greater number of bosons in comparison with fermions. It is clearly seen that fermions are stronger confined when the interaction is attractive, and their spatial distribution significantly expands when the interaction is repulsive. Similar results have been reported in Refs. [51, 52, 56].

Now, to compare the results obtained from the VA with those produced by the 3D simulations, we note that both profiles are practically identical, except for the repulsive case in which a discrepancy is observed. The inset in panel (a) of Fig. 9 shows that the difference between the two results has a magnitude of nearly three orders of magnitude lower than the density itself. We define the overall percentage error of the VA as $E_{\%,\mathrm{2D}}=\int\int\left|\rho_{\mathrm{2D}}-n_{\mathrm{2D}}\right|dxdy$ (for both species). Figure 10 shows the error for both species as a function of interspecies scattering parameter, $a_{\mathrm{BF}}$. For bosons it takes values $\sim 0.2\%$, and does not change much, as shown in the inset to panel (a) of Fig. 9. For fermions the error is greater than for bosons throughout the observed range, but it is quite small for the attractive mixture. Note that the error increases for the repulsive mixture, but remains lower that $2\%$. Thus we conclude that the 2D approximation is very accurate.

Finally, we measure the correlations of the BFM states. To this end, the spatial correlation, $C_{s}$, in the GS was calculated using the definition given in Eq. (23). Figure 11 presents the analysis of the GS synchronization of the mixture as a function of $a_{\mathrm{BF}}$, where three possible regimes are considered for the fermions: fully polarized, BCS, and unitarity. Parameters of the Lagrangian density for each fermionic regime are given in Table 1. When the interaction is attractive, there is not a large discrepancy between the correlation curves. In fact, for $a_{\mathrm{BF}}\in(-25,-15)$nm the values of $C_{s}\gtrsim 0.9$, and therefore the GS states are synchronized. In the unitarity regime, it is again observed that the correlation reaches a maximum close to $1$ at $a_{\mathrm{BF}}\approx-10$ nm, dropping to negative values when the mixture is strongly repulsive. Also, we observe that the three curves demonstrate stronger demixing when $a_{\mathrm{BF}}$ changes from negative to positive values of $a_{\mathrm{BF}}$, and for $a_{\mathrm{BF}}\gtrsim 15$ the value of $C_{s}$ tends to zero implying that the GS states are not synchronized.

3.2 The one-dimensional reduction

The 1D confinement means, as above, a cigar-shaped configuration elongated in the direction of $z$. In this case, the corresponding confining potentials trap is written as

[TABLE]

where $U_{1\mathrm{D},{\mathrm{B}/\mathrm{F}}}\left({z,t}\right)$ are the axial potentials. Assuming that the transverse trapping potential is strong enough, the dimensional reduction is carried out by means of the usual factorized ansatz for the wave functions,

[TABLE]

where $\sigma_{\mathrm{B}/\mathrm{F}}$ are the transverse GS Gaussians widths. Here, the axial functions, $f_{\mathrm{B}/\mathrm{F}}$, are normalized to $N_{\mathrm{B}/\mathrm{F}}$. For both species, we define the axial density as ${n_{1D,{\mathrm{B}/\mathrm{F}}}}\equiv{\left|f_{\mathrm{B}/\mathrm{F}}\right|^{2}}$. By means of a procedure similar to the one outlined above for the 2D reduction, we derive the Euler-Lagrange equations for the BFM in the 1D approximation:

[TABLE]

[TABLE]

In addition, the algebraic relationships between $\sigma_{\mathrm{B}/\mathrm{F}}$ and $f_{\mathrm{B}/\mathrm{F}}$ are:

[TABLE]

[TABLE]

where $\chi_{I,\mathrm{B}/\mathrm{F}}\equiv m_{\mathrm{B}/\mathrm{F}}\omega_{t,\mathrm{B}/\mathrm{F}}^{2}-2g_{\mathrm{BF}}n_{1\mathrm{D},\mathrm{F}/\mathrm{B}}/[\pi(\sigma_{\mathrm{B}}^{2}+\sigma_{\mathrm{F}}^{2})^{2}]$. Thus, Eqs. (38)-(41) constitute a system of four 1D coupled equations produced by the reduction of the underlying 3D system (21)

• (22). Simulations of the system were performed with mesh parameters $\Delta t=0.5$ $\mu$s and $\Delta z=0.25$ $\mu$m. The external potential is chosen here as the harmonic-oscillator one: $U_{1\mathrm{d},\mathrm{B/F}}=m_{\mathrm{B}/\mathrm{F}}\omega_{z,{\mathrm{B}/\mathrm{F}}}^{2}z^{2}/2$.

The effect of the magnitude and sign of the interaction parameter on the spatial profile of both species, and the accuracy of the VA compared to the 3D solution, can be analyzed by varying the scattering length, $a_{\mathrm{BF}}$. In particular, we consider a mixture with more bosons than fermions, viz., $N_{\mathrm{B}}=5\times 10^{4}$, $N_{\mathrm{F}}=2.5\times 10^{3}$. Because of this condition, the bosonic profile is mainly determined by its self-interaction and the external potential. Frames (a) and (b) of Fig. 12 show the spatial dependence of $n_{1\mathrm{D,B}}$ and $n_{1\mathrm{D},\mathrm{F}}$, respectively. These densities are calculated using both the reduced equations (38) - (41) and the full numerical simulations of Eqs. (21) and (22). In the latter case, the densities are calculated as $\bar{n}_{1\mathrm{D},\;\mathrm{j}}(z)=\int\int\left|\Psi_{\mathrm{j}}(\mathbf{r})\right|^{2}dxdy$ with $j=(F,B)$. We observe that variations of the bosonic density profile are very small in comparison to the significant changes of the inter-species scattering length. The situation is opposite for the fermionic species. As the repulsive scattering length increases, the fermions tend to be pushed to the periphery of the bosonic-gas density profile. This phenomenon is known as demixing [18, 51, 56, 19]. On the other hand, for the attraction case, fermions are, naturally, concentrated in the same region where the bosons are located. Frames (c) and (d) of Fig. 12 correspond to the profiles of $\sigma_{\mathrm{B}}$ and $\sigma_{\mathrm{F}}$. We observe that the width of the bosonic profile slightly increases while proceeding from the inter-species attraction to repulsion. A similar trend is observed for fermions, as shown in panel (d). However, the effect is amplified in the spatial zone of the interaction with the bosons, where the gas is compressed in the case of the attraction, and expands in the case of the repulsion. Note that the fermionic component expands in the confined direction much more than its bosonic counterpart, and that the fermionic width markedly varies, following changes in the density.

Further, one can see in the inset of panel (a) of Fig. 12 the difference between the density calculated by means of the VA and the full 3D simulation, $\Delta n_{\mathrm{1D,B}}=\bar{n}_{\mathrm{1D,B}}-n_{\mathrm{1D,B}}$ is really small. In fact, the difference between the bosonic profiles obtained by both methods is $\sim 2\%$ of the maximum density for all cases (the fact that the error changes very little with variations in $a_{\mathrm{BF}}$ is a consequence of the greater number of bosons). Frame (b) of Fig. 12 shows that, for the attractive mixture, the variational profile is very close to the 3D result, in particular for the case of a_{\mathrm{BF}}=-6\leavevmode\nobreak\nm. For the repulsive mixture, it is observed that the error increases, which is a consequence of the lower fermionic density at the center of the 3D harmonic-oscillator potential, which plays the dominant role for the bosons, hence a monotonously decreasing function in the transverse direction, such as the Gaussian, is not a sufficiently good approximation. We define the global error of the VA as $E_{\%,\mathrm{1D}}=\int_{-\infty}^{+\infty}\left|\bar{n}_{\mathrm{1D},\;j}-n_{\mathrm{1D}\;j}\right|dz$ (for both species). We have found that in the range of $a_{\mathrm{BF}}\in\left(-6,6\right)$ nm the global error for the bosonic species is around $2\%$ for all the values of $a_{\mathrm{BF}}$. For the fermionic species, it goes from $0.5\%$ to $5\%$ depending on $a_{\mathrm{BF}}$, such that for positive value of $a_{\mathrm{BF}}$ the error is higher than for negative ones, and the minimum error is attained at $a_{\mathrm{BF}}\approx-4$nm. This is a consequence of the fact that, for this value of $a_{\mathrm{BF}}$, the interspecies interaction practically compensates the Pauli repulsion, making the dynamics of the fully polarized Fermi gas close to that governed by the linear Schödinger equation (recall that the Gaussian is the solution for the ground state). When the mixture becomes more attractive, the fermionic dynamics is dominated by the bosons, producing a similar error for both species, while for the repulsive mixture the Gaussian approximation is not appropriate. For the non-interacting mixture, the error for the fermions is smaller than for the bosons, because the fermionic density is very low, making the self-interaction terms weak in comparison to the external potential, therefore it is appropriate to use the Gaussian ansatz to approximate the 1D dynamics. Finally, note that the error is lower in the 2D case in comparison with 1D, because the reduction to 2D case is closer to the full 3D model.

Next, we address the BFM dynamics, considering a mixture with arbitrary initial conditions for the 1D fields. To create the initial state, we start with the GS found in the absence of the inter-species interaction ($a_{\mathrm{BF}}=0$). Then, at $t=0$, we switch the interaction on, which may imply the application of the magnetic field, that gives rise to $a_{\mathrm{BF}}\neq 0$ via the FR. Figure 13 shows three cases of the temporal evolution with these initial conditions for $a_{\mathrm{BF}}=-18\leavevmode\nobreak\ \mathrm{nm}$, $a_{\mathrm{BF}}=-26\leavevmode\nobreak\ \mathrm{nm}$, and $a_{\mathrm{BF}}=-34\leavevmode\nobreak\ \mathrm{nm}$. In the first case (panels (a) and (b) of Fig. 13), it is observed that the densities converge towards the center of the potential, as may be expected, creating a pattern of oscillations around the potential minimum; in addition, the fermions are affected by bosons, as shown by the mark left by the bosons in the fermionic density. For the second case (panels (c) and (d) of Fig. 13), it is observed that the increase in the strength of the attractive interaction generates dark solitons in the fermionic density, some of which show oscillatory dynamics very close to that observed in Refs. [29, 27, 68, 69, 71]. The last case (panels (e) and (f) of Fig. 13) shows that the further increase of the strength of the interspecies interaction generates a larger number of dark solitons. In other words, we show that the attractive interaction of fermions with bosons in a state different from the GS eventually generates a gas of dark solitons.

Finally, we address the accuracy of the VA for the dynamical behavior near the GS. Figure 14 displays the spatiotemporal dynamics of the 1D density, as produced by the solution of 3D equations (21) and (22) for $a_{\mathrm{BF}}=-10$ nm. The initial conditions for the 3D dynamics are given by the ansatz based on Eq. (25), with the Gaussian profile along the $z$ axes. Panels (a) and (b) of Fig. 14 show the spatiotemporal diagrams of the bosonic and fermionic densities, making the emergence of dark solitons obvious. This result corroborates that the dark solitons emerge too in the 3D dynamics, which is approximated by the present 1D model. The other panels of Fig. 14 show a comparison of the 1D spatial profiles, as obtained from the 3D simulations, and the 1D VA, for three instants of time: $t=0$ ms ((c) and (d)), $t=25$ ms ((e) and (f)), and $t=50$ ms ((g) and (h)). The results demonstrate that the VA profiles are very similar to their counterparts produced by the 3D simulations, hence the present approximation provides good accuracy and allows one to study dynamical features of the BFM in a sufficiently simple form.

4 Conclusion

In this brief review we have summarized results produced by the VA (variational approximation) for reducing the 3D system to 1D and 2D forms for the Fermi gas and BFM (Bose-Fermi mixture) in the framework of the quasi-mean-field description [21, 58]. The method is based on the Gaussian variational ansatz, which provides very accurate results for the GSs (ground states) of the gases loaded in the disc- and cigar-shaped traps. The reduced equations are useful, in particular, for modeling systems with low atomic densities and large spatiotemporal variations of the external potential. For the 1D case, the reduced equations provide results by means of modest computational resources, allowing one to quickly explore a vast volume of the parameter space. In the 2D case, the required simulation time is still significantly lower than what is necessary for the 3D simulations. We have shown that, for the Fermi gases and BFMs alike, the VA produces results with a very good accuracy, the error, in the comparison to the 3D simulations, being lower than $5\%$, for both the GSs and dynamical states.

For the Fermi gas case in the 2D approximation, we have considered the example of the hexagonal superlattice, built as a superposition of two triangular lattices with different angles among them. This possibility may be relevant for emulating condensed-matter settings, such as graphene-like superlattices. In addition, we have presented results for dark solitons, obtained in the framework of the 1D approximation. We have verified that the interaction is repulsive [32] and strongly depends on the initial distance between the dark solitons.

Finally, for the BFM trapped in the harmonic-oscillator potential we have shown that a change in the interaction strength can generate a gas of dark solitons. The solitons oscillate under the action of the external potential.

Acknowledgents

PD acknowledges partial financial support from DIUFRO project under grant DI18-0066 and CMCC of the Universidad de La Frontera. DL acknowledge partial financial support from Centers of Excellence with BASAL/CONICYT financing, Grant FB0807, CEDENNA and CONICYT-ANILLO ACT 1410. PD and DL acknowledges financial support from FONDECYT 1180905. The work of BAM is supported, in part, by the joint program in physics between NSF and Binational (US-Israel) Science Foundation through project No. 2015616, and by the Israel Science Foundation, through grant No. 1287/17. The authors appreciate a support provided by the PAI-CONICYT program (Chile), grant No. 80160086, and hospitality of Instituto de Alta Investigación, at Universidad de Tarapacá (Arica, Chile).

Appendix A: Nonlinear Schrödinger equation for the fermionic

superfluid

Kim and Zubarev in Ref. [65] proposed an effective hydrodynamic equation for a Fermi gas, in the regime of the BCS-BEC crossover. The equation was derived from the time-dependent density-functional theory and has the form given by:

[TABLE]

where $\Psi$ is a complex field that represent the superfluid wave function, $n\left({\mathbf{r},t}\right)={\left|\Psi\left({\mathbf{r},t}\right)\right|^{2}}$ is the particle density, and $\mu$ is the chemical potential. In addition, the relationship between the chemical potential and the energy density (energy per particle), $\varepsilon\left(n\right)$, is given by:

[TABLE]

For the case of two spin states with balanced populations and a negative scattering length, ${a_{\mathrm{F}}}<0$, the BCS limit corresponds to $k_{F}{\left|{a_{\mathrm{F}}}\right|}\ll 1$, where $k_{F}=(3\pi^{2}n)^{1/3}$ is the Fermi wavenumber. In this limit $\varepsilon$ is given by [72]:

[TABLE]

where ${\varepsilon_{F}}={\hbar^{2}}k_{F}^{2}/\left({2m_{\mathrm{F}}}\right)$ is the Fermi energy. Taking the Eq. (44) into the Eq. (43) the chemical potential takes the form

[TABLE]

where the first term corresponds to the effective Pauli repulsion, and the following ones to the superfluidity due to collisions between the fermions in different spin states. Substituting the latter expression in Eq.(42), and keeping only the first collisional term, we obtain the known nonlinear Schrödinger equation for the fermionic superfluid [65, 66]

[TABLE]

where the last term is similar to one in the Gross-Pitaevskii equation for bosons, but with an extra factor of $1/2$, as the Pauli exclusion principle allows only atoms in different spin states interact via the scattering. We remark that Eq. (46) implies equal particle densities and phases of the wave functions associated with both spin states.

When we have a system with multiple atomic spin states, $\sigma_{j}$, associated with vertical projection of the spin $s$ (with $2s_{\mathrm{F}}+1$ states), we treat the atoms per state as a fully polarized Fermi gas. The term for the interactions by collisions between atoms in different spin states, with the same scattering length ($a_{\mathrm{F}}$), correspond to the scattering term in the Gross-Pitaevskii equation. The motion equation for the atoms in spin states $j$ is given by:

[TABLE]

where $\Psi_{j}$ is the wave function associated with spin projection $\sigma_{j}$, such that $n_{j}\left({\mathbf{r},t}\right)={\left|{{\Psi_{j}}\left({\mathbf{r},t}\right)}\right|^{2}}$ is the respective particle density, and $V(\mathbf{r})$ an external potential, which is assumed to be identical for all the spin states.

In the case of fully locally balanced populations, the density of particles is the same in each component, $n_{1}=n_{2}=...=n_{2s_{\mathrm{F}}+1}$, hence the total density is $n=n_{j}/(2s_{\mathrm{F}}+1)$. Assuming also equal phases of the wave-function components, we define a single wave function, $\Psi=\sqrt{2s_{\mathrm{F}}+1}\Psi_{j}$, such that the Eq. (47) take the form

[TABLE]

where $g_{\mathrm{F}}\equiv 8s_{\mathrm{F}}\pi{\hbar^{2}}{a_{\mathrm{F}}}/(2s_{\mathrm{F}}+1)m_{\mathrm{F}}$ is the scattering coefficient. This equation is the same that Eq. 3 without considered the corrections of the first principles calculations given by $\lambda_{1}$, $\lambda_{2}$ and $\beta$ [60, 61, 62, 63]. In particular, the fully polarized gas, with the interactions between identical fermions suppressed by the Pauli principle, formally corresponds to $s_{\mathrm{F}}=0$, hence $g_{\mathrm{F}}=0$, and the last term of Eq. 48 vanishes.

Finally, the equation (48) can be derived, as the Euler-Lagrange equation,

[TABLE]

from the corresponding action, $\mathcal{S}=\int{dtd{\mathbf{r}}\mathcal{L}}$, with the Lagrangian density

[TABLE]

where the asterisk stands for the complex conjugate. Similar Lagrangian formalisms have been used, in the context of the density-functional theory, in diverse settings [66, 19, 73].