Superfluids, Fluctuations and Disorder

TL;DR

This paper develops a field-theory framework to study how disorder affects ultracold bosonic atoms, revealing that disorder can destroy superfluidity even when perturbative analysis applies.

Contribution

It introduces a novel field-theory approach incorporating disorder effects on superfluid and condensate properties in ultracold bosonic systems.

Findings

Disorder modifies condensate depletion and superfluid response.

Superfluid phase can be destroyed by disorder in perturbative regimes.

Implementation of replica trick for disordered ultracold bosons.

Abstract

We present a field-theory description of ultracold bosonic atoms in presence of a disordered external potential. By means of functional integration techniques, we aim to investigate and review the interplay between disordered energy landscapes and fluctuations, both thermal and quantum ones. Within the broken-symmetry phase, up to the Gaussian level of approximation, the disorder contribution crucially modifies both the condensate depletion and the superfluid response. Remarkably, it is found that the ordered (i.e. superfluid) phase can be destroyed also in regimes where the random external potential is suitable for a perturbative analysis. We analyze the simplest case of quenched disorder and then we move to present the implementation of the replica trick for ultracold bosonic systems. In both cases, we discuss strengths and limitations of the reviewed approach, paying specific attention to possible extensions and the most recent experimental outputs.

Equations

$$\hat{H}=\sum_{i=1}^{N}\bigg{[}-\frac{\hbar^{2}\nabla_{i}^{2}}{2m}+\hat{U}_{\text{D}}(\mathbf{r}_{i})\bigg{]}+\frac{1}{2}\sum_{i\neq j}\hat{V}\big{(}|\mathbf{r}_{i}-\mathbf{r}_{j}|\big{)}\;,$$

$$Z_{N}=\sum_{\{n\}}\exp\big{(}-\beta E_{n}\big{)}$$

$$\mathcal{Z}=\int \mathcal{D}[\psi ,{\psi }^{\ast }]e^{-S[\psi ,{\psi }^{\ast }]/\hslash }.$$

$$S[\psi ,{\psi }^{\ast }]$$

$$\displaystyle\qquad\qquad\qquad\qquad\qquad+\frac{1}{2}\int d^{d}\mathbf{r}^{\prime}|\psi(\mathbf{r},\tau)|^{2}V(|\mathbf{r}-\mathbf{r}^{\prime}|)|\psi(\mathbf{r}^{\prime},\tau)|^{2}\bigg{\}}\;.$$

$$\braket{\,\ldots\,}_{\text{dis}}\equiv L^{-{Nd}}\prod_{i=1}\int d^{d}\mathbf{r}_{i}\big{(}\ldots\big{)}\;,$$

$${⟨\dots ⟩}_{\text{dis}}\equiv \int \mathcal{D}[U_{\text{D}}]\mathsf{P}[U_{\text{D}}](\dots ).$$

$$\mathsf{P}[U_{\text{D}}]=\exp\bigg{\{}-\frac{1}{2\gamma^{2}}\int d^{d+1}\mathbf{s}\int d^{d+1}\mathbf{s}^{\prime}\;U_{\text{D}}(\mathbf{s})K^{-1}(\mathbf{s}-\mathbf{s}^{\prime})U_{\text{D}}(\mathbf{s}^{\prime})\bigg{\}}$$

$$⟨U_D(\mathbf{s}){⟩}_{\text{dis}}=0,$$

$${⟨U_{\text{D}}(\mathbf{s})U_{\text{D}}({\mathbf{s}}^{\prime })⟩}_{\text{dis}}={\gamma }^{2}K(\mathbf{s}-{\mathbf{s}}^{\prime }).$$

$$\mathsf{P}[U_{\text{D}}]=\exp\bigg{\{}-\frac{1}{2\gamma^{2}}\int_{0}^{\beta\hbar}d\tau\int d^{d}\mathbf{r}\;U^{2}_{\text{D}}(\mathbf{r},\tau)\bigg{\}}\;.$$

$${⟨\mathcal{O}[\psi ,{\psi }^{\ast }]⟩}_{\text{dis+th}}=\int \mathcal{D}[U_{\text{D}}]\mathsf{P}[U_{\text{D}}]{⟨\mathcal{O}[\psi ,{\psi }^{\ast }]⟩}_{\text{th},U_{\text{D}}}$$

$${⟨\mathcal{O}[\psi ,{\psi }^{\ast }]⟩}_{\text{th},U_{\text{D}}}=\int \mathcal{D}[\psi ,{\psi }^{\ast }]\frac{1}{\mathcal{Z}(U_{\text{D}})}\mathcal{O}[\psi ,{\psi }^{\ast }]e^{-S[\psi ,{\psi }^{\ast }]/\hslash }$$

$$\mathscr{Z}[J,J^{*}]=\int\mathscr{D}[\psi,\psi^{*}]\;.\exp\bigg{\{}-\hbar^{-1}S[\psi,\psi^{*}]+\hbar^{-1}\int_{0}^{\beta\hbar}d\tau\int d^{d}\mathbf{r}\big{[}J^{*}(\mathbf{r},\tau)\psi(\mathbf{r},\tau)+J(\mathbf{r},\tau)\psi^{*}(\mathbf{r},\tau)\big{]}\bigg{\}}\;.$$

$$\braket{\mathscr{O}}_{\text{th}}=-\frac{\delta}{\delta J}\log\mathscr{Z}\bigg{|}_{J=0}\;.$$

$$\braket{\mathscr{O}}_{\text{dis+th}}=-\frac{\delta}{\delta J}\braket{\log\mathscr{Z}}_{\text{dis}}\bigg{|}_{J=0}\\ =-\int\mathscr{D}[U_{\text{D}}]\;\frac{\mathsf{P}(U_{\text{D}})}{\mathscr{Z}[J=0,U_{\text{D}}]}\frac{\delta}{\delta J}\mathscr{Z}[J,U_{\text{D}}]\bigg{|}_{J=0}\;.$$

$$\bigg{[}\hbar\frac{\partial}{\partial\tau}-\frac{\hbar^{2}\nabla^{2}}{2m}-(\mu+U_{\text{D}}(\mathbf{r},\tau))+\int d^{d}\mathbf{r}^{\prime}\,|\psi(\mathbf{r}^{\prime},\tau)|^{2}V(|\mathbf{r}-\mathbf{r}^{\prime}|)\bigg{]}\psi(\mathbf{r},\tau)=0\;.$$

$$V=g{\delta }^{(d)}(\mathbf{r})$$

$$\psi (\mathbf{r},\tau )={\psi }_0$$

