Formation of Matter-Wave Droplet Lattices in Multi-Color Periodic Confinements

Abstract

In the paper, we introduce a new model that addresses the generation of quantum droplets (QDs) in the binary Bose–Einstein condensate (BEC) mixture with mutually symmetric spinor components loaded in multi-color optical lattices (MOLs) of commensurate wavelengths and tunable intensities. The considered MOL confinement is the combination of the four-color optical lattice with an exponential periodic trap, which includes the complete set of the Fourier harmonics. Employing the one-dimensional (1D) extended Gross–Pitäevskii equation (eGPE), we calculate the exact analytical form of the wavefunction, MF/BMF nonlinearities, and MOL trap parameters. Utilizing the exact solutions, the formation of supersolid-like spatially periodic matter-wave droplet lattices and superlattices is illustrated under the space-periodic nonlinearity management. The precise positioning of the density maxima/minima of the droplet patterns at the center of the trap and tunable Anderson-like localization are observed by tuning the symmetry and amplitude of the considered MOL trap. The stability of the obtained solution is confirmed using the Vakhitov–Kolokolov (VK) criterion.

1. Introduction

Bose–Einstein condensates (BECs) and ultracold atoms confined in an optical lattice (OL) trap constitute an ideal experimental platform for the quantum simulations of emerging quantum many-body phenomena [1,2,3,4,5]. The experimental observation of the zero-temperature quantum phase transition in a strongly interacting Mott insulator [6] and weakly interacting Bose gas [7] utilizing the OL confinement provides a significant example for the same. Further, the OL is widely used to investigate the fundamental physics problems: Anderson localization [8], negative temperature [9,10,11], supersolid phase [12], etc., and it also provides the basis for the development of quantum technologies: quantum memory [13], registers [14], optical lattice clocks [15], and entanglement [16]. Experimentally, the OL trap is formed by the superposition of two counter propagating laser beams, resulting in the generation of an artificial crystal of light with spatially periodic polarization patterns, which is tunable through the power and period of the overlapping lasers [1,11]. Different engineered forms of OL geometries such as bi-periodic, kagome, hexagonal, double-well superlattices, etc., are regularly experimentally realized by interfering different sets of laser beams [17]. A number of interesting physical phenomena have been reported in the presence of the above-mentioned geometrically frustrated OLs, including frustrated quantum magnetism at negative absolute temperature [18], many-body localization [19], the exploration of the ionic Hubbard model with ultracold fermions [20], and Hund’s metal in multicomponent Fermi systems [21].

Currently, a new class of quantum liquids, ultradilute quantum droplets (QDs), has aroused a great deal of attention in the field of BECs [22,23,24]. Usually, the BEC dynamics is studied in the presence of an external trap, since it is commonly known to exist in a gaseous phase, and in the absence of a container (i.e., external trap), it expands. However, Petrov, in a pioneering theoretical proposal, pointed out that liquid-like QDs can be stabilized in a weakly interacting binary BEC mixture in free space by realizing a subtle balance between the attractive cubic mean-field (MF) interaction and quadratic repulsive beyond mean-field (BMF) interaction generated due to the quantum fluctuations [25]. Like solitons, these dilute QDs are a self-bound many-body state; however, bright solitons collapse under the influence of the attractive cubic MF interaction. Different from that, QDs can be stabilized in 3D by counter balancing MF and BMF interactions. Based on this stabilization mechanism, QDs are observed in Bose–Bose mixtures [26,27] and dipolar gases [28,29]. The idea of droplet formation without any external trapping is not new, and this is studied in classical liquid or liquid helium systems [30]. However, in comparison to classical liquids in which droplets are generated due to a balance of repulsive interactions (generated due to high density) with attractive van der Waals interactions [30], the formation of QDs is a quantum phenomena with its formation dependent on the balance of the Lee–Huang–Yang (LHY) interaction (due to quantum fluctuations), i.e., the BMF interaction with two-body MF interactions. Further, the observed equilibrium density of QDs is ${10}^8$ orders smaller than liquid helium, and due to the realization in ultracold atoms, this provides versatile control of the tunability of MF/BMF interactions and the geometry of the system.

In this work, we address the dynamics of QDs in a two-component binary BEC mixture in the presence of external multi-color optical lattice (MOL) confinement through the spatial periodic management of MF and BMF nonlinearities in 1D geometry. The QD dynamics is extensively explored in the absence of any confinement, i.e., free space [31,32], and some of us have recently reported the exact theoretical model for QDs in harmonic confinement [33]. However, the study of QDs in the presence of the OL or MOL has received less attention in the current literature. In condensed matter physics, the periodic lattice is considered as one of the fundamental problems, and even in the ultracold atoms domain, a variety of solitonic structures have been investigated in OL traps, both analytically and numerically. The study of QDs’ behavior in the presence of OL or MOL traps acts as a quantum test bed for exploring advanced solid-state physics concepts, such as topological quantum states, discrete systems, etc. Recently, Morera et al. illustrated the generation of QDs, dimerized QDs, and a variety of phases in OL confinement [34,35], whereas the supersolid-like crystallization of QDs was investigated in 1D [36] and on a periodic lattice in a quasi-2D trapped dipolar BEC [37,38]. Further, the stability of QDs is also studied in OL, and the existence of stable dipole QDs has been proposed [39]. The motivation for studying a two-component BEC in the presence of the MOL is twofold: (i) constructing a family of OL traps: the precise control of the intensity and period of overlapping laser beams results in the formation of the OL, bi-color OL (BOL), frustrated bi-color double-well superlattice, tri-color OL (TOL), and four-color OL (FOL); and (ii) a test bed for quantum simulation: optical superlattice confinements comprise a clean controllable many-body test bed, and a variety of physical phenomena [8,9,16] is observed in this trap. In principle, the multi-color beams can be used to design a variety of optical superlattice potentials necessary for supporting the existence of non-trivial QD patterns [40,41]. Therefore, theoretically, it would be important to investigate the behavior of QDs in the presence of the MOL. In this paper, we solve the 1D eGPE for the considered confinement and calculate the non-trivial exact analytical form of the wavefunction, phase, MF/BMF nonlinearities, and trap parameters. This reveals the specific form of the MOL and its correlation with the MF/BMF nonlinearities, which provide tunability for the generation of various QDs density profiles. As an illustration, we show the generation of supersolid-like density schemes in QDs: periodic lattice, bi-periodic superlattice, and bi-periodic double-well superlattice. For each of these patterns, we write the analytical solutions and identify the specific form of the multi-color OL and its parameter domain. The controllable positioning of density maxima at the center of the trap, compression, fragmentation, and Anderson-like localization of the droplet patterns are also observed by tuning the symmetry of the considered trap.

