Superfluid/Bose-glass transition in one dimension

# Superfluid/Bose-glass transition in one dimension

Zoran Ristivojevic Centre de Physique Théorique, Ecole Polytechnique, CNRS, 91128 Palaiseau, France    Aleksandra Petković Laboratoire de Physique Théorique, IRSAMC, CNRS and Université de Toulouse, UPS, 31062 Toulouse, France    Pierre Le Doussal Laboratoire de Physique Théorique–CNRS, Ecole Normale Supérieure, 24 rue Lhomond, 75005 Paris, France    Thierry Giamarchi DPMC-MaNEP, University of Geneva, 24 Quai Ernest-Ansermet, CH-1211 Geneva, Switzerland
06-07-2019  02:38
###### Abstract

We consider a one-dimensional system of interacting bosons in a random potential. At zero temperature, it can be either in the superfluid or in the insulating phase. We study the transition at weak disorder and moderate interaction. Using a systematic approach, we derive the renormalization group equations at two-loop order and discuss the phase diagram. We find the universal form of the correlation functions at the transitions and compute the logarithmic corrections to the main universal power-law behavior. In order to mimic large density fluctuations on a single site, we study a simplified model of disordered two-leg bosonic ladders with correlated disorder across the rung. Contrarily to the single-chain case, the latter system exhibits a transition between a superfluid and a localized phase where the exponents of the correlation functions at the transition do not take universal values.

###### pacs:
71.10.Pm, 64.70.pm, 64.70.Tg, 71.30.+h

## I Introduction

The combined effect of disorder and interactions is one of the most fascinating problems of quantum correlated systems. Indeed, interactions can lead to collective effects such as superconductivity or superfluidity for which the many-body function is known to be resistant to disorder. On the other hand, for single particles the quantum nature of the problem is known to strongly reinforce the effects of disorder leading to the celebrated Anderson localization, Anderson (1958), for which the system is an insulator. One can thus expect a fierce competition between these two phenomena.

One of the systems for which this competition manifests in its strongest possible way are disordered bosons. Indeed, in this case one can expect a competition between superfluidity and Anderson localization. In one dimension it was shown by a renormalization group analysis Giamarchi and Schulz (1988) that disordered interacting bosons would undergo a phase transition between a superfluid and a localized phase. This transition and this phase, nicknamed Bose glass, was also shown Fisher et al. (1989) by scaling arguments to exist in higher dimensions. The phase diagram of disordered interacting bosons in one dimension possesses several remarkable features. First, the interactions lead to a reentrant superfluid phase, leading first to a delocalization from Bose glass to superfluid, when increased, then to a second transition at stronger interactions between the superfluid and the Bose glass. Giamarchi and Schulz (1988) Second, the transition at moderate interactions is in the Berzinskii-Kosterlitz-Thouless universality class, and universal exponents exists for the various correlation functions, in particular the single-particle one, at the transition. Further analyses of the transition in terms of vortex proliferation Kashurnikov et al. (1996) confirmed and generalized this result and numerical studies of disordered bosons Prokof’ev and Svistunov (1998); Rapsch et al. (1999); Pollet et al. (2009) confirmed both the reentrant nature of the phase diagram and the nature of the moderate interaction transition.Refael and Altman (2013); Pollet (2013) Extensions of this transition to finite temperature have also recently attracted a great deal of attention. Aleiner et al. (2010)

