Abstract

We study the fate of algebraic decay of correlations in a harmonically trapped two-dimensional degenerate Bose gas. The analysis is inspired by recent experiments on ultracold atoms where power-law correlations have been observed despite the presence of the external potential. We generalize the spin wave description of phase fluctuations to the trapped case and obtain an analytical expression for the one-body density matrix within this approximation. We show that algebraic decay of the central correlation function persists to lengths of about 20% of the Thomas–Fermi radius. We establish that the trap-averaged correlation function decays algebraically with a strictly larger exponent weakly changing with trap size and find indications that the recently observed enhanced scaling exponents receive significant contributions from the normal component of the gas. We discuss radial and angular correlations and propose a local correlation approximation which captures the correlations very well. Our analysis goes beyond the usual local density approximation and the developed summation techniques constitute a powerful tool to investigate correlations in inhomogeneous systems.

pacs:

05.70.Jk, 64.60.an, 67.85.-d

With the achievement of quantum degeneracy in ultracold atom gases a new experimental platform for studying fundamental questions related to phase transitions and critical phenomena has appeared. In particular, correlation functions can be measured through interference and time-of-flight techniques Polkovnikov et al. (2006); Hadzibabic et al. (2006); Cladé et al. (2009); Tung et al. (2010); Plisson et al. (2011); Singh and Mathey (2014); Murthy et al. (2014). A crucial aspect of the experiments, however, consists in the absence of translation invariance due to the presence of the external trapping potential. Consequently, correlations between points r and r′ are not uniquely determined by the value of r−r′. This effect should be unimportant for short-range correlations, and, in fact, thermodynamic properties of ultracold gases are often very well-captured by a local density approximation. The situation changes drastically for a system at a critical point, where the correlation length diverges and thus competes with the inhomogeneity of the trapping potential.

Whereas the experimental preparation of critical systems typically requires a highly fine-tuned setup, they are rather effortlessly realized in two-dimensional (2D) systems whose low-energy excitations can be mapped onto an XY model. Examples apart from 2D ultracold quantum gases Hadzibabic and Dalibard (2011); Desbuquois et al. (2012); Ries et al. (2015); Fletcher et al. (2015); Murthy et al. (2015) are given by thin Helium films Bishop and Reppy (1978), layered magnets Dürr et al. (1989); Ballentine et al. (1990), and 2D exciton-polariton condensates Roumpos et al. (2012). To quantify correlations in these systems we introduce the one-body density matrix ρ(r,r′)=⟨^Φ†(r)^Φ(r′)⟩, where ^Φ†(r) is the creation operator for a particle at point r. For the spatially homogeneous case we then have ρhom(r,r′)=f(|r−r′|) with some function f(r) due to translation and rotation invariance. The Mermin–Wagner–Hohenberg theorem Mermin and Wagner (1966); Hohenberg (1967) forbids long-range order at any finite temperature such that limr→∞f(r)=0.
However, in the XY model an ordered low-temperature phase with infinite correlation length exists, where correlations decay with a temperature-dependent scaling exponent η as f(r)∼r−η for large r.
Since η is well below unity, this behavior of correlations is named
quasi-long-range order (QLRO). Above a transition temperature, vortices disorder the system and QLRO is lost Holzmann and Krauth (2008); Choi et al. (2013). The described picture has been developed by Berezinskii and Kosterlitz, Thouless (BKT) Berezinskii (1971, 1972); Kosterlitz and Thouless (1973); Kosterlitz (1974).

The presence of the trap raises several conceptual questions towards the validity of the BKT picture in ultracold quantum gas experiments: On which length scales is QLRO visible? If samples bigger than the state of Texas are needed to reach the thermodynamic limit in the homogeneous system Bramwell and Holdsworth (1994), do we need ultracold atomic clouds covering the whole continent for the local density approximation to hold?
Does the formula η=MkBT2πℏ2ns for a Bose gas with mass M, temperature T, and superfluid density ns imply a locally varying scaling exponent due to the inhomogeneous density? Can the 2D XY-model explain the large scaling exponent η>1 recently observed in experiment and Quantum Monte Carlo (QMC) computations Murthy et al. (2015)? Can we obtain the correlation function of the trapped system from the homogeneous one by a generalization of the local density approximation?