In the following section, we present the analytical framework for solving the 1D eGPE for a weakly interacting Bose–Bose mixture with equal masses and an equal number of atoms in the components under the MOL. The model for the calculated system variables is explained by finding the MOL potential parameters, i.e., the amplitudes and periodicity, the wavefunction, and the form of MF/BMF interactions. It is shown that with a suitable choice of the laser intensity, one can construct the following experimentally relevant forms of confinements: single-color OL, BOL, TOL, double-well superlattice BOL, and FOL. Next, we study the characteristics of QDs under the influence of the above-mentioned traps and illustrate the interesting supersolid-like periodic lattice and superlattice density patterns in QDs. Finally, the stability of the obtained solutions is confirmed using the Vakhitov–Kolokolov (VK) criterion.

2. The Model and Analytical Framework

We start by considering the two-component mass-balanced binary BEC under the influence of the BMF (LHY corrections for quantum fluctuations) in the presence of spatially varying MOL confinement. The choice of equal masses and an equal number of atoms in the BEC mixture makes the result analysis clearer and easier. In the 1D configuration, the QDs are observed under the subtle balance of the slightly repulsive MF interaction with the attractive BMF. In this geometry, the MF and BMF interactions’ contribution to the energy per particle is proportional to n and $\sqrt{n}$, where n is the density of the gas [31], and the system is described by the following equations, the 1D eGPE [25,39]:

$$\begin{array}{c}\hfill i\hslash \frac{\partial {\psi }_1}{\partial t}=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2{\psi }_1}{\partial x^2}+({\Lambda }_s\left(x\right){\left|{\psi }_1\right|}^2+{\Lambda }_c\left(x\right){\left|{\psi }_2\right|}^2){\psi }_1-\Gamma \left(x\right)({\left|{\psi }_1\right|}^2+|{\psi }_2{{|}^2)}^{1/2}{\psi }_1+v\left(x\right){\psi }_1,\end{array}$$

(1)

$$\begin{array}{c}\hfill i\hslash \frac{\partial {\psi }_2}{\partial t}=-\frac{\hslash {}^2}{2m}\frac{{\partial }^2{\psi }_2}{\partial x^2}+({\Lambda }_c\left(x\right){\left|{\psi }_1\right|}^2+{\Lambda }_s\left(x\right){\left|{\psi }_2\right|}^2){\psi }_2-\Gamma \left(x\right)({\left|{\psi }_1\right|}^2+|{\psi }_2{{|}^2)}^{1/2}{\psi }_2+v\left(x\right){\psi }_2.\end{array}$$

(2)

Here, $v\left(x\right)$ is the external MOL confinement present in the considered system with ${\psi }_1$ (${\psi }_2$) representing the wavefunctions of binary mixture components. In Equations (1) and (2), we take the interaction strengths controlling the repulsion between the atoms in each components to be equal: $g_{11}=g_{22}\equiv g=2{\hslash }^2{a}_s\left(x\right)/\left(m{a}_{\perp }^2\right)$ and $g_c=g_{12}$. Here, ${\Lambda }_s\left(x\right)=(g_c+3g)/2$ represents the self-interaction coefficients, whereas ${\Lambda }_c\left(x\right)=(g_c-g)/2$ is the cross interaction coefficients along with $\Gamma \left(x\right)=\sqrt{m}{g}^{3/2}/(\pi \hslash )$ [39]. $a_s\left(x\right)$ represent the space-dependent inter- and intra-components’ atomic scattering lengths, which are tunable through the Feshbach resonance technique [42]. Thus, the sign and strength of both (inter- and intra-) components’ atomic scattering length, i.e., MF and BMF interactions, can be experimentally modulated. Here, m is the mass of the BEC atoms and ℏ is the scaled Planck’s constant.

Next, we reduce the dynamics of the considered system to the space-dependent dimensionless single eGPE by assuming ${\psi }_1={\psi }_2={\psi }_0\psi$, i.e., mutually symmetric spinor components in the binary mixture [25,31]:

$$\begin{array}{c}\hfill i\frac{\partial \psi }{\partial t}=-\frac{1}{2}\frac{{\partial }^2\psi }{\partial x^2}-g_1\left(x\right)\left|\psi \right|\psi +g_2\left(x\right){\left|\psi \right|}^2\psi +V\left(x\right)\psi .\end{array}$$

(3)