On the experimental front, the superfluid/Bose-glass transition has regained a considerable interest thanks to cold atomic gases which have provided remarkable realization of this problem Billy et al. (2008); Roati et al. (2008); Pasienski et al. (2010) as well as realization of disordered bosons by magnetic insulators. Giamarchi et al. (2008); Hong et al. (2010); Yamada et al. (2011); Yu et al. (2012); Zheludev and Roscilde (2013); Ward et al. (2013) This has led to a hunt for the Bose glass, as well as new questions on the phase diagram. In particular new studies focused on the low interaction, strong disorder case. Lugan and Sanchez-Palencia (2011); Fontanesi et al. (2011) In this regime, a real-space renormalization group study by Altman et al. (2010) of a related disordered Josephson junction array model found again a Berzinskii-Kosterlitz-Thouless transition, but with disorder-dependent exponents at the transition. These results which are still being debated Balabanyan et al. (2005); Hrahsheh and Vojta (2012); Pollet et al. (2013); Pielawa and Altman (2013); Pollet et al. (2014) nevertheless strongly suggest the existence of a different universality class for the weak interaction Bose-glass/superfluid transition than for the moderate interaction one, with universal exponents, and thus a richer phase diagram, a possibility already discussed in Ref. Giamarchi and Schulz, 1988.

In order to analyze further the moderate interaction side of the phase diagram we had in a recent paper Ristivojevic et al. (2012a) analyzed the moderate interaction superfluid/Bose-glass transition at the next order in the renormalization group analysis. This study confirmed the universality of the exponent at the transition, in a way consistent with the previous studies, and that non-universal terms such as those appearing in slightly different disordered problems such as the Cardy-Ostlund model Cardy and Ostlund (1982) were indeed not generated. It also allowed a precise calculation of the correlation functions at the transition including their logarithmic corrections to the main power law.

In this paper we give a detailed account of the renormalization technique used in Ref. Ristivojevic et al., 2012a, since it can be useful to related problems. We also detail the calculation of the correlation functions. In addition, we generalize the study to more complex models than a single chain of bosons. In particular, we investigate the case of a disordered bosonic ladder with correlated disorder along the rung (for the uncorrelated case see Ref. Orignac and Giamarchi, 1998). This system caricatures a bosonic system where density fluctuations on a single rung can be larger than for a single chain. We show that in that case the exponent is not universal at the transition in contrast to the case of the single bosonic chain.

The plan of the paper is as follows. In Sec. II we introduce the model that we treat using the replica method. We calculate the correlation functions with respect to the harmonic part of the model. In Sec. III we calculate the effective action at weak disorder. In Sec. IV we obtain the renormalization group flow equations, and solve them finding the phase diagram. In Sec. V we calculate the generator of connected correlations that enables us to calculate the density-density correlation function in Sec. VI and the single-particle correlation function in Sec. VII. In Sec. VIII we consider a two-leg ladder problem of disordered bosons, which exhibits nonuniversal exponents at the superfluid/Bose-glass transition. Section IX contains the discussions, which are followed by conclusions. Some details of calculations are presented in Appendices A, B, C, and D.

## Ii Model

We study the system of one-dimensional interacting bosons in a disordered potential. At low energies, the clean system can be described by the Tomonaga-Luttinger Hamiltonian Haldane (1981a); Giamarchi (2003); Cazalilla et al. (2011)

 H0=ℏ2π∫dx{vK[∂xθ(x)]2+vK[∂xφ(x)]2}. (1)

Here denotes the sound velocity and is a dimensionless Luttinger liquid parameter that accounts for the interaction strength. The fields and satisfy the bosonic commutation relation .

The parameters and of the phenomenological Hamiltonian (1) can be related to the parameters of a specific microscopic Hamiltonian. For example, the microscopic Lieb-Liniger model Lieb and Liniger (1963) describes bosons with local repulsion of the strength . At low energies, the latter model can be represented by the Hamiltonian (1). The parameter of Eq. (1) is such that at very strong interaction when , we have ,Ristivojevic (2014) while in the opposite case of weak interaction , we have .Cazalilla et al. (2011) Therefore, at any strength of the repulsive contact interaction. The parameter regime can be reached only by including longer range interaction in the microscopic model. In Galilean invariant microscopic models, which include the Lieb-Liniger one, the product is fixed and determined by ratio between the mean density and the mass of physical particles forming the bosonic system.

