# Fatigue failure of amorphous alloys under cyclic shear deformation

Nikolai V. Priezjev1,2 1Department of Civil and Environmental Engineering,
Howard University, Washington, D.C. 20059  2Department of Mechanical and
Materials Engineering, Wright State University, Dayton, OH 45435

###### Abstract

The accumulation of plastic deformation and flow localization in amorphous
alloys under periodic shear are investigated using molecular dynamics
simulations. We study a well-annealed binary mixture of one million atoms
subjected to oscillatory shear deformation with strain amplitudes slightly
above a critical value. We find that upon approaching a critical strain
amplitude from above, the number of shear cycles until the yielding transition
is well described by a power-law function. Remarkably, the potential energy at
the end of each cycle as a function of the normalized number of cycles is
nearly independent of the strain amplitude, which allows for estimation of the
fatigue lifetime at a given strain amplitude. The analysis on nonaffine
displacements of atoms elucidates the process of strain localization,
including irreversible rearrangements of small clusters until the formation of
a system-spanning shear band.

Keywords: metallic glasses, fatigue, yielding transition, cyclic loading,
molecular dynamics simulations

##  I Introduction

The prediction of stability and lifetime of amorphous alloys under repeated
stress or strain deformation is important for various structural applications
Ashby06 ; Sha2020 . Although multicomponent alloys like metallic glasses
possess a number of advantageous properties, such as high strength and large
elastic strain limit, their resistance to fatigue damage is relatively poor
Egami13 ; Ritchie99 ; Menzel06 ; Tanaka22 . The failure mechanism in metallic
glasses involves the formation of nanoscale shear bands where plastic strain
becomes strongly localized, which in turn might lead to propagation of
microscale cracks Liaw05 ; Wang17 ; Hufnagel16 ; Shi19 . At the atomic level,
the elementary plastic deformation in amorphous solids consists of rapid
rearrangement of a small cluster of particles or shear transformation
Spaepen77 ; Argon79 . Notably, the results of numerical simulations of the
fibre bundle model have shown that the fatigue failure under repeated loading
of heterogeneous materials occurs after a number of cycles, and the fatigue
lifetime has a power-law dependence on the loading amplitude Herrmann07 . More
recently, using two models of elastoplastic rheology, it was demonstrated that
cyclically sheared amorphous materials initially accumulate low levels of
damage in the form of spatial strain heterogeneity, which is followed by a
sudden catastrophic material failure via shear band formation Fielding22 .
However, in spite of the considerable modeling and experimental efforts, the
precise determination of the critical loading amplitude and fatigue lifetime
remains a challenging problem.

During the last decade, the effect of cyclic loading on the yielding
transition, structural relaxation, and flow localization in amorphous
materials was extensively studied using atomistic simulations Priezjev13 ;
Reichhardt13 ; Sastry13 ; IdoNature15 ; GaoNano15 ; Priezjev16 ; Kawasaki16 ;
Priezjev16a ; Sastry17 ; Priezjev17 ; Priezjev18 ; Priezjev18a ; Sastry19band
; PriezSHALT19 ; Priez20ba ; KawBer20 ; NVP20altY ; Priez20delay ; BhaSastry21
; Priez21var ; PriezCMS21 ; Peng22 ; Procaccia22 ; PriezJNCS22 ; Maloney22 ;
Barrat23 . Interestingly, it was demonstrated that in the athermal limit,
amorphous solids under small-amplitude oscillatory shear evolve into the so-
called limit cycles, where trajectories of atoms become exactly reversible
after one or more periods, and the number of cycles to reach periodic behavior
diverges upon approaching a critical strain amplitude from below IdoNature15 ;
Reichhardt22 . On the other hand, periodic deformation at strain amplitudes
above a critical value leads to yielding and flow localization after a number
of cycles GaoNano15 ; Sastry17 ; Priezjev17 ; Sastry19band ; NVP20altY . In
general, the number of cycles until the yielding transition depends on the
degree of annealing, temperature, system size, strain amplitude and frequency.
In particular, it was found that the number of cycles to failure increases
upon increasing frequency GaoNano15  or glass stability Priez20ba ;
Priez20delay , and by decreasing strain amplitude towards a critical value
GaoNano15 ; Priezjev17 . In addition, the number of fatigue cycles can be
reduced by periodically alternating shear orientation in two or three spatial
dimensions NVP20altY  or by occasionally increasing strain amplitude above a
critical value Priez21var . Despite recent progress, however, the processes of
damage accumulation and formation of shear bands during cyclic loading near a
critical strain amplitude remain not fully understood.

In this paper, the influence of repeated shear strain on plastic deformation
and yielding transition in a disordered solid is studied via molecular
dynamics (MD) simulations. We consider a well-annealed binary glass subjected
to oscillatory shear deformation at strain amplitudes slightly above a
critical value. It will be shown that the number of shear cycles to reach the
yielding transition increases approximately as a power-law function when the
strain amplitude approaches the critical value. Moreover, we find that the
potential energy at zero strain for different strain amplitudes is well
described by a single function of the normalized number of cycles. In turn,
the appearance of local plastic events and the formation of a shear band at
the yielding transition are quantified via the fraction of atoms with large
nonaffine displacements during one shear cycle.

The rest of this paper is organized as follows. The details of molecular
dynamics simulations as well as the oscillatory shear deformation protocol are
described in the next section. The analysis of shear stress, potential energy,
and nonaffine displacements is presented in section III. A brief summary is
provided in the last section.

##  II Molecular dynamics simulations

In this study, the amorphous alloy was modeled via the standard Kob-Andersen
(KA) binary mixture composed of 80 % of atoms of type A and 20 % of type B
KobAnd95 . The total number of atoms is $10^{6}$. In this model, the
interaction between atoms of types $\alpha,\beta=A,B$ is defined via the
Lennard-Jones (LJ) potential:

$V_{\alpha\beta}(r)=4\,\varepsilon_{\alpha\beta}\,\Big{[}\Big{(}\frac{\sigma_{\alpha\beta}}{r}\Big{)}^{12}\\!-\Big{(}\frac{\sigma_{\alpha\beta}}{r}\Big{)}^{6}\,\Big{]},$
(1)

where the parameters are set to $\varepsilon_{AA}=1.0$,
$\varepsilon_{AB}=1.5$, $\varepsilon_{BB}=0.5$, $\sigma_{AA}=1.0$,
$\sigma_{AB}=0.8$, $\sigma_{BB}=0.88$, and $m_{A}=m_{B}$ KobAnd95 . A similar
parametrization was used by Weber and Stillinger to study structure and
dynamics of the amorphous metal-metalloid alloy $\text{Ni}_{80}\text{P}_{20}$
Weber85 . All physical quantities are reported in the units of length, mass,
energy, and time, as follows: $\sigma=\sigma_{AA}$, $m=m_{A}$,
$\varepsilon=\varepsilon_{AA}$, and $\tau=\sigma\sqrt{m/\varepsilon}$. The MD
simulations were carried out using the LAMMPS parallel code with the
integration time step $\triangle t_{MD}=0.005\,\tau$ and the cutoff radius
$r_{c}=2.5\,\sigma$ Allen87 ; Lammps .

The sample preparation procedure and the deformation protocol are similar to
the ones reported in the previous MD study Priez20ba . More specifically, the
binary mixture was first placed in a cubic box of linear size
$L=94.10\,\sigma$ and equilibrated at the temperature
$T_{LJ}=1.0\,\varepsilon/k_{B}$ and density
$\rho=\rho_{A}+\rho_{B}=1.2\,\sigma^{-3}$ using the Nosé-Hoover thermostat and
periodic boundary conditions Allen87 ; Lammps . For reference, the critical
temperature of the KA model at the density $\rho=1.2\,\sigma^{-3}$ is
$T_{g}=0.435\,\varepsilon/k_{B}$ KobAnd95 . Then, the sample was cooled with
computationally slow rate of $10^{-5}\varepsilon/k_{B}\tau$ from
$T_{LJ}=1.0\,\varepsilon/k_{B}$ to $0.01\,\varepsilon/k_{B}$ at constant
density $\rho=1.2\,\sigma^{-3}$. Right after cooling, the glass was subjected
to oscillatory shear deformation along the $xz$ plane, as follows:

$\gamma_{xz}(t)=\gamma_{0}\,\text{sin}(2\pi t/T),$ (2)

where $\gamma_{0}$ is the strain amplitude and $T=5000\,\tau$ is the
oscillation period. The simulations were performed for strain amplitudes
$0.069\leqslant\gamma_{0}\leqslant 0.075$ at $T_{LJ}=0.01\,\varepsilon/k_{B}$
and $\rho=1.2\,\sigma^{-3}$. The results for the potential energy, shear
stress, and nonaffine displacements of atoms are reported only for one
realization of disorder because of the considerable computational burden. As
an example, it took about 36 days to simulate 800 shear cycles at the strain
amplitude $\gamma_{0}=0.069$ using 400 processors in parallel.

##  III Results

Recent studies have shown that model glasses prepared by thermal annealing can
yield after a certain number of cycles at strain amplitudes that are smaller
than the yielding strain during uniform shear deformation Sastry17 ; Priez20ba
. The precise value of the critical strain amplitude is difficult to determine
numerically due to a large number of cycles needed to reach the yielding
transition. Within the range of about three thousand cycles, it was found that
rapidly quenched binary glasses under cyclic loading yield at the critical
strain amplitude $\gamma_{0}=0.067$, regardless of whether shear is applied
along a single plane or periodically alternated in two or three spatial
dimensions NVP20altY . In the present study, we consider a relatively large
system of one million atoms and subject a well-annealed KA glass to
oscillatory shear deformation at strain amplitudes slightly above the critical
value.

We first report the variation of shear stress along the $xz$ plane as a
function of time in Fig. 1 for two values of the strain amplitude, i.e.,
$\gamma_{0}=0.072$ and $0.075$. It can be seen that in both cases, the
amplitude of shear stress oscillations slightly decreases upon continued
loading until a sudden drop during one shear cycle. Notice that the number of
cycles until yielding becomes greater upon decreasing strain amplitude.
Specifically, the yielding transition occurs during 218-th cycle for
$\gamma_{0}=0.072$ and during 56-th cycle for $\gamma_{0}=0.075$. By contrast,
after the yielding transition, the maximum shear stress is determined by the
plastic flow within a shear band. These results are consistent with those
previously reported for a well-annealed binary glass that was periodically
deformed for only 40 shear cycles at larger strain amplitudes Priezjev17 .

Along with shear stress, we plot in Fig. 2 the time dependence of the
potential energy for the same strain amplitudes, $\gamma_{0}=0.072$ and
$0.075$, as in Fig. 1. It is evident that the yielding transition is
associated with an abrupt increase of the potential energy due to the
formation of a shear band across the system. It should be emphasized that for
each strain amplitude, the sudden change in shear stress and potential energy
occur at the same cycle number. One can further realize that before yielding,
the cyclic shear deformation results in a slow accumulation of plastic events.
This is reflected in a gradual increase of the potential energy minima when
strain is zero as a function of the cycle number.

Next, the potential energy minima at the end of each cycle are presented in
Fig. 3 for strain amplitudes in the range $0.069\leqslant\gamma_{0}\leqslant
0.075$. Note that the data at zero strain for $\gamma_{0}=0.072$ and $0.075$
are the same as in Fig. 2. It can be clearly observed in Fig. 3 that upon
reducing strain amplitude towards a critical value, the yielding transition
becomes significantly delayed. The exception to this trend is the case of
loading at the strain amplitude $\gamma_{0}=0.071$, where the number of cycles
until yielding is smaller than for $\gamma_{0}=0.072$. In turn, the maximum
number of cycles until the yielding transition is $n_{Y}=685$ for the strain
amplitude $\gamma_{0}=0.069$. We comment that simulations at smaller strain
amplitudes, $\gamma_{0}<0.069$, were not carried out due to the high
computational cost.