Equation (3) is the extended form of the Gross–Piteäevskii equation (GPE) with the external confinement added. For the case, $V\left(x\right)=0$, it becomes the 1D eGPE in free space, which is extensively explored to investigate the QDs’ dynamics [23,24]. In the equation, $\psi (x,t)$ is the condensate wave function of the QDs having mass m, and $g_1\left(x\right)=\Gamma \left(x\right)$, $g_2\left(x\right)={\Lambda }_s\left(x\right)+{\Lambda }_c\left(x\right)$ are the magnitude coupling strengths of the two-component Bose–Bose mixture representing BMF and MF interactions, respectively. Here, $\frac{\hslash g_1\left(x\right)}{\Gamma \left(x\right)}\sqrt{\frac{{\Lambda }_s+{\Lambda }_c}{2m{g}_2\left(x\right)}}$, $\frac{\hslash ({\Lambda }_s+{\Lambda }_c)g_1^2\left(x\right)}{2g_2\left(x\right){\Gamma }^2\left(x\right)}$, $\frac{\sqrt{2}{g}_2\left(x\right)\Gamma \left(x\right)}{({\Lambda }_s+{\Lambda }_c)g_1\left(x\right)}$ are the magnitudes of the scaling parameters $x_0$, $t_0$, ${\psi }_0$, respectively [31]. To investigate the structure and dynamics of QDs in the presence of MOL, we consider the form of the external trap:

$$\begin{array}{c}\hfill V\left(x\right)=\sum _{j=1}^4{V}_j{cos}^2\left(jkx\right)+V_{5}exp\left[2(p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right))\right],\end{array}$$

(4)

which is a combination of the FOL with commensurate lattice periods and an exponential periodic trap. The combination results in the generation of an MOL potential trap, which becomes an FOL for $V_5\to 0$. The choice of the exponential periodic trap ensures the presence of the complete set of the Fourier harmonics in the considered trap combination. Here, $p_1$ and $p_2$ are real constants. In Equation (4), $k=2\pi a_{\perp }/\lambda$ is the scaled lattice wave vector, which is commensurate for the four-color laser beams ($k,2k,3k,4k$) with $a_{\perp }=\sqrt{\hslash /m{\omega }_{\perp }}$, and ${\omega }_{\perp }$ is the transverse oscillator frequency. Here, $V_j$ ($j=1,2,3,4,5$) represents the potential depths of each OL and is connected to the recoil energy: $E_R=\frac{2{\pi }^2{\hslash }^2}{m{\lambda }^2}$ for the laser of wavelength (color) $\lambda$ and mass m of BEC atoms [40]. Recently, ultracold atoms were investigated in the FOL potential to obtain eightfold rotationally symmetric OLs [43] and solitonic solutions [41]. In this paper, we considered the cigar-shaped $L{i}^7$ BEC atoms trapped with transverse frequency ${\omega }_{\perp }=2\pi \times 710$ Hz, atomic scattering length $a_s=-0.21$ nm, and a CO${}_2$ laser of wavelength $\lambda =10.62$ $\mathsf{\mu }$m [44]. In the experiments, various forms of OLs can be realized by tuning the applied magnetic field, the magnitude of k, and the angle between superimposing laser beams [45].

In order to construct the analytical solution form of $\psi$ for Equation (3), based on the general similarity transformation scheme used for constructing matter-wave solitons in 1D geometry [46,47], we start with following ansatz solution:

$$\psi (x,t)=A(x,t)U\left[\eta (x,t)\right]e^{i\varphi (x,t)},$$

(5)

where $A(x,t)$, $\varphi (x,t)$, and $U\left[\eta \right(x,t\left)\right]$, being real functions, are the space- and time-modulated amplitude, phase, and similarity variables, respectively. Using the ansatz solution, our goal is to connect Equation (3) to the solvable differential equation:

$$-\frac{{\partial }^{2}U}{\partial {\eta }^2}-G_1\mid U\left(\eta \right)\mid U+G_2\mid U\left(\eta \right){\mid }^{2}U=EU,$$

(6)

such that we obtain the following consistency conditions on the amplitude and MF and BMF nonlinearities for the chosen potential of Equation (4) (see Appendix A):

$$\begin{array}{c}\hfill {\left[A^2(x,t){\eta }_x(x,t)\right]}_x=0,{\eta }_t(x,t)+{\eta }_x(x,t){\varphi }_x(x,t)=0,\end{array}$$

(7)

$$\begin{array}{c}\hfill G_1{\eta }_x^2(x,t)-2A(x,t)g_1(x,t)=0,G_2{\eta }_x^2(x,t)-2A^2(x,t)g_2(x,t)=0,\end{array}$$

(8)

$$\begin{array}{c}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{A_t(x,t)}{A(x,t)}+\frac{1}{2A^2(x,t)}{\left[A^2(x,t){\varphi }_x(x,t)\right]}_x=0,\end{array}$$

(9)

$$\begin{array}{c}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{A_{xx}(x,t)}{2A(x,t)}-\frac{{\varphi }_x^2(x,t)}{2}-{\varphi }_t(x,t)-\frac{1}{2}E{\eta }_x^2(x,t)-V\left(x\right)=0.\end{array}$$

(10)

In Equation (6), E is the eigenvalue of the equation, $G_1$, $G_2$ denote the constantBMF/MF interactions, which can take a positive or negative magnitude depending on the sign of the inter- and intra-component atomic scattering length. Here, the function with the subscript implies the partial differentiation of the function with respect to the subscripted variable. The above set of consistency conditions are simultaneously solved to obtain the amplitude, phase, and MF/BMF:

