#  Effects of quantum fluctuations on the low-energy collective modes
of two-dimensional superfluid Fermi gases from the BCS to the Bose Limit

Senne Van Loon<EMAIL_ADDRESS>School of Physics, Georgia
Institute of Technology, Atlanta, Georgia 30332, USA  TQC, Universiteit
Antwerpen, Universiteitsplein 1, B-2610 Antwerpen, België  C. A. R. Sá de Melo
School of Physics, Georgia Institute of Technology, Atlanta, Georgia 30332,
USA

###### Abstract

We investigate the effects of quantum fluctuations on the low-energy
collective modes of two-dimensional (2D) Fermi superfluids from the BCS to the
Bose limit. We compare our results to recent Bragg scattering experiments in
2D box potentials, with excellent agreement. We show that quantum fluctuations
in the phase and modulus of the superfluid order parameter are absolutely
necessary to give physically acceptable chemical potential and dispersion
relation of the low-energy collective mode throughout the BCS to Bose
evolution. Furthermore, we demonstrate that the dispersion of the collective
modes change from concave to convex as interactions are tuned from the BCS to
the Bose regime.

The study of collective modes is a fundamental component of many-particle
physics, because for every spontaneously broken continuous symmetry there are
low-energy modes that emerge as expected from Goldstone’s theorem Goldstone
(1961), and additional higher-energy excitations such as the Higgs mode Higgs
(1966); Englert  _et al._ (1966). Collective modes are essential in
understanding a variety of systems ranging from condensed matter (quantum
magnets, superconductors) Hong  _et al._ (2017); Maiti and Hirschfeld (2015),
high energy physics (standard nuclear matter, quantum chromodynamics) Greco
_et al._ (2003); Yamamoto  _et al._ (2007), and astrophysics (neutron stars,
black holes) Kobyakov and Pethick (2013); Mathur and Turton (2014) to atomic
(Bose and Fermi superfluids) physics Singh and Mathey (2021); Choudhury  _et
al._ (2020). Unfortunately, in condensed matter it is not easy to tune
parameters such as interactions, density, and dimensionality over a wide
range, in high energy physics it is very difficult, and in astrophysics it is
impossible. However, in atomic physics this is relatively easy via well
established techniques Chin  _et al._ (2010); Ibarra-García-Padilla  _et al._
(2020). This makes it possible to investigate collective modes in superfluids,
particularly important because they reveal the effects of quantum fluctuations
above the superfluid ground state.

Superfluids in 2D are inherently different from their 3D counterparts, due to
the importance of fluctuations Mermin and Wagner (1966); Hohenberg (1967)
leading to a Berezinskii-Kosterlitz-Thouless (BKT) transition Berezinskii
(1972); Kosterlitz and Thouless (1972). In the context of ultracold atoms, the
desire to study 2D Fermi superfluids is driven not only by connections to
high-temperature superconductors Keimer  _et al._ (2015); Ge  _et al._ (2015);
Cao  _et al._ (2018); Yu  _et al._ (2019), but also by the high degree of
experimental control that allows the measurement of the equation of state
Makhalov  _et al._ (2014); Fenech  _et al._ (2016); Boettcher  _et al._
(2016), the observation of the BKT transition Hadzibabic  _et al._ (2006);
Murthy  _et al._ (2015), and the examination of collective modes Ville  _et
al._ (2018); Bohlen  _et al._ (2020); Christodoulou  _et al._ (2021); Sobirey
_et al._ (2021).

Ultracold fermions with tunable interactions in nearly 2D configurations were
studied using harmonic traps and optical lattices Martiyanov  _et al._ (2010);
Feld  _et al._ (2011); Dyke  _et al._ (2011); Ong  _et al._ (2015). With the
very recent advent of box potentials, it is now possible to study
experimentally homogeneous 2D fermions Hueck  _et al._ (2018); Bohlen  _et
al._ (2020); Sobirey  _et al._ (2021). Inspired by recent measurements of
collective excitations via Bragg-spectroscopy Hoinka  _et al._ (2017); Biss
_et al._ (2021); Sobirey  _et al._ (2021), we investigate the low-energy
collective modes of 2D Fermi superfluids in box potentials. We find that mean
field (saddle point) results in 2D Randeria  _et al._ (1990) produce incorrect
values of the chemical potential and lead to the erroneous conclusion that the
speed of sound is a constant throughout the BCS to Bose evolution Marini  _et
al._ (1998); Salasnich  _et al._ (2013); Lumbeeck  _et al._ (2020). In sharp
contrast, we show that the inclusion of quantum fluctuations He  _et al._
(2015) is crucial to produce physically acceptable results for the dispersion
of collective modes and leads to a nonconstant speed of sound in the crossover
from BCS to Bose regimes at low temperatures. Furthermore, we demonstrate that
phase and modulus fluctuations of the superfluid order parameter become
increasingly more coupled with growing interaction strength, leading to a
change in curvature of the collective mode dispersion from concave in the BCS
to convex in the Bose regime. A comparison of our results to very recent
experiments Sobirey  _et al._ (2021) shows excellent agreement.

Hamiltonian: To analyze the low-energy collective modes of 2D Fermi
superfluids in box potentials, we start from the Hamiltonian density

${\cal H}=\psi^{\dagger}_{s}({\bf
r})\frac{(-i\hbar\nabla)^{2}}{2m}\psi_{s}({\bf
r})-g\psi^{\dagger}_{\uparrow}({\bf r})\psi^{\dagger}_{\downarrow}({\bf
r})\psi_{\downarrow}({\bf r})\psi_{\uparrow}({\bf r}),$ (1)

where $\psi_{s}({\bf r})$ is a fermion field operator with spin $s$ at
position ${\bf r}$. The first term is the kinetic energy and the second
represents local attractive interactions. The associated action is ${\cal
S}(\psi^{\dagger},\psi)=\int
d^{3}r\\{\psi^{\dagger}(r)\left[\hbar\partial_{\tau}-\mu\right]\psi(r)+{\cal
H}(r)\\},$ where $r=({\bf r},\tau),\,\int d^{3}r=\int_{0}^{\beta}d\tau\int
d^{2}{\bf r},~{}\,\beta=\hbar/k_{B}T$, and $\mu$ is the chemical potential.
The grand canonical partition function of the system is ${\cal Z}=\int
D\psi^{\dagger}D\psi e^{-{\cal S}/\hbar}$. We introduce the Hubbard-
Stratonovich complex pair field $\Phi(r)$ to decouple the contact interactions
and integrate out the fermionic fields to obtain an effective action $S_{\rm
eff}(\Phi^{\dagger},\Phi)$. We write
$\Phi(r)=|\Phi(r)|e^{\mathrm{i}\theta(r)}$ in terms of its modulus
$|\Phi(r)|=|\Delta|\left[1+\lambda(r)\right]$ and phase $\theta(r)$, and
expand $S_{\rm eff}(\Phi^{\dagger},\Phi)$ up to quadratic order in both the
phase $\theta(r)$ and modulus fluctuations $\lambda(r)$ around the saddle
point $|\Phi_{\rm sp}(r)|=|\Delta|$. The resulting Gaussian action is

$S_{\rm eff}=S_{\rm
sp}+\beta|\Delta|^{2}\sum_{q}\begin{pmatrix}\mathrm{i}\theta_{-q}&\lambda_{-q}\end{pmatrix}M(q)\begin{pmatrix}-\mathrm{i}\theta_{q}\\\
\lambda_{q}\end{pmatrix},$ (2)

where $q=({\bf q},\mathrm{i}\nu)$ and $\nu=2\pi n/\beta$ are bosonic Matsubara
frequencies, $S_{\rm sp}=\beta\sum_{\bf{k}}(\xi_{\bf{k}}-E_{\bf{k}})+\beta
L^{2}|\Delta|^{2}/g$ is the saddle point action, and $M(q)$ is the symmetric
Gaussian fluctuation matrix Engelbrecht  _et al._ (1997). Using the analytic
continuation $\mathrm{i}\nu\to\omega+\mathrm{i}\delta$, the matrix elements of
$M(q)$, at zero temperature, are