The similarity of the functional form for potential energy minima shown in
Fig. 3 suggests a possibility of rescaling the $\hat{x}$-coordinate by the
number of cycles, $n_{Y}$, required for the system to reach the yielding
transition at a given strain amplitude. Fig. 4 shows the same potential energy
curves as a function of the ratio $n/n_{Y}$. Remarkably, the data for
different $\gamma_{0}$ nearly collapse onto a single curve when $n<n_{Y}$. The
master curve is approximately linear in the range $0.2\lesssim n/n_{Y}\lesssim
0.8$, followed by a steep increase due to accumulation of plastic events
within a narrow region that ultimately leads to flow localization when
$n=n_{Y}$. On the other end, the initial slope of the curve is determined by
irreversible rearrangements of group of atoms that settled at relatively
shallow energy minima after thermal annealing. In practice, the function
$U(n/n_{Y})$ can be used to estimate $n_{Y}$ for a binary glass loaded for a
number of cycles ($n<n_{Y}$) at a strain amplitude in the vicinity of the
critical value.

Furthermore, the variation of $n_{Y}$ versus $\gamma_{0}$ is shown in the
inset of Fig. 4. It is readily apparent that the number of cycles until
yielding increases significantly when the strain amplitude approaches a
critical value from above. Moreover, the MD data are well described by the
power-law function, as follows:

$n_{Y}=0.024\cdot(\gamma_{0}-0.067)^{-1.66},$ (3)

where the critical strain amplitude is taken to be $0.067$. This value was
determined previously for a smaller system of $60\,000$ atoms at
$T_{LJ}=0.01\,\varepsilon/k_{B}$ and $\rho=1.2\,\sigma^{-3}$ NVP20altY . These
results imply that the number of cycles to reach the yielding transition might
further increase at lower strain amplitudes and possibly diverge in the case
of athermal systems Sollich22 .

The local plastic events in disordered solids can be accurately identified via
the analysis of nonaffine displacements of atoms Falk98 . As a reminder, the
nonaffine measure for displacement of the $i$-th atom from $\mathbf{r}_{i}(t)$
to $\mathbf{r}_{i}(t+\Delta t)$ is defined via the matrix $\mathbf{J}_{i}$
that transforms positions of its neighboring atoms and minimizes the following
expression:

$D^{2}(t,\Delta
t)=\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}\Big{\\{}\mathbf{r}_{j}(t+\Delta
t)-\mathbf{r}_{i}(t+\Delta
t)-\mathbf{J}_{i}\big{[}\mathbf{r}_{j}(t)-\mathbf{r}_{i}(t)\big{]}\Big{\\}}^{2},$
(4)

where the summation is performed over $N_{i}$ atoms that are initially located
within $1.5\,\sigma$ from $\mathbf{r}_{i}(t)$. It should be noted that the
plastic rearrangement of neighboring atoms during the time interval $\Delta t$
typically corresponds to values of the nonaffine measure $D^{2}(t,\Delta t)$
greater than the cage size, which is about $0.1\,\sigma$ for the KA binary
glass at $\rho=1.2\,\sigma^{-3}$ KobAnd95 .

In Fig. 5 we show the fraction of atoms with relatively large nonaffine
displacements during one cycle, $D^{2}[(n-1)T,T]>0.04\,\sigma^{2}$, for the
strain amplitudes $0.069\leqslant\gamma_{0}\leqslant 0.075$. We comment that
the nonaffine measure was evaluated only for selected cycles due to excessive
computational cost for the large system. It is clearly observed that the shape
of $n_{f}$ is similar to the dependence of energy minima on the number of
cycles shown in Fig. 3. Notice a small peak in $n_{f}$ during the first cycle
due to a number of atoms that become arranged in shallow energy minima upon
thermal annealing, and, as a result, these atoms are prone to plastic
rearrangement under shear deformation. As expected, the yielding transition is
clearly marked by a sharp increase in the fraction $n_{f}$, indicating
extended plastic flow.

In analogy with the potential energy minima shown in Fig. 4, we replot the
same data for $n_{f}$ as a function of the ratio $n/n_{Y}$ in Fig. 6. It is
evident that fractions $n_{f}(n/n_{Y})$ for different values of the strain
amplitude approximately follow a common curve. These results indicate that
plastic rearrangements of only about 1 % of atoms during the first $n_{Y}/2$
cycles result in the increase of the potential energy reported in Fig. 4. A
shear band forms when $n_{f}\approx 0.14$ at $n=n_{Y}$. We also note that both
$n_{f}$ and $U$ increase and level out for $n>n_{Y}$, which reflects widening
of a shear band under cyclic shear. In addition, a closer inspection of the
data in the inset to Fig. 6 reveals that, on average, the fraction $n_{f}$ is
slightly larger for cyclic loading at higher strain amplitudes.

The spatial distribution of plastic rearrangements can be visualized by
plotting positions of atoms with large nonaffine displacements during one
shear cycle, i.e., $\Delta t=T$ in Eq. (4). For example, atomic configurations
for selected number of cycles are presented in Fig. 7 for the strain amplitude
$\gamma_{0}=0.072$ and in Fig. 8 for $\gamma_{0}=0.069$. It can be seen in
Fig. 7 (a) that before yielding, atoms with $D^{2}(200\,T,T)>0.04\,\sigma^{2}$
are organized into small clusters that are homogeneously distributed. Upon
further loading, the glass yields and a shear band forms along the $xy$ plane
during the 218-th cycle, as shown in Fig. 7 (c). During the next 7 cycles, the
shear band becomes wider, which is consistent with the increase in $n_{f}$ and
$U$ after yielding reported in Figs. 3–6. Similar trends can be observed in
Fig. 8 for cyclic loading at $\gamma_{0}=0.069$, except that the orientation
of the shear band is along the $yz$ plane. Also, the sequence of snapshots for
the strain amplitude $\gamma_{0}=0.075$ during the first 100 cycles was
reported in the previous study Priez20ba . Overall, the visualization of
plastic events confirm our earlier conclusions regarding the appearance of
small clusters of atoms that rearrange irreversibly after a full cycle,
followed by the formation of a shear band at the yielding transition, and its
subsequent widening upon continued loading.

##  IV Conclusions

In summary, the effect of oscillatory shear on the damage accumulation and
yielding transition was investigated using molecular dynamics simulations. The
binary glass was prepared by cooling with a computationally slow rate deep
into the glass phase and then subjected to periodic shear deformation with
strain amplitudes slightly greater than a critical value. It was found that
the number of shear cycles until the yielding transition increases
approximately as a power-law function of the difference between the strain
amplitude and the critical value. We showed that the fatigue process proceeds
via a sequence of irreversible rearrangements of small clusters of atoms until
a sudden formation of a shear band at the yielding transition. This behavior
is reflected in the gradual increase of the potential energy at the end of
each cycle and a steep increase near the yielding point. Furthermore, the
potential energy minima for different strain amplitudes closely follow a
master curve when plotted versus the normalized number of cycles. The master
curve can be used to estimate the fatigue lifetime for a binary glass
periodically deformed for only a small number of cycles at a strain amplitude
near the critical value.

## Acknowledgments

Financial support from the National Science Foundation (CNS-1531923) is
gratefully acknowledged. Molecular dynamics simulations were carried out at
Wright State University’s Computing Facility and the Ohio Supercomputer Center
using the LAMMPS code Lammps .

Figure 1: (Color online) The shear stress, $\sigma_{xz}$ (in units
$\varepsilon\sigma^{-3}$), versus cycle number for strain amplitudes
$\gamma_{0}=0.072$ (blue curve) and $\gamma_{0}=0.075$ (black curve). The
period of oscillation is $T=5000\,\tau$.  Figure 2: (Color online) The
potential energy per atom, $U/\varepsilon$, as a function of time for strain
amplitudes $\gamma_{0}=0.072$ (blue curve) and $\gamma_{0}=0.075$ (black
curve). The period of oscillation is $T=5000\,\tau$. The binary glass was
initially prepared by cooling from $T_{LJ}=1.0\,\varepsilon/k_{B}$ to
$0.01\,\varepsilon/k_{B}$ with the rate of $10^{-5}\varepsilon/k_{B}\tau$ at
$\rho=1.2\,\sigma^{-3}$.  Figure 3: (Color online) The potential energy minima
at the end of each shear cycle for strain amplitudes $\gamma_{0}=0.069$,
$0.070$, $0.071$, $0.072$, $0.073$, and $0.075$ (from right to left). The
oscillation period is $T=5000\,\tau$.  Figure 4: (Color online) The potential
energy minima as a function of the ratio, $n/n_{Y}$, where $n=t/T$ is the
cycle number and $n_{Y}$ is the number of cycles to yield at a given strain
amplitude. The same data as in Fig. 3. The symbols are omitted for clarity.
The inset shows $n_{Y}$ versus the strain amplitude $\gamma_{0}$. The dashed
red curve is the best fit to the data given by Eq. (3).  Figure 5: (Color
online) The fraction of atoms with large nonaffine displacements,
$D^{2}[(n-1)\,T,T]>0.04\,\sigma^{2}$, as a function of the cycle number for
the indicated values of the strain amplitude $\gamma_{0}$. Here, $n=t/T$ is
the cycle number, and $T=5000\,\tau$ is the oscillation period.  Figure 6:
(Color online) The fraction of atoms with $D^{2}[(n-1)\,T,T]>0.04\,\sigma^{2}$
versus the ratio $n/n_{Y}$, where $n$ is the cycle number and $n_{Y}$ is the
number of cycles until yielding. The values of $n_{Y}$ are reported in Fig. 4.
The symbols are not shown for clarity. The inset shows an enlarged view of the
same data for $n/n_{Y}\leqslant 1$.  Figure 7: (Color online) The atomic
configurations of the binary glass loaded at the strain amplitude
$\gamma_{0}=0.072$. The nonaffine displacements are shown for atoms with the
nonaffine measure (a) $D^{2}(200\,T,T)>0.04\,\sigma^{2}$, (b)
$D^{2}(216\,T,T)>0.04\,\sigma^{2}$, (c) $D^{2}(217\,T,T)>0.04\,\sigma^{2}$,
and (d) $D^{2}(224\,T,T)>0.04\,\sigma^{2}$. The legend color denotes the
magnitude of $D^{2}$. The atoms are not drawn to scale.  Figure 8: (Color
online) Snapshots of the well-annealed glass subjected to cyclic shear with
the strain amplitude $\gamma_{0}=0.069$. The nonaffine measure is (a)
$D^{2}(600\,T,T)>0.04\,\sigma^{2}$, (b) $D^{2}(683\,T,T)>0.04\,\sigma^{2}$,
(c) $D^{2}(684\,T,T)>0.04\,\sigma^{2}$, and (d)
$D^{2}(689\,T,T)>0.04\,\sigma^{2}$. The colorcode indicates values of $D^{2}$.