In this Rapid Communication, we address and answer all of these questions within the spin wave approximation for the phase fluctuations of the trapped gas. The main results are (1) the generalization of the textbook spin wave theory to a discrete collective mode spectrum, (2) the derivation of explicit analytical expressions for the one-body density matrix and first-order correlation functions that are suited for practical applications, and (3) the demonstration that the trap-averaged correlation function decays algebraically with an increased exponent. Our calculations provide a simple picture to the large exponents observed in experiment and simulation. We limit the derivation to the key steps and invite the mathematically interested reader to consult the detailed Supplemental Material (SM) SOM ().

Spin wave theory. At sufficiently high temperatures, quantum corrections are small and the macroscopic properties of the trapped gas are described by the classical action

S[Φ]=1T∫d2r[

|∇Φ|22M+(−μ+V(r))|Φ|2+g2|Φ|4]

(1)

for the complex bosonic field Φ(r)=A(r)eiθ(r). Herein, μ and g are the chemical potential and coupling constant, and we use units such that ℏ=kB=1. We assume an isotropic and harmonic trapping potential V(r)=M2ω2r2.
The action is stationary towards small changes in A(r) if (−∇22M−μ+V+gA20+(∇θ)22M)A0=0 is satisfied.
We neglect gradient terms of the fields and obtain a Thomas–Fermi (TF) density profile given by n(r)=A20(r)=n0(1−r2/R2) with n0=μ/g and TF radius R=(2gn0/Mω2)1/2. By neglecting deviations of A(r) from the stationary configuration we then arrive at the action for the phase θ(r) given by

Sph[θ]=n02MT∫Rλd2r(1−r2R2)(∇θ)2.

(2)

We explicate the length scales that delimit the validity of our description, given by the TF radius R at large distances, and a microscopic scale λ, which is given by the thermal wavelength in cold atom experiments, or by the lattice spacing in spin models. For typical 2D ultracold quantum gases we have R≃100μm and λ≃1μm, such that R/λ≃100 is large.

After partial integration the phase-only action can be written as quadratic form Sph[θ]=n02MT(θ,DRθ) with self-adjoint operator

DR=−(1−r2R2)∇2+2r⋅∇R2

(3)

and scalar product (f,g)=∫Rd2rf(r)g(r). The normalized eigenfunctions of DR for R<∞ are given by

Since the phase-only action Sph is quadratic, though with a nontrivial discrete spectrum, correlation functions are obtained via Gaussian integration.
Within the spin wave approximation we have ρ(r,r′)=A0(r)A0(r′)⟨ei[θ(r)−θ(r′)]⟩=A0(r)A0(r′)e−12⟨Δθ(r,r′)2⟩ with Δθ(r,r′)=θ(r)−θ(%r′) and

⟨Δθ(r,r′)2⟩=MTn0∑n,m′|θnm(r)−θnm(r′)|2εnm.

(6)

In the homogeneous limit (R=∞) we have D∞=−∇2 and the eigenfunctions are plane waves θq(r)=eiq⋅r with energies εq=q2, and

⟨Δθ(r,r′)2⟩hom

=MTn0∫λ−1d2q(2π)2|θq(r)−θq(r′)|2εq

≃η0log(|r−r′|2λ2)for|r−r′|≫λ,

(7)

with η0=MT2πn0Herbut (2007); Altland and Simons (2010). Due to sm/2P(m,0)n(1−2s)∼Jm(2n√s) and εnm∼4n2 for large n and small s, we recover this formula for R→∞. The homogeneous case is recalled in the SM SOM ().
In the homogeneous system, spin wave theory breaks down for η0≃0.25Kosterlitz and Thouless (1973); Kosterlitz (1974). In a trap it is applicable to the superfluid core in the center of the atomic cloud. We estimate its radius r from n(r)λ2>4 as r2/R2<1−4η0. For η0=0.10 (0.20) we have r<0.8R (0.4R), which is sufficiently large to observe the effects studied here. Much of the elegance of BKT theory stems from replacing the integral in Eq. (7) by the logarithm, which is valid for all relevant separations, but not exact.

Evaluation of the sum. Just as the momentum integration in the homogeneous case is limited by qλ≤1, also Eq. (6) has an ultraviolet cutoff on possible values of (n,m), indicated by the prime. In fact, besides (n,m)≠(0,0), the values of n and m are limited from above according to εnm≤εmax∼(R/λ)2.
In the following we choose summation cutoffs N=e−γ(R/λ),M=e−γ(R/λ)2
with Euler’s constant γ≃0.577SOM (). We can then write

