Periodic Ultranarrow Rods as 1D Subwavelength Optical Lattices

TL;DR

This paper investigates the ground state properties of a one-dimensional Bose gas in a periodic multi-rod potential, providing semi-analytical solutions and analyzing band structures and interactions at zero temperature.

Contribution

It introduces a semi-analytical approach to solve the Gross-Pitaevskii equation for a 1D Bose gas in a periodic rod structure, including elliptic function solutions and band spectrum analysis.

Findings

Identification of energy band loops at strong interactions.

Dependence of ground state energy and chemical potential on potential height and interactions.

Prediction of subwavelength optical lattice band structures.

Abstract

We report on ground state properties of a one-dimensional, weakly-interacting Bose gas constrained by an infinite multi-rods periodic structure at zero temperature. We solve the stationary Gross-Pitaevskii equation (GPE) to obtain the Bloch wave functions from which we give a semi-analytical solution for the density profile, as well as for the phase of the wave function in terms of the Jacobi elliptic functions, and the incomplete elliptic integrals of the first, second and third kind. Then, we determine numerically the energy of the ground state, the chemical potential and the compressibility of the condensate and show their dependence on the potential height, as well as on the interaction between the bosons. We show the appearance of loops in the energy band spectrum of the system for strong enough interactions, which appear at the edges of the first Brillouin zone for odd bands and at the center for even bands. We apply our model to predict the energy band structure of the Bose gas in an optical lattice with subwavelength spatial structure. To discuss the density range of the validity of the GPE predictions, we calculate the ground state energies of the free Bose gas using the GPE, which we compare with the Lieb-Liniger exact energies.

Equations

$$i\hslash {\partial }_t\Psi (z,t)={\hat{H}}_{\mathrm{GP}}\Psi (z,t),$$

$$\Psi (z,t)=\Phi (z)e^{-i\mu t/\hslash },$$

$${\hat{H}}_{\mathrm{GP}}\Phi (z)=\mu \Phi (z),$$

$${\hat{H}}_{\mathrm{GP}}=-\frac{{\hslash }^2}{2m}{\partial }_z^2+V(z)+g|\Phi (z){|}^2,$$

$$\int |\Phi (z){|}^{2}dz=\mathcal{N},$$

$$V_{\mathrm{KP}}(z)=V_0\sum _{j=-\infty }^{\infty }\Theta [z-(j-1)l-a]\Theta [jl-z],$$

$${\int }_0^l|\Phi (z){|}^{2}dz=N,$$

$$\frac{E[\Phi ]}{N}=\frac{1}{N}{\int }_0^l{\Phi }^{\ast }(z)\left[-\frac{{\hslash }^2}{2m}{\partial }_z^2+V_{\mathrm{KP}}(z)+\frac{g}{2}|\Phi (z){|}^2\right]\Phi (z)dz.$$

$${\Phi }_k(z)=e^{ikz}{\varphi }_k(z),$$

$${\hat{H}}_{\mathrm{GP}}^k{\varphi }_k(z)={\mu }_k{\varphi }_k(z)$$

$${\hat{H}}_{\mathrm{GP}}^k=\frac{{\hslash }^2}{2m}{\left(-i{\partial }_z+k\right)}^2+V_{\mathrm{KP}}(z)+g|{\varphi }_k(z){|}^2$$

$${\varphi }_k(z)=\sqrt{n_1(z)}{e}^{iS(z)},$$

$${\partial }_{z}S(z)=-k+\frac{\alpha }{n_1(z)},$$

$$S(z)=S_0-kz+{\int }_0^z\frac{\alpha }{n_1(z^{\prime })}d{z}^{\prime },$$

$$-\frac{{\hslash }^2}{2m}$$

$$\left(\frac{{\hslash }^{2}k}{m}{\partial }_{z}S(z)+\frac{{\hslash }^2{k}^2}{2m}+V_{\mathrm{KP}}(z)\right)r(z)=\mu r(z).$$

$$-\frac{{\partial }^{2}r(z)}{\partial z^2}+\frac{{\alpha }^2}{r(z{)}^3}+\frac{2m}{{\hslash }^2}\left(V_{\mathrm{KP}}(z)-\mu \right)r(z)+\frac{2mg}{{\hslash }^2}r(z{)}^3=0.$$

$$-\frac{1}{2}{\left(\frac{\partial r(z)}{\partial z}\right)}^2-\frac{{\alpha }^2}{2r(z{)}^2}+\frac{m}{{\hslash }^2}\left(V_0-\mu \right)r(z{)}^2+\frac{mg}{2{\hslash }^2}r(z{)}^4-\sigma =0,$$

$${\left(\frac{d{n}_1}{dz}\right)}^2=\frac{4mg}{{\hslash }^2}{n}_1^3+\frac{8m}{{\hslash }^2}(V_0-\mu )n_1^2-8\sigma n_1-4{\alpha }^2.$$

$$n_1(z)=n_{\mathrm{off}}+4m_{\mathrm{j}}{\lambda }^2{\mathrm{sn}}^2\left(\sqrt{\frac{4mg}{{\hslash }^2}}\lambda (z-z_{\mathrm{off}})|m_{\mathrm{j}}\right),$$