$$\begin{array}{c}\hfill {\eta }_x(x,t)=\frac{b\left(t\right)}{A^2(x,t)},{\varphi }_z=-\frac{{\eta }_t(x,t)}{{\eta }_x(x,t)},g_1(x,t)=G_1\frac{{\eta }_x^2(x,t)}{2A(x,t)},g_2(x,t)=G_2\frac{{\eta }_x^2(x,t)}{2A^2(x,t)},\end{array}$$

(11)

where $b\left(t\right)$ is an integration constant. It is evident from Equation (11) that the form of MF/BMF nonlinearities and the phase is directly dependent on the amplitude, which will be determined by solving the consistency Equation (11). For that purpose, we substitute the trap expression from Equation (4) into the set of consistency Equations (7) and (11) and choose $\eta (x,t)=\gamma {\int }_0^{z}exp[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]\partial x$, to obtain the exact analytical form of the amplitude, phase, and nonlinearities:

$$\begin{array}{c}\hfill A(x,t)=\sqrt{\frac{b\left(t\right)}{\gamma \times exp[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]}},\end{array}$$

(12)

$$\begin{array}{c}\hfill \theta (x,t)=[2p_1^2{k}^2+8p_2^2{k}^2-8k^2(p_1+4p_2)]t,\end{array}$$

(13)

$$\begin{array}{c}\hfill g_1(x,t)=\frac{G_1{\gamma }^{3/2}}{2b\left(t\right)}exp{[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]}^{\frac{3}{2}},\end{array}$$

(14)

$$\begin{array}{c}\hfill g_2(x,t)=\frac{G_2{\gamma }^3}{2b\left(t\right)}exp{[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]}^3,\end{array}$$

(15)

with the potential depths of each overlapping OL connected in the following manner:

$$\begin{array}{c}\hfill \\ \hfill V_1=p_1{k}^2\left[1+\frac{p_2}{2}\right],V_2=k^2\left[\frac{p_1^2}{8}+4p_2\right]V_3=-\frac{p_1{p}_2{k}^2}{2}{V}_4=-\frac{p_2^2{k}^2}{2}{V}_5=\frac{E{\gamma }^2}{2}.\end{array}$$

(16)

Equation (16) reveals a non-trivial correlation in between trap parameters $p_1$ and $p_2$, and this is one of the important results of the article. In principle, by the suitable tuning of these parameters, one can realize various forms of $V\left(x\right)$: OL (Figure 1c), BOL (Figure 1d), double-well superlattice (Figure 1e), frustrated double-well optical superlattice (Figure 1f), bi-periodic frustrated double-well optical superlattice (Figure 1g), etc., and study the QDs’ profiles with precise knowledge of the MF/BMF interactions. In Equation (12), $b\left(t\right)=b$ (constant) $={\gamma }^2$.

The solution of Equation (6) can be given as: $U\left[\eta \right]=\frac{3(E/G_1)}{1+\sqrt{1-\frac{E}{{\mu }_0}\frac{G_2}{G_1^2}}cosh\left(\sqrt{-\mathrm{E}}\eta \right)}$ with ${\mu }_0=-2/9$, $E<0$, $G_1<0$, and $G_2>0$ [25,31]. Thus, the complete solution of Equation (3) can be written as:

$$\begin{array}{c}\hfill \psi (x,t)=\sqrt{\frac{b}{\gamma \times exp[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]}}\times exp\left[i(2p_1^2{k}^2+8p_2^2{k}^2-8k^2(p_1+4p_2)t)\right]\\ \hfill \frac{\frac{3E}{G_1}}{1+\sqrt{1-\frac{E}{{\mu }_0}\frac{G_2}{G_1^2}}cosh\left(\sqrt{-\mathrm{E}}\left(\gamma {\int }_0^{x}exp[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]\right)\right)}.\end{array}$$

(17)

Further, one can choose the solution of Equation (6) in terms of the Jacobi elliptic function ($cn$) as: $U\left(\eta \right)=Bcn[\beta \eta ,q]+D,$ with $B=\sqrt{\frac{2}{(2q^2-1)}}D$$>0$, $D=\frac{G_1}{3G_2}$$<0$, ${\beta }^2=-\left(\frac{6G_2}{(2q^2-1)}\right)$$D^2$, and $q^2>1/2$ [48]. For this case, the complete wavefunction form of Equation (3) becomes:

$$\begin{array}{c}\hfill \psi (x,t)=\sqrt{\frac{b}{\gamma \times exp[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]}}\times exp\left[i(2p_1^2{k}^2+8p_2^2{k}^2-8k^2(p_1+4p_2)t)\right]\\ \hfill \left[Bcn(\beta (\gamma {\int }_0^{x}exp[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]),q)+D\right],\end{array}$$

(18)

where q is the modulus parameter for the Jacobi elliptic function $cn$. In principle, the $cn$ function can also possess a family of solutions for the range of the modulus parameter, 0 $\le q\le 1/2$, and at $q=0$ will signify a periodic QD profile. Thus, it is worth indicating that we constructed a large family of exact analytical solutions of the 1D eGPE for the considered MOL trap configuration.

3. Results

The exact expressions of the distributed confinement parameters are already solved in Equation (16) where the parameters ($p_1,p_2$) can be widely tuned. We are now in a position to analyze the relevant QDs’ profiles along with the corresponding confinement form. In the following, we would like to perform a detailed study of the tuning of the trap under consideration and of the density patterns of the QDs.

3.1. Potential Profiles and Corresponding Trap Parameters

From the constructed model, the general form of the potential can be written as:

$$\begin{array}{c}\hfill V\left(x\right)=p_1{k}^2\left[1+\frac{p_2}{2}\right]{cos}^2\left(kx\right)+k^2\left[\frac{p_1^2}{8}+4p_2\right]{cos}^2\left(2kx\right)-\frac{p_1{p}_2{k}^2}{2}{cos}^2\left(3kx\right)\\ \hfill -\frac{p_2^2{k}^2}{2}{cos}^2\left(4kx\right)+\frac{E{\gamma }^2}{2}exp\left[2(p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right))\right],\end{array}$$

(19)

where $-2/9<E<0$ and $\gamma >0$. With suitable tuning of the physical parameters $p_1$ and $p_2$, one can construct various potential profiles and obtain the exact analytical expression of the wavefunction. Physically, $p_1$ and $p_2$ are connected to the power of overlapping laser intensities. Figure 1 illustrates the various trap patterns by tuning of the trap parameters ($p_1$, $p_2$). In Figure 1a,b, the variation of the potential profile is depicted for $\gamma =0.05$, $k=0.84$, $E=-2/9$ in the interval $[-2\pi /k,2\pi /k]$: (a) for fixed $p_2=1$, with $p_1$ varying from $-5$ to $+5$; and in (b) for fixed $p_1=1$, $p_2$ varying from $-5$ to $+5$. This clearly indicates the various shapes of the resultant potential by changing the magnitude of $p_1$ and $p_2$. For $p_1=p_2=0$, $V\left(x\right)$ becomes constant, representing the free space scenario. Further, the trap takes the form of the two-color BOL with periodicity ($k,2k$) for $p_1>0$ and $p_2=0$, whereas if $p_2\ne 0$, then the MOL becomes a disordered double-well superlattice with frustration present both in the inter- and intra-well separations [41]. Additionally, we observe interesting trap forms for $p_1,p_2<0$. For $p_1<0$ with $p_2$ positive constant, this leads to the formation of the triple-well superlattice, which is evident from Figure 1a.

Similarly, we illustrate a half-wavelength shift of this double-well superlattice by tuning $p_2<0$ with $p_1$ positive constant. For better insight into the tunability of the generated MOL trap form, we plot the shape of the trap at some specific points and illustrate the experimentally realizable trap configuration: (a) OL: for $p_1=0$, $p_2=0.1$ (black solid line) and the flipping of this OL by the half-wavelength for $p_1=0$, $p_2=-0.1$ (dashed blue line) (Figure 1c); (b) BOL: overlapping of two OLs of frequency $3k$ and $4k$ for $p_1=0$, $p_2=3.5$ (black solid line) (Figure 1d); (c) triple-well optical superlattice: $p_1=-5$, $p_2=1$ (red solid line) (Figure 1e); (d) frustrated double-well optical superlattice: $p_1=1$, $p_2=-5$ (blue solid line) (Figure 1f); (e) bi-periodic frustrated double-well optical superlattice: $p_1=1$, $p_2=4$ (red line) (Figure 1g). Thus, potential Equation (19) consists of the family of OLs, and in principle, the complete set of the Fourier harmonics can be generated from it.

For ${\gamma }^2<1$ or $p_1,p_2<0$, the higher-order terms of the exponential trap tend towards zero, and the trap behaves as an FOL confinement. Utilizing Equation (19), we identify the points in

Table 1. Various shapes of the MOL potential by tuning the magnitude of the power of the laser beam, i.e., $p_1$ and $p_2$.

Multi-Color OL (for $\mathit{\mu }=0$)
$p_1$	$p_2$	Trap form
0	0	Free space
<1	0	OL (k)
0	≠0	BOL ($2k$, $4k$)
>1	0	BOL (k, $2k$)
8	−2	BOL ($3k$, $4k$)
8	2	BOL ($3k$, $4k$)
≠8	−2	TOL ($2k$, $3k$, $4k$)
≠8	2	TOL ($2k$, $3k$, $4k$)
Other points	Other points	FOL (k, $2k$, $3k$, $4k$)

at which the shape of the MOL becomes the OL, BOL, or TOL in the $(p_1$, $p_2)$ space for $\mu =0$ [41]. On the contrary, the FOL is obtained in the entire space, excluding the points indicated in the table. It needs to be emphasized here that the constructed MOL trap configuration provides a large variety of experimentally useful potential profiles with exact analytical solutions, which may find applications for quantum information processing and simulations [5,16,17]. We further illustrate the results by displaying various QD patterns in the presence of the above-mentioned trap configurations.

3.2. Periodic Lattice Density Patterns in QDs

In this section, we correlate the above-mentioned understanding of the confinement engineering with the atomic condensate density for a variety of experimentally relevant forms of the MOL confinement to investigate the droplet characteristics in it. As discussed earlier, a family of OL traps can be constructed from (19), but we begin by investigating the generation of QDs in the presence of the two-color BOL($k,2k$) with competing MF and BMF nonlinearities. The superposition of two-color OLs results in the formation of quasi-periodic optical superlattice confinement, which is used to study various interesting physical phenomena: Anderson localization, frustrated quantum magnetism, negative absolute temperature, etc. [10,18,19]. Motivated by that, we construct the BOL with commensurate period $(k,2k)$ from potential Equation (20) by taking $p_2=0$. Further, the form of the MF and BMF nonlinearities is: $g_2\left(x\right)=(G_2{\gamma }^3/2b)exp{\left[p_1{cos}^2\left(kx\right)\right]}^3$ and $g_1\left(x\right)=(G_1{\gamma }^{3/2}/2b)exp{\left[p_1{cos}^2\left(kx\right)\right]}^{\frac{3}{2}}$, respectively, with $b\left(t\right)=b$ (constant). Thus, the resulting form of the potential from Equation (19) takes the form:

$$\begin{array}{c}\hfill V\left(x\right)=p_1{k}^2{cos}^2\left(kx\right)+\frac{{\left(p_{1}k\right)}^2}{8}{cos}^2\left(2kx\right)+\frac{E{\gamma }^2}{2}exp\left[2\left(p_1{cos}^2\left(kx\right)\right)\right],\end{array}$$

(20)

for which the wavefunction solution of Equation (3) can be given as:

$$\begin{array}{c}\hfill \psi \left(x\right)=\sqrt{\frac{b}{\gamma \times exp\left[p_1{cos}^2\left(kx\right)\right]}}\frac{\frac{3E}{G_1}}{1+\sqrt{1-\frac{E}{{\mu }_0}\frac{G_2}{G_1^2}}cosh(\sqrt{-\mathrm{E}}\left(\gamma {\int }_0^{x}exp\left[p_1{cos}^2\left(kx\right)\right)\right)}\\ \hfill \times exp\left[i(2p_1^2{k}^2-8k^2{p}_1)t)\right].\end{array}$$

(21)

Using the above wavefunction, we illustrate the density profiles of QDs for $\mp p_1$ for the BOL potential Equation (20) in Figure 2 and Figure 3, respectively. In these figures, each plot (a–d) has two panels: the upper panel shows the density plot, and the lower panel consists of a 2D plot of the density. In Figure 2a, initially, we take $p_1=0$, making $V\left(x\right)=$ constant, i.e., the free space potential and with the MF interaction ($G_2=0.999999999$) and the BMF interaction ($G_1=-1$). As there is a subtle balance between MF and BMF, thus we observe the flat top density profile, which is the signature feature of QDs.

Similarly, condensate density has previously been illustrated in the free space potential, and the observed pattern is in conformity with the physical situations reported in the literature [32]. Here, the magnitude of other physical parameters: $b=1$, $\gamma =1$, $k=0.84$, $E=-2/9$. Next, we investigate the impact of $p_1<0$ on the droplet density profile for the same physical parameters. We took $p_1=-0.5$, $-1$, $-1.5$ and depict their corresponding condensate density patterns in Figure 2b–d. It is evident from the figure that with $p_1$ tending from $0\to -1.5$, this leads to the expansion of the QDs, and due to the increase in the lattice depth of the trap, it forms the lattice patterns in the QDs. Importantly, the signature of the potential is superimposed over the flat density profile, and correspondingly, we observe periodic lattice density pattern on the QDs. We observe that for $p_1<0$, the density periodic lattice maxima are located at $x=0$. However, to establish supersolidity in the system, we need to show the spontaneous breaking of the translational symmetry and the spontaneous breaking of the gauge symmetry in the considered system. However, in absence of that, we term this supersolid-like periodic density patterns, as done in previous studies on spinor BECs [38,49]. Further, we note that the decreasing magnitude of $p_1$ also results in the increase of the width of the QDs, i.e., leading to its expansion. Thus, in principle, the depth and width of these formed periodic lattice density patterns in QDs are connected to the magnitude of $p_1$, and from the constructed analytical model, we reveal a non-trivial correlation in between them.

In Figure 3, we illustrate the density variation of the wavefunction solution (21) for $p_1>0$ with its magnitude changing from $0\to 1.5$ for the identical physical parameter as taken in the case of $p_1<0$. Here, we take the parameter values: $b=1$, $\gamma =1$, $k=0.84$, $E=-2/9$, $G_2=0.999999999$, $G_1=-1$. We begin with $p_1=0$ in Figure 3a, i.e., the free space situation. Next, the magnitude of $p_1$ increased to $0.5$, resulting in the formation of multiple peaks in the density profile, which are due to the presence of the OL trap (Figure 3b). However, in comparison to $p_1=-0.5$ shown in Figure 2b, here, the number of droplet lattice peaks is less and the width is compressed. Subsequently, we observe a reduction in the width of the droplet and the number of density lattice wells with the increase of $p_1$ (Figure 3c,d). The condensate atomic density illustrates an Anderson-like localization with the increasing magnitude of the laser intensity $p_1$. This behavior was confirmed by comparing the maximum amplitude of the density for $\pm p_1$. In comparison to $p_1<0$, in which the maximum amplitude remains the same (Figure 2), here, its magnitude increases with $p_1$ tending from 0 to $1.5$, showing the localization of condensate atoms. Further, in comparison to $p_1<0$, we note here that the density periodic lattice minima are located at $x=0$ for $p_1>0$.

In order to physically understand the reasons for the observation of the QDs’ density patterns in Figure 2 and Figure 3, we plot the profile of MF/BMF nonlinearities with respect to the variation of the MOL trap parameters with $p_2=0$ and $p_1$ varying from $-1.5$ to $+1.5$ in Figure 4a. It is evident from the figure that the magnitude of MF and BMF interactions remains comparable till $p_1$ changes in the interval [$-1.5,0$]. This is due to the negative magnitude of $p_1$, which is present in the exponential term of both interactions, and it ensures the observation of the flat top density profile in Figure 2a–d. Furthermore, the increase in the depth of the QDs’ periodic lattice patterns is due to the potential depth of the trap with changing the value of $p_1$ from 0 to $-1.5$, which leads to the expansion of the droplet profile. For $p_1>0$ and its increase from 0 to $+1.5$, this leads to the amplification of both MF/BMF amplitudes (due to the presence of the exponential factor); however, in this region, the MF term dominates the BMF due to the $exp{\left[p_1{cos}^2\left(kx\right)\right]}^3$ term in it. This leads to an imbalance of the MF and BMF interaction strengths, as is evident from Figure 4a, and simultaneously, the system tends towards the soliton region with a dominant repulsive MF. Since the balance of these nonlinearities is essential for the observation of the droplet profile, the increase of $p_1$ results in the gradual decrease of the flat top density profile (see Figure 3) and the phenomena of the Anderson-like localization observed, which is a characteristic feature of disordered optical lattices [40]. Thus, we illustrate an interesting and non-trivial tunability of QDs’ density patterns by changing the sign (i.e., symmetry) and potential depth of the BOL ($k,2k$) confinement.