⟨Δθ(r,r′)2⟩

=η02[F(N)0(s,s)−2F(N)0(s,s′)+F(N)0(s′,s′)]

+η0M∑m=1[F(Nm)m(s,s)−2cos(mΔϕ)F(Nm)m(s,s′)

+F(Nm)m(s′,s′)]

(8)

with Nm=−m+12+√2M+m2+12 and

F(N)0(s,s′)

=N∑n=12n+1n(n+1)Pn(1−2s)Pn(1−2s′),

F(Nm)m(s,s′)

=Nm∑n=02n+m+1m2+n(n+m+1)(ss′)m/2

×P(m,0)n(1−2s)P(m,0)n(1−2s′).

(9)

Here and in the following we denote

r=(rcosϕrsinϕ),s=r2R2,Δϕ=ϕ−ϕ′,

s>=max(s,s′),s<=min(s,s′),^λ=λR.

(10)

Equation (8) can be implemented numerically and thus allows the computation of the phase fluctuations in a trap for arbitrary values of R/λ. Note in particular that ⟨Δθ(r,r)2⟩=0. If R/λ is large, however, analytical approximation formulas can be derived from the N,M→∞ limits.

Evaluating Eq. (9) for N,Nm=∞ relies on the observation that for a self-adjoint operator ^O with eigenvalues εi and eigenfunctions fi(s), the function g(s,s′)=−∑ifi(s)fi(s′)εi is the Green function of the operator due to the completeness of the fi. Accordingly, if we know the Green function by some other means, we also know the result of the summation Englis and Peetre (1998); Gustavsson (2001). Now, for fixed m the operator DR is of Sturm–Liouville-type and thus its Green function is given by u(s<)v(s>), where u,v are zero modes of the operator. With this approach we find

F(∞)0(s,s′)

=−1−log[s>(1−s<)],

(11)

F(∞)m(s,s′)

=1m(s<s>)m/2um(s<)vm(s>)

(12)

for 0<s,s′<1. We have um(s)=2F1(am,bm,m+1,s) and vm(s)=Γ(am)Γ(bm)Γ(m)2F1(1−am,1−bm,1,1−s)
with hypergeometric function 2F1(a,b,c,z) and am=m+12+√m2+12, bm=m+12−√m2+12 such that am+bm=m+1 and ambm=m/2Abramowitz and Stegun (1964); SOM (). The key properties of the functions um and vm are um(0)=vm(0)=1 and

limm→∞um(s)=limm→∞vm(s)=1√1−s.

(13)

We define the remainder function

R(s,s′,Δϕ)

=M0∑m=11m(s<s>)m/2cos(mΔϕ)

(14)

×(um(s<)vm(s>)−1√(1−s)(1−s′)),

which converges quickly in M0 and is well-approximated by summing up the first M0=10 terms. The careful limit N,M→∞ in Eq. (8) is presented in the SM SOM ().

Correlation functions. The phase fluctuations due to Eq. (8) for R/λ≫10 are given by

⟨Δθ(r,r′)2⟩

≃η02(1+11−s>)log(¯s>^λ2(1−s>))−η02(1−11−s<)log(¯s<^λ2(1−s<))

+η0√(1−s)(1−s′)log(Δ¯r2¯s>)+η0(R(s,s,0)+R(s′,s′,0)−2R(s,s′,Δϕ))

(15)

with ¯s=max{s,^λ2},Δ¯r2=max{|r−r′|2R2,^λ2(1−s)}.
The latter two definitions give a proper way to treat the singularities that occur when r,r′,|r−r′|→0. Their origin is already visible in the homogeneous limit (7): Whereas the integral is well-defined for all r,r′, the logarithm expression requires a meaningful procedure of how to treat the ultraviolet singularities. The remainder term involving the R-functions in Eq. (15) can be neglected for most cases of interest. The homogeneous limit η0log(|r−r′|2/λ2) is recovered for ^λ2≪s,s′≪1, i.e., λ≪r,r′≪R.