$\displaystyle\frac{M_{\\!{}_{\pm\pm}}}{L^{2}}=$
$\displaystyle\int\\!\frac{d^{2}{\bf
k}}{4\pi^{2}}\\!\left[\frac{E_{+}+E_{-}}{2E_{+}E_{-}}\frac{E_{+}E_{-}+\xi_{+}\xi_{-}\pm|\Delta|^{2}}{\hbar^{2}\omega^{2}-(E_{+}+E_{-})^{2}}\\!+\\!\frac{1}{2E_{\bf
k}}\right],$ $\displaystyle\frac{M_{\\!{}_{+-}}}{L^{2}}=$
$\displaystyle\int\\!\frac{d^{2}{\bf
k}}{4\pi^{2}}\\!\left[\frac{\hbar\omega}{2E_{+}E_{-}}\frac{E_{+}\xi_{-}+\xi_{+}E_{-}}{\hbar^{2}\omega^{2}-(E_{+}+E_{-})^{2}}\right],$

where $E_{\bf k}=\sqrt{\xi_{\bf k}^{2}+|\Delta|^{2}}$ is the quasiparticle
energy, $\xi_{\bf k}=\epsilon_{\bf k}-\mu$ is the energy of the free fermions
of mass $m$ $(\epsilon_{\bf k}=\hbar^{2}|{\bf k}|^{2}/2m)$ with respect to
$\mu$, and $|\Delta|$ is the modulus of the order parameter. We have used the
shorthand notations $E_{\pm}=E_{{\bf{k}}\pm{\bf{q}}/2}$ and
$\xi_{\pm}=\xi_{{\bf{k}}\pm{\bf{q}}/2}$, $M_{++}=M_{\theta\theta}$,
$M_{--}=M_{\lambda\lambda}$, and $M_{+-}=M_{\theta\lambda}$.

Equation of state: In Eq. (2), the action $S_{\rm eff}$ is fully characterized
by $|\Delta|$ and $\mu$ or by the dimensionless parameters $x=\mu/|\Delta|$
and $|\Delta|/\epsilon_{F}$, where $\epsilon_{F}$ is the Fermi energy for a
specified density $n=k_{F}^{2}/2\pi$, with $k_{F}$ the Fermi momentum.
However, to study the evolution from the BCS to the Bose limit, it is
experimentally more relevant to relate $|\Delta|$ and $\mu$ to the 2D
scattering length $a$ and the density $n$. The order parameter is found from
$\left[\partial\Omega_{\rm sp}/\partial|\Delta|\right]_{T,V}=0$, where
$\Omega_{\rm sp}=S_{\rm sp}/\beta$ is the saddle-point thermodynamic
potential. Replacing the interaction strength $g$ in favor of the (positive)
two-body binding energy $\epsilon_{\rm b}$, using the Lippmann-Schwinger
relation $L^{2}/g=\sum_{\bf k}1/(2\epsilon_{\bf k}+\epsilon_{b})$ Botelho and
Sá de Melo (2006), one finds
$|\Delta|=\sqrt{\epsilon_{b}(2\mu+\epsilon_{b})}\Theta(2\mu+\epsilon_{b})$,
which is explicitly only a function of $\mu$ and $\epsilon_{b}$. The chemical
potential can be found by solving the saddle point number equation $n_{\rm
sp}=-\left[\partial\Omega_{\rm sp}/\partial\mu\right]_{T,V}/L^{2}$, while
fixing the density $n_{\rm sp}=n=k_{F}^{2}/2\pi$, resulting in $\mu_{\rm
sp}=\epsilon_{F}-\epsilon_{b}/2$, which substituted in the order parameter
relation leads to $|\Delta|_{\rm sp}=\sqrt{2\epsilon_{F}\epsilon_{b}}$. These
expressions are connected to the 2D scattering length $a$ via the relation
$\epsilon_{b}=8\epsilon_{F}/\exp(2\gamma_{\mathrm{E}}+2\ln k_{F}a)$, with
$\gamma_{\mathrm{E}}\approx 0.577$ the Euler–Mascheroni constant Petrov and
Shlyapnikov (2001). However, as discussed below, this analysis leads to the
unphysical result of a constant speed of sound $c=v_{F}/\sqrt{2}$, where
$v_{F}$ is the Fermi velocity, over the entire BCS-to-Bose evolution, because
$\mu$ is inaccurately calculated.

A more precise determination of $\mu$ for fixed density $n$ requires not only
$\Omega_{\rm sp}$, but also the Gaussian contribution $\Omega_{\rm
g}=(k_{B}T/2)\sum_{q}\ln\det[|\Delta|M(q)]$, found by integrating the bosonic
fields in the second term of Eq. (2). The full number equation $n=n_{\rm
sp}+n_{g}$, including fluctuations of the order parameter at the Gaussian
level, is necessary to determine $\mu$. This is important both in 3D Nozières
and Schmitt-Rink (1985); Sá de Melo  _et al._ (1993) and 2D Botelho and Sá de
Melo (2006); He  _et al._ (2015). Here,
$n_{g}=-\left[\partial\Omega_{\mathrm{g}}/\partial\mu\right]_{T,V}/L^{2}$ must
be always positive and $n(\mu)=k_{F}^{2}/2\pi$ must be solved numerically,
leading to a significant reduction of $\mu$ in the Bose regime 111Note that
one has to introduce convergence factors to regularize the Matsubara sum in
the thermodynamic potential $\Omega_{\mathrm{g}}$. This is most easily done in
the basis of the real and imaginary part of the pair field $\Phi$, instead of
the phase-modulus basis used here Engelbrecht  _et al._ (1997); He  _et al._
(2015)., see Fig. 1.

Figure 1:  The ratio $x=\mu/|\Delta|$ for different approximations of the
equation of state: the saddle point approximation (dashed blue line) and
including Gaussian fluctuations (solid red line). The BCS and Bose limits are
shown as dotted black and dot-dashed green lines.

Collective Modes: As seen from Eq. (2), the fluctuation matrix $M(q)$ acts as
the inverse propagator of modulus and phase fluctuations. The collective mode
frequency $\omega_{\bf q}$ is found from the poles of $\left[M(q)\right]^{-1}$
or equivalently from $\det M({\bf{q}},\omega_{\bf q})=0$ Engelbrecht  _et al._
(1997). A real solution is found below the two-particle continuum
$\epsilon_{c}({\bf{q}})=\min_{\bf{k}}(E_{+}+E_{-})$, above which it is
energetically more favorable to break pairs. In a 3D Fermi gas, $\omega_{\bf
q}$ can hit the two-particle continuum at some finite value of ${\bf q}$
causing damping Combescot  _et al._ (2006). However, for a 2D Fermi gas, a
real $\omega_{\bf q}<\epsilon_{c}({\bf q})$ is found for all values of ${\bf
q}$, given that a two-body bound state always exists for a 2D contact
potential Landau and Lifshitz (1977). This physics arises from
$M_{++}=M_{\theta\theta}$, which always diverges when
$\hbar\omega\rightarrow\epsilon_{c}$ making $\det M({\bf{q}},\epsilon_{c})=0$
impossible, and thus there is no damping of the mode at the Gaussian level.
For arbitrary ${\bf q}$ there is no general analytical solution for
$\omega_{\bf q}$, but we obtain numerical results shown as solid red lines in
Fig. 2.

Figure 2:  Frequency $\omega_{\bf q}$ vs. $|{\bf q}|$, using Fermi units
$\omega_{F}=E_{F}/\hbar$ and $k_{F}$, for different interactions. The solid
red line represents the numerical solution for $\omega_{\bf q}$, while
analytical approximations are shown up to linear (dotted blue line), cubic
(dashed green line), and quintic (dot-dashed magenta line) order in $|{\bf
q}|$. The parameters for each panel are: (a) $\ln
k_{F}a=-0.1,\,(\mu/|\Delta|\simeq-1.0)$, (b) $\ln
k_{F}a=0.7,\,(\mu/|\Delta|\simeq 0.2)$, (c) $\ln
k_{F}a=1.1,\,(\mu/|\Delta|\simeq 0.7)$, and (d) $\ln
k_{F}a=1.4,\,(\mu/|\Delta|\simeq 1.2)$, using the Gaussian value of $\mu$. The
double-dotted-dashed black lines represent $\omega_{\bf q}^{\mathrm{PO}}$
using the saddle point $\mu$. The gray area indicates the two-particle
continuum $\varepsilon\geq\epsilon_{c}({\bf{q}})$ at the Gaussian level. The
interactions are the same as used in experiments Sobirey  _et al._ (2021).