The density of bosonic particles can be expressed in terms of the displacement fields as Haldane (1981b)

 ρ(x)=ρ0−1π∂xφ(x)+2ρ2cos[2φ(x)−2πρ0x]. (2)

In Eq. (2), denotes the mean density while the remaining terms account for the density fluctuations. The second term in Eq. (2) describes long-wavelength fluctuations, while the oscillatory term describes the density fluctuations around the wave vector . The constant is nonuniversal and depends on microscopic details. Cazalilla et al. (2011)

In the following, we consider the system of disordered bosons and account for the effects of the disorder via its coupling to the particle density (2). Thus, we include additional terms in the model that are of the form

 Hd=∫dx{−1πη(x)∂xφ+ρ2[ξ∗(x)ei2φ+h.c.]}. (3)

Here we distinguish the so-called forward and backward scattering caused by the disorder potential. They respectively have Fourier components around the wavevectors and .Giamarchi (2003) We consider Gaussian disorder where the fields and have mean values zero and correlations

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯η(x)η(x′)=ℏ2Dfδ(x−x′), (4) ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ξ(x)ξ∗(x′)=ℏ2Dbδ(x−x′), (5) ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯η(x)ξ(x′)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ξ(x)ξ(x′)=0. (6)

We therefore assume that and characterize the disorder strength. Here and in the following denotes the disorder average. The total Hamiltonian is

 H=H0+Hd. (7)

In order to treat the effects of disorder we employ the replica method.Giamarchi (2003) We reexpress the free energy , where denotes the partition function, as . The latter form enables us to introduce the replicated action that contains the disorder average. It is defined via the relation

 ¯¯¯¯¯¯¯Zn=∫(n∏α=1DθαDφα)e−Srep/ℏ, (8)

where and are the replicated fields. A technical difficulty of the replica method is that the replicated action depends on fields, in contrast to two fields present in the initial Hamiltonian.

For the model (7), we conveniently split the replicated action as

 Srep=S0+Sf+Sb. (9)

The first term in Eq. (9) is quadratic as it arises from . It has the form

 S0ℏ= 12π∑α∫dxdτ[vK(∂xθα)2+vK(∂xφα)2 +2i(∂xθα)(∂τφα)]. (10)

The second term in Eq. (9) originates from -dependent part of the disordered potential (3). It is also quadratic, but off-diagonal in replica space,

 Sfℏ=−Df2π2∑αβ∫dxdτdτ′[∂xφα(x,τ)][∂xφβ(x,τ′)]. (11)

The third term in Eq. (9) is anharmonic and originates from the disordered potential expressed through the field in Eq. (3). It describes the effects of backward scattering on the disordered potential and reads

 Sbℏ=−ρ22Db∑αβ∫dxdτdτ′cos[2φα(x,τ)−2φβ(x,τ′)]. (12)

We therefore see that Gaussian disorder enters into the replicated action through the parameters describing its variance, and , see Eqs. (4) and (5).

The disorder potential does not couple to the phase field . It is therefore sometimes useful to integrate out the fields from the replicated action. As it involves Gaussian integrations, the resulting quadratic part in the action can be found exactly,

 S0ℏ= v2πK∑α∫dxdτ[(∂xφα)2+1v2(∂τφα)2+m2φ2α]. (13)

In Eq. (13) we have introduced a small mass that represents an infrared cutoff. On some occasions, we perform calculations at finite , taking the limit at the end. An example is the evaluation of the effective action, see Sec. III. However, for calculation of some quantities we can use from the beginning. We encounter this situation when calculating the density-density correlation function and the single-particle correlation function. We emphasize that in the most part of this article we use the expression (13) rather than (II). The latter expression we only employ to calculate the single-particle correlation function. This is expected, as fields are involved in its definition.

In this work we study the case of small anharmonic term (12) that enables us to use perturbation theory. We therefore need the correlation function

 Gαβ(x,τ)=⟨φα(x,τ)φβ(0,0)⟩, (14)