Figure 1: Selection of correlation functions for R/λ=100 and η0=0.15 in terms of ρ(r,r′,Δϕ)=⟨Φ∗(r)Φ(r′)⟩.
Black data represents the sum from Eq. (6) (points) and the analytic result from Eq. (15) (solid line). The LCA from Eq. (19) is indicated by the dashed orange lines. Upper panel. Central correlation function g1(r)=ρ(r,0,0) (upper curve) decaying algebraically with exponent η0, and trap-averaged correlation function ¯g1(r) (lower curve) decaying with strictly larger exponent ηtrap>η0. Long-dashed red lines correspond to (r/λ)−η with η=η0 and η=1.6η0, respectively. The inset shows the same data on a linear scale. Lower panel. Radial correlation function grad(r,r′)=ρ(r,r′,0) for fixed r/R=0.1,0.5,0.8 (from left to right). Note that the peaks heights are (1−r2/R2). The inset shows the angular correlation function gang(r,Δϕ)=ρ(r,r,Δϕ) for the same values of r/R (from top to bottom). The long-dashed blue lines show (1−r2/R2)(|r−r′|/λ)−η0, which deviates considerably.

The central correlation function g1(r)=ρ(r,0) derived from Eq. (15) is given by

g1(r)=n0(1−s)12[1+η(s)](rλ)−η(s)e−η02R(s,s,0)

(16)

with η(s)=η02(1+11−s).
The factor e−η02R is close to unity and may be neglected. The algebraic decay with η0 for small s≪1 agrees with the central correlation functions obtained in Refs. Petrov et al. (2000); Crecchi and Vicari (2011); Ceccarelli et al. (2013). The angular correlations for r=r′ read

⟨Δθ(r,r,Δϕ)2⟩

=η01−slog(2s(1−cosΔϕ)^λ2(1−s))

+2η0[R(s,s,0)−R(s,s,Δϕ)]

(17)

for Δϕ2≥^λ2(1−s)/s. (This corresponds to |r−r′|2≥λ2(1−s). For Δϕ→0 the correlations vanish.) In Fig. 1 we plot a selection of correlation functions.

Figure 2: Scaling exponent ηtrap extracted from ¯g1(r)∼(r/λ)−ηtrap for R/λ=20,100,1000 (blue, dark red, orange). The solid lines give the result of applying Eq. (15) to the whole cloud. As a first estimate of the influence of the superfluid transition we replace the phase correlations in the outer regions of the cloud by an exponential decay with correlation length λ. This yields an increase of ηtrap for small R/λ (dashed lines). The red triangles correspond to the QMC data from Ref. Murthy et al. (2015) for the quasi-2D gas with ~g=0.60 and R/λ≈50, with η0 determined from the scaling of g1(r). The error bars estimate the fitting error of η0, which is substantial due to its small value.

Applications. To obtain a useful estimate of the trapped correlation function from the homogeneous result, we write Eq. (7) as

⟨Δθ(r,r′)2⟩hom=2G(r,r′)−G(r,r)−G(%r′,r′)

(18)

with G(r,r′)=η02log(|r−r% ′|2/R2),G(r,r)=η02log(^λ2), the Green function of −MTn0∇2. In each G separately we apply a local density approximation η0→η(r,r′)=MT2π√n(r)n(r′)=η0√(1−s)(1−s′) and obtain

⟨Δθ(r,r′)2⟩LCA

=η02[11−s+11−s′]log(1^λ2)

+η0√(1−s)(1−s′)log(|r−r′|2R2).

(19)

We name this particular procedure local correlation approximation (LCA), and observe from Fig. 1 that it works sufficiently well to approximate the actual result. In contrast, the approximation ρ(r,r′)≈A(r)A(r′)(|r−r′|/λ)−η0 fails considerably.

In a translation invariant setting we have ¯g1(r)∝g1(r), whereas both functions differ in the trapped case. As shown in Figs. 1 and 2, ¯g1(r) decays algebraically for λ≤r≪R with an increased exponent ηtrap≃(1.5-2)η0, which weakly depends on R/λ. In the SM SOM () we analytically estimate ηtrap, which supports the exceptionally weak R/λ-dependence. Although ηtrap→η0 for infinitely large systems, it would require absurd values of R/λ.

The spin wave description only applies to the superfluid core of the cloud. Hence, applying Eq. (15) to the whole cloud becomes less accurate for increasing values of η0 and underestimates ηtrap. To estimate the effect of the normal gas on the trap-averaging we set ⟨ei(θ(r)−θ(r′))⟩=e−|r−r′|/λ for (r,r′)∉[0,rs]×[0,rs] with approximate core radius rs=(1−4η0)1/2RSOM (). Assuming the correlation length in the normal gas to coincide with λ should give the right order of magnitude of the effect. Computing ¯g1(r) with this modification of ρ(r,r′), we find an increased scaling exponent ηtrap which is slightly above the spin wave prediction for η0≲0.15. Increasing η0 further, we obtain exponents up to ηtrap≈1.5 for small R/λ≲100, whereas the algebraic decay disappears for larger R/λ in the typical fitting range r≲0.1R. Since most experiments are in the first regime, an enhanced scaling exponent is likely to be measured. Neglecting the spin wave part and only keeping exponential correlations, no algebraic decay of ¯g1(r) is observed.