Although numerical solutions are useful for comparison to recent experiments
Bohlen  _et al._ (2020); Sobirey  _et al._ (2021), analytical insight is
essential to understand the underlying physics. To reveal the interplay
between modulus and phase fluctuations, we show next that their coupling
increases dramatically as the system evolves from the BCS to the Bose regime.
We invert $M(q)$ and obtain the propagators in Fourier space, which in the
long wavelength limit become

$\displaystyle\begin{pmatrix}\langle\theta_{-q}\theta_{q}\rangle&\langle\theta_{-q}\lambda_{q}\rangle\\\
\langle\lambda_{-q}\theta_{q}\rangle&\langle\lambda_{-q}\lambda_{q}\rangle\end{pmatrix}$
(3)
$\displaystyle\quad=\frac{|\Delta|^{2}}{\hbar^{2}c^{2}q^{2}-\hbar^{2}\omega^{2}}\begin{pmatrix}\chi_{\theta\theta}&\mathrm{i}\hbar\omega\,\chi_{\theta\lambda}\\\
-\mathrm{i}\hbar\omega\,\chi_{\theta\lambda}&{\bf{q}}^{2}\,\chi^{\bf{q}}_{\lambda\lambda}-\hbar^{2}\omega^{2}\,\chi^{\omega}_{\lambda\lambda}\end{pmatrix}.$

The different $\chi$ coefficients defined in this equation are shown in Fig.
3(a) as a function of interaction parametrized by the ratio $x=\mu/|\Delta|$.
They are explicitly given by
$\chi_{\theta\theta}=4\pi,\,\chi_{\theta\lambda}=2\pi(-x+\sqrt{1+x^{2}})/|\Delta|,\,\chi^{\bf{q}}_{\lambda\lambda}=\hbar^{2}\pi\sqrt{1+x^{2}}/2m|\Delta|$,
and $\chi^{\omega}_{\lambda\lambda}=\pi/|\Delta|^{2}$. Most notably, the solid
red line in Fig. 3(a) shows $\chi_{\theta\lambda}$, which controls the
coupling between phase and modulus. Notice that $\chi_{\theta\lambda}$ is
large in the Bose regime $(\mu/|\Delta|\ll-1)$, indicating that phase and
modulus are strongly mixed, while it is negligible in the BCS regime
$(\mu/|\Delta|\gg 1)$, showing that phase and modulus are essentially
decoupled. For $\omega>0$, the pole of $\left[M(|{\bf q}|,\omega)\right]^{-1}$
occurs at $\omega_{\bf q}=c|{\bf q}|$, with

$2mc^{2}=\mu+\sqrt{\mu^{2}+|\Delta|^{2}},$ (4)

where $c$ is the speed of sound. This expression includes modulus and phase
fluctuations of the order parameter, and explicitly illustrates the dependence
of $c$ on $\mu$. In Fig. 3(b), we show the behavior of $c/v_{F}$ at three
levels of approximation: the dotted green line includes only phase
fluctuations, the dashed blue line includes phase and modulus fluctuations
with the saddle point value of $\mu$, while the solid red line includes phase
and modulus fluctuations with a self-consistent Gaussian fluctuation value of
$\mu$. It is clear that the first two levels of approximation Marini  _et al._
(1998); Salasnich  _et al._ (2013) give completely unphysical results, in
particular predicting the constant value $c=v_{F}/\sqrt{2}$ for any coupling
in the saddle point approximation, as quoted in the literature Marini  _et
al._ (1998); Salasnich  _et al._ (2013). However, we show here that the
inclusion of phase and modulus fluctuations with the correct $\mu$ leads to
the appropriate behavior of $c$ both in the Bose and BCS limits, leading to
physical results in excellent agreement to recent experiments of 2D fermions
in a box potential Bohlen  _et al._ (2020); Sobirey  _et al._ (2021).

Figure 3:  (a) Different $\chi$ coefficients vs. $\mu/|\Delta|$ or $\ln
k_{F}a$ as defined in Eq. (3). (b) Speed of sound $c/v_{F}$ vs. $\ln k_{F}a$:
the dotted green line includes phase-only fluctuations with saddle point
$\mu$, diverging in the Bose limit; blue dashed line includes phase and
modulus fluctuations with saddle point $\mu$, giving always a constant value;
solid red line combines phase and modulus fluctuations with Gaussian $\mu$
self-consistently. The dot-dashed magenta and black lines show the results in
the BCS and Bose limits, respectively.

Further insight is gained by studying the long wavelength limit of $S_{\rm
eff}$ in Eq. (2) by expanding the fluctuation matrix $M({\bf q},\omega)$ for
small ${\bf{q}}$ and $\omega$. Performing an inverse Fourier transform back to
real space, the second term in Eq. (2) reduces to

$S_{g}=\frac{1}{2}\int
d^{3}r\left[\rho_{s}(\nabla\theta)^{2}+A(\hbar\partial_{\tau}\theta)^{2}+\mathrm{i}D\lambda\hbar\partial_{\tau}\theta+C\lambda^{2}\right].$
(5)

The first two coefficents are the $T=0$ superfluid density
$2m\rho_{s}/\hbar^{2}=n_{\mathrm{sp}}/2$ and the compressibility
$A=(1+x/\sqrt{1+x^{2}})m/8\pi\hbar^{2}$, while
$D=m|\Delta|/\sqrt{1+x^{2}}/2\pi\hbar^{2}$ controls the phase-modulus
coupling, and $C=|\Delta|^{2}(1+x/\sqrt{1+x^{2}})m/2\pi\hbar^{2}$ describes
the mass term for the modulus fluctuations. Neglecting modulus fluctuations
$(\lambda=0)$ leads to a sound velocity
$c=\sqrt{\rho_{s}/A}/\hbar=(\mu^{2}+|\Delta|^{2})^{1/4}/\sqrt{m}$, shown as
the dotted green line in Fig. 3(b), which diverges in the Bose limit. However,
when $\lambda\neq 0$, the coupling between modulus and phase renormalizes $A$.
Integrating out $\lambda$ leads to a renormalized phase-only action with
unchanged superfluid density $\rho_{sR}=\rho_{s}$, but renormalized
compressibility $A_{R}=A+D^{2}/4C=m/4\pi\hbar^{2}$. This renormalization leads
to the corrected speed $c=\sqrt{\rho_{sR}/A_{R}}/\hbar$ given by Eq. (4), as
expected.

Change in concavity: To investigate the low-energy collective modes beyond
linear dispersion, it is necessary to expand $M_{\pm\pm}$ and $M_{+-}$ up to
sixth order in $\omega$ and $|{\bf q}|$. In this case, the condition $\det
M({\bf q},\omega_{\bf q})=0$ leads to

$\omega_{\bf q}=c|{\bf q}|\left[1+\frac{\gamma}{8}\left(\frac{\hbar|{\bf
q}|}{mc}\right)^{2}\\!\\!+\frac{\eta}{16}\left(\frac{\hbar|{\bf
q}|}{mc}\right)^{4}\\!\\!+\mathcal{O}\\!\left(\frac{\hbar|{\bf
q}|}{mc}\right)^{6}\right]\\!,$ (6)

where the coefficients of the cubic and quintic order corrections have
analytic expressions

$\displaystyle\gamma=$
$\displaystyle\,\frac{1}{24}\left(1-4x^{4}-x\frac{7+4x^{2}}{\sqrt{1+x^{2}}}\right),$
(7) $\displaystyle\eta=$
$\displaystyle\,-\frac{365+2802x^{2}+2048x^{4}+160x^{6}}{23040(1+x^{2})}$
$\displaystyle-x\frac{685+1813x^{2}+1064x^{4}+80x^{6}}{11520(1+x^{2})^{3/2}}.$
(8)