where denotes an average with respect to the quadratic action . Introducing the Fourier transform as

 φα(x,τ)=∫dkdω(2π)2ei(kx+ωτ)φα(k,ω), (15)

we easily diagonalize Eqs. (13) and (11), and obtain

 S0+Sfℏ=12∑αβ∫dkdω(2π)2φα(k,ω)φβ(−k,−ω)Gαβ(k,ω). (16)

The propagator in the previous equation is given by 111We used the inversion formula for the matrix which reads .

 Gαβ(k,ω)= [δαβvπK(k2+ω2v2+m2)−2Dfπk2δ(ω)]−1 = πK/vk2+ω2v2+m2δαβ+2πK2Dfv2 ×k2δ(ω)(k2+ω2v2+m2)2+O(n). (17)

After performing inverse Fourier transform and setting the number of replicas to zero, for the correlation function (14) we obtain

 Gαβ(x,τ)=δαβG(x,τ)+G0(x), (18)

where

 G(x,τ)=K2K0(m√x2+v2τ2+a2), (19) G0(x)=K2Df4v2e−m|x|m(1−m|x|). (20)

In Eq. (19) we have introduced the parameter as an ultraviolet cutoff, while denotes the modified Bessel function of the second kind.Amit et al. (1980) In the limit of small distances , the correlation functions (19) and (20) become

 G(x,τ)=−3+2δ8ln[c2m2(x2+v2τ2+a2)], (21) G0(0)−G0(x)=K2Df2v2|x|,|x|≪m−1. (22)

Here is the constant, where denotes the Euler constant. In Eq. (21) we have conveniently introduced

 δ=K−3/2. (23)

The parameter of Eq. (23) is introduced as its renormalized value measures the distance from the phase transition between the superfluid and the disordered phase that occurs at moderate interaction correspondingGiamarchi and Schulz (1988) to the renormalized Luttinger liquid parameter . In the superfluid phase, the renormalized parameter corresponding to is positive. The strength of the anharmonic disordered potential (12) and are small parameters in our study.

## Iii Effective action

In this section we derive the renormalization group scaling equations at two-loop order. We use the standard field-theoretical method, which is particularly suitable for calculations beyond lowest order.Zinn-Justin (2002); Amit (1984); Amit et al. (1980) While being somewhat different from the standard Kadanoff-Wilson approach popular in condensed matter physics, it yields equivalent results.Domb and Green (1976); Benettin et al. (1976) The field-theoretical technique is used in a number of earlier studies. To mention the most relevant to our study, Ref. Amit et al., 1980 studied the sine-Gordon model at two-loop order, while Ref. Ristivojevic et al., 2012b considered the random-phase sine-Gordon model at two-loop order. Here we study the quantum version of the random-phase sine-Gordon model.

The main quantity needed to derive the scaling equations is the effective action. For an action , it is defined as , where is the generator of connected correlations defined as Using , we obtain an implicit equation for the effective action Zinn-Justin (2002)

 e−Γ(φ)=∫Dχe−S(φ+χ)ℏ+∫dxχ(x)δΓδφ(x). (24)

Equation (24) applies for an arbitrary action .

We consider the model (9) at weak disorder that enables us to solve the implicit equation (24) using a perturbation theory. It is controlled by the small parameter that determines the strength of the anharmonic term (12). Up to an additive constant, the effective action for the model (9) is given by the expressionRistivojevic et al. (2012b)

 Γ=S0+Sfℏ+Γ1+Γ2+O(D3b), (25)

where

 Γ1= 1ℏ⟨Sb(φ+χ)⟩χ, (26) Γ2= −⟨S2b(φ+χ)⟩χ2ℏ2+12Γ21 +12∑αβ∫dxdτdx′dτ′Gαβ(x−x′,τ−τ′) ×δΓ1δφα(x,τ)δΓ1δφβ(x′,τ′). (27)