In order to relate to the enhanced exponents found in Ref. Murthy et al. (2015), we compare to the QMC results where η0 is determined by the decay of the central correlation function Holzmann et al. (2010)2. As shown in Fig 2, we find good qualitative agreement with the QMC results when including the exponential decay of correlations in the normal component. Hence we identify the superfluid core, a significant normal component, and a ratio of R/λ∼102 as a possible explanation for the large universal transition exponent ηtrap,c≃1.4.
Although it is possible to include more quantitative knowledge of the quasi-2D equation of state Holzmann et al. (2008, 2010); Prokof’ev and Svistunov (2002); Boettcher et al. (2016); Wu et al. (2015) into the LCA approximation, it is clear that our qualitative findings are robust against these refinements. A full quantitative analysis of the BKT transition in the trapped gas is left for future work.

In conclusion, we have generalized the spin wave theory of phase fluctuations towards the trapped 2D Bose gas, analytically computed correlation functions, and applied this to obtain the trap-averaged correlation function. Our analysis provides a natural explanation for the unexpected large exponents observed in Ref. Murthy et al. (2015). It should be possible to apply the method to other dimensions and polytropic equations of state, just as it is the case for superfluid hydrodynamic modes Fuchs et al. (2003); Heiselberg (2004); Boettcher et al. (2011). We proposed an LCA prescription which can be further developed by considering these other mode spectra. An improved understanding of correlations in inhomogeneous systems is also mandatory for conducting experiments on universality in disordered systems Allard et al. (2012); Beeler et al. (2012); Carleo et al. (2013) and situations far from equilibrium Mathey et al. (2011); Altman et al. (2015); Dagvadorj et al. (2015).

Acknowledgments. We gratefully acknowledge inspiring discussions with A. Böttcher, R. Ganesh, I. F. Herbut, P. A. Murthy, and M. M. Scherer. IB acknowledges funding by the European Research Council (ERC) under Advanced Grant No. 290623 and the German Research Foundation (DFG) under Grant No. BO 4640/1-1. MH was supported by ANR-12-BS04-0022-01, ANR-13-JS01-0005-01 and DIM Nano’K from Région Île-de-France.

Appendix A Computation of the correlation function

In this section we derive the formula for the phase fluctuations in detail. The presentation provides all the intermediate steps that have been left out in the main text for reasons of brevity.

We start by solving the eigenvalue problem for the operator

DR=−(1−r2R2)∇2+2r⋅∇R2

(S1)

appearing in the phase-only action

Sph[θ]=n02MT∫Rλd2r(1−r2R2)(∇θ)2=n02MT(θ,DRθ).

(S2)

We introduce dimensionless variables ^r=r/R and ^D=R2DR to obtain the eigenvalue problem

^Dθε=εθε.

(S3)

We employ a polar coordinate Ansatz θε(r)=∑∞m=−∞θεm(^r)eimϕ to arrive at

−(1−^r2)(d2d^r2+1^rdd^r−m2^r2)θεm+2^rdd^rθεm=εθεm.

(S4)

Note that since ^D is invariant under r→−r, the eigenvalues only depends on |m|. With the new coordinate

s=^r2=r2/R2

(S5)

we eventually arrive at ^D(m)θεm=εθεm with

^D(m)=−dds[4s(1−s)dds]+m21−ss.

(S6)

Note that this operator is of Sturm–Liouville form, see Eq. (S114). The eigenfunctions of ^D(m) are easily seen to be of the form s|m|/2P(|m|,0)n(1−2s) with P(α,β)n being the Jacobi polynomials, see Eq. (S98). The corresponding eigenvalues are given by

εnm=2|m|+4n(n+|m|+1)

(S7)

with n=0,1,2,… such that (n,m)≠(0,0). For m=0 the solution reduces to the Legendre polynomial Pn(1−2s). In App. D.1 we collect relevant properties of the Legendre and Jacobi polynomials.

We define the normalized eigenfunctions