## References

  * (1) M. F. Ashby and A. L. Greer, Metallic glasses as structural materials, Scr. Mater. 54, 321 (2006).
  * (2) Z. Sha, W. Lin, L. H. Poh, G. Xing, Z. Liu, T. Wang, and H. Gao, Fatigue of metallic glasses, Appl. Mech. Rev. 72, 050801 (2020).
  * (3) T. Egami, T. Iwashita, and W. Dmowski, Mechanical properties of metallic glasses, Metals 3, 77 (2013).
  * (4) C. J. Gilbert, V. Schroeder, and R. O. Ritchie, Mechanisms for fracture and fatigue-crack propagation in a bulk metallic glass, Metall. Mater. Trans. A 30, 1739 (1999).
  * (5) B. C. Menzel and R. H. Dauskardt, Stress-life fatigue behavior of a Zr-based bulk metallic glass, Acta Mater. 54, 935 (2006).
  * (6) Y. Kurotani and H. Tanaka, Fatigue fracture mechanism of amorphous materials from a density-based coarse-grained model, Commun. Mater. 3, 67 (2022).
  * (7) B. Yang, M. L. Morrison, P. K. Liaw, R. A. Buchanan, G. Wang, C. T. Liu, and M. Denda, Dynamic evolution of nanoscale shear bands in a bulk-metallic glass, Appl. Phys. Lett. 86, 141904 (2005).
  * (8) X. D. Wang, R. T. Qu, Z. Q. Liu, and Z. F. Zhang, Shear band propagation and plastic softening of metallic glass under cyclic compression, J. Alloys Compd. 695, 2016 (2017).
  * (9) T. C. Hufnagel, C. A. Schuh, and M. L. Falk, Deformation of metallic glasses: Recent developments in theory, simulations, and experiments, Acta Mater. 109, 375 (2016).
  * (10) Y. Shi, Size-dependent mechanical responses of metallic glasses, Int. Mater. Rev. 64, 163 (2019).
  * (11) F. Spaepen, A microscopic mechanism for steady state inhomogeneous flow in metallic glasses, Acta Metall. 25, 407 (1977).
  * (12) A. S. Argon, Plastic deformation in metallic glasses, Acta Metall. 27, 47 (1979).
  * (13) F. Kun, M. H. Costa, R. N. Costa Filho, J. S. Andrade Jr, J. B. Soares, S. Zapperi, and H. J. Herrmann, Fatigue failure of disordered materials, J. Stat. Mech. P02003 (2007).
  * (14) J. O. Cochran, G. L. Callaghan, and S. M. Fielding, Slow fatigue and highly delayed yielding via shear banding in oscillatory shear (2022). arXiv:2211.11677
  * (15) N. V. Priezjev, Heterogeneous relaxation dynamics in amorphous materials under cyclic loading, Phys. Rev. E 87, 052302 (2013).
  * (16) I. Regev, T. Lookman, and C. Reichhardt, Onset of irreversibility and chaos in amorphous solids under periodic shear, Phys. Rev. E 88, 062401 (2013).
  * (17) D. Fiocco, G. Foffi, and S. Sastry, Oscillatory athermal quasistatic deformation of a model glass, Phys. Rev. E 88, 020301(R) (2013).
  * (18) I. Regev, J. Weber, C. Reichhardt, K. A. Dahmen, and T. Lookman, Reversibility and criticality in amorphous solids, Nat. Commun. 6, 8805 (2015).
  * (19) Z. D. Sha, S. X. Qu, Z. S. Liu, T. J. Wang, and H. Gao, Cyclic deformation in metallic glasses, Nano Lett. 15, 7010 (2015).
  * (20) N. V. Priezjev, Reversible plastic events during oscillatory deformation of amorphous solids, Phys. Rev. E 93, 013001 (2016).
  * (21) T. Kawasaki and L. Berthier, Macroscopic yielding in jammed solids is accompanied by a non-equilibrium first-order transition in particle trajectories, Phys. Rev. E 94, 022615 (2016).
  * (22) N. V. Priezjev, Nonaffine rearrangements of atoms in deformed and quiescent binary glasses, Phys. Rev. E 94, 023004 (2016).
  * (23) P. Leishangthem, A. D. S. Parmar, and S. Sastry, The yielding transition in amorphous solids under oscillatory shear deformation, Nat. Commun. 8, 14653 (2017).
  * (24) N. V. Priezjev, Collective nonaffine displacements in amorphous materials during large-amplitude oscillatory shear, Phys. Rev. E 95, 023002 (2017).
  * (25) N. V. Priezjev, Molecular dynamics simulations of the mechanical annealing process in metallic glasses: Effects of strain amplitude and temperature, J. Non-Cryst. Solids 479, 42 (2018).
  * (26) N. V. Priezjev, The yielding transition in periodically sheared binary glasses at finite temperature, Comput. Mater. Sci. 150, 162 (2018).
  * (27) A. D. S. Parmar, S. Kumar, and S. Sastry, Strain localization above the yielding point in cyclically deformed glasses, Phys. Rev. X 9, 021018 (2019).
  * (28) N. V. Priezjev, Accelerated relaxation in disordered solids under cyclic loading with alternating shear orientation, J. Non-Cryst. Solids 525, 119683 (2019).
  * (29) N. V. Priezjev, Shear band formation in amorphous materials under oscillatory shear deformation, Metals 10, 300 (2020).
  * (30) W. T. Yeh, M. Ozawa, K. Miyazaki, T. Kawasaki, and L. Berthier, Glass stability changes the nature of yielding under oscillatory shear, Phys. Rev. Lett. 124, 225502 (2020).
  * (31) N. V. Priezjev, Alternating shear orientation during cyclic loading facilitates yielding in amorphous materials, J. Mater. Eng. Perform. 29, 7328 (2020).
  * (32) N. V. Priezjev, A delayed yielding transition in mechanically annealed binary glasses at finite temperature, J. Non-Cryst. Solids 548, 120324 (2020).
  * (33) H. Bhaumik, G. Foffi, and S. Sastry, The role of annealing in determining the yielding behavior of glasses under cyclic shear deformation, PNAS 118, 2100227118 (2021).
  * (34) N. V. Priezjev, Accessing a broader range of energy states in metallic glasses by variable-amplitude oscillatory shear, J. Non-Cryst. Solids 560, 120746 (2021).
  * (35) N. V. Priezjev, Yielding transition in stable glasses periodically deformed at finite temperature, Comput. Mater. Sci. 200, 110831 (2021).
  * (36) S. Cui, H. Liu, and H. Peng, Anisotropic correlations of plasticity on the yielding of metallic glasses, Phys. Rev. E 106, 014607 (2022).
  * (37) B. P. Bhowmik, H. G. E. Hentchel, and I. Procaccia, Fatigue and collapse of cyclically bent strip of amorphous solid, EPL 137, 46002 (2022).
  * (38) N. V. Priezjev, Mechanical annealing and yielding transition in cyclically sheared binary glasses, J. Non-Cryst. Solids 590, 121697 (2022).
  * (39) A. Elgailani, D. Vandembroucq, and C. E. Maloney, Anomalous softness in amorphous matter in the reversible plastic regime (2022). arXiv:2212.10472
  * (40) B. Shang, N. Jakse, P. Guan, W. Wang, and J.-L. Barrat, Influence of oscillatory shear on nucleation in metallic glasses: A molecular dynamics study, Acta Mater. 246, 118668 (2023).
  * (41) C. Reichhardt, I. Regev, K. Dahmen, S. Okuma, and C. J. O. Reichhardt, Perspective on reversible to irreversible transitions in periodic driven many body systems and future directions for classical and quantum systems (2022). arXiv:2211.03775
  * (42) W. Kob and H. C. Andersen, Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture: The van Hove correlation function, Phys. Rev. E 51, 4626 (1995).
  * (43) T. A. Weber and F. H. Stillinger, Local order and structural transitions in amorphous metal-metalloid alloys, Phys. Rev. B 31, 1954 (1985).
  * (44) S. J. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comp. Phys. 117, 1 (1995).
  * (45) M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon, Oxford, 1987).
  * (46) J. T. Parley, S. Sastry, and P. Sollich, Mean-field theory of yielding under oscillatory shear, Phys. Rev. Lett. 128, 198001 (2022).
  * (47) M. L. Falk and J. S. Langer, Dynamics of viscoplastic deformation in amorphous solids, Phys. Rev. E 57, 7192 (1998).


# Confinement-induced resonance from the generalized Gross-Pitaevskii
equations

Alexander Yu. Cherny<EMAIL_ADDRESS>Bogoliubov Laboratory of
Theoretical Physics, Joint Institute for Nuclear Research, 141980, Dubna,
Moscow region, Russia

###### Abstract

The confinement-induced resonances for trapped bosons in the cigar-shaped and
pancake geometries are studied within the generalized Gross-Pitaevskii
equations, which are a simplified version of the Hartree-Fock-Bogoliubov
approximation. Although the Hartree-Fock-Bogoliubov method is considered
applicable only for small interparticle interactions, the resonance
denominators for the chemical potential are obtained in both quasi-one and
quasi-two dimensions. A useful integral representation of the one-particle
Green’s function are found for the cylindrical confinement. We find the
position of a smoothed resonance for the chemical potential in the pancake
geometry at positive scattering length.

##  I Introduction

The confinement-induced resonance (CIR) is an efficient tool to tune effective
interparticle interactions in low-dimensional systems of cold gases [1, 2, 3,
4, 5]. The resonance was predicted in the cigar-shaped [1, 2] and pancake [3,
4] geometries and confirmed experimentally [5]. The cases of very tight and
anisotropic confinements were studied in Refs. [6] and [7, 8, 9],
respectively.

The tight cylindrical or flat atomic waveguide creates, in effect, a low-
dimensional system of cold atoms, and the problem is how to describe the
corresponding interactions in low dimensions. The standard approach [1, 2, 3,
4] is to obtain interparticle interactions in low dimensions from the two-body
scattering problem, thus reducing the effective dimension of the many-body
system.

For identical bosons, it seems more attractive to develop a scheme that
directly describes the wavefunction of the Bose-Einstein condensate (BEC) for
trapping potentials of arbitrary shape. This is all the more interesting
because various optical traps are currently being produced experimentally,
including an optical-box trap [10, *Navon21]. The well-known Gross-Pitaevskii
(GP) equation [12, 13, *gross63] is not suitabe for this purpose, since it
does not reproduce the resonance denominators for effective interactions in
low dimensions. For instance, in the cigar-shaped geometry, the wavefunction
of BEC coincides up to a factor with the ground state of the trapping harmonic
oscillator in the $x-y$ plane provided its frequency is sufficiently high.
Then one can integrate out the $x$ and $y$ coordinates in the interaction term
of the GP functional and obtain the effective coupling constant in one
dimension [15]: $g_{\text{1D}}={2\hbar^{2}a}/{(m\,l^{2})}$. This result is
valid only away from the resonance when $a\ll l$. Here $a$, $m$, and $l$ are
the 3D scattering length in free space, boson mass, and length of the two-
dimensional oscillator, respectively. On the other hand, the correct coupling
constant contains the resonance denominator [1], see Eq. (34) below.

In the paper [16], a system of non-linear equations in both stationary and
time-dependent cases was proposed and called the generalized Gross-Pitaevskii
(GGP) equations. Apart from the standard one-body wavefunction of BEC, the
equations also contain the two-body wavefunction of BEC. This approach is
nothing else but a simplified version of the full Hartree-Fock-Bogoliubov
(HFB) approximation (see Appendix A below). The GGP equations are applicable
for arbitrary shape of the trapping potential. In particular, the resonance
denominators were obtained [16] in both cases of cigar-shaped and pancake
geometries. However, the resulting denominators of Ref. [16], which determine
the position of the resonances, contained incorrect coefficients, since the
solutions were derived with an additional quite rough approximation. As is
shown in this paper, the GGP equations indeed yield the exact coefficients,
which are consistent with the two-body scattering problem [1, 2, 3, 4].