$$\big{[}-(\mu+\braket{U_{\text{D}}}_{\text{dis}})+g\psi_{0}^{2}\big{]}\psi_{0}=0\Longrightarrow\begin{cases}\psi_{0}^{2}=0\qquad\text{for $\mu+\braket{U_{\text{D}}}_{\text{dis}}\leq 0$}\\ \\ \psi_{0}^{2}=(\mu+\braket{U_{\text{D}}}_{\text{dis}})/g\qquad\text{for $\mu+\braket{U_{\text{D}}}_{\text{dis}}>0$}\;.\end{cases}$$

$$\psi (\mathbf{r},\tau )={\psi }_0+\eta (\mathbf{r},\tau ),$$

$$\mathscr{L}_{\text{dis}}[\psi,\psi^{*},U_{\text{D}}]=-U_{\text{D}}(\mathbf{r})\big{|}\psi(\mathbf{r},\tau)\big{|}^{2}=-U_{\text{D}}(\mathbf{r})\big{[}\psi_{0}^{2}+|\eta(\mathbf{r},\tau)|^{2}+v\eta(\mathbf{r},\tau)+v\eta^{*}(\mathbf{r},\tau)\big{]}$$

$${\rm Y}={\psi }_0\left(\begin{array}{c}1\\ 1\end{array}\right)\text{and}\chi =\left(\begin{array}{c}\eta \\ {\eta }^{\ast }\end{array}\right)$$

$${\mathcal{L}}_{\text{dis}}[\eta ,{\eta }^{\ast },U_{\text{D}}]=-U_{\text{D}}(\mathbf{r}){{\rm Y}}^T\chi .$$

$$\mathcal{Z}(U_{\text{D}})$$

$$\displaystyle=e^{\beta L^{d}\mu^{2}/(2g)}\int\mathscr{D}[\chi]\exp\bigg{\{}-\frac{1}{\hbar}\int d^{d+1}\mathbf{s}\bigg{[}\frac{1}{2}\chi^{\dagger}\mathscr{G}^{-1}(\partial_{\tau},\nabla)\chi-U_{\text{D}}(\mathbf{r})\Upsilon^{T}\chi\bigg{]}\bigg{\}}\;,$$

$$\mathscr{Z}(U_{\text{D}})=e^{-\beta\Omega_{\text{mf}}}\;\big{(}\det\mathscr{G}^{-1}\big{)}^{-1}\exp\bigg{\{}\frac{1}{2}\int d^{d+1}\mathbf{s}\int d^{d+1}\mathbf{s}^{\prime}U_{\text{D}}(\mathbf{r})\Upsilon^{T}\mathscr{G}(\mathbf{s}-\mathbf{s}^{\prime})\Upsilon U_{\text{D}}(\mathbf{r}^{\prime})\bigg{\}}\;.$$

$${\mathcal{G}}^{-1}(\mathbf{q},\omega )=\left(\begin{array}{cc}-i\hslash {\omega }_n+{ϵ}_q-\mu +2g{\psi }_0^2& g{\psi }_0^2\\ & \\ g{\psi }_0^2& +i\hslash {\omega }_n+{ϵ}_q-\mu +2g{\psi }_0^2\end{array}\right)$$

$$\big{(}\det\mathscr{G}^{-1}\big{)}^{-1}=\exp\big{(}-\text{tr}\log\mathscr{G}^{-1}\big{)}=\exp\big{(}-\log\det\mathscr{G}^{-1}\big{)}\;.$$

$$\mathscr{Z}(U_{\text{D}})=e^{-\beta\Omega_{\text{mf}}}\exp\bigg{\{}\underbrace{-\frac{1}{2}\sum_{\mathbf{q}}\sum_{n}\log\det\mathscr{G}^{-1}(\mathbf{q},\omega_{n})\bigg{\}}}_{\text{Pure Gaussian term}}\exp\bigg{\{}\underbrace{\frac{1}{2}\sum_{\mathbf{q}}\big{|}U_{\text{D}}(\mathbf{q})\big{|}^{2}\Upsilon^{T}\mathscr{G}(\mathbf{q},0)\Upsilon}_{\text{disorder Gaussian term }}\bigg{\}}\;.$$

$$\frac{1}{2\beta}\sum_{n}\log\det\mathscr{G}^{-1}(\mathbf{q},\omega_{n})=\frac{1}{2\beta}\sum_{n}\log\big{[}\beta^{2}\big{(}\hbar^{2}\omega_{n}^{2}+E_{q}^{2}\big{)}\big{]}=\frac{1}{2}\,E_{q}+\frac{1}{\beta}\log\big{(}1-e^{-\beta E_{q}}\big{)}\;,$$

Full text

cold atoms, superfluidity, disorder

I Introduction

Ultracold atomic gases have been the cornerstone of atomic physics since 1995 cornell-1995 ; ketterle-1995 , when Bose-Einstein condensation was experimentally achieved for the first time. This groundbreaking result was the starting point of an ongoing and intense research effort to improve atoms trapping and cooling procedures pethick-book . Nowadays, modern laboratories have reached an exquisite control on the relevant physical parameters for atomic gases, such as the number density or the strength and range of atom-atom interaction. At the same time, it is also possible to tune the coupling with the external environement, enabling the observation of quantum dynamics without spurious smearing effects langen-2015 .

It is then natural to think of ultracold atomic gases as an extremely effective platform to probe quantum theories up to the macroscopic scale. One of the most striking example is given by the Bose-Hubbard model and its superfluid-insulator transition originally predicted in fisher-1989 . This quantum phase transition was first observed in 2002 by loading a repulsive condensate in an optical lattice confining potential greiner-2002 . Similar achievements have been reached for low-dimensional quantum gases and their remarkable phenomenology giamarchi-book ; dalibard-review-2d . Among a wide and interesting literature, it is worth mentioning the first direct observation of the Tonks-Girardeau regime in an optical lattice filled with alkali atoms paredes-2004 ; kinoshita-2004 or the vortex proliferation associated to the Berezinskii-Kosterlitz-Thouless physics in two-dimensional setups hadzibabic-2006 ; cornell-2007 . More recently, a growing interest has aroused on cold atoms as a promising quantum simulation platform. Together with the refined experimental expertise, it has been shown that, within peculiar setups, the superfluid dynamics of atoms obeys to equations which can be mapped to other ones appearing in high-energy physics or cosmology. For instance, Rabi-coupled bosonic mixtures seem to be a burgeoning candidate to explore inflationary dynamics and space-time expansion in table top experiments marquardt-2015 ; brandt-2017 ; braden-2018 . Other proposals involve cold-atoms analog of quantum gravity liberati-2006 or the implementation of system mimicking QCD features sademelo-2018 .