In Fig. 2, we show $\omega_{\bf q}$ for various interaction regimes and
different levels of approximation. The solid red curves represent the full
numerical solutions, while the other curves represent the linear (dotted
blue), cubic (dashed green), and quintic (dot-dashed magenta) approximations
in $|{\bf q}|$ with the Gaussian corrected $\mu$. The double-dotted-dashed
line represents the numerical phase-only $(\lambda=0)$ dispersion $\omega_{\bf
q}^{PO}$ using the saddle point value of $\mu$. One of the important messages
from the panels in Fig. 2 is that there is a change in curvature in the red
solid lines, also found in 3D Fermi gases Engelbrecht  _et al._ (1997);
Combescot  _et al._ (2006); Kurkjian  _et al._ (2016); Biss  _et al._ (2021),
where the dispersion $\omega_{\bf q}$ is supersonic $(\gamma>0)$ in the Bose
regime shown in panel (a), and subsonic in the BCS regime $(\gamma<0)$ shown
in panel (d), where it bends downwards due to the pair-breaking continuum. A
second important observation is that the dispersion $\omega_{\bf q}^{PO}$
severely overestimates the correct $\omega_{\bf q}$ in the Bose regime, that
is, $\omega_{\bf q}^{PO}\gg\omega_{\bf q}$, while in the BCS limit
$\omega_{\bf q}^{PO}\approx\omega_{\bf q}$, because the modulus and phase
fluctuations are nearly decoupled.

The coefficients $\gamma$ and $\eta$ are presented in Fig. 4 as function of
$\ln k_{F}a$. In all panels, $\gamma$ and $\eta$ are evaluated at different
levels of approximation: the dotted green lines describe phase-only
fluctuations using the saddle point $\mu$, the dashed blue lines include
modulus and phase fluctuations using the saddle point $\mu$, and the solid red
lines include modulus and phase fluctuations using the Gaussian $\mu$. Panels
(a)-(b) ((c)-(d)) show $\gamma$ and $\eta$ in units that elucidate their
limiting behavior in the Bose (BCS)regime given by the dot-dashed black
(magenta) lines.

Figure 4:  Parameters $\gamma$ and $\eta$ vs. $\ln k_{F}a$ at different levels
of approximation: Phase-only fluctuations with saddle point $\mu$ (dotted
green lines), phase and modulus fluctuations with saddle point $\mu$ (dashed
blue lines) and with Gaussian corrected $\mu$ (solid red lines). The dot-
dashed black and magenta lines describe the Bose and BCS limits. In (c) and
(d), $\gamma$ and $\eta$ are scaled to reveal their behavior in the BCS limit.

The concavity of $\omega_{\bf q}$ is controlled by $\gamma$, which changes
from $\gamma>0$ (convex) to $\gamma<0$ (concave) in the Bose and BCS regimes,
respectively. The parameter $\gamma$ changes sign at
$x=\sqrt{(2\sqrt{13}-7)/12}\approx 0.133$, corresponding to $\ln k_{F}a\simeq
0.65$ using the Gaussian $\mu$. In this case $(\gamma=0)$, the first
correction to the linear spectrum is a quintic $(|{\bf q}|^{5})$ term
controlled by $\eta$. In contrast to 3D results Kurkjian  _et al._ (2016),
$\eta$ is always negative in 2D for the entire crossover.

In the Bose limit $(x\ll-1)$, expanding the matrix elements of $M({\bf
q},\omega)$ to order $(\Delta/|\mu|)^{2}$ and to lowest order in $\hbar|{\bf
q}|/\sqrt{2m|\mu|}$ leads to

$\hbar^{2}\left(\omega_{\bf q}^{\rm
B}\right)^{2}=\frac{\hbar^{2}q^{2}}{2m_{B}}\left(\frac{\hbar^{2}q^{2}}{2m_{B}}+2m_{B}c^{2}_{B}\right),$
(9)

where $c_{B}=|\Delta|/\sqrt{2m_{B}|\mu|}$ is the Bogoliubov speed of sound,
and $m_{B}=2m$ is the boson mass. Using the saddle point $\mu$ leads to the
incorrect value $c_{B}=v_{F}/\sqrt{2}$, while using the Gaussian corrected
$\mu$ leads to $c_{B}=v_{F}/\sqrt{8|\ln k_{F}a|}$ in the Bose limit. The
Bogoliubov-Popov interaction energy $2n_{B}V_{B}(0)=2m_{B}c^{2}=E_{F}/|\ln
k_{F}a|$, with boson density $n_{B}=n/2$ leads to the boson-boson interaction
parameter $V_{B}(0)=(E_{F}/n)/|\ln k_{F}a|$. The values of $\gamma$ and $\eta$
from Eq. (9) are equal to limiting results obtained from Eqs. (7–8), that is,
$\gamma\rightarrow 1/4$, and $\eta\rightarrow-1/128$, as seen in Figs.
4(a)-(b).

In the BCS limit $(x\gg 1)$, the saddle point and Gaussian correction tend to
the same results, as fluctuations are less important. Rescaling energies by
$|\Delta|$, such that $\hbar\omega_{\bf q}/|\Delta|$ tends to a universal
function of $\hbar c|{\bf q}|/|\Delta|$, leads to
$(|\Delta|/mc^{2})^{2}\gamma\rightarrow-1/3$ and
$(|\Delta|/mc^{2})^{4}\eta\rightarrow-1/72$, as revealed in Figs. 4(c)-(d).
This is a consequence of the two-particle continuum pushing down the
collective mode branch Combescot  _et al._ (2006); Kurkjian  _et al._ (2016).
In this case, the expansion in $|{\bf q}|$ is limited to $\hbar c|{\bf q}|\leq
2|\Delta|$ and $c\to v_{F}/\sqrt{2}$.

Conclusions: We analyzed low-energy collective modes of 2D Fermi superfluids
from the BCS to the Bose regime giving excellent results when compared to
Bragg spectroscopy experiments in 2D box potentials. We showed that quantum
fluctuations in the phase and modulus of the superfluid order parameter are
absolutely necessary to give physically acceptable chemical potential and
speed of sound. Lastly, we showed that the dispersion of the collective mode
changes from concave to convex as interactions are changed from the BCS to the
Bose regime, a feature also seen experimentally.

###### Acknowledgements.

Acknowledgments: We thank the Belgian American Educational Foundation for
support.