The GGP equations are nonlinear, but the effective coupling constants in low
dimensions are determined by the short-range correlations, where the
nonlinearity is not important. The problem can be reduced to the two-body
problem and the corresponding Green’s functions (GFs). We find and examine the
integral representations for the one-particle GFs at low energies for the flat
and cylindrical harmonic confinements. In the latter case, the equation for
the GF differs from the previously found integral representation [17,
*Idziaszek15, 19], which is actually the Fourier transform of the Feynman
propagator for a harmonic oscillator with contour rotation in the complex
plane. We believe that the found representation is more convenient for
obtaining the low-energy expansions of the GF.

This paper is organized as following. In the next section, the GGP equations
are written and discussed in the case of the tight confinements. In Sec. III,
the integral representations and expansions of the GFs are obtained for the
cylindrical and flat confinements. The CIRs for the effective coupling
constant and chemical potential are found within the suggested scheme in Secs.
IV and V. The main results and prospects are discussed in Conclusions. In
Appendix A, the GGP are derived through the HFB approximation.

##  II Generalized Gross-Pitaevskii equations

We consider a dilute many-body system of identical bosons, interacting with a
pairwise short-range111A potential is of the short-range type if it decreases
at infinity as $1/r^{\alpha}$ with the exponent $\alpha>3$ or faster (see,
e.g., Ref. [20]). For such a potential, only the $s$-wave scattering amplitude
survives in the limit of zero scattering energy. It is called the scattering
length [see Eq. (27) below]. potential $V(r)$, in an external field
$V_{\mathrm{ext}}(\bm{r})$. For simplicity, the bosons are assumed to be
spinless. The equilibrium properties of such system at zero temperature are
described by Eqs. (16) and (17) of Ref. [16], which can be written in the form

$\displaystyle\big{[}\hat{H}(\bm{r}_{1})\\!-\\!\mu\big{]}\phi(\bm{r}_{1})\\!+\\!\int\\!d\bm{r}_{2}\phi(\bm{r}_{2})V(|\bm{r}_{1}\\!-\\!\bm{r}_{2}|)\varphi(\bm{r}_{1},\bm{r}_{2})\\!=\\!0,$
(1)
$\displaystyle\big{[}\hat{H}(\bm{r}_{1})\\!+\\!\hat{H}(\bm{r}_{2})\\!-\\!2\mu\big{]}\psi(\bm{r}_{1},\bm{r}_{2})\\!+\\!V(|\bm{r}_{1}\\!-\\!\bm{r}_{2}|)\varphi(\bm{r}_{1},\bm{r}_{2})\\!=\\!0,$
(2)

where $\mu$ is the chemical potential. Here we set by definition

$\displaystyle\psi(\bm{r}_{1},\bm{r}_{2})=\varphi(\bm{r}_{1},\bm{r}_{2})-\phi(\bm{r}_{1})\phi(\bm{r}_{2}),$
(3)

and the one-body Hamiltonian is given by

$\displaystyle\hat{H}(\bm{r})=-\frac{\hslash^{2}\nabla_{\bm{r}}^{2}}{2m}+V_{\mathrm{ext}}(\bm{r})-E_{0}$
(4)

with $E_{0}$ being its ground-state energy. By definition,
$\phi(\bm{r})\equiv\langle\hat{\Psi}({\bm{r}})\rangle$ and
$\varphi(\bm{r}_{1},\bm{r}_{2})\equiv\langle\hat{\Psi}({\bm{r}_{1}})\hat{\Psi}({\bm{r}}_{2})\rangle$
are the anomalous averages of the field Bose operators. The derivation of
these equations and their relation to the HFB approximation are given in
Appendix A.

The function $\phi(\bm{r})$ is the one-body function of inhomogeneous BEC and
can be considered as the order parameter. The function
$\varphi(\bm{r}_{1},\bm{r}_{2})$ is called pair wavefunction. It describes two
bosons in the Bose-Einstein condensate, and this function is _not_ reduced in
general to the product of the one-body functions. The exact definitions can be
done through the one-body and two-body matrices [21, 22, 23]. The property of
decay of correlations dictates the boundary condition

$\displaystyle\varphi(\bm{r}_{1},\bm{r}_{2})\simeq\phi(\bm{r}_{1})\phi(\bm{r}_{2})\quad\text{when}|\bm{r}_{1}-\bm{r}_{2}|\gg\xi_{\mathrm{h}},$
(5)

where $\xi_{\mathrm{h}}$ is the healing (or coherence) length defined through
the chemical potential: $\xi_{\mathrm{h}}=\hbar/\sqrt{\mu m}$. The system of
equations (1) and (2) is nonlinear due to the presence of the condensate
wavefunction $\phi(\bm{r}_{2})$ under the integral in Eq. (1). The GGP
equations (1) and (2), being a particular case of HFB approximation, are
applicable to weak interactions when the mean distance between bosons is much
smaller than the healing length.

Various solutions for different external potentials were considered in Ref.
[16]. In particular, the 1D and 2D harmonic CIRs were obtained [see Eqs. (48)
and (57) of Ref. [16]]. However, the used method did not enable us to obtain
the precise coefficients in the resonance denominators, which determine the
CIR positions. In this paper, we show that the correct coefficients can be
derived from Eqs. (1) and (2).

We assume that the radius $r_{0}$ of the short-range interaction potential
$V(r)$ is much smaller than the characteristic length of the external
potential $V_{\mathrm{ext}}(\bm{r})$ and the mean distance between bosons.
This implies that at low energies, physical processes are governed by the
$s$-wave scattering length $a$ of the interaction potential, which is obtained
from the symmetric solution of the two-body Scrödinger equation at zero energy
in free space.

Thus, our task is to get the chemical potential of $N$ strongly confined
bosons at zero temperature as a function of the 3D scattering length $a$ and
oscillator frequency $\omega$ in the case of cylindrical and flat harmonic
confinements.

##  III Green’s functions of the one-body Schrödinger equation at low energies
for cylindrical and flat harmonic confinements

To obtain the chemical potentials of the Bose-Einstein condensates for the
cigar-shaped and pancake geometries, we need to calculate the corresponding
GFs for the one-body Hamiltonian (4), which is defined as

$\displaystyle\hat{G}_{E}=-(\hat{H}-E)^{-1}$ (6)

with the scattering energy $E$, which is assumed to be sufficiently small (see
the discussions below).

In this section, we find useful integral representations for the GFs in the
cases of cylindrical and flat harmonic confinements. These representations are
very convenient for obtaining the low-energy expansions of the GFs.

###  III.1 Cylindrical confinement

The cylinder harmonic confinement implies that a particle is confined in the
$x$-$y$ plane but can move freely in $z$-direction with a finite momentum,
whose absolute value is denoted by $\hbar p$. The one-particle Hamiltonian (4)
is given by

$\displaystyle\hat{H}=-\frac{\hslash^{2}\nabla^{2}}{2m_{0}}+\frac{m_{0}\omega^{2}(x^{2}+y^{2})}{2}-E_{0},$
(7)

where the shift $E_{0}=\hbar\omega$ is introduced to get zero ground-state
energy.

It is convenient to use the following notation for radius vector:
$\bm{r}=(\bm{\rho},z)$ with $\bm{\rho}=(x,y)$. We calculate the GF (6) in
coordinate representation with the radial coordinate $\rho=\sqrt{x^{2}+y^{2}}$
being equal to zero and $E=\frac{\hbar^{2}p^{2}}{2m_{0}}$:

$\displaystyle G_{p}(z-z^{\prime})\equiv\langle
z,\rho=0|\hat{G}_{E}|\rho^{\prime}=0,z^{\prime}\rangle$
$\displaystyle=-\frac{1}{2\pi^{2}l^{2}}\sum_{n=0}^{\infty}\text{p.v.}\\!\\!\int_{-\infty}^{\infty}dq\frac{e^{iq(z-z^{\prime})}}{2\hbar\omega
n+\frac{\hbar^{2}}{2m_{0}}(q^{2}-p^{2})}.$ (8)

Here $l_{0}=\sqrt{\hbar/(m_{0}\omega)}$ is the harmonic oscillator length, and
p.v. denotes the Cauchy principal value. This kind of regularization ensures
the reality of the GF in coordinate representation [24].

Here we use the expression for the GF in Dirac notation
$\langle\bm{r}|\hat{G}_{E}|\bm{r}^{\prime}\rangle=-\sum_{k}\langle\bm{r}|k\rangle\langle
k|\bm{r}^{\prime}\rangle/(E_{k}-E)$ with $|k\rangle$ and $E_{k}$ being the
eigenfunctions and eigenvalues of a Hamiltonian, respectively. For the
Hamiltonian (7), $k=(n,m,q)$ is a multi-index, where $n=0,1,2,\ldots$ is the
main quantum number, $m=0,\pm 1,\pm 2,\ldots$ is the orbital momentum in 2D,
and $q$ is the $z$-component of the wavevector. The eigenvalues and normalized
eigenfunctions are given by $E_{k}=\hbar\omega(2n+|m|)+\frac{q^{2}}{2m_{0}}$
and
$\varphi_{k}(\rho,\alpha,z)=\frac{e^{iqz}}{\sqrt{2\pi}}\phi_{n,m}(\rho,\alpha)$,
respectively. Here $\phi_{n,m}(\rho,\alpha)$ are the eigenvalues for the 2D
oscillator in the polar coordinates $x=\rho\cos\alpha$, $y=\rho\sin\alpha$.
One can show that for an arbitrary main quantum number,
$\phi_{n,0}(\rho=0,\alpha)=1/(\sqrt{\pi}l_{0})$ and
$\phi_{n,m}(\rho=0,\alpha)=0$ for $m\not=0$. Since we are interested at the GF
of the Hamiltonian (7) only at $\rho=\rho^{\prime}=0$, we are left with one
sum over the main quantum number and thus arrive at Eq. (8).

When the scattering momentum is sufficiently small, $p<2/l_{0}$, only the
first term in the sum in Eq. (8) is singular. Separating this term and using
relations
$\int_{-\infty}^{\infty}\frac{dq}{2\pi}\frac{e^{iqz}}{(q^{2}+\varkappa^{2})}=\frac{e^{-\varkappa|z|}}{2\varkappa}$
and
$\text{p.v.}\\!\\!\int_{-\infty}^{\infty}\frac{dq}{2\pi}\frac{e^{iqz}}{(p^{2}-q^{2})}=\frac{\sin
p|z|}{2p}$ yield

$\displaystyle\frac{\hbar^{2}}{2m_{0}}G_{p}(z)=$ $\displaystyle\frac{1}{\pi
l_{0}^{2}}\left[\frac{\sin
p|z|}{2p}-\frac{l_{0}}{4}\Lambda\left(\frac{2|z|}{l_{0}},-\frac{p^{2}l_{0}^{2}}{4}\right)\right],$
(9) $\displaystyle\Lambda(\xi,\epsilon)=$
$\displaystyle\sum_{n=1}^{\infty}\frac{e^{-\xi\sqrt{n+\epsilon}}}{\sqrt{n+\epsilon}},$
(10)

where dimensionless variables $\xi=2|z|/l_{0}$ and
$\epsilon=-p^{2}l_{0}^{2}/4$ are introduced.

The first term of the GF (9) is proportional to the one-dimensional GF of a
free particle, which obeys the relation $\left(\frac{\partial^{2}}{\partial
z^{2}}+p^{2}\right)\frac{\sin p|z|}{2p}=\delta(z)$. At large distances $|z|\gg
l_{0}$, the second term of the GF falls off exponentially as it follows
directly from the expansion (10). Then the asymptotics of the scattering part
of the wavefunction, proportional to the GF, is essentially one-dimensional,
as expected.