Restraining ourselves to the condensed matter domain, ultracold atoms offer the interesting possibility to investigate the intriguing issue of quantum transport in a disordered environment kamenev-review . Indeed, from the dawn of solid state theory it was clear that a well-founded investigation has to consider the multifaceted role played by microscopic impurities and disordered external potential. For instance, common incandescent lamp would not exist without a certain degree of impurity scattering. At the same time, it is reasonable to expect a smearing of experimental data due to randomly wrinkled substrate or the destruction of supertransport properties. Examples of this parasitic role can be traced back to dirty superconductors or 4He diffusion in porous medium ma-1985 ; ma-1986 . The complexity of this problem lies in the interplay between superfluidity (i.e. the occuring of some kind of ordered phase) and the localization phenomenon first theorized by Anderson in anderson-1959 .

Given the broad interest and the cross-disciplinary relevance of this problem, it is reasonable to turn towards cold atoms, since they offer reliable experimental protocols to engineer disordered external potential and control the corresponding relevant parameters. The dirty bosons problem then gathers all the investigations addressing the properties of ultracold bosons subject to an external random potential. The interplay between superfluidity and localization may lead to the observation of a new phase, the so-called Bose glass fisher-1989 . Similarly to the Bose-Hubbard model in absence of disorder, optical lattices have been proved again to be a reliable setup to explore these exotic features damski-2003 . While interfering counterpropagating beams create the first regular pattern, randomness is added by superimposing a second optical lattice, providing a viable strategy to control the degree of disorder schulte-2005 ; billy-2008 ; roati-2008 . An alternative is given by the speckle optical potential, generated by the interference of waves with different phases and amplitudes but with the same frequency dainty-1980 ; goodman-book . The resulting wave varies randomly in space and can be used to study the transport properties of ultracold atoms in presence of static external disorder lye-2005 ; clement-2005 ; ghabour-2014 .

In order to tackle down the crucial aspects of the dirty bosons problem, the relevant theoretical tools can be divided in two different subsets. At first, one can adopt the point of view of Bogoliubov theory, treating disorder in a perturbative way, similarly to quantum and thermal fluctuations in the pure system stoof-book ; salasnich-2016 . Since the first seminal work on this topic huang-1992 ; giorgini-1994 ; tauber-1997 , this approach has led us to a better understanding of the role played by disorder, no matter how small, in condensed bosonic systems. The additional depletion to condensate and superfluid was computed in huang-1992 , while in giorgini-1994 the authors investigated the modifications induced by a random external potential to the sound propagation. At the same time, critical properties are investigated by focusing, for instance, on the condensation temperature shift driven by the interplay of disorder and interaction lopatin-2002 ; falco-2007 .

A non-perturbative approach paves the way to the exploration of the superfluid-glass (quantum) transition in bosonic ensembles. Moving from the seminal renormalization-group analysis for one-dimensional systems giamarchi-1988 ; giamarchi-book , different strategies and methods have been used to understand the phase diagram of dirty bosons in three and lower dimensions, such as the random-phase approximation navez-2007 , stochastic field-theory techniques and the replica formalism yukalov-2007 ; falco-2009 ; khellil-2016 ; khellil-2016-2 .

The ambition of this review is obviously more limited, since it is not possible face in the main text all the topics we have mentioned during this introduction. Nevertheless, we aim to present, in a clear and pedagogical way, the crucial ingredients of a field-theory analysis of superfluid bosons. Throughtout the proceeding of the paper, we are going to adopt the functional integration framework, which is flexible enough to enable the implementation of different approximation schemes altland-book . Our discussion is organized in the following way: in Sec. II we specify how one can build up a quantum field theory for superfluid bosons subject to a random confinement. In particular, one needs to clarify the way in which disorder is introduced and characterized. In Sec. III we compute, up to the Gaussian level of approximation, the additional contribution to the condensate and the superfluid component due to the presence of a quenched disordered potential. Moving to Sec. IV, we present the core ideas at the foundation of the replica formalism, which allows to generalize our analysis to more complex disorder realization and does not rely, in principle, upon the assumption of a perturbatively small disorder. Within this scheme, we show how to compute the meaningful correlation functions and describe the superfluid response of the system in Sec. V. Further comments and future perspective are addressed in the Conclusions.

II Disorder in Field Theories

Let us begin by considering $N$ interacting bosonic particles of mass $m$ enclosed in a $L^{d}$ volume and subject to an external disordered potential $U_{\text{D}}(\mathbf{r})$. Thus, the Hamiltonian is simply given by

[TABLE]

with \hat{V}\big{(}|\mathbf{r}_{i}-\mathbf{r}_{j}|\big{)} the two-body central interaction potential. The thermodynamic description of the system can be derived from the partition function. Within the canonical ensemble, it is defined as

[TABLE]

where $\beta\equiv(k_{B}T)^{-1}$ and $\{n\}$ denoting the bosonic states (i.e. totally symmetric with respect to particle exchange) with eigenenergy $E_{n}$ of the operator $\hat{H}$ in Eq. (1). In order to adopt the functional integration formalism altland-book , Eq. (1) has to be expressed in second quantization. By adding $-\mu\hat{N}$, with $\mu$ the chemical potential and $\hat{N}$ the number operator, we can also move to the grand canonical ensemble. Here, a basis of coherent bosonic states can be chosen, leading to the imaginary-time representation of the (grand canonical) partition function

[TABLE]

In the equation above, the Bose statistics is encoded in the periodic boundary condition for the fields, namely $\psi(0)=\psi(\beta\hbar)$ and $\psi^{*}(0)=\psi^{*}(\beta\hbar)$. The Eucliden action appearing in Eq. (3) is given by

[TABLE]

We have now to clarify what is meant by a (quantum) field theory for the disordered Bose gas. Also in the simplest case of non-interacting bosons, given a certain realization of the disordered confinement, an exact diagonalization of Eq. (1) does not appear as a feasible task.

A statistical approach is certainly more reasonable, since it allows, in principle, to characterize the system through a limited range of variables as, for instance, the strength of the impurity potential leading to the disorder, or the scale of space-time potential fluctuations. This obviously implies that it is possible to average (in a way to define) over different microscopic realization of the disorder potential.