## References

  * Goldstone (1961) J. Goldstone, Field theories with superconductor solutions, Nuovo Cimento 19, 154 (1961).
  * Higgs (1966) P. W. Higgs, Spontaneous symmetry breakdown without massless bosons, Phys. Rev. 145, 1156 (1966).
  * Englert  _et al._ (1966) F. Englert, R. Brout, and M. F. Thiry, Vector mesons in presence of broken symmetry, Nuovo Cimento A43, 244 (1966).
  * Hong  _et al._ (2017) T. Hong, M. Matsumoto, Q. Y., W. Chen, T. R. Gentile, S. Watson, F. F. Awwadi, M. M. Turnbull, D. S. E., H. Agrawal, R. Toft-Petersen, B. Klemke, K. Coester, K. P. Schmidt, and D. A. Tennant, Higgs amplitude mode in a two-dimensional quantum antiferromagnet near the quantum critical point, Nature Physics 13, 638 (2017).
  * Maiti and Hirschfeld (2015) S. Maiti and P. J. Hirschfeld, Collective modes in superconductors with competing s- and d-wave interactions, Phys. Rev. B 92, 094506 (2015).
  * Greco  _et al._ (2003) V. Greco, M. Colonna, M. Di Toro, and F. Matera, Collective modes of asymmetric nuclear matter in quantum hadrodynamics, Phys. Rev. C 67, 015203 (2003).
  * Yamamoto  _et al._ (2007) N. Yamamoto, M. Tachibana, T. Hatsuda, and G. Baym, Phase structure, collective modes, and the axial anomaly in dense QCD, Phys. Rev. D 76, 074001 (2007).
  * Kobyakov and Pethick (2013) D. Kobyakov and C. J. Pethick, Dynamics of the inner crust of neutron stars: Hydrodynamics, elasticity, and collective modes, Phys. Rev. C 87, 055803 (2013).
  * Mathur and Turton (2014) S. D. Mathur and D. Turton, Comments on black holes I: The possibility of complementarity, J. High Energ. Phys. 2014, 34.
  * Singh and Mathey (2021) V. P. Singh and L. Mathey, Collective modes and superfluidity of a two-dimensional ultracold Bose gas, Phys. Rev. Research 3, 023112 (2021).
  * Choudhury  _et al._ (2020) S. Choudhury, K. R. Islam, Y. Hou, J. A. Aman, T. C. Killian, and K. R. A. Hazzard, Collective modes of ultracold fermionic alkaline-earth-metal gases with SU(N) symmetry, Phys. Rev. A 101, 053612 (2020).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010).
  * Ibarra-García-Padilla  _et al._ (2020) E. Ibarra-García-Padilla, R. Mukherjee, R. G. Hulet, K. R. A. Hazzard, T. Paiva, and R. T. Scalettar, Thermodynamics and magnetism in the two-dimensional to three-dimensional crossover of the Hubbard model, Phys. Rev. A 102, 033340 (2020).
  * Mermin and Wagner (1966) N. D. Mermin and H. Wagner, Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic Heisenberg models, Phys. Rev. Lett. 17, 1133 (1966).
  * Hohenberg (1967) P. C. Hohenberg, Existence of long-range order in one and two dimensions, Phys. Rev. 158, 383 (1967).
  * Berezinskii (1972) V. L. Berezinskii, Destruction of long-range order in one-dimensional and two-dimensional systems possessing a continuous symmetry group. II. Quantum systems, Sov. Phys. JETP 34, 610 (1972).
  * Kosterlitz and Thouless (1972) J. M. Kosterlitz and D. J. Thouless, Long range order and metastability in two dimensional solids and superfluids.(Application of dislocation theory), Journal of Physics C: Solid State Physics 5, L124 (1972).
  * Keimer  _et al._ (2015) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, From quantum matter to high-temperature superconductivity in copper oxides, Nature 518, 179 (2015).
  * Ge  _et al._ (2015) J.-F. Ge, Z.-L. Liu, C. Liu, C.-L. Gao, D. Qian, Q.-K. Xue, Y. Liu, and J.-F. Jia, Superconductivity above 100 K in single-layer FeSe films on doped SrTiO3, Nature Materials 14, 285 (2015).
  * Cao  _et al._ (2018) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018).
  * Yu  _et al._ (2019) Y. Yu, L. Ma, P. Cai, R. Zhong, C. Ye, J. Shen, G. D. Gu, X. H. Chen, and Y. Zhang, High-temperature superconductivity in monolayer Bi2Sr2CaCu2O8+δ, Nature 575, 156 (2019).
  * Makhalov  _et al._ (2014) V. Makhalov, K. Martiyanov, and A. Turlapov, Ground-state pressure of quasi-2D Fermi and Bose gases, Phys. Rev. Lett. 112, 045301 (2014).
  * Fenech  _et al._ (2016) K. Fenech, P. Dyke, T. Peppler, M. G. Lingham, S. Hoinka, H. Hu, and C. J. Vale, Thermodynamics of an attractive 2D Fermi gas, Phys. Rev. Lett. 116, 045302 (2016).
  * Boettcher  _et al._ (2016) I. Boettcher, L. Bayha, D. Kedar, P. A. Murthy, M. Neidig, M. G. Ries, A. N. Wenz, G. Zürn, S. Jochim, and T. Enss, Equation of state of ultracold fermions in the 2D BEC-BCS crossover region, Phys. Rev. Lett. 116, 045303 (2016).
  * Hadzibabic  _et al._ (2006) Z. Hadzibabic, P. Krüger, M. Cheneau, B. Battelier, and J. Dalibard, Berezinskii–Kosterlitz–Thouless crossover in a trapped atomic gas, Nature 441, 1118 (2006).
  * Murthy  _et al._ (2015) P. A. Murthy, I. Boettcher, L. Bayha, M. Holzmann, D. Kedar, M. Neidig, M. G. Ries, A. N. Wenz, G. Zürn, and S. Jochim, Observation of the Berezinskii-Kosterlitz-Thouless phase transition in an ultracold Fermi gas, Phys. Rev. Lett. 115, 010401 (2015).
  * Ville  _et al._ (2018) J. L. Ville, R. Saint-Jalm, É. Le Cerf, M. Aidelsburger, S. Nascimbène, J. Dalibard, and J. Beugnon, Sound propagation in a uniform superfluid two-dimensional Bose gas, Phys. Rev. Lett. 121, 145301 (2018).
  * Bohlen  _et al._ (2020) M. Bohlen, L. Sobirey, N. Luick, H. Biss, T. Enss, T. Lompe, and H. Moritz, Sound propagation and quantum-limited damping in a two-dimensional Fermi gas, Phys. Rev. Lett. 124, 240403 (2020).
  * Christodoulou  _et al._ (2021) P. Christodoulou, M. Gałka, N. Dogra, R. Lopes, J. Schmitt, and Z. Hadzibabic, Observation of first and second sound in a BKT superfluid, Nature 594, 191 (2021).
  * Sobirey  _et al._ (2021) L. Sobirey, H. Biss, N. Luick, M. Bohlen, H. Moritz, and T. Lompe, Comparing fermionic superfluids in two and three dimensions, arXiv:2106.11893 (2021).
  * Martiyanov  _et al._ (2010) K. Martiyanov, V. Makhalov, and A. Turlapov, Observation of a two-dimensional Fermi gas of atoms, Phys. Rev. Lett. 105, 030404 (2010).
  * Feld  _et al._ (2011) M. Feld, B. Fröhlich, E. Vogt, M. Koschorreck, and M. Köhl, Observation of a pairing pseudogap in a two-dimensional Fermi gas, Nature 480, 75 (2011).
  * Dyke  _et al._ (2011) P. Dyke, E. D. Kuhnle, S. Whitlock, H. Hu, M. Mark, S. Hoinka, M. Lingham, P. Hannaford, and C. J. Vale, Crossover from 2D to 3D in a weakly interacting Fermi gas, Phys. Rev. Lett. 106, 105304 (2011).
  * Ong  _et al._ (2015) W. Ong, C. Cheng, I. Arakelyan, and J. E. Thomas, Spin-imbalanced quasi-two-dimensional Fermi gases, Phys. Rev. Lett. 114, 110403 (2015).
  * Hueck  _et al._ (2018) K. Hueck, N. Luick, L. Sobirey, J. Siegl, T. Lompe, and H. Moritz, Two-dimensional homogeneous Fermi gases, Phys. Rev. Lett. 120, 060402 (2018).
  * Hoinka  _et al._ (2017) S. Hoinka, P. Dyke, M. G. Lingham, J. J. Kinnunen, G. M. Bruun, and C. J. Vale, Goldstone mode and pair-breaking excitations in atomic Fermi superfluids, Nature Physics 13, 943 (2017).
  * Biss  _et al._ (2021) H. Biss, L. Sobirey, N. Luick, M. Bohlen, J. J. Kinnunen, G. M. Bruun, T. Lompe, and H. Moritz, Excitation spectrum and superfluid gap of an ultracold Fermi gas, arXiv:2105.09820 (2021).
  * Randeria  _et al._ (1990) M. Randeria, J.-M. Duan, and L.-Y. Shieh, Superconductivity in a two-dimensional Fermi gas: Evolution from Cooper pairing to Bose condensation, Phys. Rev. B 41, 327 (1990).
  * Marini  _et al._ (1998) M. Marini, F. Pistolesi, and G. Strinati, Evolution from BCS superconductivity to Bose condensation: analytic results for the crossover in three dimensions, The European Physical Journal B - Condensed Matter and Complex Systems 1, 151 (1998).
  * Salasnich  _et al._ (2013) L. Salasnich, P. A. Marchetti, and F. Toigo, Superfluidity, sound velocity, and quasicondensation in the two-dimensional BCS-BEC crossover, Phys. Rev. A 88, 053612 (2013).
  * Lumbeeck  _et al._ (2020) L.-P. Lumbeeck, J. Tempere, and S. Klimin, Dispersion and damping of phononic excitations in Fermi superfluid gases in 2D, Condensed Matter 5, 13 (2020).
  * He  _et al._ (2015) L. He, H. Lü, G. Cao, H. Hu, and X.-J. Liu, Quantum fluctuations in the BCS-BEC crossover of two-dimensional Fermi gases, Phys. Rev. A 92, 023620 (2015).
  * Engelbrecht  _et al._ (1997) J. R. Engelbrecht, M. Randeria, and C. A. R. Sá de Melo, BCS to Bose crossover: Broken-symmetry state, Phys. Rev. B 55, 15153 (1997).
  * Botelho and Sá de Melo (2006) S. S. Botelho and C. A. R. Sá de Melo, Vortex-antivortex lattice in ultracold fermionic gases, Phys. Rev. Lett. 96, 040404 (2006).
  * 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).
  * Nozières and Schmitt-Rink (1985) P. Nozières and S. Schmitt-Rink, Bose condensation in an attractive fermion gas: From weak to strong coupling superconductivity, J. Low Temp. Phys. 59, 195 (1985).
  * Sá de Melo  _et al._ (1993) C. A. R. Sá de Melo, M. Randeria, and J. R. Engelbrecht, Crossover from BCS to Bose superconductivity: Transition temperature and time-dependent Ginzburg-Landau theory, Phys. Rev. Lett. 71, 3202 (1993).
  * Note (1) Note that one has to introduce convergence factors to regularize the Matsubara sum in the thermodynamic potential $\Omega_{\mathrm{g}}$. This is most easily done in the basis of the real and imaginary part of the pair field $\Phi$, instead of the phase-modulus basis used here Engelbrecht  _et al._ (1997); He  _et al._ (2015).
  * Combescot  _et al._ (2006) R. Combescot, M. Y. Kagan, and S. Stringari, Collective mode of homogeneous superfluid Fermi gases in the BEC-BCS crossover, Physical Review A 74, 042717 (2006).
  * Landau and Lifshitz (1977) L. D. Landau and E. M. Lifshitz,  _Quantum Mechanics - $3^{\rm rd}edition$,_ (Pergamon Press, 1977) p. 63.
  * Kurkjian  _et al._ (2016) H. Kurkjian, Y. Castin, and A. Sinatra, Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover, Phys. Rev. A 93, 013623 (2016).