Let us calculate the GF at low momenta $|\epsilon|\ll 1$ and small coordinates
$\xi\ll 1$. To expand $\Lambda(\xi,\epsilon)$ in $\xi$ and $\epsilon$ around
zero, we use the identity for exponential
$e^{-x}=\frac{1}{\sqrt{\pi}}\int_{0}^{\infty}\frac{dt}{t^{3/2}}e^{-1/t}e^{-x^{2}t/4}$
and find the integral representation

$\displaystyle\Lambda(\xi,\epsilon)=\frac{1}{\sqrt{\pi}}\int_{0}^{\infty}\frac{dt}{t^{3/2}}e^{-\frac{1}{t}}e^{-\frac{\xi^{2}(1+\epsilon)t}{4}}\Phi\left(e^{-\frac{\xi^{2}t}{4}},\frac{1}{2},1+\epsilon\right),$
(11)

where $\Phi\left(x,s,\alpha\right)=\sum_{n=0}^{\infty}x^{n}/(n+\alpha)^{s}$ is
the Lerch transcendent, or Lerch zeta function [25]. Upon expanding the
integrand

$\displaystyle
e^{-\frac{\xi^{2}(1+\epsilon)t}{4}}\Phi\left(e^{-\frac{\xi^{2}t}{4}},\frac{1}{2},1+\epsilon\right)=$
$\displaystyle\frac{2\sqrt{\pi}}{\sqrt{t}\xi}+\zeta({1}/{2})-\frac{\epsilon}{2}\zeta({3}/{2})$
$\displaystyle+O(\xi^{2})+O(\epsilon^{2})$

and substituting it into Eq. (11), we obtain

$\displaystyle\Lambda(\xi,\epsilon)=$
$\displaystyle\frac{2}{\xi}+\zeta({1}/{2})+\frac{\xi}{2}-\frac{\epsilon}{2}\zeta({3}/{2})+O(\xi^{2})+O(\epsilon^{2})$
$\displaystyle+O(\xi\epsilon)$ (12)

with $\zeta(x)$ being the Riemann zeta function. The symbol $O(x)$ is the big
$O$ notation, which designates terms of the order of $x$ or smaller. When
$|z|\ll l_{0}$ and $p\ll 1/l_{0}$, the expansion of the GF (9) takes the form

$\displaystyle\frac{\hbar^{2}}{2m_{0}}G_{p}(z)=$
$\displaystyle-\frac{1}{4\pi|z|}\\!-\\!\frac{1}{4\pi
l_{0}}\zeta({1}/{2})\\!+\\!\frac{|z|}{4\pi
l_{0}^{2}}\\!-\\!\frac{p^{2}l_{0}}{32\pi}\zeta({3}/{2})$
$\displaystyle+O(z^{2}/l_{0}^{2})+O(|z|\,p^{2}l_{0})+O(p^{4}l_{0}^{4}).$ (13)

###  III.2 Flat confinement

The flat harmonic confinement along the $z$-axis is considered by close
analogy with the previous subsection. A particle, being confined in
$z$-direction, can move freely in the $x$-$y$ plane. The Hamiltonian reads

$\displaystyle\hat{H}=-\frac{\hslash^{2}\nabla^{2}}{2m_{0}}+\frac{m_{0}\omega^{2}z^{2}}{2}-E_{0}$
(14)

with the shift $E_{0}=\hbar\omega/2$. A complete set of its eigenfunctions is
given by
$\varphi_{\bm{q},n}(\bm{\rho},z)=\phi_{n}(z)e^{i\bm{q}\cdot\bm{\rho}}/(2\pi)$
with $\phi_{n}(z)$ being the normalized eigenfunctions of the linear harmonic
oscillator. We need only their absolute values at $z=0$: when $n$ is odd then
$\phi_{n}(0)=0$, while for even $n$

$\displaystyle|\phi_{2k}(0)|=\frac{1}{\pi^{1/4}l_{0}^{1/2}}\frac{\sqrt{(2k)!}}{2^{k}k!},\quad
k=0,1,2,\ldots,$ (15)

see, e.g., the textbook [26].

Using the same method as in the previous subsection, we find the expression
for the GF at $z=z^{\prime}=0$ and $E=\frac{\hbar^{2}p^{2}}{2m_{0}}$

$\displaystyle
G_{p}(|\bm{\rho}-\bm{\rho}^{\prime}|)\equiv\langle\bm{\rho},z=0|\hat{G}_{E}|z^{\prime}=0,\bm{\rho}^{\prime}\rangle$
$\displaystyle=-\sum_{k=0}^{\infty}|\phi_{2k}(0)|^{2}\text{p.v.}\\!\\!\int\frac{d^{2}q}{(2\pi)^{2}}\frac{e^{i\bm{q}\cdot(\bm{\rho}-\bm{\rho}^{\prime})}}{2\hbar\omega
k+\frac{\hbar^{2}}{2m_{0}}(q^{2}-p^{2})},$ (16)

where $|\phi_{2k}(0)|$ is given by Eq. (15). When the momentum is small,
$p<2/l_{0}$, only the first term in the sum is singular. Separating it and
using the identities
$\text{p.v.}\\!\\!\int\frac{d^{2}q}{(2\pi)^{2}}\frac{e^{i\bm{q}\cdot\bm{\rho}}}{p^{2}-q^{2}}=\frac{Y_{0}(p\rho)}{4}$,
$\int\frac{d^{2}q}{(2\pi)^{2}}\frac{e^{i\bm{q}\cdot\bm{\rho}}}{\varkappa^{2}+q^{2}}=\frac{K_{0}(\varkappa\rho)}{2\pi}$
for positive $p$ and $\varkappa$ yield

$\displaystyle\frac{\hbar^{2}}{2m_{0}}G_{p}(\rho)=$
$\displaystyle\frac{1}{\sqrt{\pi}l_{0}}\left[\frac{Y_{0}(p\rho)}{4}-\frac{1}{2\pi}S\left(\frac{2\rho}{l_{0}},-\frac{p^{2}l_{0}^{2}}{4}\right)\right],$
(17) $\displaystyle S(\xi,\epsilon)=$
$\displaystyle\sum_{k=1}^{\infty}K_{0}(\xi\sqrt{k+\epsilon})\frac{(2k)!}{4^{k}(k!)^{2}}$
(18)

with the dimensionless variables $\xi=2\rho/l_{0}$ and
$\epsilon=-p^{2}l_{0}^{2}/4$. Here $Y_{0}(x)$ and $K_{0}(x)$ are the Bessel
and modified Bessel functions of the second kind, respectively (see, e.g.
[27], Ch. 9).

Similar to the one-dimensional case considered in Sec. III.1, the first term
is proportional to the GF of a free particle in two dimensions:
$\left(\frac{\partial^{2}}{\partial\rho^{2}}+\frac{1}{\rho}\frac{\partial}{\partial\rho}+p^{2}\right)\frac{Y_{0}(p\rho)}{4}=\delta(\bm{\rho})$.
When $\rho\gg l_{0}$, the second term in Eq. (17) is exponentially small,
which actually implies two-dimensional scattering.

The relation $K_{0}(x)=\int_{0}^{\infty}\frac{dt}{2t}e^{-t}e^{-x^{2}/(4t)}$
enables us to sum up the series (18) under the integral and find its integral
representation in elementary functions

$\displaystyle
S(\xi,\epsilon)=\int_{0}^{\infty}\frac{dt}{2t}\exp\left(-t-\frac{\epsilon\,\xi^{2}}{4t}\right)\left[\frac{1}{\sqrt{1-e^{-\xi^{2}/(4t)}}}-1\right].$
(19)

This relation is a particular case of a more general equation, which was
obtained by another method in Refs. [17, 19]. This general equation is
actually the Fourier transform of the well-known Feynman propagator for a
harmonic oscillator.

It is easy to obtain the expansion

$\displaystyle S(\xi,\epsilon)=$
$\displaystyle\frac{\sqrt{\pi}}{\xi}+\ln\xi+\gamma-\sqrt{\pi}-\ln
2+C-\xi\frac{\sqrt{\pi}}{8}-\epsilon\ln 2$
$\displaystyle+O(\xi^{2})+O(\epsilon^{2})+O(\xi\epsilon),$ (20)

where $\gamma=0.5772\ldots$ is Euler’s constant and

$\displaystyle
C\\!=\\!\int_{0}^{\infty}\\!\frac{dy\,e^{-y}}{2y}\left[\frac{1}{\sqrt{1-e^{-y)}}(1+\sqrt{1-e^{-y)}})}\\!+\\!1-\\!\frac{1}{\sqrt{y}}\right].$

Its numerical value is given by $C=0.8035589\ldots$.

To find the GF at $\rho\ll l_{0}$ and $p\ll 1/l_{0}$, we substitute the
expansions (20) and $Y_{0}(x)=2\ln(xe^{\gamma}/2)/\pi+O(x^{2}\ln x)$ into Eq.
(17) and finally arrive at

$\displaystyle\frac{\hbar^{2}}{2m_{0}}G_{p}(\rho)=$
$\displaystyle-\frac{1}{4\pi\rho}+\frac{1}{4\pi^{3/2}l_{0}}\ln\left(\frac{\pi
p^{2}l_{0}^{2}}{2B}\right)$ $\displaystyle+\frac{1}{8\pi
l_{0}}\left[\frac{\rho}{l_{0}}-\frac{\ln
2}{\sqrt{\pi}}p^{2}l_{0}^{2}\right]+O(\rho^{2}p^{2}\ln p\rho)$
$\displaystyle+O(\rho^{2}/l_{0}^{2})+O(\rho p^{2}l_{0})+O(p^{4}l_{0}^{4}),$
(21)

where we put by definition

$\displaystyle B\equiv 2\pi\exp[2(C-\sqrt{\pi})]=0.904916\ldots.$ (22)

The first two terms of GF (21) were first calculated in [3] with the constant
$B$ being equal to one. In the subsequent paper [4], its numerical value was
estimated to be $0.915$. Further, the value of $B$ being equal to $0.9049$ was
obtained in Ref. [17, *Idziaszek15], which is consistent with Eq. (22). We
believe that this is the most accurate estimation.

##  IV Confinement-induced resonance in one dimension: cylindrical geometry

###  IV.1 The tight confinement

When the bosons are trapped in the $x$-$y$ transverse direction with the
axially symmetric harmonic potential
$V_{\mathrm{ext}}(\rho)=m\omega^{2}\rho^{2}/2$, the order parameter depends
only on the transverse coordinates, but it is homogeneous along the $z$-axis.

The tight confinement implies that the chemical potential is small in scale of
$\hbar\omega$, that is, $l\ll\xi_{\mathrm{h}}$, in terms of the oscillator
length $l=\sqrt{\hbar/(m\omega)}$. Then with a good accuracy one can put for
the wavefuncton of BEC [28]
$\phi=\phi(\rho_{1})=\sqrt{n_{\mathrm{1D}}}\phi_{0}(\rho_{1})$ with
$\phi_{0}(\rho_{1})=\frac{\exp[-\rho^{2}_{1}/(2l^{2})]}{l\sqrt{\pi}}$ being
the ground state wavefunction of the one-particle Schrödinger equation with
Hamiltonian (4), and $n_{\mathrm{1D}}$ is the linear density of the Bose-
Einstein condensate. Thus, the nonlinear term in Eq. (1) is treated as a
perturbation, and the chemical potential is obtained by multiplying Eq. (1) by
$\phi_{0}(\rho_{1})$ and integrating out over $\bm{r}_{1}$:

$\displaystyle\mu=\frac{1}{L}\iint
d\bm{r}_{1}d\bm{r}_{2}\phi_{0}(\rho_{1})\phi_{0}(\rho_{2})V(|\bm{r}_{1}\\!-\\!\bm{r}_{2}|)\varphi(\bm{r}_{1},\bm{r}_{2}),$
(23)

where $L$ is the length of the system along $z$-axis.

To take the integral (23), we need to examine the properties of the pair
wavefunction. In the case of transverse confinement, Eqs. (1) and (2) are
translationally invariant under a spatial translation $z_{1}\to z_{1}+b$,
$z_{2}\to z_{2}+b$ for any $b$, which yields
$\varphi=\varphi(\bm{\rho}_{1},\bm{\rho}_{2},|z_{1}-z_{2}|)$. Moreover, for
the harmonic confinement, the spatial variables in Eq. (2) are separated into
the center-of-mass position and relative coordinates by analogy with the two-
body problem. We can look for a solution in the form

$\displaystyle\varphi(\bm{r}_{1},\bm{r}_{2})=\varphi(\rho,z)\exp\left[-\big{(}\bm{\rho}_{1}+\bm{\rho}_{2}\big{)}^{2}/(4l^{2})\right]$
(24)

and the same relation for $\psi(\bm{\rho}_{1},\bm{\rho}_{2},z)$, see Eq. (3).
Here $\bm{\rho}=\bm{\rho}_{1}-\bm{\rho}_{2}$ and $z=z_{1}-z_{2}$ are the
relative coordinates, and $\rho=|\bm{\rho}|$. The boundary conditions (5) read
$\varphi(\rho,z)\simeq n_{\rm 1D}\exp[-\rho^{2}/(4l^{2})]/(\pi l^{2})$ when
$z\gg\xi_{\mathrm{h}}$.

In the range $r\ll l$, one can use the approximation

$\displaystyle\varphi(\rho,z)=\eta\varphi^{(0)}_{\rm 3D}(r)\frac{n_{\rm
1D}}{\pi l^{2}},$ (25)

since $r_{0}\ll l$. Here $r=\sqrt{\rho^{2}+z^{2}}$, and $\varphi^{(0)}_{\rm
3D}(r)$ is the 3D solution of the two-body Schrödinger equation in the center-
of-mass system at zero energy

$-(\hbar^{2}/m)\nabla^{2}\varphi^{(0)}_{\rm 3D}(r)+V(r)\varphi^{(0)}_{\rm
3D}(r)=0.$ (26)

The solution is assumed to be regular in the vicinity of the origin and obeys
the boundary condition when $r\gg r_{0}$

$\displaystyle\varphi^{(0)}_{\rm 3D}(r)\simeq$ $\displaystyle\,1-a/r,$ (27)
$\displaystyle a=$ $\displaystyle\,\frac{m}{4\pi\hbar^{2}}\int
d\bm{r}V(r)\varphi^{(0)}_{\rm 3D}(r).$ (28)

Substituting Eqs. (24) and (25) into Eq. (23) and using the difference of the
scales $r_{0}\ll l$ and Eq. (28), we are left with

$\displaystyle\mu=\eta\frac{2\hbar^{2}a}{m\,l^{2}}n_{\rm 1D}.$ (29)

Here the point is that the constant $\eta\not=1$, which is crucial for the
appearance of CIR. This constant can be recovered from Eq. (2).

###  IV.2 The application of the Green’s function

Substituting Eq. (24) into Eq. (2), we arrive at the equation for the
wavefunction $\varphi(\rho,z)$:

$\displaystyle\left[-\frac{\hbar^{2}\nabla^{2}}{2m_{0}}+\frac{m_{0}\omega^{2}\rho^{2}}{2}-E_{0}-2\mu\right]\psi(\rho,z)=-V(r)\varphi(\rho,z)$
(30)

with $m_{0}=m/2$ being the reduced mass. It is sufficient to find the
wavefunction $\varphi(\rho,z)$ at $\rho=0$. According to Eqs. (25) and (27),
it has the asymptotics

$\displaystyle\varphi(\rho=0,z)\simeq\eta\frac{n_{\rm 1D}}{\pi
l^{2}}\left(1-\frac{a}{|z|}\right)$ (31)

in the short-range regime $r_{0}\ll|z|\ll l$.

The solution of Eq. (30) can be found with the Green’s function (GF) at
$E=2\mu=\frac{\hbar^{2}p^{2}}{2m_{0}}$ as $\psi(\rho,z)=\int
d\bm{r}^{\prime}\langle
z,\bm{\rho}|\hat{G}_{E}|\bm{\rho}^{\prime},z^{\prime}\rangle
V(r^{\prime})\varphi(\rho^{\prime},z^{\prime})$. In the short-range regime
$|z|\ll l$ and at low energies $p\ll 1/l$, the GF at $\rho=\rho^{\prime}=0$
was obtained in Sec. III.1 and given by Eq. (13). At $m_{0}=m/2$ and
$l_{0}=\sqrt{\hbar/(m_{0}\omega)}=\sqrt{2}l$, we get

$\displaystyle\varphi(\rho=0,z)=\frac{n_{\rm 1D}}{\pi l^{2}}+G_{p}(z)\int
d\bm{r}^{\prime}V(r^{\prime})\varphi(\rho^{\prime},z^{\prime}),$

where Eq. (3) is used: $\varphi(\rho=0,z)=n_{\rm 1D}/(\pi
l^{2})+\psi(\rho=0,z)$.

The potential $V(r)$ is localized within $r\lesssim r_{0}$, where the
wavefunction $\varphi(\rho,z)$ takes the form (25). Using the relation (28)
and the first two terms of Eq. (13) for the Green’s function, we find

$\displaystyle\varphi(\rho=0,z)\simeq\frac{n_{\rm 1D}}{\pi
l^{2}}\left[1-\eta\,a\left(\frac{1}{|z|}+\frac{\zeta(1/2)}{\sqrt{2}l}\right)\right].$

Comparing this asymptotics with Eq. (31) finally yields

$\displaystyle\eta=\left(1+\frac{a\,\zeta(1/2)}{\sqrt{2}l}\right)^{-1},$ (32)

where $\zeta(1/2)=-1.46035\ldots$.

###  IV.3 The chemical potential

The chemical potential is given by Eqs. (29) and (32)

$\displaystyle\mu$ $\displaystyle=g_{\text{1D}}n_{\rm 1D},$ (33)
$\displaystyle g_{\text{1D}}$
$\displaystyle=\frac{2\hbar^{2}a}{m\,l^{2}}\frac{1}{1+\frac{a\,\zeta(1/2)}{\sqrt{2}l}}.$
(34)

Taking into account the relation
$g_{\text{1D}}=-\frac{2\hbar^{2}}{m\,a_{\mathrm{1D}}}$ for the coupling
constant in 1D, we arrive at the expression for the effective 1D scattering
length $a_{\mathrm{1D}}$ in the two-body problem (see, e.g., Ref. [20]), first
obtained in Ref. [1]

$\displaystyle{a_{\mathrm{1D}}}=-\frac{l^{2}}{a}\left(1+\frac{a\,\zeta(1/2)}{\sqrt{2}l}\right).$

##  V Confinement-induced resonance in two dimensions: flat geometry

###  V.1 The tight confinement

The flat confinement of identical bosons is provided by harmonic forces along
the $z$-axis: $V_{\mathrm{ext}}(z)=m\omega^{2}z^{2}/2$. The order parameter is
homogeneous within the $x$-$y$ plane, and in the case of the tight confinement
$l\ll\xi_{\mathrm{h}}$, one can put with a good accuracy
$\phi=\phi(z_{1})=\sqrt{n_{\mathrm{2D}}}\phi_{0}(z_{1})$ with the oscillator
ground-state wavefunction
$\phi_{0}(z_{1})=\exp[-z_{1}^{2}/(2l^{2})]/(l^{1/2}\pi^{1/4})$. Here
$n_{\mathrm{2D}}$ is the 2D density of the Bose-Einstein condensate. The
chemical potential can be obtained from Eq. (1) by analogy with the previous
section

$\displaystyle\mu=\frac{1}{S}\iint
d\bm{r}_{1}d\bm{r}_{2}\phi_{0}(z_{1})\phi_{0}(z_{2})V(|\bm{r}_{1}\\!-\\!\bm{r}_{2}|)\varphi(\bm{r}_{1},\bm{r}_{2}),$
(35)

where $S$ is the area of the system occupied in the $x$-$y$-plane.

The pair wavefunction is invariant under translations along the $x$-$y$ plane:
$\varphi(\bm{r}_{1},\bm{r}_{2})=\varphi(z_{1},z_{2},\rho)$. It can be written
for the harmonic confinement as
$\varphi(\bm{r}_{1},\bm{r}_{2})=\varphi(z,\rho)\exp\left[-(z_{1}+z_{2})^{2}/(4l^{2})\right]$.

The approximation for the pair wavefunction for $r\ll l$

$\displaystyle\varphi(z,\rho)=\eta\varphi^{(0)}_{\rm 3D}(r)\frac{n_{\rm
2D}}{l\sqrt{\pi}},$ (36)

is analogous to Eq. (25). In the same manner, we find from Eqs. (28), (35) and
(36)

$\displaystyle\mu=\frac{4\pi\hbar^{2}}{m}n_{\rm 2D}u,$ (37)

where we introduce a more convenient dimensionless parameter

$\displaystyle u=\frac{\eta a}{\sqrt{2\pi}l}.$ (38)

The prefactor $\eta$ can also be found by analogy with Sec. IV.

###  V.2 The application of the Green’s function

On the one hand, in the shot-range regime $r_{0}\ll\rho\ll l$, Eqs. (27) and
(36) yield

$\displaystyle\varphi(z=0,\rho)\simeq\eta\frac{n_{\rm
2D}}{l\sqrt{\pi}}\left(1-\frac{a}{\rho}\right).$ (39)

One the other hand, one can use the GF for Eq. (2) (see Sec. III.2) and obtain
in the same regime

$\displaystyle\varphi(z=0,\rho)\simeq$ $\displaystyle\,\frac{n_{\rm
2D}}{\sqrt{\pi}l}+G_{p}(\rho)\int
d\bm{r}^{\prime}V(r^{\prime})\varphi(z^{\prime},\rho^{\prime})$
$\displaystyle=$ $\displaystyle\,\frac{n_{\rm
2D}}{\sqrt{\pi}l}+G_{p}(\rho)\eta\frac{n_{\rm
2D}}{l\sqrt{\pi}}\frac{4\pi\hbar^{2}a}{m}$ (40)

with the GF (21), where $2\mu=\frac{\hbar^{2}p^{2}}{2m_{0}}$, $m_{0}=m/2$,
$l_{0}=\sqrt{2}l=\sqrt{2\hbar/(m\omega)}$.

###  V.3 The chemical potential

In the quasi-1D condensate, considered in the previous section V, the GF
remains finite in the limit $p\to 0$. By contrast, in the quasi-2D condensate,
the GF (21) exhibits a logarithmic divergence as $p\to 0$. This implies that
the chemical potential should be determined self-consistently. Indeed,
comparing Eq. (39) with Eq. (40) and using the definition (38) yield the
algebraic equation for $u$

$\displaystyle\frac{1}{u}+\ln u=\delta,$ (41)
$\displaystyle\delta=\frac{\sqrt{2\pi}l}{a}-\ln\left[\frac{8\pi^{2}n_{\rm
2D}l^{2}}{B}\right],$ (42)