θnm(r)=√2n+|m|+1πs|m|/2P(|m|,0)|n(1−2s)eimϕ

(S8)

such that

^Dθnm=R2DRθnm=εnmθnm,

(S9)

and, due to Eq. (S95), (θnm,θn′m′)=R2δnn′δm,−m′. Due to the completeness of the Jacobi polynomials, every regular phase configuration θ(r) can be expanded according to

θ(r)=∑n,manmθnm(r).

(S10)

Since θ(r) is real we have a∗nm=an,−m. We arrive at the regular phase-only action

Missing dimension or its units for \hskip

(S11)

The smallest and largest length scales λ and R introduced in Eq. (S2) manifest in lower and upper boundaries on the summation (indicated by a prime). The implications of this are addressed below in detail.

We now consider the one-body density matrix

ρ(r,r′)=⟨Φ∗(r)Φ(r′)⟩.

(S12)

Note that since ⟨Φ(r)⟩=0 this correlation function coincides with the connected correlation function ⟨Φ∗(r)Φ(r′)⟩−⟨Φ∗(r)⟩⟨Φ(r′)⟩. We write Φ(r)=A(r)eiθ(r) and assume θ(r) to be regular. Approximating A(r) by the stationary solution A0(r)=√n(r)=√n0(1−s) we arrive at

ρ(r,r′)

=√n(r)n(r′)⟨e−i(θ(r)−θ(r′))⟩

(S13)

with

⟨O⟩=∫DθOe−Sph[θ]∫Dθe−Sph[θ].

(S14)

Since Sph[θ] is quadratic and ⟨θ(r)⟩=0 we have

⟨e−i(θ(r)−θ(r′))⟩=e−12⟨[θ(r)−θ(r′)]2⟩

(S15)

according to the rules of Gaussian integration.

We proceed by computing

⟨Δθ(r,r′)2⟩=⟨[θ(r)−θ(r′)]2⟩

(S16)

with the Gaussian action Sph. We write the functional measure as Dθ=∏ν,μ∫daνμ and find

⟨[θ(r)−θ(r′)]2⟩

=∫Dθ(θ(r)−θ(%r′))(θ(r)−θ(r′))e−Sph[θ]∫Dθe−Sph[θ]

=∑n,m∑n′,m′(θnm(r% )−θnm(r′))(θn′m′(r)−θn′m′(r′))

×∏ν,μ∫daνμanman′m′e−12γνμaνμaν,−μ∫daνμe−12γνμaνμaν,−μ

(S17)

with γnm=n0εnm/MT. In the product, all terms except one yield unity. The nontrivial one is δnn′δm,−m′1γnm due to

∫daa2e−12aγa∫dae−12aγa=1γ

(S18)

for a real variable a. This eventually yields

⟨Δθ(r,r′)2⟩=MTn0∑n,m′|θnm(r)−θnm(r′)|2εnm.

(S19)

The macroscopic description of the trapped Bose gas in terms of the phase-only action is restricted to length scales r≫λ. This implies an ultraviolet energy cutoff according to Emax∼λ−2. The energies εnm=2|m|+4n(n+|m|+1) are thus restricted according to

0<εnm<εmax∼(Rλ)2

(S20)

with εmax=R2Emax. For m=0 the largest term in the n-sum is

N=N0=−1+√1+εmax2≃√εmax2∼(Rλ).

(S21)

The largest value of m, denoted by M, is attained for n=0 and is found to be

M=εmax2∼(Rλ)2.

(S22)

If 0<m<M, the possible values of n are limited from above by Nm<N0. Given m>0 we have

2m+4Nm(Nm+m+1)<εmax

(S23)

for the largest value of n=Nm. Thus

Nm

=−m+12+√εmax+m2+12

=−m+12+√2M+m2+12.

(S24)

In the following we choose

N=e−γ(Rλ),M=e−γ(Rλ)2

(S25)

with Euler’s constant γ≃0.577. For large N and M this implies

2HN=log(1^λ2),HM=log(1^λ2)

(S26)

with harmonic number Hn=∑nν=11ν. We will see below that fixing N and M according to Eqs. (S26) leads to convenient expressions for the phase fluctuations by trading N,M for ^λ. Note, however, that there is no εmax such that Eqs. (S25) can be true at the same time. On the other hand, due to e−γ≃0.56≃12, Eqs. (S25) are close approximations to εmax=(R/λ)2, i.e., with a proportionality factor of unity in Emax∼λ−2.