This outline may seem very vague: how can we implement a theoretical description in statistical terms? Apparently, it is simple to imagine a set of $N$ impurities at $\{\mathbf{r}_{i}\}_{i=1,\ldots,N}$ positions. When they are assumed equivalent, the resulting disorder (i.e. impurity) potential takes the obvious form $U_{\text{D}}(\mathbf{r})=\sum_{i}U(\mathbf{r}-\mathbf{r}_{i})$. Within this recipe the disorder average is equivalent to the integration over the whole configuration space of impurities, namely

[TABLE]

a technical task not easy to implement in a functional formalism.

As pointed out throughtout a vast literature (hertz-1985 ; nelson-1990 ; tauber-1997 ; schakel-1997 ; lopatin-2002 ; falco-2007 ; lubensky-1975 ; grinstein-1976 , just to mention a few), a more fruitful approach relies upon the statistical features of the disorder potential or, in other words, a description in terms of probability density function (PDF) $\mathsf{P}(U_{\text{D}})$. It follows then that

[TABLE]

For instance, we can surely imagine a Gaussian distributed disorder, i.e.

[TABLE]

where $\mathbf{s}=(\mathbf{r},\tau)$ and $d^{d+1}{\bf s}=d\tau\,d^{d}{\bf r}$. The expectation value of a central Gaussian distributed variable is zero, i.e.

[TABLE]

while its second (central) momentum corresponds to

[TABLE]

For a $\delta$-correlated disorder , where $K({\bf s}-{\bf s}^{\prime})=\delta^{(d+1)}({\bf s}-{\bf s}^{\prime})$ with $\delta^{(d+1)}({\bf s})$ the Dirac delta function in $d+1$ dimensions, the above equation is further simplified in

[TABLE]

If we deal with atomic bosons in the superfluid phase the atomic wavelength is much larger than the scattering range of impurities generating the disorder potential $U_{\text{D}}(\mathbf{s})$. This physical situation clarifies the simplifying assumption of a $\delta$-correlated disorder as in $\eqref{pdf per disordine delta correlato}$. Moreover, in this section we explicitly consider the case of a quenched disorder, where the time scale of $U_{\text{D}}$ is much longer than the thermodynamic one. From a technical point of view this implies that $U_{\text{D}}$ retains only a spatial dependence on $\mathbf{r}$ and that the disorder average in Eq. (6) has to be performed after the thermal one falco-2007 . Thus, the average over a quenched disorder has to be intended as

[TABLE]

where

[TABLE]

is the usual thermal average over the grand canonical ensemble for a given impurity configuration. Within the scheme outlined by Eqs. (11) and (12), $U_{\text{D}}(\mathbf{r})$ and $\psi(\mathbf{r},\tau)$ are not being treated on the same footing, since the random field $U_{\text{D}}$ is fixed when we average over the bosonic field $\psi$. More precisely, Eq. (12) is telling us to compute the partition function over a precise microscopic realization of the system (i.e. with a given $U_{\text{D}}$). Then, according to Eq. (11), we have to average all the partition functions obtained in this way over the disorder field.

The quenched average presented above is certainly a viable scheme to adopt when the disorder can be treated in a perturbative way. However, this is not always the case. It can occur a physical situation where one has to perform the disorder average early in the calculation and/or the random field $U_{\text{D}}$ has to be considered at the same level of $\psi$. In order to highlight the eventual technical obstacle, we make use of the source method to compute expectation values of physical observable. According to this field-theory technique, the expectation value of a given observable $\mathscr{O}$ can be extracted by means of the source method altland-book ; schakel-book . Here, the source is intended as a term linear in the fields and proportional to $J(\mathbf{r},\tau)$, namely

[TABLE]

A simple differentiation leads us to the pursued expectation value. The most immediate example is the chemical potential, indeed $\braket{N}=\beta^{-1}\partial_{\mu}\log\mathscr{Z}$. The thermal average of $\mathscr{O}$ is defined as

[TABLE]

As a consequence of Eq. (6), the joint thermal and disorder average reads

[TABLE]

Since one has both to differentiate with respect to $J$ and integrate over $U_{\text{D}}$, the disorder field appears both at the numerator and denominator. Consequently, Eq. (14) can be challenging.

In the following sections we approach both the situations: first we consider the case of a quenched disorder potential, then we move to a more general framework by presenting the implementation of the replica trick for superfluid atomic bosons.

III Perturbative approach to quenched disorder

III.1 Thermodynamic picture and disorder-driven condensate depletion

Let us consider the Euclidean action in Eq. (4) and apply the saddle-point method by stationarizing it, i.e.

[TABLE]

In the case of a two-body contact interaction, we have

[TABLE]

where $g=\tilde{V}(\mathbf{q}=0)$ and $\tilde{V}(\mathbf{q})$ is the Fourier transform of the two-body potential. For a uniform configuration of the field, with

[TABLE]

not depending on $\mathbf{r}$ and $\tau$, we get

[TABLE]

This is approximation must not be underrated: actually, in presence of an external potential $\psi_{0}$ should not be uniform. However, at low temperatures it seems reasonable to assume that the wavelengths of atoms are much larger than the spatial variations caused by the impurity potential.

By comparing the equations above with the ones for a pure system salasnich-2016 , it is immediate to realize that, in principle, the disorder may shift the critical point of the superfluid transition. However, if $U_{\text{D}}$ follows a centered Gaussian distribution like the one in Eq. (7), this shift equates to zero. From the saddle-point result it descends that the total density $n=\psi^{2}_{0}$. By identifying $\psi_{0}^{2}$ with the condensate density $n_{0}$, this signals that, up to this (very rough) level of approximation, all the particles take part in the condensate within the broken-symmetry phase.

In this phase, where a $\mathsf{U}(1)$ symmetry is spontaneously broken, fluctuations can be taken into account by considering the following splitting of the field

[TABLE]

with $\eta(\mathbf{r},\tau)\in\mathbb{C}$. Replacing the equation above in the Euclidean action of Eq. (4), the disorder contribution reads

[TABLE]

where $\mathscr{L}$ denotes the (Euclidean) Lagrangian density, i.e. $S[\psi]=\int d\tau\int d^{d}\mathbf{r}\;\mathscr{L}[\psi]$. The first two terms give rise to an unobservable shift of the chemical potential, so we will neglect them. Moreover, let us remark that we are considering fluctuations over the saddle-point non-trivial configuration $\psi_{0}^{2}=\mu/g$ with $\braket{U_{\text{D}}}_{\text{dis}}=0$, but this is not the true stationary point. Consequently, we are not allowed to neglect the linear terms in Eq. (20), not even up to the Gaussian order in the fluctuation fields. By introducing the column vectors

[TABLE]

Eq. (20) is brought to its more compact form

[TABLE]