3.3. Double-Well Superlattice Density Patterns in QDs

In this section, we illustrate the formation, expansion, and compression of double-well superlattice density patterns in QDs under the MOL confinement. For that purpose, we take $p_2\ne 0$ in the potential Equation (20) with $\gamma =1$ and $k=0.84$. Utilizing the corresponding wavefunction solution from Equation (20), we plot the atomic condensate density in Figure 5a–d for $p_2=0$; $p_2=-0.5$, $p_2=-1$, and $p_2=-1.5$, respectively, for $p_1=0.5$, $G_1=-1$, $G_2=0.999999999$, $E=-2/9$. It is evident from Figure 5a that the periodic lattice density patterns are formed for $p_2=0$ and $p_1=0.5$. However, as $p_2$ changes from $0\to -0.5$, the double-well superlattice density patterns become visible in the QDs’ profile (Figure 5b). The decreasing magnitude of $p_2$ from $-0.5\to -1.5$ results in the expansion of the QDs along with the increase in the depths of these double-well superlattice density patterns. Further, in Figure 6, we depict the impact of $p_2>0$ and $p_1=0.5$ on the QDs’ profile for the same physical parameter values. Here, we plot for (a) $p_2=0$, (b) $p_2=0.5$, (c) $p_2=1$, and (d) $p_2=1.5$, respectively. In comparison to the double-well superlattice density pattern of QDs for $p_2=-0.5$, we observe the formation of a bi-periodic density lattice for the case of $p_2=0.5$ (Figure 6b). The depth of this bi-periodic density pattern increases, and the width of the droplet decreases with the increasing magnitude of $p_2$ from $0.5\to 1.5$ (Figure 6c,d).

In order to physically understand the reasons for the observation of the QDs’ density superlattice patterns in Figure 5 and Figure 6, we plot the profile of MF/BMF nonlinearities with respect to the variation of the MOL trap parameters with $p_1=0.5$ and $p_2$ varying from $-1.5$ to $+1.5$ in Figure 4b. In comparison to Figure 4a, in which $p_2=0$, here, the double-well superlattice trap is formed due to $p_2\ne 0$. As discussed earlier, in this regime also, the magnitude of MF and BMF interactions remains comparable till $p_2$ changes from [$-1.5,0$] (shown in Figure 4b). This results in the observation of the flat top density profile as shown in Figure 6a–d. The superlattice QDs’ density patterns are due to the superlattice periodicity of MF/BMF interactions and the corresponding resulting trap. For $p_2>0$ and tending from 0 to $+1.5$, then the MF term dominates the BMF due to the $exp{[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]}^3$ factor. Similar to the previous regime, we observe a decrease in the flat top density profile with the increasing magnitude of $p_2$. However, different from the previous case, here, the droplet density patterns are bi-periodic. Therefore, we observe that the compression and expansion of the droplet density patterns are due to the spatial periodic management of MF and BMF nonlinearities: $g_1(x,t)=[G_1{\gamma }^{3/2}/2b\left(t\right)]exp{[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]}^{\frac{3}{2}}$ and $g_2(x,t)=[G_2{\gamma }^3/2b\left(t\right)]exp{[p_1{cos}^2\left(kx\right)+p_2{cos}^2\left(2kx\right)]}^3,$ respectively, with respect to the changing sign of $\pm p_2$. Thus, we illustrate an interesting transition of the supersolid-like double-well superlattice and bi-periodic density patterns in droplets from the starting periodic density profile by changing the magnitude of $p_2$ and MF/BMF nonlinearities.

3.4. Stability of QDs in MOL Confinement

In the above sections, we illustrated the generation of a variety of droplet density patterns in the MOL confinement by tuning the magnitude of ($p_1$, $p_2$) in Equation (17). In this section, we evaluate the stability of the obtained wavefunction solution (17), and for that purpose, we employed the VK criterion, which is extensively utilized to determine the stability of nonlinear Schrödinger equation solutions [50]. According to the VK criterion, a solution is found to be stable if the slope of the number of atoms with respect to the chemical potential, i.e., $N_E=dN/dE$, is positive. Here, N is the normalization, and E is the chemical potential of the system. For the case $N_E<0$, the solution is unstable, whereas $N_E=0$ provides the instability threshold of the obtained solution. Now, using Equation (17) and $N={\int }_{-\infty }^{+\infty }{\left|\psi \right|}^2\partial x$, one can estimate the correlation between normalization N and E as:

$$N=\frac{4}{3}\left[ln\left(\frac{1+\sqrt{\frac{E}{{\mu }_0}}}{\sqrt{1-\frac{E}{{\mu }_0}}}\right)-\sqrt{\frac{E}{{\mu }_0}}\right],$$

(22)

where $\frac{G_2}{G_1^2}\approx 1$ and $G_1=-1$. Equation (22) estimates the magnitude of N in the presence of the MOL and is equal to the N reported for the free space [31,32]. Thus, even in the presence of the MOL confinement, N is conserved, and the considered system shows a continuous symmetry property according to Noether’s theorem [51]. In Figure 7, using Equation (22), we plot $N_E$ with respect to E, where $G_1=-1$. It is evident from Figure 7a that the magnitude of $N_E$ is positive with respect to its variation E, which indicates the stable nature of the obtained solution.