Here denotes an average with respect to the quadratic action . For example, the average of the functional is defined as

 ⟨F(φ+χ)⟩χ=∫DχF(φ+χ)e−S0(χ)+Sf(χ)ℏ∫Dχe−S0(χ)+Sf(χ)ℏ. (28)

Now we evaluate the effective action. Using Eq. (26), we easily obtain

 Γ1= −B∑αβ∫dxdτdτ′{[e4G(0,τ−τ′)−1]δαβ+1} ×cos[2φα(x,τ)−2φβ(x,τ′)], (29)

where . We should note that from the correlation function (18) only the diagonal part in replica indices (19) enters the result (III). This is because the off-diagonal part (20) does not depend on the imaginary time. As will become obvious below in this and in the next section, the lowest order term contains all the information necessary to obtain the scaling equations of Giamarchi and Schulz. Giamarchi and Schulz (1987)

The second order term in the expansion of the effective action is given by Eq. (III). After a careful calculation, we express the final result as a sum , where denotes the number of different replica indices contained in , which read

 Γ(1)2= −12B2∑α∫x,τ,τ′,x1,τ1,τ′1f1(x−x1,τ,τ′,τ1,τ′1)cos[2φα(x,τ)−2φα(x,τ′)+2φα(x1,τ1)−2φα(x1,τ′1)], (30) Γ(2)2= −B2∑αβ∑s=±1∫x,τ,τ′,x1,τ1,τ′1{f2(x−x1,τ−τ′1,τ−τ1)cos[2φα(x,τ)+2φα(x1,τ1)−2φα(x1,τ′1)−2φβ(x,τ′)] +12f3(x−x1,τ−τ1,τ′−τ′1,s)cos[2φα(x,τ)−2sφα(x1,τ1)−2φβ(x,τ′)+2sφβ(x1,τ′1)]}, (31) Γ(3)2= −B2∑αβγ∑s=±1∫x,τ,τ′,x1,τ1,τ′1f4(x−x1,τ−τ1,s)cos[2φα(x,τ′)−2φβ(x,τ)+2sφβ(x1,τ1)−2sφγ(x1,τ′1)]. (32)

Here we have used the shorthand notation to denote . The functions , , and in the previous equations are

 f2(x,τ,τ′)= e4G(x,τ)−4G(x,τ′)+4G(0,τ−τ′)−e4G(x,τ) −e−4G(x,τ′)−e4G(0,τ−τ′)+2 +4[e4G(0,τ−τ′)−1][G(x,τ′)−G(x,τ)], (33) f3(x,τ,τ′,s)= [e4sG(x,τ)−1][e4sG(x,τ′)−1], (34) f4(x,τ,s)= e4sG(x,τ)−4sG(x,τ)−1. (35)

The function in [see Eq. (30)] is cumbersome and not necessary for the present study. Namely, the term leads to the correction proportional to in the scaling equation for [see the last term in Eqs. (55)]. This correction is important at three-loop order calculation and thus is not needed here.

The effective action of the model [Eqs. (13), (11), and (12)] up to second order in the anharmonic coupling is given by Eq. (25), where is the sum of the terms in Eqs. (III)-(32). It contains all the information about critical properties of our model at two-loop order. In order to derive the scaling equations one should find the most relevant operators in the effective action. In the limit when the ultraviolet cutoff goes to zero, those operators contain divergent prefactors.

The most relevant contributions from Eq. (III) are

 Γ1= 2Ba1∑α∫dxdτ[∂τφα(x,τ)]2 −B∑αβ∫x,τ,τ′cos[2φα(x,τ)−2φβ(x,τ′)]+… (36)

where

 a1= ∫dττ2[e4G(0,τ)−1], (37)

while the ellipsis denotes many irrelevant operators.