At this point, we can replace the field $\psi(\mathbf{r},\tau)$ with Eq. (19). By retaining terms up to the quadratic order in $\eta$ and $\eta^{*}$, the grand canonical partition function reads

[TABLE]

with $\mathbf{s}=(\mathbf{r},\tau)$. In the equation above, the integral in $\chi$ is Gaussian and, consequently, can be computed exactly, leading us to the following result

[TABLE]

Let us also recall that, from Eq. (18) it is easy to realize that the mean-field grand potential is simply $\Omega_{\text{mf}}/L^{d}=-\mu^{2}/2g$. The remaining steps of the calculation are easier in the Fourier space, where the (inverse) Gaussian propagator is the following $2\times 2$ matrix

[TABLE]

with $\omega_{n}=2\pi n/\beta$ the bosonic Matsubara frequencies and $\epsilon_{q}=\hbar q^{2}/(2m)$. The (functional) determinant in Eq. (24) can be handled more easily thanks to the identities

[TABLE]

By putting all these pieces together, Eq. (24) is built by three blocks, namely

[TABLE]

The disorder contribution in the equation above does not depend on the frequency, since the time scales of the quenched disordered are frozen, implying that the corresponding frequencies goes to zero. Concerning the pure Gaussian contribution, we have salasnich-2016

[TABLE]

where a crucial role is played by the excitation spectrum of collective excitations above the uniform saddle-point configuration, i.e.

[TABLE]

Focusing on the zero-temperature case, we find

[TABLE]

with

[TABLE]

the zero-temperature Gaussian grand potential. At this point, by recalling the relation $\Omega=\beta^{-1}\log\mathscr{Z}$, we are ready compute the quenched disorder average average of the grand potential. It is immediate to realize the $\Omega_{\text{mf}}$ and $\Omega_{g}^{(0)}$ are transparent to this operation. On the contrary, from Eq. (30) we easily derive that

[TABLE]

whose quenched average reads

[TABLE]

thanks to Eq. (10) (with only the spatial integration). Thus, we have now isolated the disorder contribution to the zero-temperature grand potential, namely

[TABLE]

Moving to the continuum limit, i.e. $\sum_{\mathbf{q}}\rightarrow\frac{L^{d}}{(2\pi)^{d}}\int d^{d}\mathbf{q}$, and performing the resulting integral with the aid of dimensional regularization, we finally obtain, in agreement with schakel-1997 ; schakel-book , that

[TABLE]

with $\Gamma(x)$ being the Euler Gamma function. From the equation above we can compute the corresponding disorder contribution to the total number density $n$ as a function of the condensate one $n_{0}$. Indeed, from $\braket{\Omega_{\text{dis}}}_{\text{dis}}$ in Eq. (35) we have

[TABLE]

resulting in

[TABLE]

Concerning physical values of $d$, for $d=3$ we have $g=4\pi\hbar^{2}a_{s}/m$ with $a_{s}$ the s-wave scattering length. Consequently, one can verify that

[TABLE]

The two-dimensional case is much more complicated, since the zero-range coupling constant displays a peculiar logarithmic dependence on the number density boronat-2009 ; salasnich-2017 ; tononi-2018 ; tononi-2019 . For two-spatial dimensions a separate investigation is then needed, in order to investigate the eventual divergences in the free energy and the relation between the disorder potential and the prediction of the Mermin-Wagner-Hohenberg theorem lewenstein-2006 ; boudjema-2013 ; boudjema-2015 . At $T=0$ the total number density is then given by

[TABLE]

with $n_{g}^{(0)}$ being the pure Gaussian contribution, derived from Eq. (31) as stated by Eq. (39). According to tononi-2018 we have

[TABLE]

It is crucial to remark that Eqs. (37) and (39) imply that, for a given number of total particles $n$, a smaller fraction of them resides in the condensate. Thus, we can conclude that a disorder potential will cause an additional depletion of the condensate. In order to extract the condensate fraction $n_{0}/n$, Eq. (39) has to be numerically solved. Thus, for $d=3$, by replacing Eq. (40) and Eq. (38) in Eq. (39) we obtain

[TABLE]

The equation above represents an implicit expression for the condensate fraction. At a given number density $n$, we can numerically solve it in $n_{0}$ to compute $n_{0}/n$. Concerning Eq. (41) we have to underline another important detail. It seems that the disorder contribution to $n$ (cfr. Eq. (38)) diverges in the very weakly interacting limit $a_{s}\rightarrow 0$. However, the presence of a random external potential implies the presence of another characteristic length scale. Following falco-2007 , we define it in terms of the disorder correlator $\gamma$ given by Eq. (10), i.e.

[TABLE]

Now, we recall that our perturbative Gaussian scheme requires that $n_{g}^{(0)}$ in Eq. (40) and $n_{\gamma}$ from Eq. (37) are both small compared to $n_{0}\simeq n$. Thus, by imposing $n_{g}^{(0)}\ll n_{0}\simeq n$ we get the condition $n^{1/3}a_{s}\ll 1$, the usual condition on the interaction strength. At the same time, we have to require that $n_{\gamma}\ll n_{0}\simeq n$, leading us $(n^{1/3}a_{s})^{1/2}\,2(n^{1/3}l_{\text{dis}})/\sqrt{\pi}\ll 1$. The latter highlights the fact that we have to require $n^{1/3}l_{\text{dis}}\gg 1$ in addition to the usual condition on the gas parameter. Provided that both the constrains hold, Eq. (41) can be perturbatively inverted, reading

[TABLE]

It is worth remembering at the end of this section that, in the case of quenched disorder, the additional contribution to the grand potential is effectively a zero-temperature one (or, in other words, it does not depend on $\beta$). The reason is that, as we have stressed in the previous section, we are considering an external random confinement whose time scale is much larger than the thermodynamic one.

III.2 Superfluid response and quenched disorder

Now, what can we say about the disorder influence on the superfluid motion of the system? Actually, within the Landau-Khalatnikov two-fluid formulation landau-book ; khalatnikov-book , the calculation proceeds in a similar fashion, provided that the theory is modified accordingly to the following points schakel-book ; tononi-2018 ; tononi-2019 :

• •

within the Landau two-fluid model, we assume that the normal part of the system is in motion with velocity $\mathbf{v}$, so $\partial_{\tau}\rightarrow\partial_{\tau}-i\mathbf{v}\cdot\nabla$;

• •

Superfluid and normal parts do not exchange momentum. This means that a superflow can be imposed by means of phase twist of the field $\psi(\mathbf{r},\tau)$, i.e. \psi\rightarrow\exp\big{(}im\mathbf{v}_{s}\cdot\mathbf{r}/\hbar\big{)}\psi with $\mathbf{v}_{s}$ the superfluid velocity fisher-1973 ; taylor-2006 ;