where the constant $B$ is defined by Eq. (22). Thus, the chemical potential is
given by Eqs. (37), (41), and (42).

Equation (41) has no real solutions when $\delta<1$ and two solutions when
$\delta>1$: one is greater than 1, and the other is less than 1. It follows
from the thermodynamic stability condition $\frac{\partial\mu}{\partial n_{\rm
2D}}=\frac{4\pi\hbar^{2}}{m}\frac{u}{1-u}>0$ that only the solution $u<1$ is
physical. This fits well with the condition of weak coupling, which tells us
that the healing length should be at least smaller than the mean distance
between bosons on the plane: $\frac{\hbar}{\sqrt{\mu m}}>\frac{1}{\sqrt{n_{\rm
2D}}}$, see Eq. (37). When the dimensionless parameter $\delta$ is big then
$u\simeq 1/\delta$ and we arrive at the weak-coupling relation for the
chemical potential

$\displaystyle\mu\simeq\frac{4\pi\hbar^{2}n_{\rm 2D}}{m\delta},$ (43)

which allows us to interpret $\delta$, given by Eq. (42), as the resonance
denominator. Note, however, that this approximation overestimates the values
of the chemical potential, see Fig. 1a.

Figure 1:  (a) The solution of Eq. (41), $u(\delta)$ (blue solid line). The
parameter $u$ is the chemical potential in units of $4\pi\hbar^{2}n_{\rm
2D}/m$, see Eq. (37). At large $\delta$, the function $u(\delta)$
asymptotically approaches $1/\delta$ (dashed red line). (b) The chemical
potential in the pancake geometry [in units of $4\pi\hbar^{2}n_{\rm 2D}/m$],
obtained from Eqs. (41) and (42), vs. the oscillator length $l$ of the
trapping potential [in units of the 3D scattering length $a$] (blue solid
line). The same diagram in the approximation $u\simeq 1/\delta$ is shown in
dashed red line. This approximation overestimates the values of the chemical
potential. The absolute value of $a$ and 2D density $n_{\rm 2D}$ are supposed
to be fixed, while the sign of $a$ can change. When $a$ is negative, the CIR
is observed: the chemical potential grows rapidly with increasing of the
oscillator length, because the parameter $\delta$ [Eq. (42)] decreases. When
the scattering length is positive, one can also observe a kind of smoothed
resonance behaviour of the chemical potential with the maximum at
$l/a=0.7978\ldots$, see Eq. (44). The position of the maximum is independent
of the 2D density.

The CIR arises (see Fig. 1b) when the 3D scattering length is negative while
the 2D density is sufficiently small for $\delta$ to be positive and close to
1, which is its minimum possible value within our approach. At positive
scattering length, a reminiscent of the resonance is also observed, when the
dimensionless chemical potential reaches the local maximum at

$\displaystyle\frac{l}{a}=\sqrt{\frac{2}{\pi}}=0.7978\ldots$ (44)

as it follows from Eq. (42). This differs from experimentally found [5] value
${l}/{a}=1.19(3)$, at which the CIR is realized at positive scattering
lengthes.

Relations (37), (41), and (42) should be used with some caveats. First, the
quasi 2D regime imposes the restriction on the oscillator length: the
characteristic wavevector of bosons $p\simeq\sqrt{n_{\rm 2D}}$ should be
smaller than $2/l$ as discussed in Sec. III.2. This yields
$\frac{l}{a}\lesssim\frac{2}{\sqrt{n_{\rm 2D}a^{2}}}$, which for the chosen
parameter $n_{\rm 2D}a^{2}=0.05$ in Fig. 1b is of order of 9. Second, we
assume that the range of interaction potential is much smaller than the
oscillator length: $r_{0}\ll l$. Thus, in the vicinity of $l=0$, the
approximation of zero-range interactions ($r_{0}=0$) is used. We believe that
the extremely tight confinement $r_{0}\gtrsim l$ also can be described by
means of the GGP equations with methods similar to those considered in Ref.
[6].

##  VI Conclsusions

Using the GGP equations (1) and (2), we obtained the correct resonance
denominators, which determine the positions of CIRs in the cigar-shaped and
pancake geometries, see Eqs. (33), (34), and (42), (43).

The GGP equations, being a particular case of the HFB approximation (see
Appendix A), are applicable for a weakly interacting Bose gas, strictly
speaking. Nevertheless, the GGP equations correctly reproduce the effective
coupling in the low-dimensional geometries. The reason is that these equations
describe well the short-range correlations even inside the radius of
interaction potential. This confirms that the GGP equations can be considered
as a universal tool for describing BECs in various regimes. For instance, the
method has good prospects for studying the extremely tight confinements,
crossover between 3D and 2D or 1D regimes, and the dipolar BEC in various
regimes, including quantum droplets [29].

In addition, the integral representation for GF was found for the cylindrical
confinement [Eqs. (9), (11)], which is valid for small scattering wavevectors
$p<2/l_{0}$ and arbitrary values of coordinates. It is convenient for
obtaining the low-density expansion (13) of the GF.

Within our approach, we found the position (44) of the smoothed resonance at
positive scattering length in the pancake geometry. The position is
independent of the density. This is nearly two-third of the experimentally
observed value [5]. The discrepancy might be attributed to condensate
depletion at finite temperatures or to many-body effects that are beyond the
used weak-coupling approximation.

###### Acknowledgements.

The author is grateful to Azis Embareck for useful and fruitful discussions
and Vladimir Melezhik for careful reading of the manuscript and useful
comments.

##  Appendix A Hartree-Fock-Bogoliubov approximation

Let make some remarks about the Hartree-Fock-Bogoliubov (HFB) method. We
consider only the stationary HFB theory at zero temperature. The HFB
approximation for fermions was suggested by Bogoliubov [30] as a variational
scheme formulated in terms of the normal and anomalous averages of two Fermi-
field operators. These scheme can be simplified [31, 32, 33] if we note that
these averages are given by bilinear forms of $u$ and $v$ parameters of the
Bogoliubov transformation. Then the variational procedure can be reformulated
directly in terms of $u$ and $v$ parameters, and the resulting equations are
called Bogoliubov-de Gennes equations.

For bosons, the relations analogous to the Bogoliubov-de Gennes equations can
be written down [12, 34]. They are called the Bogoliubov equations. In
addition to the $u$ and $v$ parameters, these equations contain the
wavefunction of the inhomogeneous Bose-Einstein condensate.

On the other hand, one can also formulate the HFB method for bosons as a
variational scheme with three variational functions: the condensate
wavefunction and normal and anomalous averages of two Bose-field operators.
Below we follow this scheme [16] and derive the generalized GP equations by
further simplifying the HFB scheme.

###  A.1 Hartree-Fock-Bogoliubov approximation for bosons in the variational
formulation

The ground-state energy and mean number of interacting bosons in an external
field is given by (see, e.g., textbooks [35, 36])

$\displaystyle E=$ $\displaystyle-\frac{\hslash^{2}}{2m}\int
dx\left.\nabla^{2}_{x}\langle{\hat{\Psi}}^{\dagger}(x^{\prime}){\hat{\Psi}}(x)\rangle\right|_{x^{\prime}=x}$
$\displaystyle+\int dx\,V_{\rm
ext}(x)\langle{\hat{\Psi}}^{\dagger}(x){\hat{\Psi}}(x)\rangle$
$\displaystyle+\frac{1}{2}\int\\!d{x}_{1}d{x}_{2}\,V({x}_{1},{x}_{2})\langle{\hat{\Psi}}^{\dagger}({x}_{1}){\hat{\Psi}}^{\dagger}({x}_{2}){\hat{\Psi}}({x}_{2}){\hat{\Psi}}({x}_{1})\rangle,$
(45) $\displaystyle N=$ $\displaystyle\int
dx\,\langle{\hat{\Psi}}^{\dagger}(x){\hat{\Psi}}(x)\rangle,$ (46)

respectively. Here $x=({\bf r},\sigma)$ are the coordinate and spin of a
particle, respectively, and $\int dx\cdots=\sum_{\sigma}\int d{\bf r}\cdots$.
The brackets $\langle\cdots\rangle$ stand for the ground-state average. The
Bose field operators are denoted as ${\hat{\Psi}}(x)$ and
${\hat{\Psi}}^{\dagger}(x)$.

The Bose-Einstein condensation implies that the global gauge symmetry is
broken. Then one can separate the condensate contribution to the field
operators: $\hat{\Psi}(x)=\phi(x)+\hat{\vartheta}(x)$ and
$\hat{\Psi}^{{\dagger}}(x)=\phi^{*}(x)+\hat{\vartheta}^{{\dagger}}(x)$. Here
$\phi(x)$ and $\hat{\vartheta}(x)$ are the numerical and operator parts,
respectively, for which we have $\phi(x)=\langle\hat{\Psi}(x)\rangle$ and
$\langle\hat{\vartheta}(x)\rangle=0$. Then the two-boson correlators read

$\displaystyle\langle{\hat{\Psi}}^{\dagger}(x){\hat{\Psi}}(x^{\prime})\rangle=\,$
$\displaystyle\phi^{*}(x)\phi(x^{\prime})+\langle\hat{\vartheta}^{\dagger}(x)\hat{\vartheta}(x^{\prime})\rangle,$
(47) $\displaystyle\langle{\hat{\Psi}}(x){\hat{\Psi}}(x^{\prime})\rangle=\,$
$\displaystyle\phi(x)\phi(x^{\prime})+\langle\hat{\vartheta}(x)\hat{\vartheta}(x^{\prime})\rangle.$
(48)

The normal two-boson correlator (47) is nothing else but the one-body density
matrix, and the anomalous correlator (48) can be treated as the two-body
wavefunction of the Bose-Einsten condensate [21, 22].

The diagonal part of the two-body density matrix, which determines the
interaction term, can also be found by separating the condensate part and
using Wick’s theorem for the correlators of $\hat{\vartheta}$ and
$\hat{\vartheta}^{\dagger}$ operators

$\displaystyle\langle{\hat{\Psi}}^{\dagger}({x}_{1}){\hat{\Psi}}^{\dagger}({x}_{2}){\hat{\Psi}}({x}_{2}){\hat{\Psi}}({x}_{1})\rangle=\phi^{*}(x_{1})\phi(x_{2})^{*}\phi(x_{2})\phi(x_{1})$
$\displaystyle+\phi(x_{1})\phi(x_{2})\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}^{\dagger}(x_{2})\rangle+\phi^{*}(x_{1})\phi^{*}(x_{2})\langle\hat{\vartheta}(x_{2})\hat{\vartheta}(x_{1})\rangle$
$\displaystyle+\phi(x_{1})\phi^{*}(x_{2})\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}(x_{2})\rangle+\phi^{*}(x_{1})\phi(x_{2})\langle\hat{\vartheta}^{\dagger}(x_{2})\hat{\vartheta}(x_{1})\rangle$
$\displaystyle+\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}^{\dagger}(x_{2})\rangle\langle\hat{\vartheta}(x_{2})\hat{\vartheta}(x_{1})\rangle$
$\displaystyle+\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}(x_{1})\rangle\langle\hat{\vartheta}^{\dagger}(x_{2})\hat{\vartheta}(x_{2})\rangle$
$\displaystyle+\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}(x_{2})\rangle\langle\hat{\vartheta}^{\dagger}(x_{2})\hat{\vartheta}(x_{1})\rangle.$
(49)

Wick’s theorem is applicable, since the effective HFB Hamiltonian contains
only a bilinear form of creation and annihilation operators [34].