&=-\frac{1}{2}\bigg{(}1-\frac{3052L^{2}+9792(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{4}}\lambda^{2}\eta^{2}\bigg{)}\sum_{k=0}^{K-2}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]\\\
&+\frac{34\lambda^{2}}{N(1-\lambda^{2})^{3}}\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}+\frac{34\lambda^{2}\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}|\mathcal{B}|}\bigg{(}1+\frac{10\beta}{1-\lambda^{2}}\bigg{)}+\frac{136\lambda^{2}\beta^{2}K\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}}\bigg{(}1+\frac{5\beta}{1-\lambda^{2}}\bigg{)}.\end{split}$
(75)

With the condition that
$\eta\leq\frac{(1-\lambda^{2})^{2}}{\lambda\sqrt{12844L^{2}+9792G^{2}C_{\upsilon}}}$,
it is immediately that

$1-\frac{3052L^{2}+9792(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{4}}\lambda^{2}\eta^{2}\geq
0$

such that the first term on the right hand side of Eq. 75 can be removed.
Moreover, since
$\eta\leq\frac{\sqrt{N(1-\lambda^{2})}\lambda}{31\sqrt{L^{2}+G^{2}C_{\upsilon}}}$,
which leads to $\beta\leq\frac{(1-\lambda^{2})\lambda^{2}}{10}$ such that

$1+\frac{10\beta}{1-\lambda^{2}}<2,\;1+\frac{5\beta}{1-\lambda^{2}}<\frac{3}{2}.$

By now we can conclude that

$\begin{split}-\frac{1}{2}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{17}{16\eta^{2}N}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]&\leq\frac{34\lambda^{2}}{N(1-\lambda^{2})^{3}}\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}+\frac{68\lambda^{2}\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}|\mathcal{B}|}\\\
&+\frac{204\lambda^{2}\beta^{2}K\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}}.\end{split}$
(76)

Substituting Eq. 76 into Eq. 74 and dividing both sides by $\frac{1}{K+1}$,
and with $\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|^{2}]=\frac{1}{N(K+1)}\sum_{i=1}^{N}\sum_{k=0}^{K}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k})\|^{2}]$ completes the proof. ∎

### Proof of Corollary 1

###### Proof.

Substituting $\eta,\beta,\textnormal{and}|\mathcal{B}|$ into the Eq. 9 and
conducting some simple mathematical manipulations can easily attain the
desirable result. ∎

### Analysis for MDPGT with Single Trajectory Initialization

######  Theorem 2.

Let Assumptions 1,2 and 3 hold. Let the momentum coefficient
$\beta=\frac{96L^{2}+96G^{2}C_{\upsilon}}{N}\eta^{2}$. If MDPGT is initialized
by a single trajectory and the step size satisfies the following condition

$\begin{split}0<\eta\leq\textnormal{min}&\bigg{\\{}\frac{(1-\lambda^{2})^{2}}{\lambda\sqrt{12844L^{2}+9792G^{2}C_{\upsilon}}},\frac{\sqrt{N(1-\lambda^{2})}\lambda}{31\sqrt{L^{2}+G^{2}C_{\upsilon}}},\\\
&\frac{1}{6\sqrt{6(L^{2}+G^{2}C_{\upsilon})}}\bigg{\\}},\end{split}$ (77)

then the output $\tilde{\mathbf{x}}_{K}$ satisfies: for all $K\geq 2$:

$\begin{split}&\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|^{2}]\leq\frac{4(J^{*}-J(\bar{\mathbf{x}}_{0}))}{\eta
K}+\frac{4\bar{\sigma}^{2}}{N\beta K}+\frac{8\beta\bar{\sigma}^{2}}{N}\\\
&+\frac{34\lambda^{2}}{KN(1-\lambda^{2})^{3}}\|\nabla\mathbf{J}(\bar{\mathbf{x}}_{0})\|^{2}+\frac{68\lambda^{2}\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}K}\\\
&+\frac{204\lambda^{2}\beta^{2}\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}},\end{split}$
(78)

where $J^{*}$ is the upper bound of $J(\mathbf{x})$ and
$\|\nabla\mathbf{J}(\bar{\mathbf{x}}_{0})\|^{2}\triangleq\sum^{N}_{i=1}\|\nabla
J_{i}(\bar{\mathbf{x}}_{0})\|^{2}$.

###### Proof.

As the initialization is single trajectory, $|\mathcal{B}|=1$. Following the
similar proof techniques as already shown in Theorem 1 yields the desirable
result. The proof process is not repeated in this context. ∎

Similarly, in view of Theorem 2, the asymptotic behavior is the same as in
Theorem 1, having the same steady-state error. Due to an infinite time
horizon, regardless of the initialization strategy, diverse agents are able to
learn effectively in a collaborative manner. However, in terms of the non-
asymptotic property, the single trajectory initialization strategy makes a
difference in the sampling complexity, which is reflected by the following
result.

######  Corollary 3.

Let $\eta=\frac{N^{3/4}}{8LK^{1/4}},\beta=\frac{DN^{1/2}}{64L^{2}K^{1/2}}$, in
Theorem 2. We have

$\begin{split}&\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|^{2}]\leq\frac{32L(J^{*}-J(\bar{\mathbf{x}}_{0}))}{(NK)^{3/4}}+\frac{2048L^{4}\bar{\sigma}^{2}+D^{2}\bar{\sigma}^{2}N}{8L^{2}DN^{3/2}K^{1/2}}\\\
&+\frac{\lambda^{2}}{(1-\lambda^{2})^{3}K}\bigg{(}\frac{34\|\nabla\mathbf{J}(\bar{\mathbf{x}}_{0})\|^{2}}{N}+68\bar{\sigma}^{2}+\frac{51\bar{\sigma}^{2}D^{2}N}{1024L^{4}}\bigg{)},\end{split}$
(79)

for all