• •

as a consequence of the previous points, we can redefine the chemical potential as

[TABLE]

The resulting Gaussian partition function (or grand potential) is formally similar to the one in Eq. (27), with the replacement $\mu\rightarrow\mu_{\text{eff}}$ and the (inverse) Gaussian propagator given by

[TABLE]

whose determinant is given by

[TABLE]

where $E_{q}$ as in Eq. (29) with $\mu\rightarrow\mu_{\text{eff}}$.

By following this prescription, the phase-twisted grand potential can be derived from Eq. (33), reading

[TABLE]

We are interested in the linear-response regime, where the density current

[TABLE]

is taken up to the first order in $(\mathbf{v}-\mathbf{v}_{s})$. Therefore, we begin with an expansion of Eq. (47) up to the quadratic order in $\hbar\mathbf{q}\cdot(\mathbf{v}-\mathbf{v}_{s})$. For symmetry reasons the linear term of this expansion is identically zero, such that we have to deal with the following equation,

[TABLE]

where we have made use of the identity

[TABLE]

Since the first term in the equation above depends only on $\mu_{\text{eff}}$ and it holds $\frac{\partial\mu_{\text{eff}}}{\partial\mathbf{v}}=m\mathbf{v}_{s}$, it will be the disorder contribution to the $n\mathbf{v}_{s}$ term in $\mathbf{g}$. More precisely, by considering also the thermal depletion due to the quasiparticle excitations, Eq. (48) reads

[TABLE]

At this point, we finally replace Eq. (18) for the broken-symmetry case in the equation above. A simple comparison between the resulting equation and the the two-fluid assumption $\mathbf{g}=n_{s}\mathbf{v}_{s}+n_{n}\mathbf{v}$ leads us to

[TABLE]

with $n_{\gamma}(n_{0})$ as in Eq. (37).

In order to understand the relevance of Eq. (52), one has to remember that, in the pure system (i.e. $\gamma=0$), all the particles belong to the superflow. On the contrary, a disorder potential induce a superfluid depletion also at $T=0$ and, in addition, the parameter $\gamma$ may act as knob to destroy the superflow since

[TABLE]

In a similar way to the condensate fraction, we are not allowed to expand $n_{s}/n$ in the limit $g\rightarrow 0$ because Eq. (37) is divergent in that regime. Obviously, the superfluid density can be computed from the two-fluid assumption $n=n_{s}+n_{n}$ where $n_{n}$ is given by Eq. (52).

In fig. 1 we report the behaviour of the superfluid fraction in presence of a disordered but uncorrelated external potential with $\gamma^{2}=0.001$ for two different values of the gas parameter $na_{s}^{3}$. In order to highlight the crucial contribution due to disorder, we also consider the case where $\gamma^{2}=0$. It is immediate to realize that in pure system ($\gamma^{2}=0$) the whole system is superfluid at $T=0$, while the presence of a disorder potential lowers the curves, since there is an additional depletion at $T=0$. It is also interesting to underline that the more dilute is the system, the more evident is the detachment from the pure system.

IV The replicated action for superfluid bosons

The previous section has been focused on the extension of usual field-theory techniques superfluid bosons subject to a random external confinement. Moreover, the whole calculation has been carried out in the quenched regime, where we effectively average over different microscopic realization of the system with a given disorder configuration.

We have already made clear that this is not the only viable strategy and that it does not handle physical realizations where, for instance, the disorder potential also has a thermal component. This means that, looking at Eq. (14), it is not clear in which order integration and differentiation have to be performed. In order to overcome this ambiguity, the replica trick represents one of the most ingenious strategy. It was first developed in the context of spin glasses anderson-1975 ; parisi-book ; parisi-review but, in the following, we are going to detail how it can be applied to describe the role of disorder in ultracold bosonic gases.

Let us begin by considering again Eq. (14). The following three hypothesis represent the starting of the replica formalism. Thus,

• •

We assume that $\mathscr{Z}[J=0]=1$, establishing a proper normalization and implying that Eq. (14) has no dependence on $V_{\text{D}}$ in the denominator;

• •

it holds $\braket{\mathscr{O}}_{\text{th+dis}}=-\delta\braket{\log\mathscr{Z}}_{\text{dis}}/\delta_{J}|_{J=0}$;

• •

for $J\neq 0$ the partition function retains its algebraic properties.

Actually, the last point is a somewhat obscure way to state that the average procedure involves functionals linear in the disorder potential, since $V_{\text{D}}$ has been linearly added to the action. The great advantage provided by the replica trick can be understood from the formal relations reported below,

[TABLE]

with $R\in\mathbb{N}$. The replacement of $\log\mathscr{Z}$ with an algebraic function of $\mathscr{Z}$ strongly reduces the complexity. Indeed, concerning the calculation of the expectation values, we have

[TABLE]

Actually, within the replicated formalism, expectation values can be obtained by simply considering $R$ identical copies of the system under consideration. Concerning the functional integration, it does not alter the linear dependence of the action on the disorder field, because $[\exp(\int dV_{\text{D}})]^{R}=\exp(\sum^{R}\int dV_{\text{D}})$. Unfortunately all this comes with a price. Looking at Eq. (55) we realize that, at the end, we should perform an analytical continuation $R\rightarrow 0$. The point is that there is no guarantee at all that the function $\delta_{J}\braket{\mathscr{Z}^{R}}_{\text{dis}}$ is analytic all the way down to $R=0$ altland-book . This means that, from a mathematical standpoint, the replica trick is not a well-founded technique. Surprisingly, it is also true that it seems to rarely fail, at least when the disorder can be treated in a perturbative way parisi-book ; parisi-review .

Let us consider in the following the generalization of Eq. (10), namely given by

[TABLE]

where, actually, the disorder is assumed translationally invariant in time and space. In order to implement the replica trick within the context of superfluid bosonic systems, as an alternative to Eq. (19) we move to the phase-density representation, where

[TABLE]

with $n_{0}=\psi_{0}^{2}$, $\psi_{0}^{2}$ being given by Eq. (18). From the equation above we deduce that $n(\mathbf{r},\tau)=n_{0}+\pi(\mathbf{r},\tau)$ and a perturbative strategy is easily implemented by assuming the smallness of density fluctuation. Thus, Eq. (57) becomes

[TABLE]