The most relevant terms from are

 Γ2= −2b2B2∑αβ∫dxdτdτ′[∂xφα(x,τ)][∂xφβ(x,τ′)] −2b1B2∑αβ∫dxdτdτ′cos[2φα(x,τ)−2φβ(x,τ′)] +O(B2)∑α∫dxdτ[∂τφα(x,τ)]2+…, (38)

As will be discussed below, the last term Eq. (III) corresponds to a higher order contribution in the scaling equation for . The coefficients and in Eq. (III) read

 b1=∫dxdτdτ′f2(x,τ,τ+τ′), (39) b2=∫dxdτdτ′x2f3(x,τ,τ′,1). (40)

Considering and performing the derivative expansion, one would naively expect to find the operator with two replica indices of the form

 ∑αβ∫dxdx′dτ[∂τφα(x,τ)][∂τφβ(x′,τ)]. (41)

However, the latter operator is not generated in the expanded effective action as a consequence of the equality . The absence of the operator (41) has important consequences for the universality of the exponent of the singe particle correlation function, as we discuss in the following.

The three-replica part (32) contains two contributions that turn out to be unimportant. Namely, for after a gradient expansion a free sum over index delivers a factor , that vanishes in the replica limit. For , one finds the operator of the form where for simplicity we omitted the arguments in replica fields. Such term is irrelevant close to the critical point and has a nondivergent prefactor when .

Collecting the most relevant terms given by Eqs. (III) and (III), the effective action (25) becomes

 Γ= ∑α∫dxdτ{v2πK[(∂xφα)2+m2(φα)2]+[12πKv+2Ba1+O(B2)](∂τφα)2}−[Df4π2+2B2b2] ×∑αβ∫dxdτdτ′[∂xφα(x,τ)][∂xφβ(x,τ′)]−(B+2b1B2)∑αβ∫dxdτdτ′cos[2φα(x,τ)−2φβ(x,τ′)]. (42)

The parameters , , and in Eq. (III) can be evaluated in a series expansion in the small parameter that is defined in Eq. (23). The details of calculation are presented in Appendices A,B, and C, while here we state the final results. In the limit we find

 a1=−2λ+δλ2−4(ln2−1)δλ+O(δ2)+2c12(cmv)3, (43)

where and is the constant defined earlier. In Eq. (43), is a constant. The other two terms are

 b1=π[9λ2+2(51−54ln2)λ+O(δ)+4c2]8v2(cm)3, (44) b2=2πv2(cm)61a[1+O(δ)]+c3v2m5. (45)

The constants , , and are nonuniversal since they depend the choice of the infrared cutoff function in Eq. (13). We will see below that only appears in the renormalization group equations. In general, the appearance of nonuniversal terms in renormalization group equations is not surprising, since other low-dimensional models also contain them. Some examples are the sine-Gordon model Amit et al. (1980) and the Cardy-Ostlund model.Ristivojevic et al. (2012b) However, those nonuniversal terms do not appear in physical observables at the transition. In the present case an example is the single particle correlation function that does not depend on , as shown in Sec. VII.

## Iv Renormalization group equations

Having calculated the effective action (III), we can derive the renormalization group equations. The effective action contains divergent terms in the limit . We absorb them by introducing the renormalized coupling constants that are in the following denoted by the subscript . Characterizing the disorder strength by the dimensionless parameters

 D=πρ22a3Db/v2, (46) Df=aK2Df/v2, (47)

we define the renormalized coupling constants as

 D=ZbDR,δ=Z(3/2+δR)−3/2, (48) Df=ZfDfR,v=ZvR, (49) K=Z(3/2+δR),m=mR. (50)

Here , and are the functions that depend on renormalized parameters and on . Those functions are determined in such a way that the effective action expressed in terms of renormalized parameters is finite order by order in an expansion at small and , in the limit . We find

 Zb= 1−12[(39−54ln2+9c1)DR+2δR]λ +14(9DR+2δ2R)λ2−DR(6c1+c2)+O(DRδR), (51) Z= 1−[3DR+(6ln2−4)DRδR]λ+32DRδRλ2 +3c1DR+2c1DRδR+O(DRδ2R), (52) Zf= (53)