$\begin{split}&K\geq\textnormal{max}\bigg{\\{}\frac{N^{3}D^{2}}{4096L^{4}},\frac{923521N(L^{2}+G^{2}C_{\upsilon})^{2}}{4096L^{4}(1-\lambda^{2})^{2}\lambda^{4}},\\\
&\frac{(12844L^{2}+9792G^{2}C_{\upsilon})^{2}\lambda^{4}N^{3}}{4096L^{4}(1-\lambda^{2})^{8}}\bigg{\\}},\end{split}$
(80)

where $D=96L^{2}+96G^{2}C_{\upsilon}$.

###### Proof.

Substituting $\eta,\beta$ into the Eq. 78 and conducting some simple
mathematical manipulations can easily attain the desirable result. ∎

Corollary 3 implies that with only a single trajectory initialization for
MDPGT, when $K$ is sufficiently large, the mean-squared convergence rate is

$\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|^{2}]\leq\mathcal{O}\bigg{(}\frac{1}{(NK)^{1/2}}\bigg{)},$

which is slower than that obtained in Corollary 1. With similar mathematical
manipulation, the eventual sampling complexity is
$\mathcal{O}(N^{-1}\epsilon^{-4})$. Though variance reduction techniques has
not reduced the order of $\epsilon^{-1}$, compared to the SOTA approaches, the
linear speed up still enables the complexity to be $N$ times smaller than that
in (Xu, Gao, and Gu 2019; Huang et al. 2020). Additionally, different from
traditional decentralized learning problems, MARL has more significant
variances in the optimization procedure due to the non-oblivious
characteristic. Using just a single trajectory for each agent to initialize is
a quite poor scheme, but the adopted variance reduction techniques can
successfully maintain the SOTA sampling complexity in a decentralized setting.

### Implication for Gaussian Policy

In this section, we study the sample complexity when the policy function
$\pi^{i}(a^{i}|s)$ of each agent is explicitly a Gaussian distribution. For a
bounded action space $\mathcal{A}^{i}\subset\mathbb{R}$, a Gaussian policy
parameterized by $\mathbf{x}_{i}$ is defined as

$\pi^{i}(a^{i}|s)=\frac{1}{\sqrt{2\pi}}\textnormal{exp}\bigg{(}-\frac{((\mathbf{x}^{i})^{\top}\phi_{i}(s)-a^{i})^{2}}{2\xi^{2}}\bigg{)},$
(81)

where $\xi^{2}$ is a constant standard deviation parameter and
$\phi_{i}(s):\mathcal{S}\to\mathbb{R}^{d_{i}}$ is mapping from the state space
to the feature space. Note that the standard deviation parameter can be
varying in terms of different agents, while for simplicity in this context, we
assume that it is the same for each agent. Next, we verify the assumptions on
the Gaussian policy. We first impose a mild assumption that the action space
and feature space are bounded, i.e., there exist constants $C_{a}>0$ and
$C_{f}>0$ such that for all $i\in\mathcal{V}$, $|a^{i}|\leq C_{a},\forall
a^{i}\in\mathcal{A}^{i}$ and $\|\phi_{i}(s)\|\leq C_{f},\forall
s\in\mathcal{S}$. Then we can show that Assumption 1 hold. Due to the limit of
space, we defer the derivation to the supplementary materials. While
Assumption 3 does not always hold for all Gaussian distributions, but based on
a result from (Cortes, Mansour, and Mohri 2010), it has been shown that for
two Gaussian distributions,
$\pi^{i}_{\mathbf{x}_{1}}(a^{i}|s)=\mathcal{N}(\mu_{1},\xi_{1}^{2})$ and
$\pi^{i}_{\mathbf{x}_{2}}(a^{i}|s)=\mathcal{N}(\mu_{2},\xi_{2}^{2})$, if
$\xi_{2}>\frac{\sqrt{2}}{2}\xi_{1}$, then the variance of the importance
sampling weight $\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})$ is
bounded. Due to a constant standard deviation parameter defined for the
Gaussian distributions, it is easily verified that for any time step $k\geq
0$,
$\mathbb{V}(\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}_{k}^{i}))$
should be bounded by some constant $\mathcal{M}>0$. We recall Corollary 1 that
can apply to any general policy models. Therefore, based on the above
discussion, it also applies to the Gaussian policy function scenario. We
present the result in the following corollary to explicitly state the
relationship between $\gamma$ and $\epsilon$ in the sampling complexity. We
first present an established result that bounds the variance of policy
gradient estimator $\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})$ in Eq. 6 when the
policy function $\pi^{i}$ follows a Gaussian distribution.

######  Lemma 11.

(Lemma 5.5 in (Pirotta, Restelli, and Bascetta 2013)) Given a Gaussian policy
$\pi^{i}_{\mathbf{x}^{i}}=\mathcal{N}((\mathbf{x}^{i})^{\top}\phi_{i}(s),\xi^{2})$,
for all $i\in\mathcal{V}$, if $|r_{i}(s,a^{i})|\leq R$ and
$\|\phi^{i}(s)\|\leq C_{f}$ for all $s\in\mathcal{S},a^{i}\in\mathcal{A}^{i}$
and $R>0,C_{f}>0$ are constants, then the variance of policy gradient
estimator $\mathbb{V}(\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i}))$ can be bounded
as

$\mathbb{V}(\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i}))\leq\frac{R^{2}C_{f}^{2}}{(1-\gamma)^{2}\xi^{2}}\bigg{(}\frac{1-\gamma^{2H}}{1-\gamma^{2}}-H\gamma^{2H}-2\gamma^{H}\frac{1-\gamma^{H}}{1-\gamma}\bigg{)}.$
(82)

In this context, we have defined a constant standard deviation parameter for
the policy of each agent to simplify the analysis. However, one can still
define separately this parameter for each agent, i.e., $\xi_{i}$. Then in the
upper bound of last inequality,
$\xi^{2}=\textnormal{min}\\{\xi^{2}_{1},\xi^{2}_{2},...,\xi^{2}_{N}\\}$. We
show the proof for Corollary 2 in the following.

###### Proof.

According to Corollary 1, when $K$ is sufficiently large, $T_{1}$ dominates
the convergence such that we investigate it. According to Lemma 11, we can
immediately obtain that

$\bar{\sigma}^{2}\leq\frac{R^{2}C_{f}^{2}}{(1-\gamma)^{2}\xi^{2}}\bigg{(}\frac{1-\gamma^{2H}}{1-\gamma^{2}}-H\gamma^{2H}-2\gamma^{H}\frac{1-\gamma^{H}}{1-\gamma}\bigg{)}=\mathcal{O}\bigg{(}\frac{1}{(1-\gamma)^{3}}\bigg{)}.$

As $L=\frac{C_{h}R}{(1-\gamma)^{2}}$ and $G=\frac{C_{g}R}{(1-\gamma)^{2}}$, we
have
$D=96L^{2}+96G^{2}C_{\upsilon}=\mathcal{O}\bigg{(}\frac{1}{(1-\gamma)^{4}}\bigg{)}$.
Thus,

$\frac{256L^{3}D(J^{*}-J(\bar{\mathbf{x}}_{0}))+2048L^{4}\bar{\sigma}^{2}+D^{2}\bar{\sigma}^{2}}{8L^{2}D(NK)^{2/3}}\leq\mathcal{O}\bigg{(}\frac{1}{(1-\gamma)^{3}(NK)^{2/3}}\bigg{)}.$

We can easily obtain that the sampling complexity is
$\mathcal{O}\bigg{(}\frac{1}{(1-\gamma)^{4.5}N\epsilon^{3}}\bigg{)}$ when the
policy is parameterized by a Gaussian distribution. ∎

######  Remark 2.

When the policy function is a Gaussian distribution for each agent, Corollary
2 shows that the sampling complexity is inherently in the order of
$\mathcal{O}(\epsilon^{-3})$, which matches the generalized result. Moreover,
the linear speed up is still retained. Additionally, we also show that the
complexity does not rely on the horizon. The dependence on $(1-\gamma)^{-4.5}$
stems from the variance of the policy gradient estimator, which has been known
in the centralized counterpart (Xu, Gao, and Gu 2019). When the initialization
is with a single trajectory, it can be inferred from Corollary 3 that the
sampling complexity is $\mathcal{O}((1-\gamma)^{-6}N^{-1}\epsilon^{-4})$.