In order to understand the limitations of our approach, it is crucial to recall that $\braket{n(\mathbf{r},\tau)}_{\text{th}}=\braket{}{\psi(\mathbf{r},\tau)}{{}^{2}}_{\text{th}}$, while we often set $n\simeq n_{0}$ with n_{0}=\big{|}\braket{\psi(\mathbf{r},\tau)}_{\text{th}}\big{|}^{2}. Up to the Gaussian level this is an acceptable approximation since $\braket{\pi(\mathbf{r},\tau)}_{\text{th}}\simeq 0$. In the regime of small fluctuations, we can expand the mean-field analysis by retaining only the quadratic terms in $\pi(\mathbf{r},\tau)$ and $\varphi(\mathbf{r},\tau)$. The change of variables $\{\psi,\psi^{*}\}\rightarrow\{\pi,\varphi\}$ has a constant Jacobian determinant, therefore Eq. (3) reads

[TABLE]

with the periodic boundary conditions $\pi(\mathbf{r},\beta\hbar)=\pi(\mathbf{r},0)$ and $\varphi(\mathbf{r},\beta\hbar)=\varphi(\mathbf{r},0)$. The smallness of fluctuations allows to extend the integration range of the phase variable from $[-\pi,+\pi]$ to $[-\infty,+\infty]$. This decompactification has obviously to be addressed with particular attention when dealing with $d\leq 2$, where the enhancement of fluctuations prevents the occurring of a true long-range order and a more refined analysis is required altland-book ; mermin-1966 ; hohenberg-1967 ; villain-1975 ; kadanoff-1977 . The Gaussian action in Eq. (59) is given by

[TABLE]

where we have imposed Eq. (18) to eliminate the dependence on the chemical potential and $U_{\text{D}}$ is distributed according to Eq. (56). The pertubative expansion up to the terms contained in Eq. (60) is equivalent to a one-loop expansion within a diagrammatic representation. Equivalently, in the second-quantization formalism this corresponds to the Bogoliubov approximation andersen-2004 .

The implementation of the replica trick requires the calculation of $\braket{\mathscr{Z}^{R}}_{\text{dis}}$. The (not-averaged) partition function is simply

[TABLE]

where $S_{g}[\pi_{\alpha},\varphi_{\alpha}]$ is provided by Eq. (60) and $\mathscr{D}[\pi,\varphi]=\prod_{\alpha=1}^{R}\mathscr{D}[\pi_{\alpha},\varphi_{\alpha}]$. In case of a Gaussian disorder obeying to Eq. (56), Eq. (6) reads

[TABLE]

where $S^{\text{(pure)}}_{g}[\pi_{\alpha},\varphi_{\alpha}]$ is Eq. (60) with $U_{\text{D}}=0$ and $S^{\text{(coup)}}_{g}[U_{\text{D}},\pi_{\alpha}]$ is given by

[TABLE]

The second line of Eq. (62) is a Gaussian integral in $U_{\text{D}}$ such that

[TABLE]

In this way, the replicated partition in Eq. (62) results in

[TABLE]

For the sake of clarity we report below $S_{g}^{\text{(pure)}}$ and $S_{g}^{\text{(dis)}}$, namely

[TABLE]

with $\mathbf{s}=(\mathbf{r},\tau)$. It is crucial to underline the relevance of Eq. (65). Indeed, it shows that also in presence of the simplest disorder potential ($\delta$-correlated both in time and space), the outlined formalism generates an effective quartic interaction between different replica indices. This fact can be understood by recalling the physical meaning of Feynman’s path. Despite its complexity, and indifferent to different replicas, the path integrals prefer to evolve towards the lowest possible potential. In other words they have a tendency to populate the same regions of the energy landscape.

Moving to the calculation of expectation values, the replicated formalism replaces the logarithm of the partition function with, basically, Eq. (65), where the Euclidean action of Eq. (60) is splitted into Eq. (66) and Eq. (67). According to Eq. (55), the joint thermal and disorder average is obtained through

[TABLE]

where $\braket{\bullet}_{\text{rep}}$ has to be intended as a thermal average over the replicated partition function in Eq. (65) or, equivalently, over the action given by the sum of Eq. (66) and Eq. (67), i.e.

[TABLE]

It is useful to expande the fields in Fourier series such that Eq. (66) reads

[TABLE]

where we have defined the vector $\chi_{\alpha}$ as \chi_{\alpha}(\mathbf{q},\omega_{n})^{T}\equiv\big{(}\varphi_{\alpha}(\mathbf{q},\omega_{n});\pi_{\alpha}(\mathbf{q},\omega_{n})\big{)} The pure inverse propagator in Eq. (70) is given by

[TABLE]

On the other hand, for Eq. (63) we have now

[TABLE]

The equation above can be decomposed in its diagonal and off-diagonal contributions, leading us to the final equation for the Euclidean action of Gaussian fluctuations. By renaming $Q=(\mathbf{q},\omega_{n})$ and using the fact that $x(-Q)=x^{*}(Q)$, we have

[TABLE]

In the equation above, the sum over the replicas is encoded in the $(R\times 2)$-dimensional column vector

[TABLE]

The $R\times 2$ matrix $\mathbb{M}^{-1}(Q)$ in Eq. (73) is built as

[TABLE]

with the generating blocks

[TABLE]

V Correlations functions in the replicated formalism

As detailed in App. A, the matrix $\mathbb{M}$ in Eq. (75) is diagonalized by a well-known unitary transformation, reading

[TABLE]

with

[TABLE]

Up to this level of approximation, correlation functions are specified by the entries of the fluctuation propagator. For the pure system, by assuming translational invariance both in time and space, we have

[TABLE]

On the other hand, according to Eq. (68), the replicated formalism provides both the disorder and the thermal average at the same time. For instance, concerning the density-density correlation function, we have

[TABLE]

where the final result is obtained by expanding the second line up to the first order in $R$. Moreover, we have

[TABLE]

nothing else than Eq. (29) with the saddle-point replacement $\mu=\tilde{V}(0)\,n_{0}$. The remaining correlations are computed in the same way, leading to

[TABLE]

Besides the phase and density correlators, one of the most important correlation functions is the current one, defined as

[TABLE]

where the notation $\braket{\ldots}$ represents, for now, a general average procedure. In the phase-density representation, from Eq. (57) it is easy to verify that the current density $\mathbf{g}=\hbar(n_{0}+\pi)\nabla\varphi/m$. Thus, by making use of the Wick theorem fetter-book to disentagle the four-field correlator in Eq. (84), in the Fourier space one obtains tauber-1997

[TABLE]

Thanks to the replica trick, the disorder contribution is taken into account by considering the correlators over the replicated Gaussian action, namely Eqs. (81) and (83).