The renormalization group equations are obtained by requiring that the derivatives of the bare coupling constants with respect to the scale nullify. By differentiating the expressions (48) we obtain a system of two coupled linear equations that can be solved in terms of and , yielding

 dDRdℓ=−2DRδR+AD2R+O(D2RδR), (54) dδRdℓ=−9DR+BDRδR+O(D2R). (55)

The constants and are nonuniversal. They depend on that is determined by the choice of the infrared cutoff function [see Eq. (13) for one possibility]. Since this choice is not unique, the form of the correlation function (19) beyond distances of the order of is nonuniversal. Hence the constant in Eq. (43) is also nonuniversal, see Appendix A. However, the sum

 A+B=36ln2−33 (56)

is a universal number that characterizes our model. We discuss below its role.

It is interesting to note that similar equations were derived for the action [see Eqs. (II) and (12)] in a different context, namely how a higher dimension affects the transition.Herbut (2000) This work employs expansion. The dimension plays the role of a regulator for the theory. However, in that case a different universal combination of the parameters appears.

In renormalization group equations (54) and (55) we have determined the first subleading terms, being proportional to and . The term in Eq. (55) is controlled by the last term in Eq. (III) that is beyond two-loop. Setting in Eqs. (54) and (55) and using , we obtain the lowest order scaling equations that correspond to those first obtained by Giamarchi and Schulz. Giamarchi and Schulz (1987, 1988) We note that those references use the parameter instead of our .

The other two renormalization group equations are

 dDfRdℓ=0+O(D2RδR), (57) ddℓ(vRKR)=0. (58)

The last equation turns out to be valid beyond our perturbative calculation. The absence of renormalization to the parameter in our model (7) is an exact result, due to the statistical symmetry , where is an arbitrary function, of the disordered part of the action.Schulz et al. (1988); Hwa and Fisher (1994) Therefore, Eq. (58) is an exact result valid at all orders. The renormalization of forward scattering part of disorder has no consequences for the localization properties of the system.222From one can obtain the equation . It determines very small depletion of the bare value that is proportional to at the transition. We neglect this effect.

Let us analyze the renormalization group equations. We see that the possible fixed points in Eqs. (54) and (55) are (i) , (ii) , and (iii) . The fixed points (iii) are not in the domain of applicability of our calculation as we now show. We assumed that both and are much smaller than unity and thus . Using the condition (56), we obtain . Since measures the disorder strength that must be non-negative, see Eq. (5), the fixed point (iii) is unphysical.

Now we consider the fixed points (i) and (ii). The solution of the flow equations (54) and (55) up to third order in can be expressed as

 δ2R−9DR−AδRDR+2(A+B)δ3R/27=C, (59)

where is an arbitrary constant. The case corresponds to the insulating phase and to the superfluid phase. Thus, the fixed point (i) marks the transition between the superfluid and the insulating phase. The superfluid is characterized by a line of fixed points (ii) where renormalized disorder strength is zero and . The insulating phase cannot be described using the perturbative calculation since the strength of disorder grows with increasing the length scale, becoming higher than unity at a finite scale, when our perturbation theory becomes unapplicable.

Using Eq. (59) we obtain the connection between the parameters and at the critical line determined by the condition :

 DR=(δR/3)2+(2B−A)δ3R/243+O(δ4R). (60)

Substituting this expression in Eq. (55), at large scales we obtain

 δR(ℓ)=1ℓ+(A+B)27lnℓℓ2+O(1ℓ2). (61)