### Algorithmic Framework for MDPG

We next present a new decentralized algorithm adapted from the centralized
MBPG (Huang et al. 2020), which serves as a new baseline for empirical
comparison. Different from MBPG, the updates are simplified with constant step
size and $\beta$ value. This simplification is to enable a fair comparison
among all algorithms.

Result: $\tilde{\mathbf{x}}_{K}$ chosen uniformly random from
$\\{\mathbf{x}^{i}_{k},i\in\mathcal{V}\\}^{K}_{k=1}$

Input:
$\mathbf{x}^{i}_{1}=\bar{\mathbf{x}}_{1}\in\mathbb{R}^{d},\eta\in\mathbb{R}^{+},\beta\in(0,1),\mathbf{W}\in\mathbb{R}^{N\times
N},K,\mathcal{B}\in\mathbb{Z}^{+},k=1$;

while  _$k <K$_ do

for  _each agent_ do

if $k=1$: compute the local policy gradient surrogate by sampling a trajectory
$\tau^{i}_{1}$ from
$p_{i}(\tau^{i}|\mathbf{x}^{i}_{1}):\mathbf{u}^{i}_{1}=\mathbf{g}_{i}(\tau^{i}_{1}|\mathbf{x}^{i}_{1})$,
or by sampling a mini-batch of trajectories
$\\{\tau_{1}^{i,m}\\}_{m=1}^{|\mathcal{B}|}$ from
$p_{i}(\tau^{i}|\mathbf{x}^{i}_{1}):\mathbf{u}^{i}_{1}=\frac{1}{|\mathcal{B}|}\sum_{m=1}^{|\mathcal{B}|}\mathbf{g}_{i}(\tau^{i,m}_{1}|\mathbf{x}^{i}_{1})$;

if $k>1$: sample a trajectory $\tau^{i}_{k}$ from
$p_{i}(\tau^{i}|\mathbf{x}^{i}_{k})$ and compute the local policy gradient
surrogate using Eq. 8;

Update the local estimate of the policy network parameters
$\mathbf{x}^{i}_{k+1}=\sum_{j\in
Nb(i)}\omega_{ij}(\mathbf{x}^{j}_{k}+\eta\mathbf{u}^{j}_{k})$;

end for

$k=k+1$;

end while

Algorithm 2 Momentum-based Decentralized Policy Gradient (MDPG)

### Computing Resource Details

All numerical experiments presented in the main article and supplementary
section were performed on a computing cluster with 2 x 16 Core Intel Xeon CPU
E5-2620 v4 and a memory of 256 GB. All codes were implemented using PyTorch
version 1.8.1.

### Additional Experimental Results

In the following section, we include additional results for comparing MDPGT,
MDPG and DPG in a gridworld environment with two and three agents.
Additionally, we also provide results for five, ten, twenty and thirty agents
in a simplified lineworld environment.

#### Gridworld Environment

Figure 4: Experimental results comparing MDPGT, MDPG and DPG in the gridworld
environment with 2 agents using fully-connected topology. Figure 5:
Experimental results comparing MDPGT, MDPG and DPG in the gridworld
environment with 3 agents using fully-connected topology. Figure 6:
Experimental results comparing MDPGT, MDPG and DPG in the gridworld
environment with 10 agents using fully-connected topology.

Figures 4, 5, and 6 illustrate the average reward obtained by the three
algorithms, MDPGT, MDPG and DPG in with two, three and ten agents respectively
using a fully-connected network topology with $\beta$ = 0.5. It should be
noted that in Figure 6 the number of trajectories is only 25,000 such that
MDPGT and MDPG perform similarly. However, when the number is 50,000, it
suffices to show that the result should resemble the plot in Figure 1, with
even a larger gap between MDPGT and MDPG. In all scenarios, we observed that
both MDPGT and MDPG still outperforms DPG. Interestingly, with fewer agents,
the performance of MDPGT and MDPG are similar, in contrast to the five agent
scenario presented in the main article. Fewer agents may reduce the impact of
linear speed up on the error bound, which is w.r.t $\mathcal{O}(N^{-1})$. By
observing carefully, when the number of agents becomes larger, MDPGT is more
advantageous in a complex graph than MDPG as the extra tracking step
contributes to correcting the gradient bias caused by diverse agents. This
interesting finding will be validated in the following simplified lineworld
environment where more agents are incorporated.

#### Lineworld Environment

As an initial proof of concept and to demonstrate our algorithm’s ability to
scale to a larger number of agents, we created a simplification of the
gridworld environment called the Lineworld environment. In this environment,
all agents are randomly initialized with a 1-D coordinate. The goal of all the
agents are to cooperatively arrive at the coordinate 0. In this setting, the
agents are not allowed to collide and have a smaller action set of either
moving up, down or stay stationary. The reward function we used for this
environment is also defined as the Euclidean distance of individual agent to
its respective goal. All agent’s policy is represented by a 3-layer neural
network with 64 hidden units with $tanh$ activation functions. The agents were
trained for 10,000 episodes with a horizon of 500 steps and discount factor of
0.99. A learning rate of 3E-4 was used in all our experiments and we average
the results conducted over 5 random seeds.

Figures 7 and 8 shows the rewards of our proposed methods and the baseline DPG
in the Lineworld environment with five, ten, twenty and thirty agents
respectively. In summary, the findings from the performance in the Lineworld
environment strengthens the conclusion attained from the Gridworld setting
that MDPGT outperforms MDPG in a more sophisticated graphs involving more
agents. With five agents, we observed that MDPGT performs better than DPG but
worse than MDPG. However, as the number of agent increases, we see that the
performance of MDPGT improves to be on-par or slightly better than MDPG and
significantly better than DPG. Thus, a practical guideline from an application
point-of-view is that one can choose MDPG for simple graphs and expect MDPGT
to perform much better in complex graphs.

Figure 7: Experimental results comparing MDPGT, MDPG and DPG in the lineworld
environment with 5 agents (left) and 10 agents (right) using fully-connected
topology.

Figure 8: Experimental results comparing MDPGT, MDPG and DPG in the lineworld
environment with 20 agents (left) and 30 agents (right) using fully-connected
topology. Figure 9: Experimental results comparing MDPGT in the gridworld
environment with 5 agents and $\beta=0.5$ using fully-connected topology with
different mini-batch initializations.

#### Minibatch Initialization

To validate the theoretical findings obtained in this work for the difference
between two different versions of MDPGT, we implement several experiments with
5 agents using fully-connected topology with different mini-batch
initializations. In this context, mini-batch size equal to 1 corresponds to
MDPGT with only single trajectory initialization. From Figure 9, one can
observe that the reward curves due to different mini-batch initializations are
depicted in three phases. Before approximately 20,000 trajectories, regardless
of whichever initialization, MDPGT yields similar performance. After 20,000
trajectories, the difference becomes appealing, which is reflected by the
outperforming of MPDGT with mini-batch size equal to 4. This supports the
claim of Corollary 1 that $K$ has to satisfy a certain condition to start
showing the non-asymptotic behavior. Additionally, when the number of
trajectories is larger than 60,000, another interesting phenomenon is that all
reward curves start converging together again, which can be explained by using
Theorem 1 as MDPGT enables the asymptotic results once $K$ is sufficiently
larger than a certain number. We have pointed out in the analysis that they
all converge to the same steady-state error defined in Eq. 10. Intuitively,
practitioners can set a sufficiently large number of trajectories to get the
same optimal solution regardless of what initialization is taken. While
different approximate solutions can also be obtained when different specific
mini-batch sizes for initialization are defined. One may argue how to set the
mini-batch size for a specific problem. For example, we can observe from
Figure 9 that to get better approximate solutions, the mini-batch size could
be selected in between 4 and 16. This is also suggested from Corollary 1 that
the mini-batch size is w.r.t both the number of agents and the number of
trajectories. Thus, either mini-batch size equal to 64 or 256 results in
slightly worse result than only one trajectory initialization.