The current response function in Eq. (85) is crucial since it encodes the superfluid character of the system. This point is clarified by noticing that $\mathscr{C}_{ij}$ is a 2-rank tensor, so it can be decomposed in its longitudinal and transverse components by means of the projectors $\mathscr{P}^{L}_{ij}=q_{i}q_{j}/q^{2}$ and $P^{T}_{ij}=\delta_{ij}-\mathscr{P}^{L}_{ij}$. It has been shown (see baym-1967 and ueda-book for a pedagogical review), that the static limit of transverse component of Eq. (85) $\mathscr{C}_{ij}^{T}(\mathbf{q},\omega_{n})\equiv(d-1)^{-1}\sum_{ij}\mathscr{P}_{ij}^{T}(\mathbf{q},\omega_{n})\mathscr{C}_{ij}(\mathbf{q},\omega_{n})/\hbar$ corresponds to the normal component of the system. By replacing the replicated correlator given by Eqs. (81) and (83) into Eq. (85), after burdensome (but standard) algebra, we find

[TABLE]

When the time scale of the disorder potential is frozen, i.e. for

[TABLE]

Eq. (86) becomes

[TABLE]

With the further simplification of a constant disorder correlator $\gamma^{2}$, the equation above equates to Eq. (52). So, Eq. (88) reduces corretly to the quenched disorder regime and agrees with results reported in tauber-1997 ; schakel-1997 ; giorgini-1994 .

In order to compute the condensate depletion through the replicated correlators in Eqs. (81) and (83), we rely upon the equation fetter-book

[TABLE]

Following the phase-density representation in Eq. (57), we expand the fields $\psi$ and $\psi^{*}$ up to the linear order in the fluctuations $\psi(\mathbf{r},\tau)\simeq\sqrt{n_{0}}[1+\pi(\mathbf{r},\tau)/(2n_{0})+i\varphi(\mathbf{r},\tau)]$, such that the resulting correlators in Eq. (89) are Gaussian. In the Fourier space we have

[TABLE]

where the prime signal the dependence on $Q^{\prime}=(\mathbf{q}^{\prime},\tau^{\prime})$ and, similarly to Eq. (84), $\braket{\ldots}$ corresponds to a general average procedure. In order to include the disorder contribution one simply has to consider the replicated correlators in Eqs. (81) and (83), where space-time translational invariance restored by assuming Eq. (80).

We then have

[TABLE]

The first line represents the contribution coming from the pure system, reading Eq. (40) at $T=0K$. On the other hand, the second line accounts for the disorder contribution. For a point-like frozen disorder as in Eq. (87), one can verify that

[TABLE]

By comparing the equation above and Eq. (88), at zero-temperature the disorder contribution to the condensate depletion and the normal component of the system are deeply related, namely

[TABLE]

VI Conclusions and future perspectives

In this review, we have reviewed the crucial aspects of a field-theory approach to superfluid bosons moving in a random environment. Within the powerful framework of functional integration, we have recovered the important result on the superfluid and condensate depletion in presence of an external disorder huang-1992 ; giorgini-1994 ; schakel-1997 . Moreover, we have shown how the theoretical description can be generalized beyond the assumption of a static external disorder by using the replicated formalism navez-2007 . Obviously, our pedagogical overview has to left out other interesting and crucial issues. For instance, we have pointed out in Sec. III that our perturbative approach lacks self-consistency, since we are expanding above a uniform background also in presence of an external potential. In order to overcome this point, in khellil-2016 ; khellil-2016-2 the authors present a self-consistent implementation relying upon the Hartee-Fock mean-field theory within the replica formalism. Moreover, our review does not include any discussion about the eventual occuring of the superfluid-glass transition in three or lower dimensions. This is a crucial topic, deserving a separate investigation. Numerical simulations based on quantum Monte Carlo methods are likewise important, since they are ab-initio calculation and may shed light on unexpected effects. For instance, in giorgini-2002 the non-trivial relation between superfluidity and condensation is numerically investigated, while in ng-2015 the worm algorithm provides an estimation for critical exponents of the superfluid-glass transition in two-dimensions.

Another promising line of research concerns the possibility that disorder does not always act in a parasitic way, generating, on the contrary, surprising dynamical effects. In the context of condensed matter theory, one of the most striking example is given by the quantum Hall transition pruisken-1987 , where impurity scattering plays a crucial role. On the other hand, very recently it was shown that a certain degree of random fluctuations, such as thermal noise, may enhance the transport properties of a particles ensemble. While Anderson localization anderson-1959 acts to halt the flow of a certain quantity, in caruso-2015 ; potocnik-2018 the authors observe a boosting of transport through optical fibers and superconducting circuits. The reason is the so-called environment-assisted quantum trasport, which relies crucially upon the presence of a certain degree of disorder. Indeed, while the coherence of the system is reduced by propagating in a disordered medium, this also reduces the possibility of destructive interference responsible for Anderson localization.

This phenomenon can pave the way to novel interesting protocol to engineer more efficient quantum devices, by tuning the coupling with the environment to enhance the trasport of a desired quantity. Recent experimental confirmations have been produced by using a one-dimensional array of trapped ions maier-2019 . Thus, it would be extremely interesting to review to effect in the context of ultracold atoms, both from an experimental and theoretical point of view.

Acknowlegdements

LS acknowledges for partial support the FFABR grant of Italian Ministry of Education, University and Research. AC thanks Andrea Tononi for useful discussions and meaningful insight on the manuscript.

Appendix A Diagonalization of a block circulant matrix

It is possible to perform a block diagonalization of $\mathbb{M}^{-1}$ in Eq. (75) by noticing its circulant structure. Indeed, $\mathbb{G}^{-1}$ is generated by cyclically permuting the elements (in this case $2\times 2$ matrices) of its first row. In this sense, we can formally write down

[TABLE]

It has been shown davis-book ; olson-2014 that, given $R$ the number of the generating matrices and $M$ their dimension, every block circulant matrix $\mathbb{A}$ is block diagonalized by the same unitary transformation. Indeed, one can verify that

[TABLE]

where the $\Lambda_{j}$ are the $2\times 2$ eigenblocks, specified by $\Lambda_{j}=\mathscr{M}^{-1}+\sum_{\alpha=1}^{R-1}\mathbb{B}\;w_{R}^{\alpha j}$. The $\otimes$ symbol denotes the direct product of the two matrices. In Eq. (95) the key ingredient is the Fourier matrix, defined as

[TABLE]

with $w_{R}=\exp(2\pi i/R)$ the fundamental root of the unity.