The small-scale fluctuations in the droplet under the MOL confinement can be estimated by linearizing the 1D eGPE (3) around the ground state given by Equation (17) [32]. We take $\psi (x,t)={\psi }_0\left(x\right)+\delta \psi (x,t)$, where ${\psi }_0\left(x\right)$ is the stationary solution (17) and $\delta \psi (x,t)=\left[\begin{array}{c}\delta {\psi }_R\\ \delta {\psi }_I\end{array}\right]<<1$ is the small perturbation with $\delta {\psi }_R$; $\delta {\psi }_I$ are the real and imaginary parts of $\delta \psi$. With the substitution of $\delta \psi$ in Equation (3), this leads to the well-known Bogoliubov–de Gennes (BdG) equation:

$$\begin{array}{c}\hfill \left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]\frac{\partial }{\partial t}\left[\begin{array}{c}\delta {\psi }_R\\ \delta {\psi }_I\end{array}\right]=\left[\begin{array}{cc}T& 0\\ 0& T^{\prime }\end{array}\right]\left[\begin{array}{c}\delta {\psi }_R\\ \delta {\psi }_I\end{array}\right],\end{array}$$

(23)

with $T=-\frac{1}{2}\frac{{\partial }^2}{\partial x^2}+3n{g}_2\left(x\right)-\frac{1}{2}{g}_1\left(x\right)n^{1/2}+v\left(x\right)$, $T^{\prime }=-\frac{1}{2}\frac{{\partial }^2}{\partial x^2}+g_2\left(x\right)n-\frac{1}{2}{g}_1\left(x\right)n^{1/2}+v\left(x\right)$, and $n=|{\psi }_0{\left(x\right)|}^2$. Here, we consider $\delta \psi =exp\left[i\right(lx-\omega t\left)\right]$, and applying it in Equation (23), then this yields the perturbation eigenmodes, where l denotes the wave number and $\Omega$ stands for the frequency. The resulting dispersion relation can be given as:

$${\Omega }^2=\frac{l^4}{4}+l^2\left(2V\left(x\right)-4g_1\left(x\right)n^{1/2}+4n{g}_2\left(x\right)\right),$$

(24)

by neglecting the l independent terms. From Equation (24), $\Omega$ is imaginary for $4g_1\left(x\right)n^{1/2}>2v\left(x\right)+4n{g}_2\left(x\right)+l^2/4$, and these are the instability region. In Figure 7b, $Im(\Omega )$ is depicted as a function of $p_1$ and the wavenumber (l), keeping $p_2=0$, representing the two-color BOL ($(k,2k)$) case, whereas in Figure 7c, we locate $Im(\Omega )$ with respect to the changing magnitude of $p_2$ and the wavenumber (l) with $p_1=0.5$ for the MOL confinement. The magnitudes of the other physical parameters: $b=1$, $\gamma =0.05$, $k=0.84$, $G_1=-1$, $G_2=0.999999999$, $E=-2/9$.

4. Conclusions

In conclusion, we obtained a family of exact analytical solutions of the 1D eGPE for the generation of QDs in the binary BEC mixture with mutually symmetric spinor components in the presence of the MOL confinement. In the constructed model, the choice of the MOL confinement, which is a combination of the FOL and exponential periodic trap, results in the generation of various experimentally relevant trap profiles: OL, BOL, TOL, FOL, symmetric and asymmetric double-well superlattice, etc., and the corresponding exact wavefunction solution. From the constructed model, we identified the two parameters $p_1$ and $p_2$ for tuning the shape of the MOL confinement and revealed interesting potential symmetry with its tuning. Further, by taking the spatial periodic variation of the competing repulsive cubic MF and attractive quadratic BMF interactions, we illustrated the generation of interesting supersolid-like periodic, bi-periodic, and double-well superlattice density patterns in QDs under the BOL and MOL confinements. By tuning the symmetry of the BOL trap, we demonstrated two possible types of periodic lattice density patterns of the droplets: (a) at the central potential site (i.e., at $x=0$), the maxima of the periodic density lattice for $p_1<0$ and $p_2=0$ and (b) for $p_1>0$ and $p_2=0$, its minima $x=0$. Interestingly, the strength of ($p_1,p_2$) was identified as a key parameter for the fragmentation, compression, and inter-well transport of droplets. We observed an Anderson-like localization for $p_1,$$p_2>0$, i.e., the compression of the droplet, whereas the expansion of the droplet width was noted for $p_1>0$, $p_2<0$, forming the MOL, as well as in two-color ($p_2=0$) OL confinements. This can be attributed to the change in the subtle balance of the MF and BMF interactions, leading to a decrease in the flat top density profile and the observation of localization. Finally, the stability of the obtained droplet solutions was confirmed using the VK criterion and by the linear stability analysis.

The results of this paper can be tested experimentally by taking binary Bose gases with mass-balanced two-mixture components under the influence of BMF interactions in the presence of spatially varying MOL confinement with the present know-how. The shape of the MOL confinement is dependent on the power of the laser intensity and the wavelength of the laser light, which were connected with ($p_1,$$p_2$) in the present work. In particular, for typical values of the physical parameters, one can estimate the potential depths: $E_R=\frac{2{\pi }^2{\hslash }^2}{m{\lambda }^2}$. As the constructed model is based on a 1D geometry, it may be interesting to explore the impact in a 2D/3D confinement setup and droplet crossover to the 2D/3D configuration as future extensions of the present work.