$${\varphi }_k(z)=\{\begin{array}{ll}{\varphi }_k^{\mathrm{w}}(z)=\sqrt{n_1^{\mathrm{w}}(z)}{e}^{i{S}^{\mathrm{w}}(z)},& V_{\mathrm{KP}}(z)=0\\ {\varphi }_k^{\mathrm{b}}(z)=\sqrt{n_1^{\mathrm{b}}(z)}{e}^{i{S}^{\mathrm{b}}(z)},& V_{\mathrm{KP}}(z)=V_0\end{array}$$

$$n_1^{\mathrm{b}}(0)$$

$$S_0^{\mathrm{w}}-S_0^{\mathrm{b}}$$

$$n_1^{\mathrm{b}}(-b)$$

$$k(a+b)$$

$${\partial }_z{n}_1^{\mathrm{b}}(0^{-})$$

$${\alpha }^{\mathrm{b}}$$

$${\partial }_z{n}_1^{\mathrm{b}}(-b)={\partial }_z{n}_1^{\mathrm{w}}(a).$$

$$N(z)$$

$$\displaystyle=(n_{\mathrm{off}}+4\lambda^{2})z\;-\frac{4\lambda}{\sqrt{4m{g}/\hbar^{2}}}\Big{(}\mathcal{E}\left(u(z)|m_{\mathrm{j}}\right)-\mathcal{E}\left(u(0)|m_{\mathrm{j}}\right)\Big{)},$$

$$N={\int }_{-b}^0{n}_1^{\mathrm{b}}(z^{\prime })d{z}^{\prime }+{\int }_0^a{n}_1^{\mathrm{w}}(z^{\prime })d{z}^{\prime },$$

$$\int_{0}^{z}\frac{dz^{\prime}}{n_{1}(z^{\prime})}=\frac{1}{\sqrt{4m{g}/\hbar^{2}}\lambda n_{\mathrm{off}}}\;\times\Big{(}\Pi(n_{\mathrm{j}};u(z)|m_{\mathrm{j}})-\Pi(n_{\mathrm{j}};u(0)|m_{\mathrm{j}})\Big{)},$$

$${\mu }_k=\frac{1}{N}{\int }_0^l{\varphi }_k^{\ast }(z){\hat{H}}_{\mathrm{GP}}^k{\varphi }_k(z)dz.$$

Full text

Quantum Gases; Bose-Einstein Condensation; Bloch States; Gross-Pitaevskii Equation

I Introduction

Since the realization of a Bose-Einstein condensate (BEC) in 1995 bib:anderson_science.1995 ; bib:davis_prl.1995 , many new types of experiments have been proposed and realized involving BEC. In particular, we know that a successful research line has been to study BEC within periodic optical lattices created by superposition of two opposing laser lights in one, two or three mutually perpendicular directions bib:jaksch_prl-81.3108.1998 ; bib:bloch_nat-phys-1.23.2005 , such that the atomic gases are trapped in 3D multilayers, multitubes, or in a simple cubic array of dots, respectively. Atoms are trapped in one direction by the standing wave formed by two opposite laser lights whose effective potential acting on atoms has the generic sinusoidal form $V(x)=A\sin^{2}(2\pi\,x/\lambda)$, where $\lambda$ is the wavelength of the laser light and $A$ is the lattice potential height given in energy units. This ability to generate optical lattices of various types has become a fundamental tool to study the physics of bosonic or/and fermionic many atoms systems bib:gross_science.357.995.2017 . For instance, the superfluid-Mott insulator quantum phase transitions were experimentally probed “in a Bose–Einstein condensate with repulsive interactions, held in a three-dimensional optical lattice potential” bib:jaksch_prl-81.3108.1998 ; bib:greiner_nature.2002 , as well as in a 1D optical lattice bib:haller_nature.2010 . However, these kind of optical lattices are limited by its spatial resolution, which is of the order $\lambda$, to manipulate atoms. Fortunately, there has recently been a notorious interest and advances in developing tools to overcome the diffraction limit, arriving to the physical realization of subwavelength optical lattices (SWOL) of nearly $\delta$-function potential with ultranarrow barriers of width below $\lambda/50$ bib:lacki_prl-117.233001.2016 ; bib:wang_prl-120.083601.2018 ; reported results include, among others, the energy band structure. These SWOL can be seen as a very close experimental realization of the Dirac comb (DC) potential bib:merzbacher_qm.jws.2nd.1969 , as well as attractive setups to achieve $p$-wave superfluidity in a gas of fermions in 2D optical lattices bib:fedorov_pra.2017 .

From the theoretical point of view, the effect of a point-like, Dirac delta potential on the weakly-interacting Bose gas has been studied bib:seaman-pra-3.033609.2005 . Also, there has been substantial research on the properties of the Bose gas within a Dirac comb potential in the mean-field approximation, for instance, on the behavior of density profile of the condensate (both analytically and numerically), and how the interactions can have a profound impact on the energy spectrum, like the appearance of “swallow tails” in the band structure bib:theodorakis_jpa-13.4835.1997 ; bib:machholm_pra-67.053613.2003 ; bib:weidong-li_phys-rev-e-70.016605.2004 ; bib:seaman-pra-3.033622.2005 , as well as analysis on the influence of the periodic structure on the stability of superfluid currents bib:dong_lp-2.190.2007 ; bib:danshita_pra.2007 . Also, the DC potential has been extended to model interactions between atoms and ions forming a lattice bib:negretti_prb.2014 where both $s$-wave and $p$-wave scattering is present. However, to our knowledge, research on the properties of the Bose gas within the true Kronig-Penney (KP) potential bib:kronig-prsa-130.499.1931 with barriers and wells is missing, since all of the referred previous works employ the KP potential in the limit when barriers become exact Dirac deltas, i.e., the so-called DC potential.

In this work we study a one-dimensional weakly-interacting Bose gas within an infinite permeable multi-rods periodic structure which we use to discuss the ultranarrow rods limit as a 1D optical lattice with subwavelength spatial structure. The structure is modeled by a KP potential, which we analytically solve in the weak interaction regime where the GPE bib:gross_inc.1961 ; bib:pitaevskii_spjetp.1961 is applicable. The KP potential, i.e., periodic structure of well plus barrier, has the advantage that it is closer to the sinusoidal optical lattice potential than the Dirac comb, but at the same time it retains the simplicity to be solved analytically. We analyze the effects of the height and width of the barriers, as well as the interaction strength between bosons, on the ground state properties such as the density profile, the chemical potential, and the energy spectrum. Then we use our model to represent SWOL by means of very high, narrow barriers which retain their finite extent. In other words, although the barriers are very narrow compared to the optical lattice period, they have a nonzero width. The finite extent of the barriers holds even for future prospects of SWOL with period $\lambda/4$ bib:subhankar_arxiv.2019 .

On the other hand, the Bose gas within a multirods periodic structure is the same as the beautiful exact soluble Lieb-Liniger (LL) model bib:lieb_phys-rev-130.1605.1963 but within an external KP potential, which we analytically solve in the weakly interacting regime where the Gross-Pitaevskii equation is applicable. We show the effects of the KP potential on the ground state energy of the LL Bose gas, recovering the LL results in the weak interaction regime when we delete the KP potential. In order to establish the density regime of the GPE applicability, given an interaction magnitude, we use both ground state energies of the free Bose gas calculated exactly (LL case) and approximately (GP case), to fix the density regime where both energies are approximately equal. Correspondingly, for the case of the trapped Bose gas we give a lower estimate for the average linear density of $250$ bosons per potential period within the interaction interval used.

This work is developed in the following way. In Sec. II we present our model of a 1D Bose gas within permeable multi-rods; we establish the boundary conditions to obtain the constants on which the solutions depend. In Sec. III we give the ground state density profiles and we calculate the chemical potential and compressibility. In Sec. IV we calculate and plot the nonlinear energy band structure where the most remarkable thing is the appearance of energy loops, also known as swallow tails. In Sec. V we show the behavior of the density profile and the energy spectrum in the limit of very narrow barriers. We employ our model to predict the energy spectrum of a interacting Bose gas within an optical lattice with subwavelength spatial structure. Finally in the Sec. VI we give our conclusions.

II Bose gas within permeable multi-rods

We study a one-dimensional, weakly-interacting Bose gas constrained by a periodical structure composed of an infinite sequence of permeable rods of length $b$, separated a distance $a$; the rods repeat along the $z$ direction. We consider that the interactions between bosons are weak enough so that the physical properties of our system can be correctly described by the GPE bib:pitaevskii.oup.2016

[TABLE]

where $\Psi(z,t)$ is the wave function of the condensate. Since we are interested in the stationary states, i.e., those that evolve in time like

[TABLE]

where $\mu$ is the chemical potential of the system, Eq. (1) becomes the stationary Gross-Pitaevskii equation,

[TABLE]

where the l.h.s. operator $\hat{H}_{\mathrm{GP}}$ is the Gross-Pitaevskii time independent hamiltonian

[TABLE]

$m$ is the mass of the bosons, ${g}$ is the parameter that measures the strength of the interaction between particles, and $V(z)$ is the external potential. The stationary wave function ${\Phi}(z)$ is subject to the normalization condition

[TABLE]

where $\mathcal{N}$ is the number of bosons in the condensate.

The multi-rods structure is generated via an external Kronig-Penney (KP) potential bib:kronig-prsa-130.499.1931 $V(z)\equiv V_{\mathrm{KP}}(z)$. This periodic potential is an array of barriers of width $b$ separated by a distance $a$, each one with height $V_{0}$, see Fig. 1. The KP potential can be written as

[TABLE]

where $\Theta(z)$ is the Heaviside step function and $l\equiv a+b$ the potential period. For an infinite system like this, which repeats over and over, the normalization condition can be defined within a single period $l$, in the following way

[TABLE]

where $N$ is the average number of bosons in the condensate over a length equal to the potential period, such that the average linear density of the system becomes ${n}=N/l$. This condition fixes the value of the chemical potential of the system, since the average number of bosons remains constant. The energy per particle of the condensate can be defined in a similar way by

[TABLE]

Given the geometry of multi-rods, we can identify the potential period $l$ as a characteristic length which is either the distance between the midpoints of any two consecutive barriers, or the distance between the midpoints of two consecutive wells. This period $l$ is equal to that of an optical lattice produced by two counter-propagating lasers with wavelength $\lambda_{\mathrm{OL}}$ and wave number $k_{\mathrm{OL}}=2\pi/\lambda_{\mathrm{OL}}$, which is represented by an external potential $V_{\mathrm{OL}}(z)=sE_{\mathrm{R}}\sin^{2}(k_{\mathrm{OL}}z)$, with $s$ being the lattice height in recoil energy $E_{R}$ units, where $E_{\mathrm{R}}\equiv\hbar^{2}k_{\mathrm{OL}}^{2}/2m$. Since the period of the optical potential, i.e., the distance between two consecutive maximums, is $l_{\mathrm{OL}}=\lambda_{\mathrm{OL}}/2$, and hence $k_{\mathrm{OL}}=\pi/l_{\mathrm{OL}}$, the recoil energy of the optical lattice becomes $E_{\mathrm{R}}=\hbar^{2}\pi^{2}/2ml_{\mathrm{OL}}^{2}$ bib:bloch_nat-phys-1.23.2005 . Doing an analogy with the optical lattice, we can identify the recoil energy of our multi-rods system as $E_{\mathrm{R}}\equiv\hbar^{2}\pi^{2}/2ml^{2}$, that corresponds to the recoil energy of an optical lattice with the same period of the Kronig-Penney potential.

Since the density of our multi-rods system $|{\Phi}(z)|^{2}$ is a physically periodic function with the period of the KP potential, $\hat{H}_{\mathrm{GP}}=(-\hbar^{2}/2m)\partial^{2}_{z}+V_{\mathrm{KP}}(z)+g|{\Phi}(z)|^{2}$ should be invariant under translations by a distance $l$, i.e., $\hat{D}_{l}\hat{H}_{\mathrm{GP}}{\Phi}(z)=\hat{H}_{\mathrm{GP}}\hat{D}_{l}{\Phi}(z)=\mu\hat{D}_{l}{\Phi}(z)$, where $\hat{D}_{l}$ is the translation operator whose action is $\hat{D}_{l}f(z)=f(z+l)$ for any function $f(z)$. Therefore, the eigenstates of $\hat{D}_{l}$ are solutions of the GPE with chemical potential $\mu$, i.e., they can be written as Bloch waves:

[TABLE]

where the function $\phi_{k}(z)$ has the same periodicity as the potential, i.e., $\phi_{k}(z+l)=\phi_{k}(z)$, and $\hbar k$ is the quasi-momentum of bosons in the condensate. Introducing (9) in (3) we obtain

[TABLE]

with

[TABLE]

the “shifted” hamiltonian bib:kramer.epjd-3.247.2003 .

Each value of the lattice wave number $k$ fixes a solution for Eq. (10). To solve it with the corresponding boundary conditions we express the function $\phi_{k}(z)$ in its complex form

[TABLE]

where the function $S(z)$ represents the phase and $n_{1}(z)\equiv|\phi_{k}(z)|^{2}=|{\Phi}_{k}(z)|^{2}$ is the particle number density as a function of $z$.

Substituting (12) in Eq. (10), we arrive to a pair of coupled differential equations for the real and imaginary parts in terms of $n_{1}(z)$ and ${S}(z)$. The equation for the phase is a first order differential equation,

[TABLE]

where $\alpha$ is a constant of integration. Equation (13) is easily integrated through separation of variables, such that the phase is given by

[TABLE]

with $S_{0}$ a constant of integration. In order to find the phase we require the solution for the density $n_{1}(z)$, which comes from the real part of the GPE (10):

[TABLE]

Here, $r(z)=\sqrt{n_{1}{(z)}}$. We substitute Eq. (13) in Eq. (II), and after rearranging terms we obtain

[TABLE]

To proceed further, we take into account that the Kronig-Penney is a piecewise potential with a constant magnitude in each barrier or well region, i.e., $V_{\mathrm{KP}}(z)=V_{0}$ within the barrier and $V_{\mathrm{KP}}(z)=0$ inside the well. Let’s focus on the barriers region. We multiply Eq. (16) by $\partial_{z}r(z)$ and integrate the equation, arriving to

[TABLE]

where $\sigma$ is a second constant of integration. Finally, we multiply this equation by $-8r(z)^{2}$; after some algebraic steps, we obtain the corresponding differential equation for the density:

[TABLE]

An analog procedure for the wells region gives a similar equation, but with $V_{0}=0$, so the differential equation (18) changes only in the quadratic term.

The ODE (18) has a set of analytical solutions given by the Jacobi elliptic functions bib:abramowitz-stegun.1964 . The explicit form of the density is

[TABLE]

where the function $\mathrm{sn}(u|m_{\mathrm{j}})$ is the Jacobi elliptic sine. The factor $m_{\mathrm{j}}$ is a real number known as the elliptic modulus, $n_{\mathrm{off}}$ is a constant offset on the value of $n_{1}(z)$, while $\lambda$ is a parameter that fixes the amplitude of spatial density variations. Equation (19) defines a whole family of functions whose properties are deeply linked to the value of $m_{\mathrm{j}}$, the value of ${g}$, and $\lambda$.

II.1 Boundary conditions

The solution for the density (19) assumes that the potential magnitude remains constant over the interval of $z$ being evaluated. Hence, we have two solutions for $\phi_{k}(z)$: one within the “well” (w) regions and another within the “barrier” (b) regions. Then,

[TABLE]

We consider also that: both functions must match in a smooth way at the interface of each potential barrier, the periodic nature of the potential and that the system is infinite. In Fig. 2 we show the behavior of the function $\phi_{k}(z)$ in both well and barrier regions. Within the barriers $\phi_{k}^{\tiny\mathrm{b}}(z)$ has a depletion which is complemented by the accretion in $\phi_{k}^{\tiny\mathrm{w}}(z)$ within the wells.

A direct consequence of the nature of the potential is that we have a density function and a set of parameters $n_{\mathrm{off}}$, $m_{\mathrm{j}}$, $\lambda$ and $z_{\mathrm{off}}$, as well as a phase function and parameters $S_{0}$ and $\alpha$, for each region. These sets are related by the boundary conditions imposed on the system. We can exploit the periodicity of $\phi_{k}(z)$ to focus our analysis to a single period of the system fixed at the origin $z=0$, which extends from $z=-b$ to $z=a$. In this picture the edge of the barrier is located at the origin, therefore the boundary condition for continuity is $\phi_{k}^{\tiny\mathrm{b}}(0)=\phi_{k}^{\tiny\mathrm{w}}(0)$. This equality results in the conditions

[TABLE]

The first of these equations forces the density to be continuous at $z=0$, while the second states that the difference of phase between regions is discrete and equal to an integer multiple of $2\pi$. The periodicity of $\phi_{k}(z)$ can be stated as $\phi_{k}^{\tiny\mathrm{b}}(-b)=\phi_{k}^{\tiny\mathrm{w}}(a)$, which in turn implies that

[TABLE]

with $\alpha^{\mathrm{w}}$ in the wells and $\alpha^{\mathrm{b}}$ in the barriers. In addition, the derivative of $\phi_{k}(z)$ must be continuous at $z=0$, i.e., $\partial_{z}\phi_{k}^{\tiny\mathrm{b}}(0^{-})=\partial_{z}\phi_{k}^{\tiny\mathrm{w}}(0^{+})$, which is equivalent to

[TABLE]

Finally, the derivative of $\phi_{k}(z)$ must be also periodic, hence

[TABLE]

Conditions (21) to (27), along with the normalization (7) define the complete set of solutions for the wave function of the condensate $\phi_{k}(z)$.

The definition and properties of the Jacobi elliptic functions permit us to express the normalization (7) in closed-form. For this, we can define the integral

[TABLE]

as the average number of particles contained in the interval $[0,z]$. The function $\mathcal{E}(u(z)|m_{\mathrm{j}})$ is the incomplete elliptic integral of the second kind with argument $u(z)=\sqrt{4m{g}/\hbar^{2}}\lambda(z-z_{\mathrm{off}})$ expressed in canonical form accordingly to bib:abramowitz-stegun.1964 . Then the normalization condition becomes

[TABLE]

where each of the integrals can be evaluated using the equation (II.1). Analogously, the integral in (14) becomes

[TABLE]

where $\Pi(n_{\mathrm{j}};u(z)|m_{\mathrm{j}})$ is the incomplete elliptic integral of the third kind of order $n_{\mathrm{j}}=-4m_{\mathrm{j}}\lambda^{2}/n_{\mathrm{off}}$. Then the quasi-momentum $k$ in (24) can be expressed in terms of (30).

III Ground state density profile and chemical potential

The periodic structure as well as the interaction between particles have notorious effects on the density of the condensate, even in the regime where there is no relative velocity of the gas with respect to the potential frame.

In Figs. 3a and 3b we show the behavior of the density as a function of the position $z$ for two square lattices with potential barrier height $V_{0}=5E_{\mathrm{R}}$ and $V_{0}=25E_{\mathrm{R}}$, respectively. The geometric ratio of the potential, defined as $r=b/a$, is equal to unity since $b=a$. In both plots we have calculated the density profile for several values of a repulsive two-body interaction strength ${g}{n}$, including the interactionless Bose gas (${g}{n}=0$) affected only by the periodic potential. We show the density profile for ${g}{n}=0.5,1,2,3$ times $V_{0}$. For all cases the density has a maximum value at the midpoint of the well region, and a minimum at the midpoint of the barrier region. As the repulsive interaction between particles increases, we observe that the density variations in a spatial period diminish in such a way that the average value of the density approaches to unity. This occurs because the interaction between particles dominate the repulsive effect of the potential barriers, reducing the particle localization in the wells. An opposite effect appears when we increase the barrier height, keeping constant the interaction between particles, then the density profile raises in the well regions while diminish inside the barriers. Also, the geometry of the lattice has a significant influence on the density profile, as it is shown in the Figs. 3c and 3d, where we have calculated $n_{1}(z)$ for two different nonsquare lattices with geometric ratios $b=0.5a$ and $b=2a$, respectively, while keeping constant the potential height $V_{0}=25E_{\mathrm{R}}$. The results show a greater particle localization in the wells when the potential barriers becomes wider, since the barrier repulsion dominate over the repulsive interactions between particles. Eventually, as the repulsive interaction increases, the density profile becomes flat.

We can obtain a relationship between the energy (8) and the chemical potential directly from the Gross-Pitaevskii equation. Multiplying both sides of (10) by $\phi_{k}^{*}(z)$ and integrating over a potential period, we arrive to

[TABLE]

We recognize that the only difference between the total energy and chemical potential arises from the nonlinear term. It follows that $E_{k}$ and $\mu_{k}$ are related by

[TABLE]

where $\tilde{n}_{1k}(\tilde{z})=\left|\phi_{k}(\tilde{z})\right|^{2}/n$, with $\tilde{z}=z/(a+b)$. The second term of the r.h.s. of (32) accounts for the two-body interactions in the gas. It vanishes for ${g}=0$, therefore, $\mu_{k}=E_{k}/N$, which is the exact result for the energy per particle of a noninteracting Bose gas at zero temperature. In this case, the chemical potential can be obtained using the dispersion relation of an ideal Bose gas subject to a Kronig-Penney potential bib:rodriguez.jltp-3.144.2016 .

We calculated the chemical potential and the compressibility for the state $k=0$ which are shown in Fig. 4. We observe that the chemical potential is a monotonic, increasing function of the interaction parameter ${g}{n}$ (Fig. 4a). In the limit when the interaction goes to zero, $\mu$ tends correctly to the value of the ideal Bose gas, which in general is nonzero due to the presence of the lattice: the repulsive effect of the barriers raises the chemical potential, since we require more energy to add a single particle to the system. Obviously, larger potentials raise $\mu$ even further. When the interaction is strong enough so that the kinetic energy becomes small compared to the potential energy, we can obtain a closed-form formula for $\mu$. Neglecting the kinetic energy term in Eq. (11) and using the normalization condition (7), it follows that

[TABLE]

which is the chemical potential in the so-called Thomas-Fermi limit. It basically implies that, when the interaction strength is strong enough, the chemical potential increases as the corresponding one of the free, but interacting, Bose gas, plus a shift due to the external potential. When the potential $V_{0}$ vanishes, Eq. (33) reduces to the case of the free, interacting Bose gas. However, Eq. (33) is not valid anymore when $V_{0}$ approximates to ${g}{n}$, since the spatial variations of $\phi_{k}(z)$ grow and the kinetic term (proportional to $|\partial_{z}\phi_{k}(z)|$) of the Gross-Pitaevskii equation becomes significant. We can see this change in the Fig. 4c, where the linear dependence of $\mu$ on $V_{0}$ is lost as ${g}{n}$ becomes smaller.

The compressibility of the gas $\kappa$ is related to the chemical potential by the relation

[TABLE]

where we implicitly assume that $\mu$ is calculated for a specific Bloch state with a fixed momentum $\hbar k$. For a free, homogeneous interacting Bose gas, the chemical potential is $\mu={g}{n}$ bib:pitaevskii.oup.2016 , so the inverse of the compressibility is $\kappa^{-1}={g}{n}$. In general, the inverse compressibility will grow at the same rate, i.e., $\kappa^{-1}\approx{g}{n}$, in both ${g}{n}\to 0$ and ${g}{n}\gg V_{0}$ limits. In the latter case, the compressibility resembles the one of the free gas because the relatively large interactions screen out the effects of the external potential (Fig. 4b). For a nonzero lattice height in an intermediate range of ${g}{n}$ the compressibility will deviate from the free Bose gas behavior, since the presence of the external potential reduces the compressibility of the gas due to the repulsive nature of the barriers (Fig. 4d).

IV Nonlinear energy band structure

The presence of the nonlinear term in the Gross-Pitaevskii equation (10) alters the band spectrum of the non-interacting case in striking ways. The most notable phenomena is the appearance of energy loops the so-called “swallow tails” bib:machholm_pra-67.053613.2003 . Figure 5 shows the nonlinear band spectrum for a square lattice $b=a$ with potential height $V_{0}=2E_{\mathrm{R}}$. In addition, the lattice contain a strongly repulsive condensate where ${g}{n}=4E_{\mathrm{R}}$. The characteristic shape of swallow tails is readily visible.

The energy loops belong to a specific band, and appear in a regular way: at the end of the first Brillouin zone for odd bands, or at the center for even bands. They become larger as the repulsive interaction magnitude increases respect to the lattice height. The swallow tails emerge because, as the interaction increases, two states appear that share the same crystal momentum but different energies. Both of these states $\phi_{k}(z)$ are minimizers of the energy functional (8) subject to the normalization condition (5) for a fixed chemical potential bib:rogel-ejp.2013 . The origin of the swallow tails is a consequence of the change of the the energy landscape of the Bose gas, i.e., the shape of $E[{\Phi}]$ as a function of ${\Phi}(z)$. The nonlinear term in the Gross-Pitaevskii results in a situation where the system has more than one state that minimize the energy functional (8) bib:mueller-pra.6.063603.2002 . The appearance of two local minima implies (only by physical considerations) the existence of a local maximum in the landscape between both minima. The state with maximum energy will lie in the upper portion of the swallow tail (for a fixed momentum $k$), while the remaining state will lie in the lower portion as it has less energy.

The appearance of the swallow tails depends mainly on the ratio $V_{0}/{g}{n}$, and, at a lesser extent, on the geometric ratio $b/a$. For optical lattices the swallow tail for the first band appears when the interaction factor ${g}{n}$ becomes equal (or greater) than the lattice height bib:machholm_pra-67.053613.2003 . For upper bands there is not a similar analytical relation.

Figure 6a and  6b show the density profiles of some states in the first an second energy bands (within the first Brillouin zone), respectively, for a square lattice ($b=a$) with height $V_{0}=2E_{\mathrm{R}}$ and interaction strength ${g}{n}=4E_{\mathrm{R}}$. The energy spectrum for this system is the one shown in Fig. 5. The states plotted in Fig. 6a lie in the first band, in the lower part of it. The solid line represents the ground state $k=0$, while the rest of the curves have increasing values of the momentum, up to $k(a+b)=\pi$. As we can see, the maximum of the density at the origin decreases as the momentum increases, while the density reduces in the midpoint of the potential barriers. This corresponds to a greater kinetic energy as the density spatial variations are magnified. The curves in the Fig. 6b show the density variations for a set of states that lie in the lower part of the second energy band. In this case the variations of the density get reduced as the momentum grows, from an initial state at the left edge of the Brillouin zone, with $k(a+b)=-\pi$. This state has the peculiarity that it becomes zero at the origin, in the midpoint between two potential barriers. The periodicity of the Bloch states implies that there are zero density surfaces at $z=j(a+b)$, with $j$ an integer. This states form an array of so-called “dark solitons” bib:machholm_pra-67.053613.2003 ; bib:tsuzuki_jltp.4.441.1971 .

V Subwavelength optical lattices as experimental KP potential

realizations

In recent years, the experimental realization of optical lattices with subwavelength spatial structure bib:wang_prl-120.083601.2018 has open a way to study the physics of quantum manybody fluids subject to periodic potentials that closely resemble the well known Kronig-Penney potential in the Dirac $\delta$-function limit bib:merzbacher_qm.jws.2nd.1969 . In this approximation, the width of the potential barriers $b$ goes to zero, and the potential magnitude $V_{0}$ goes to infinity, but the product $V_{0}b$ remains constant. Then, the KP potential becomes a succession of $\delta$-functions centered in the positions $jl$, being $j$ an integer and $l$ the KP potential period, as well as the separation between two contiguous deltas. The expression (6) for $V_{\mathrm{KP}}(z)$ becomes the Dirac-comb potential,

[TABLE]

where the finite, constant value $V_{0}b$ is the area below a single barrier of the KP potential. In the context of the Dirac comb potential, it can be interpreted as a measure of the impermeability (the strength) of a single Dirac delta barrier. When $V_{0}$ becomes zero we recover the homogeneous, free interacting Bose gas. On the opposite side, as the delta strength becomes larger the system resembles more a succession of independent wells of infinite walls, each one having a width $l$. Since $V_{0}b$ is an energy times a length, we can redefine the delta strength as $V_{0}b=sE_{\mathrm{R}}l$. Although this definition seems somewhat arbitrary, it is very useful when we write Eq. (35) in terms of the dimensionless length $z^{\prime}=z/l$,

[TABLE]

where we used the scaling property of the delta function, $\delta(lz)=\delta(z)/l$. Then, in Eq. (36) the factor $s$ represents the dimensionless strength of the potential in $E_{\mathrm{R}}$ units, and relates the parameters of our multi-rods model with the Dirac-comb potential parameters. As an ubiquitous potential, the stationary states of mean-field BEC subject to an external, Dirac-comb potential (36) have been previously studied bib:theodorakis_jpa-13.4835.1997 ; bib:weidong-li_phys-rev-e-70.016605.2004 ; bib:seaman-pra-3.033609.2005 ; bib:dong_lp-2.190.2007 .

In Fig. 7a we present the density profile of the Bose gas as the Kronig-Penney potential approaches to the Dirac-comb potential. We can see that the shape of the density profile becomes flatter in the middle of the lattice cell as the interaction strength increases. There is a sharp change around the edged of the cell, where the Dirac deltas are located. Results are very similar for both systems, even when they have very different values of $V_{0}$ and $b$. Our results confirm that, for potentials equal or greater than $V_{0}=20E_{\mathrm{R}}$, and for ratios of the order or smaller than $b/a=0.05$, the multi-rods potential is a very close representation of the Dirac-comb potential with strength $s=1$. The calculated density profiles are in very good agreement with those in Fig. 7 of bib:seaman-pra-3.033622.2005 . Naturally, for a stronger Dirac comb potential, i.e., for a greater value of $s$, the established thresholds for $V_{0}$ and $b$, such that a potential barrier represents a delta barrier, will change. In Fig. 7b we show the energy spectrum for a relatively large interaction parameter ${g}{n}=10E_{\mathrm{R}}$. The swallow tails are notorious, and appear in all the plotted bands. The difference with the spectrum of the ideal gas is complete.

Based on the previous analysis of the Dirac comb potential, we attemp to model the optical potential bib:wang_prl-120.083601.2018

[TABLE]

in the limit when $\epsilon\ll 1$. Under this condition, $V_{\mathrm{OL}}(z)$ becomes a lattice of narrow barriers spaced by (a period) $\lambda_{\mathrm{OL}}/2$, with a peak value of $E_{\mathrm{R}}/\epsilon^{2}$ and width at half maximum scaling of $\Delta=\epsilon\lambda_{\mathrm{OL}}/2\pi\ll\lambda_{\mathrm{OL}}$ bib:lacki_prl-117.233001.2016 . For $\epsilon\ll 1$, this potential has a subwavelength spatial structure that is a very close approximation of the Dirac comb potential (35) with strength $s=1/(2\epsilon)$. We realize our analysis by fixing the KP barrier width as $b=0.05a$, and equating the KP potential period with the optical lattice period, i.e., $l=\lambda_{\mathrm{OL}}/2$. Then, to get the corresponding barrier height $V_{0}$, we equate the barrier area $V_{0}b$ with the area of $V_{\mathrm{OL}}(z)$ over a potential period for $\epsilon=0.14$. This procedure results in a barrier height of $V_{0}=74.3E_{\mathrm{R}}$, somewhat larger than the peak of $V_{\mathrm{OL}}(z)$, which is $51.02E_{\mathrm{R}}$.

First, we have calculated the predicted nonlinear band structure of the condensate for the system with ${g}{n}=10E_{\mathrm{R}}$, whose results are plotted in Fig. 8. This is a relatively strong condensate, with an energy spectrum that significantly deviates from the ideal one. Swallow tails are not significant in the spectrum, however.

V.1 Ground state energy of the Lieb-Liniger Bose gas within a subwavelength lattice

The extent to which a 1D Bose gas with interactions of the order of $gn=20E_{R}$ or less, is well modeled by the GPE will depend not only on the value of $g$, but also on the average linear density of the gas ${n}$. Unlike the 3D Bose gas, the weakly interacting regime for a 1D Bose gas corresponds to a high average density of the gas, while low densities correspond to the strongly interacting regime bib:pitaevskii.oup.2016 ; bib:dunjko_prl.86.5413.2001 . The exact description of a 1D Bose gas with contact-like, repulsive interactions is given by Lieb-Liniger (LL) model bib:lieb_phys-rev-130.1605.1963 and the LL parameter $\gamma=m{g}/\hbar^{2}{n}$, which must satisfy $\gamma\ll 1$ in order to the GPE picture to be valid. We have $\gamma={(nl/\pi)}^{-2}/2\times{g}{n}/E_{\mathrm{R}}$. Figure 9 shows the ground state energy of the gas predicted by the LL theory as a function of ${g}{n}$, for average densities of ${n}=500$, $250$, $100$ and $10$ times $l^{-1}$. We compare these results with the energy per particle ${g}{n}/2$ predicted by the GPE for the homogeneous Bose gas. Results show that only in the high density case ${n}=500l^{-1}$ the GPE gives accurate results over the full interval of interaction ${g}{n}$. Higher densities should provide results even more accurate. For ${n}=100l^{-1}$ we can see small discrepancies between the GPE predictions and the LL theory. Numerical results indicate that for ${n}=250l^{-1}$ the GPE and LL exact results are almost identical, and it can be taken as a lower limit for suitable densities; greater values of ${n}$ are within the range of typical experimental densities for Bose gases in 1D regime bib:mewes_prl.77.416.1996 ; bib:kruger_prl.105.265302.2010 . In the presence of the subwavelength optical potential the ground state energy of the Bose gas is shown in Fig. 9 by the solid, red line. The lattice raises the energy with respect to the free Bose gas, and its dependence is not linear with ${g}{n}$. Then, the observed energy band structure of the gas for ${g}{n}=10E_{\mathrm{R}}$ should be similar to the one shown in Fig. 8.

VI Conclusions

We have studied a 1D interacting Bose gas at zero temperature subject to a periodic, multi-rods potential by quasi-analytically solving the Gross-Pitaevskii equation. We model the periodic structure using a Kronig-Penney potential, which has the remarkable property of having analytical solutions for the wave function of the condensate. In this work, we focused on the Bloch state type solutions of the GP equation.

We were able to find analytical expressions for the wave function of the condensate. The density profile, the normalization condition, and the complex phase can be expressed in terms of the Jacobi elliptic functions, as well as in terms of the incomplete elliptic integrals of the first, second and third kind. We have obtained the density profile, the chemical potential and the energy spectrum of an interacting trapped Bose gas. The energy spectrum consists of bands separated by prohibited regions, like the spectrum of the trapped ideal gas. However, we found that the nonlinear spectrum may strongly differ from the ideal one since the first one shows loops, or “swallow tails”, at the edges of the first Brillouin zone for odd bands and at the center for even bands. No analytical expression for a threshold of the appearance of swallowtails was obtained so that further research in this subject is required. We obtained the chemical potential as well as the compressibility, numerically, as functions of the potential magnitude $V_{0}$ and the interaction parameter ${g}{n}$. When the interaction between particles is relatively large compared to the lattice height, our results agree satisfactorily with those predicted by the Thomas-Fermi approximation, where we have obtained closed-form expressions for the chemical potential and the compressibility.

Our periodic lattice becomes the Dirac comb potential when the potential barriers become very high and very thin, but the area below them remains finite. Although we have not used a mathematical strict expression for a Dirac $\delta$-potential, we found that for barriers as thin as $b=0.05a$ the KP potential is a very good approximation of the Dirac comb potential. In this case, we can reproduce the results for the density profile and the energy spectrum done in previous studies for the Bose gas in Dirac comb potentials. Moreover, because we have full control of the width and height of the potential barriers, we have employed our periodic ultranarrow rods model to predict the energy band structure of an interacting 1D BEC in optical lattices with subwavelength spatial structure. By comparison of the ground state energies of the free Bose gas calculated exactly (LL case) and approximately (GP case), for the trapped Bose gas, we give a lower estimate for the average linear density of $250$ bosons per potential period, such that the Gross-Pitaevskii equation results should be considered accurate in the full interaction interval shown in Fig. 9.

We acknowledge partial support from grants PAPIIT-DGAPA-UNAM IN-107616 and IN-110319, and CONACyT 221030.