At the critical line, the flow of is universal at large scales at two leading orders, see Eq. (61). The term of the order of depends on the bare parameter . Despite the scaling equations (54) and (55) contain nonuniversal parameters and , their universal combination (56) determines . This has important consequences for the universality of correlation functions, as we discuss below. The flow of the disorder strength is nonuniversal beyond the lowest order, see Eq. (60).

The correlation length close to the critical line from the insulating side () takes the form . Therefore, higher-order corrections do not affect in an essential way and the transition is of Berezinskii-Kosterlitz-Thouless type. Giamarchi and Schulz (1988, 1987) In the superfluid phase, the correlation length is infinite.

## V Generator of connected correlations

In this section we calculate another useful functional that is called the generator of connected correlations, which will be used in the following sections to calculate the correlation functions. Let us consider an action

 S(φ)=S0(φ)+gV(φ), (62)

where the action is expressed as a sum of the quadratic part and the anharmonic part . For the action (62), the generator of connected correlations is defined by the relation Zinn-Justin (2002)

 eW(J)=∫Dφe−S(φ)/ℏ+Jφ, (63)

where denotes an arbitrary source field. The expression (63) can be transformed by introducing a new field , such that . We then obtain , where is the propagator such that we have /2 in the shorthand notation. Equation (63) thus becomes

 eW(J)= ∫Dχexp[−S0(χ)ℏ+S0(GJ)ℏ−gV(GJ+χ)ℏ] = Z0exp[S0(GJ)ℏ]⟨exp[−gV(GJ+χ)ℏ]⟩χ, (64)

where denotes an average with respect to the quadratic action . At small we can use the cumulant expansion, yielding to

 W(J)= lnZ0+S0(GJ)ℏ−gℏ⟨V(GJ+χ)⟩χ+O(g2). (65)

Another expression for that directly follows from the definition (63) is

 W(J)=lnZ+ln⟨exp(Jφ)⟩S, (66)

where is the partition function of the model and denotes the average with respect to the action . Combining the last two expressions for we obtain expressed in powers of the anharmonic coupling .

Now we use the expressions (65) and (66) for our model (9). We should have in mind that the above derivation is written in shorthand notation where we suppressed all internal indices (e.g., replica indices) and is determined by the structure of the theory. For our model (9), the abbreviation denotes

 ∑β∫dx′dτ′Gαβ(x−x′,τ−τ′)Jβ(x′,τ′), (67)

where is defined in Eq. (II). The generator of connected correlations then has the form

 W(J)=lnZ0+W0+W1+O(D2b) (68)

where

 W0(J)= 1ℏ[S0(GJ)+Sf(GJ)], (69) W1(J)= −1ℏ⟨Sb(GJ+χ)⟩χ. (70)

The final expression for the correlation function with an arbitrary source field is

 ⟨⟨e∑α∫dxdτJα(x,τ)φα(x,τ)⟩⟩= Z0Zexp[W0(J)+W1(J)+O(D2b)], (71)

where the average is with respect to the replicated action , see Eq. (9). Note that in the replica limit , where and . In the following sections we consider special choices of the source field in Eq. (V) to evaluate the density-density and the single-particle correlation functions.

Comparing Eqs. (70) and (26) we note similarities between the perturbative expansions of and the effective action .Le Doussal et al. (2013) Either directly calculating or using the this correspondence, we obtain

 W0(J)= 12∑αβ∫dxdτdx′dτ′Jα(x,τ) ×Gαβ(x−x′,τ−τ′)Jβ(x′,τ′), (72) W1(J)= ρ22Dbe−4G(0,0)∑αβ∫dxdτdτ′cos(2Θαβ) ×{[e4G(0,τ−τ′)−1]δαβ+1}, (73) Θαβ= ∫dx′′dτ′′[G(x−x′′,τ−τ′′)Jα(x′′,τ′′) −G(x−x′′,τ′−τ′′)Jβ(x′′,τ′′)]. (74)

## Vi Density-density correlation function

Let us consider the oscillatory part of density-density correlation function. It can be expressed as , where

 R1(