The normal
$\langle\hat{\vartheta}^{\dagger}(x)\hat{\vartheta}(x^{\prime})\rangle$ and
anomalous $\langle\hat{\vartheta}(x)\hat{\vartheta}(x^{\prime})\rangle$
averages are not independent quantities. The relations between them at zero
temperature can be obtained [16] using the fact that the HFB ground state is
the quasiparticle vacuum. The relations do not depend explicitly on parameters
of the Hamiltonian and take the form

$\displaystyle\int
dx\,\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}^{\dagger}(x)\rangle\langle\hat{\vartheta}(x)\hat{\vartheta}(x_{2})\rangle=\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}(x_{2})\rangle$
$\displaystyle\phantom{\int dx\,\langle\hat{\vartheta}^{\dagger}}+\int
dx\,\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}(x)\rangle\langle\hat{\vartheta}^{\dagger}(x)\hat{\vartheta}(x_{2})\rangle,$
(50) $\displaystyle\int
dx\,\langle\hat{\vartheta}^{\dagger}(x)\hat{\vartheta}(x_{1})\rangle\langle\hat{\vartheta}(x)\hat{\vartheta}(x_{2})\rangle=\int
dx\,\langle\hat{\vartheta}(x_{1})\hat{\vartheta}(x)\rangle$
$\displaystyle\phantom{\int
dx\,\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}(x)\rangle\langle\hat{\vartheta}(x)}\times\langle\hat{\vartheta}^{\dagger}(x)\hat{\vartheta}(x_{2})\rangle.$
(51)

Finally, the stationary HFB equations can be obtained by means of variation of
$E-\mu N$, given by Eqs. (45)-(49), under the conditions (50) and (51). The
variational functions are $\phi(x)$, $\phi^{*}(x)$,
$\langle\hat{\vartheta}(x)\hat{\vartheta}(x^{\prime})\rangle$,
$\langle\hat{\vartheta}^{\dagger}(x)\hat{\vartheta}^{\dagger}(x^{\prime})\rangle$,
and $\langle\hat{\vartheta}^{\dagger}(x)\hat{\vartheta}(x^{\prime})\rangle$.
We do not write explicitly the resulting equations.

The full HBF scheme is rather cumbersome, and additional simplifications of
Eq. (49) are often used [37]. The simplest approximation is to neglect all the
quantum correlators and leave only the wavefunction of the condensate. Then we
arrive at the well-known GP functional with the help of the pseudopotential.

###  A.2 Generalized Gross-Pitaevskii equations as a simplified version of the
Hartee-Fock-Bogoliubov approximation

The GGP equations can be obtained [16] from the full HFB theory when the
normal two-boson terms are completely neglected in Eq. (49)

$\displaystyle\langle{\hat{\Psi}}^{\dagger}({x}_{1}){\hat{\Psi}}^{\dagger}({x}_{2}){\hat{\Psi}}({x}_{2}){\hat{\Psi}}({x}_{1})\rangle\simeq\phi^{*}(x_{1})\phi(x_{2})^{*}\phi(x_{2})\phi(x_{1})$
$\displaystyle+\phi(x_{1})\phi(x_{2})\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}^{\dagger}(x_{2})\rangle+\phi^{*}(x_{1})\phi^{*}(x_{2})\langle\hat{\vartheta}(x_{2})\hat{\vartheta}(x_{1})\rangle$
$\displaystyle+\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}^{\dagger}(x_{2})\rangle\langle\hat{\vartheta}(x_{2})\hat{\vartheta}(x_{1})\rangle.$

Besides, if the condensate depletion is small, one can neglect the second term
in the r.h.s of Eq. (50), which is of the next order of the condensate
depletion. Thus, we are left with the relation

$\displaystyle\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}(x_{2})\rangle\simeq\int
dx\,\langle\hat{\vartheta}^{\dagger}(x_{1})\hat{\vartheta}^{\dagger}(x)\rangle\langle\hat{\vartheta}(x)\hat{\vartheta}(x_{2})\rangle.$
(52)

In the approximation (52), Eq. (51) turns into identity.

Then the variational procedure yields [16] the GGP equations (1)-(4).

## References

  * Olshanii [1998] M. Olshanii, Atomic scattering in the presence of an external confinement and a gas of impenetrable bosons, Phys. Rev. Lett. 81, 938 (1998).
  * Bergeman  _et al._ [2003] T. Bergeman, M. G. Moore, and M. Olshanii, Atom-atom scattering under cylindrical harmonic confinement: Numerical and analytic studies of the confinement induced resonance, Phys. Rev. Lett. 91, 163201 (2003).
  * Petrov  _et al._ [2000a] D. S. Petrov, M. Holzmann, and G. V. Shlyapnikov, Bose-Einstein condensation in quasi-2D trapped gases, Phys. Rev. Lett. 84, 2551 (2000a).
  * Petrov and Shlyapnikov [2001] D. S. Petrov and G. V. Shlyapnikov, Interatomic collisions in a tightly confined Bose gas, Phys. Rev. A 64, 012706 (2001).
  * Haller  _et al._ [2010] E. Haller, M. J. Mark, R. Hart, J. G. Danzl, L. Reichsöllner, V. Melezhik, P. Schmelcher, and H.-C. Nägerl, Confinement-induced resonances in low-dimensional quantum systems, Phys. Rev. Lett. 104, 153203 (2010).
  * Naidon  _et al._ [2007] P. Naidon, E. Tiesinga, W. F. Mitchell, and P. S. Julienne, Effective-range description of a bose gas under strong one- or two-dimensional confinement, New J. Phys. 9, 19 (2007).
  * Peng  _et al._ [2010] S.-G. Peng, S. S. Bohloul, X.-J. Liu, H. Hu, and P. D. Drummond, Confinement-induced resonance in quasi-one-dimensional systems under transversely anisotropic confinement, Phys. Rev. A 82, 063633 (2010).
  * Zhang and Zhang [2011] W. Zhang and P. Zhang, Confinement-induced resonances in quasi-one-dimensional traps with transverse anisotropy, Phys. Rev. A 83, 053615 (2011).
  * Melezhik and Schmelcher [2011] V. S. Melezhik and P. Schmelcher, Multichannel effects near confinement-induced resonances in harmonic waveguides, Phys. Rev. A 84, 042712 (2011).
  * Lopes  _et al._ [2017] R. Lopes, C. Eigen, N. Navon, D. Clément, R. P. Smith, and Z. Hadzibabic, Quantum depletion of a homogeneous bose-einstein condensate, Phys. Rev. Lett. 119, 190404 (2017).
  * Navon  _et al._ [2021] N. Navon, R. P. Smith, and Z. Hadzibabic, Quantum gases in optical boxes, Nat. Phys. 17, 1334 (2021).
  * Pitaevskii [1961] L. P. Pitaevskii, Vortex lines in an imperfect Bose gas, Zh. Eksp. Teor. Fiz. 40, 646 (1961), [Sov. Phys. JETP 13, 451 (1961)].
  * Gross [1961] E. P. Gross, Structure of a quantized vortex in boson systems, Nuovo Cim. 20, 454 (1961).
  * Gross [1963] E. P. Gross, Hydrodynamics of a superfluid condensate, J. Math. Phys. 4, 195 (1963).
  * Pitaevskii and Stringari [2016] L. Pitaevskii and S. Stringari,  _Bose-Einstein Condensation and Superfluidity_ (Oxford University, Oxford, 2016).
  * Cherny and Brand [2004] A. Yu. Cherny and J. Brand, Self-consistent calculation of the coupling constant in the Gross-Pitaevskii equation, Phys. Rev. A 70, 043622 (2004).
  * Idziaszek and Calarco [2006] Z. Idziaszek and T. Calarco, Analytical solutions for the dynamics of two trapped interacting ultracold atoms, Phys. Rev. A 74, 022712 (2006).
  * Idziaszek  _et al._ [2015] Z. Idziaszek, K. Jachymski, and P. S. Julienne, Reactive collisions in confined geometries, New J. Phys. 17, 035007 (2015).
  * Liang and Zhang [2008] J.-J. Liang and C. Zhang, Two ultracold atoms in a completely anisotropic trap, Phys. Scr. 77, 025302 (2008).
  * Jeszenszki  _et al._ [2018] P. Jeszenszki, A. Yu. Cherny, and J. Brand, $s$-wave scattering length of a Gaussian potential, Phys. Rev. A 97, 042708 (2018).
  * Cherny [2000] A. Yu. Cherny, Pair wave functions in a Bose liquid (2000), arXiv:cond-mat/9807120 [cond-mat.stat-mech] .
  * Cherny and Shanenko [2000] A. Yu. Cherny and A. A. Shanenko, Short-range particle correlations in a dilute Bose gas, Phys. Rev. E 62, 1646 (2000).
  * Cherny and Shanenko [2001a] A. Yu. Cherny and A. A. Shanenko, Dilute bose gas: short-range particle correlations and ultraviolet divergence, Eur. Phys. J. B 19, 555 (2001a).
  * Cherny and Shanenko [2001b] A. Yu. Cherny and A. A. Shanenko, Dilute Bose gas in two dimensions: Density expansions and the Gross-Pitaevskii equation, Phys. Rev. E 64, 027105 (2001b).
  * Erdélyi  _et al._ [1981] A. Erdélyi, W. Magnus, F. Oberhettinger, and F. G. Tricomi,  _Higher Transcendental Functions_ , Vol. 1 (Krieger, N.Y., 1981) pp. 27-31.
  * Landau and Lifshitz [1977] L. D. Landau and E. M. Lifshitz,  _Quantum Mechanics (Non-relativistic Theory)_ , 3rd ed., Course of Theoretical Physics, Vol. 3 (Pergamon, N.Y., 1977).
  * Abramowitz and Stegun [1974] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables_ (Dover Publications, Incorporated, 1974).
  * Petrov  _et al._ [2000b] D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Walraven, Regimes of quantum degeneracy in trapped 1D gases, Phys. Rev. Lett. 85, 3745 (2000b).
  * Böttcher  _et al._ [2020] F. Böttcher, J.-N. Schmidt, J. Hertkorn, K. S. H. Ng, S. D. Graham, M. Guo, T. Langen, and T. Pfau, New states of matter with fine-tuned interactions: quantum droplets and dipolar supersolids, Rep. Prog. Phys. 84, 012403 (2020).
  * Bogoliubov [1959] N. N. Bogoliubov, The compensation principle and the self-consistent field method, Sov. Phys. Usp. 2, 236 (1959).
  * de Gennes [1964] P. G. de Gennes, Boundary effects in superconductors, Rev. Mod. Phys. 36, 225 (1964).
  * de Gennes [1966] P. G. de Gennes,  _Superconductivity of Metals and Alloys_ (Benjamin, N. Y., 1966).
  * Baranger [1963] M. Baranger, Reformulation of the theory of pairing correlations, Phys. Rev. 130, 1244 (1963).
  * Fetter [1972] A. L. Fetter, Nonuniform states of an imperfect Bose gas, Ann. Phys. 70, 67 (1972).
  * Feynman [1972] R. P. Feynman,  _Statistical mechanics: a set of lectures_ , Frontiers in physics: a lecture note and reprint series (W.A. Benjamin, Reading, Mass., 1972) Sec. 6.8.
  * Bogoliubov and Bogoliubov [1982] N. N. Bogoliubov and N. N. Bogoliubov, Jr.,  _Introduction to Quantum Statistical Mechanics_ (World Scientific, Singapore, 1982) p. 219.
  * Griffin [1996] A. Griffin, Conserving and gapless approximations for an inhomogeneous Bose gas at finite temperatures, Phys. Rev. B 53, 9341 (1996).