An exactly solvable model of randomly pinned charge density waves in two dimensions

The nature of the interplay between fluctuations and quenched random disorder is a long-standing open problem, particularly in systems with a continuous order parameter. This lack of a full theoretical treatment has been underscored by recent advances in experiments on charge density wave materials. To address this problem, we formulate an exactly solvable model of a two-dimensional randomly pinned incommensurate charge density wave, and use the large-N technique to map out the phase diagram and order parameter correlations. Our approach captures the physics of the Berezinskii–Kosterlitz–Thouless phase transition in the clean limit at large N. We pay particular attention to the roles of thermal fluctuations and quenched random field disorder in destroying long-range order, finding a novel crossover between weakly- and strongly-disordered regimes.


I. INTRODUCTION
It has been well-understood for many decades that below four dimensions random impurities in a material will "pin" the U (1) phase of an incommensurate charge density wave (ICDW) [1].The pinning leads to a frustration that prevents a well-defined CDW state with a uniform phase θ from forming [2].In fact, this destruction of long-range symmetry breaking order by random fields (when not forbidden by gauge invariance) is a general feature in systems with a global U (1) symmetry.The special case of two spatial dimensions is particularly subtle since a state with long-range order is already forbidden at finite temperature by the Mermin-Wagner theorem [3].In the absence of quenched random field disorder, the Berezinskii-Kosterlitz-Thouless (BKT) transition, which involves the proliferation of topological defects (vortices in magnets with planar anisotropy, and dislocations in the case of 2D ICDWs), is special since it separates a (critical) phase with power law correlations and a disordered phase with unbroken U (1) symmetry [4,5].The question of how the critical phase is destroyed by a quenched random field has been understood only qualitatively.
Over the last four decades, a significant amount of work has been done to attempt to answer this question.Initially, advances were made using perturbative field-theoretic renormalization group (RG) techniques [6][7][8].Those authors studied the random field XY model as an idealization of, among other systems, two-dimensional CDWs in systems with charge disorder.They found that the original BKT transition separating the high temperature neutral vortex plasma and low temperature power law correlated phases was destroyed by disorder due to impurities becoming (RG) relevant at an intermediate temperature.However, the nature of the low temperature impurity-dominated state remained undetermined, spurring the application of non-perturbative techniques.For example, the functional RG was used to disprove the long-standing dimensional reduction hypothesis [9][10][11].To this day, the nature of the low temperature state remains hotly debated, with some proposing the ground state to be a so-called "Bragg glass" of dislocations [12,13], while others claim that this state is featureless and lacks any form of order [14].Despite the disagreements, it is now evident that non-perturbative techniques are necessary for shedding light on the physics of the disordered system [15][16][17][18][19].
Interest in disordered ICDWs has not purely been driven by the theoretical challenges described above.The cuprate high-temperature superconductors, which have focused the attention of much of the condensed matter physics community, are strongly layered compounds which display quasi-two-dimensional charge order proximate to the superconducting phase [20][21][22][23][24].The possibility of intertwined [25,26] charge and superconducting orders makes understanding the role of CDW order crucial to the physics of high-T c superconductors.Recently, significant experimental advances have been made in x-ray scattering [27][28][29], scanning electron microscopy and spectroscopy [24,30] and momentum-resolved electron energy loss spectroscopy (M-EELS) [31,32].The increased resolution of measurements has allowed for precise determination of dynamic charge correlations in CDW materials.
Nevertheless, despite the intense theoretical and experimental interest, a full theoretical treatment of the interplay of thermal and quantum fluctuations in the presence of random disorder is an open problem.In this paper, we introduce a model which is exactly solvable in a suitable large-N limit even in the presence of a quenched random field.Our results apply to the specific case of an ICDW in two dimensions which we will take to be unidirectional.Although our main motivation stems from the experimentally observed ICDWs in the lanthanum cuprates [33][34][35], the problem is of much broader interest since such CDWs are seen in many systems, including the dichalcogenides and many other quasi-2D materials.
In this paper we will focus on the classical theory, leaving the quantum theory to another publication.
The order parameter of a unidirectional CDW is the Fourier component of the local charge density ρ(x) at the ordering wave vector Q: ρ(x) = ρ 0 (x) + ρ Q (x) e iQ•x + ρ −Q (x) e −iQ•x + higher harmonics, (1.1) where ρ 0 (x) is the slowly-varying uniform component and the also slowly-varying field is the CDW order parameter.Here we will assume that the ordering wave vector is incommensurate with the underlying lattice and that we have an ICDW.A CDW, commensurate or not and regardless of which microscopic mechanism is responsible for it, is a phase of an electronic system in which ⟨ρ Q (x)⟩ ̸ = 0.This state breaks translation invariance and the point group symmetry of the underlying lattice.States of this type are inherently fragile to disorder since a local (effectively random) electrostatic potential V imp (x) due to charged impurities couples linearly to the local charge density ρ(x) and, consequently, disorder couples linearly to the CDW order parameter.In what follows, we will denote the complex field of the ICDW order parameter as ρ Q (x) ≡ σ(x) and the ordering wave vector Q will be left implicit.
In this work, we present a non-perturbative approach to treating the physics of a U (1) order parameter coupled to quenched disorder.Our method is based on the well-known large-N technique (see, for example, Refs.[36][37][38][39]).The large-N approach has been used to investigate the O(N ) model in a random field by Nie, Tarjus and Kivelson [16], who applied their results to the case of randomly pinned ICDWs.Taken literally, their results apply to dimensions d > 2, since in the absence of a random field the 2D O(N ) model does not have long-range order.For this reason, we were motivated to explore other large-N models that are, in principle, able to capture unique aspects of 2D physics, including the BKT transition (in the absence of disorder), even in the regime where N is large.
To address this problem, in this paper we consider a theory of a two-component generalization of the CP N model with a global U (1) symmetry between the two components.This global symmetry characterizes the order parameter manifold of interest.Our generalized model includes an interaction term between the two components which is solvable in the large-N limit.Here we use the fact that the large-N limit of the CP N model is well understood [40,41], and show that this coupled theory is also solvable in the large-N limit.In the absence of disorder, the N = ∞ ground state we find appears to spontaneously break the U (1) symmetry, but we show that the inclusion of the leading order fluctuations about the N = ∞ ground state is sufficient for a BKT phase transition to emerge.In fact, we show explicitly that the BKT transition appears as a 1/N correction to the N = ∞ transition.
Having confirmed that the clean theory behaves as expected, we then include quenched disorder and solve this model exactly in the N = ∞ limit as a function of disorder and of the coupling strength between the two components.Similarly to the clean theory, the disordered theory appears to have a broken symmetry phase.However, we again demonstrate that allowing for fluctuations about the large-N state necessarily restores the symmetry.The large-N limit of the model predicts a complex phase diagram in the presence of quenched random fields.While in the strong disorder regime we find a phase with short-range order (as expected), in the weak disorder regime the naive large-N limit is the same as in the clean theory.We also derive explicit results for the order parameter correlator in the strong disorder regime (to leading order in the 1/N expansion) and show that it has the same functional form as found in earlier theories of the random field Ising model [42] and O(N ) model with a random field [16,17].Finally, an analysis of the 1/N corrections reveals that the correlation functions have essentially the same behavior in the weak and strong disorder regimes, resolving the naive disagreement with the Imry-Ma argument [2], and implying the existence of a subtle crossover between the two regimes.
The rest of this paper is structured as follows: In Sec.II we present the model we will be using in this work.We demonstrate our method for studying a U where a µ is a fluctuating U (1) gauge field and D µ [a] = ∂ µ + ia µ .Importantly, we take both z and w to be coupled to a µ , as opposed to two independent fluctuating gauge fields.
This model has a formal similarity with the chiral Gross-Neveu model which has a wellunderstood large-N limit [43,44].It will be convenient to work in a representation where the U (1) × U (1) ⊂ G subgroup is generated by two types of transformations: Note that the quartic interaction term is invariant under the global relative U (1) symmetry (ii) which can be spontaneously broken (with the usual caveats on the restrictions imposed by the Mermin-Wagner theorem).

B. Large-N Solution
We will now demonstrate that this model has the desired properties.The partition function for the theory is where λ 1 and λ 2 are Lagrange multiplier fields imposing the constraints |z| 2 = |w| 2 = 1.
Additionally, the quartic interaction term |z * • w| 2 can be decoupled using a Hubbard-Stratonovich (HS) transformation: Since the Gaussian integral over σ is peaked around the value σ max = (K/g)z • w * , we identify the complex field σ as an emergent U (1) order parameter.In fact, invariance of the right hand side of Eq. (2.4) under a relative U (1) transformation imposes the transformation rule σ → e 2iϕ σ.After the HS transformation, the functional integrals over z and w become Gaussian, and hence, can be performed exactly to obtain where we have defined g = g 0 /N and K = K 0 /N .As N → ∞, the partition function can be evaluated exactly within mean field theory.We then take the spatially uniform ansatz σ(x) = σ 0 , λ 1 (x) = λ 2 (x) = m 2 (since the theory also has a Z 2 symmetry which exchanges the z and w fields) and a µ (x) = 0 (by gauge invariance).Then, where V is the volume of d-dimensional space and Λ denotes some regulator; e.g., a UV cutoff.
We will now specialize to the case d = 2.The effective potential (2.6b) is logarithmically ultraviolet divergent, but renormalization of the coupling constant g 0 suffices to cure this.
for some renormalization scale µ, the regulator Λ can then be removed from Eq. (2.6b), yielding a renormalized effective potential: The ground state values of m 2 and |σ 0 | are then obtained by minimizing the potential.First, Next, it is straightforward to check that for K 0 < K c = 4πµ 2 e −4π/g R , the potential is minimized by a vanishing order parameter σ 0 = 0.However, for K 0 > K c , the magnitude of σ 0 is determined by a Curie-Weiss type transcendental equation where |σ| = |σ 0 |/(µ 2 e −4π/g R ), which is well-known to yield a β = 1/2 critical exponent in the limit |σ| ≃ 0. Since the effective potential Eq. (2.8) is invariant under the relative U (1) symmetry, a mean-field solution σ 0 spontaneously breaks this symmetry.However, it is evident that in d = 2 this spontaneous symmetry breaking (SSB) is an artifact of the large-N limit, since the Mermin-Wagner theorem forbids SSB [3].This signals the need to include fluctuations around the mean-field limit, for which the large-N technique provides a systematic method.

C. Fluctuations & BKT Transition
Since we are primarily interested in the infrared physics of this model, we can focus on the Goldstone manifold of degenerate ground states parameterized by the U (1) phase θ of the complex order parameter field; σ(x) ≡ ρ(x)e iθ(x) , where ρ(x) is the (real) amplitude of the order parameter (not to be confused with the charge density).By freezing ρ(x) = ρ 0 and allowing θ(x) to vary slowly in space, the leading corrections to the effective action will be determined by a gradient expansion.However, θ ∈ [0, 2π) is a compactified boson, and its periodicity allows the existence of vortices [4,5].
Following the usual analysis (see, for example, Ref. [45]), the periodicity is imposed by a non-fluctuating source of vorticity A µ which satisfies for some configuration {m j } of vortices with topological charge m j ∈ Z.This source is then minimally coupled to the relative U (1) symmetry, so that the partition function Eq. ( becomes where the sum is over all possible configurations of vortices.Expanding to quadratic order in θ, A µ , and their derivatives, yields where is the effective coupling constant for the "electrodynamic" response of A µ , γ θ is the phase stiffness and ρ = ρ 0 /(µ 2 e −4π/g R ); see Appendix A for the derivation of these quantities.Note that σ has U (1) charge q = 2 (not to be confused with the topological charge of a vortex) since it is a charged composite of q = ±1 fields.
FIG. 1. Phase diagram corresponding to the action Eq.(2.1).For K 0 > K BKT dislocations are confined, leading to quasi-long-range order and power law correlations with exponent α > 0.
For K c < K 0 < K BKT there is a condensate but dislocations proliferate; correlations of σ decay exponentially with length ξ 1 .For K 0 < K c no U (1) condensate forms, so correlations of σ decay exponentially with different length ξ 2 .K BKT is split from K c to order 1/N .
This result implies a critical value of the phase stiffness γ BKT = 2/(πq 2 ) = 1/2π such that vortices/dislocations proliferate for γ θ < γ BKT (see Appendix A).Since γ BKT is of order 1/N relative to generic values of γ θ (ρ), for large enough N we can safely approximate , where ρBKT is the value of ρ which solves the equation γ BKT = 1/2π.Substituting ρBKT ≃ 6/ √ N into the saddle point equation Eq. (2.10) and expanding to first order in 1/N yields one of our first main results: That is, we have explicitly shown how the BKT transition emerges at order 1/N relative to the mean-field transition; K c remains the point at which the amplitude of the condensate forms, but only exponentially-short range order exists for K c < K < K BKT .The fluctuationdominated regime between K c and K BKT is narrow in the large-N regime and it is expected to become broad for smaller values of N .The existence of this fluctuational regime is a feature of two-dimensional physics; this behavior is analogous to the situation in (quasi-)two-dimensional superconductors, where the condensate develops at a higher temperature than the onset of zero resistivity.The results of this section are summarized in the phase diagram in Fig. 1.Finally, we note that the apparent contradiction with the Mermin-Wagner theorem is resolved at this order in the fluctuations, since the critical Goldstone phase K > K BKT also has no long-range order, though with power-law decaying correlations.This is unsurprising given that the Mermin-Wagner theorem is fundamentally a statement of the importance of fluctuations in low dimensions.This behavior is closely analogous to what was found long ago in the chiral Gross-Neveu model, though, importantly, the phase stiffness in that model is a pure number that cannot be tuned [44].

III. THE ROLE OF DISORDER A. Coupling to Quenched Disorder
Having established that our model reproduces the salient features of a U (1) order parameter within the large-N limit, we turn to the primary focus of this work: disorder.Any quenched disorder must couple only to gauge invariant combinations of fluctuating fields.
With our model, we are spoilt for choice with possibilities, not all of which are physically interesting.For example, the simplest option, a complex scalar disorder field coupled linearly to the emergent U (1) order parameter z * • w, can be easily shown to produce a trivial large-N limit since the dimension of the disorder field does not scale with N .As such, it is useful to draw inspiration from the conventional CP N model, the large-N solution of which is presented in Appendix C. We find that it is natural to consider "adjoint disorder" where, for simplicity and without loss of generality, we take the disorder z a (x) and w a (x) to be N 2 component real vectors in the adjoint representation of U (N ), with generators τ a satisfying τ a αβ τ a γδ = N δ αδ δ βγ (implied summation over repeated indices), and distributed with variance η 2 1 according to where overlines denote averaging over disorder configurations.Note that cross-correlations between z a (x) and w a (x) are generally allowed by symmetry, but we do not observe any qualitative impact as a result of this added model complexity.A general approach, however, should also include the new gauge invariant bilinear where h a (x) is a N 2 component complex random vector with variance η 2 2 :  4).In this way, this second form of disorder has the same U (1) symmetry properties as a disorder field coupled linearly to z * • w.We will see in a later section that symmetry-breaking disorder is vital to the physics of fluctuations, while the neutral disorder is important for stabilizing the ground state of the theory.Other forms of disorder are certainly possible, but the adjoint disorder presented here provides useful mathematical simplifications, and, we believe, is the most physically motivated and natural choice leading to a non-trivial large-N limit.
Since we are interested in thermodynamic observables that are independent of any specific realization of the disorder, we use the replica trick: where S j are the replicas of the original action Eq.(2.1).We are excluding the terms with i = j from the sum in the last line because the unit vector constraints them trivial constants.Note that the novelty of quenched disorder in CP N -type models is the generation of quartic inter-replica interaction terms, reminiscent of the Sherrington-Kirkpatrick model of spin glasses [46].We will see that this produces fundamentally different large-N solutions compared to O(N ) models (see Appendix B).

B. Large-N Solution
Next, we must take some care with the Hubbard-Stratonovich transformation to avoid over-counting the effective degrees of freedom in the theory.We start by introducing the collective coordinates ζ ij and ω ij for i ̸ = j and fix ζ ii = ω ii = 0 for the reason discussed just above.Since this expression is quadratic in ζ ij and ω ij , the Lagrange multipliers κ ij and ψ ij can be eliminated exactly by the saddle point equations Therefore, the interaction term can be expressed as after rescaling ζ ij and ω ij by 1/N .The z and w fields can then be integrated out of the partition function to obtain the effective action where we have rescaled (ζ, ω) → g 0 (ζ, ω), and defined g = g 0 /N , K = K 0 /N , and η 1,2 = η1,2 /g 0 to obtain a well-defined large-N limit.Tr( • ) includes the functional operator trace as well as the trace over replica indices, diag( • ) denotes a matrix which is diagonal in replica indices, and ζ and ω are the matrices with elements ζ ij and ω ij , respectively.We note that this effective action is positive definite as long as η1 > η2 .While the disorder averaged action Eq.(3.5) is bounded from below regardless of the disorder strengths, maintaining η1 > η2 is necessary if one wishes to avoid spontaneously breaking the replica permutation symmetry, which is beyond the scope of this work.Then, observe that the HS disorder fields inherit the transformation rules show that the replica-symmetric configuration is stable in the regime of interest η1 > η2 .
We then take the spatially uniform and replica-diagonal ansatz 2 ) (suppressing tildes from hereon for notational clarity).Next, one can simplify the block matrix structure in the effective action to obtain the replicated effective potential where Î is the n × n identity matrix and M is the n × n matrix of ones.After the same renormalization procedure as in all the previous cases, it follows that the disorder-averaged physical effective potential is (3.12) In d = 2 this evaluates to The full multivariate structure is complicated, so we first consider the case σ 0 = 0.This reduces the problem to two decoupled CP N models with disorder strength η 2 1 +η 2 2 ≡ η 2 tot (see Appendix C); the lack of any dependence on the relative magnitudes of η 1 and η 2 is a result of the replica-symmetric ansatz.The saddle point equations for ω 0 and m 2 are, respectively, However, there is an apparent problem with these equations, as they seem to imply that ω 0 < 0 for η 2 tot < η 2 c,0 = 4πµ 2 e −4π/g R ; if this were the case, the scalar field propagator in Eq. (3.9) would not be positive-definite.However, in this problem regime, the configuration with ω 0 = 0 is more energetically favorable, being the global extremum (maximum, because of the effectively negative number of degrees of freedom of ζ ij and ω ij when n → 0) [47].This implies that η c,0 is a crossover scale which divides a "clean" or weakly-disordered regime from a strongly-disordered regime: Observe that m 2 and ω 0 are continuous at η c,0 , and that for η 0 < η c,0 , m = µe −2π/g R has the same value as in the theory without disorder.
For σ 0 ̸ = 0, the full saddle point equations cannot be solved explicitly, so we simply state them here for completeness: ) ) The phase boundary for the onset of a condensate density |σ 0 | follows from the non-trivial solution to Eq. (3.16c): where K c,0 = η 2 c,0 = 4πµ 2 e −4π/g R is the critical coupling for the mean-field U (1) transition in the clean limit (η tot = 0).Therefore, at the mean-field level, the disorder only shifts the boundary for the onset of the condensate; a larger coupling K 0 is required to overcome increasing amounts of disorder.Importantly, however, the disorder does not eliminate the condensate (at least at the mean-field level).Similarly, the boundary for the disorder-driven FIG. 2. N = ∞ phase diagram of the theory with quenched disorder.ω 0 = 0 everywhere to the left of the dashed line.σ 0 = 0 everywhere below the solid blue line.K c,0 is the critical coupling in the clean limit and η c,0 is the critical disorder in the absence of condensate.
crossover is shifted at finite condensate density: where σ = σ 0 /(µ 2 e −4π/g R ), and since ω 0 = 0 on the phase boundary, σ(K 0 ) is determined by the clean saddle point equation Eq. (2.10).The above results are summarized in the phase diagram in Fig. 2.
Just as the large-N solution in the absence of disorder appeared to violate the Mermin-Wagner theorem, our present solution appears to violate the Imry-Ma condition [2]; in dimensions d < 4, any infinitesimal amount of disorder should not only destroy long-range order, but correlations must decay at least exponentially; i.e., the existence of a Goldstone phase with power law correlations must also be eliminated.Note that a condensate of ω 0 does not break any physical symmetries of the non-replicated theory, and hence, is not problematic.In the following section, we will show that in the dirty regime ω 0 > 0 the contradiction with Imry-Ma is directly resolved by the inclusion of fluctuations to leading order in N .This will mirror the resolution of the Mermin-Wagner theorem in the clean theory.
Next, we will demonstrate that sub-leading order fluctuations in the clean regime ω 0 = 0 must be included, and generate qualitatively the same interactions which are responsible for the destruction of long-range order in the dirty regime.

Correlations in the Symmetric Phase
We begin with the case σ 0 = 0 and ω 0 > 0. The physical observable of most interest in this regime is the correlation function of the order parameter.The natural parameterization here is the Cartesian form σ j (x) = σ R j (x) + iσ I j (x).Expanding the effective action Eq.(3.9) to quadratic order in the σ R and σ I yields j (y), (3.19)where the full integral expressions for the kernels Π (1) (x−y; n) and Π (2) (x−y; n) are given in Appendix D. It follows directly from the matrix structure of this action that the propagator and hence, that the disorder-averaged propagator is σ (p; 0) + Π (2) (3.21) Immediately, we see that the Imry-Ma condition holds: if Π σ (0; 0) ̸ = 0, the replica-diagonal kernel Π (1) σ (p; 0) must be gapped for all d < 4, barring any pathological momentum dependence.Specifically, in d = 2 we find for small momentum Π (1)  σ where the effective mass and (static) disorder strength are Therefore, the propagator has the well-known [42,48] double-Lorentzian form; any approximation of the double-Lorentzian term in powers of momentum should at least respect the double-pole structure.A potentially surprising feature of this result, compared to the momentum-independent numerator of the propagator in the O(N ) model, is the "disorder kernel" Π (2) (p).We note that this momentum dependence is simply a consequence of the fact that σ i is a bound state of the z i and w i , and not of the composite nature of the disorder fields ζ ij and ω ij , since to this order in N fluctuations of σ do not couple to fluctuations of the disorder fields.

Nature of the Symmetry-Breaking Phase
Here we consider the regime σ 0 ̸ = 0 and ω 0 > 0. In the absence of explicit symmetry breaking (η 2 disorder), the low energy degrees of freedom of the effective action Eq.(3.9) are the n gapless Goldstone modes corresponding to long-wavelength distortions of the replicated relative symmetry phases θ i (x) where the phases θ i are scaled so as to correspond to the replicated order parameter phase.
With the U (1) n global symmetry explicitly broken down to the replica-diagonal U (1) subgroup, n−1 of these modes will become gapped.However, as long as the symmetry breaking disorder is sufficiently weak, other gapped modes will remain frozen out at comparatively higher energies.Given these considerations, the infrared sector of the effective action in where the three parameters γ θ (n), Γ θ (n) and m θ (n) depend explicitly on the number n of replicas.The derivation of this effective action and the full (rather uninstructive) expressions for the three parameters are given in Appendix D. Note that the m θ term also naively appears to break the periodicity of the θ j .This is forbidden, since Eq.(3.9) is invariant under shifts θ j → θ j + 2π, so the constant term must actually be the total contribution to quadratic order from terms in the "true" IR effective theory of the form where p ∈ N are the degrees of p-fold anisotropy and L 0 is the part of the theory which remains fully U (1) n invariant [the gradient terms in Eq. (3.26)].Our conventional 1/N expansion does not give us direct access to the coefficients α p , but it is clear that the full effective theory should respect the periodicity of the θ i .
We will now unpack this result.The spectrum of the inverse propagator is ) where the φ i are orthogonal eigenvectors of the kernel in the basis of the θ i , and their normalization is chosen to preserve the compactification radius (L 1 norm).The single remaining gapless mode φ 0 corresponds to in-phase fluctuations of all the replicas; the gapped modes correspond to the n − 1 linearly independent out-of-phase motions, where, to leading order in the replica number n, and ).Since m 2 θ > 0 for η 1 > η 2 , this result confirms that the ground state obtained in the large-N limit is stable.
Next, we must address the nature of physical observables in this dirty regime and compatibility with Imry-Ma.It is useful to consider an effective disorder-averaged phase stiffness.
To this end, we first project Eq. (3.26) into the "deep IR" scale below the mass gap m θ , We again emphasize that the factor of n is needed to preserve the compactification radius on changing basis from θ i to φ 0 .Next, we impose a uniform static twist φ 0 (x) = Q • x.Finally, the phase stiffness should be identified with the usual thermodynamic helicity modulus, determined from the disorder-averaged effective potential in the presence of the twist: We could repeat our analysis from the clean theory to determine what this expression predicts for the evolution of the BKT transition as a function of disorder.However, it is well-understood that γ θ controls the relevance (in the renormalization group sense) of the disorder-induced p-fold anisotropy Eq. (3.27) [6,8].Since the theory can be tuned to make ω 0 small, we can treat disorder as a perturbation to the sine-Gordon representation of the clean theory [see Eq. (A.9)].Then, the results of Ref. [6] imply that order-p anisotropy is a relevant perturbation when γ θ > p 2 /16π.Therefore, the BKT transition at γ θ = 1/2π is actually preempted by random field (p = 1) and random bond (p = 2) anisotropy.Unfortunately, our large-N analysis does not reveal the nature of the disorder-dominated phase.
However, it is clear that the conventional Goldstone phase with power law correlations cannot survive, and the BKT vortex plasma is known to have exponentially-decaying correlations.
To determine the boundary of the disorder dominated phase K dis (η tot ), we begin by using the saddle point equations Eqs.(3.16), to write γ θ as a function of only the dimensionless parameters ρ = ρ 0 /(µ 2 e −4π/g R ) and η tot /η c,0 .Expanding the resulting expression to quadratic order in ρ yields Similarly, the saddle point equation Eq. (3.16c) can be expressed entirely in terms of ρ, η tot /η c,0 and K 0 /K c,0 .Expanding the result to quadratic order in ρ yields Solving γ dis = p 2 /16π for ρdis and substituting into the saddle point equation yields Predictably, for η tot ≫ η c,0 the vortex-dominated BKT phase eventually vanishes, though only logarithmically slowly.The surprising feature of this result is its non-monotonicity; the vortex phase is actually slightly extended to a maximum value K dis (η max )/K c (η max ) ≈ (1 + 0.17p 2 /N ) at η max /η c,0 ≈ 1.2; see Fig. 4. Note that this result is not an artifact of the expansion in powers of ρ as it can be confirmed by solving the equation numerically.It may seem counter-intuitive for small amounts of disorder to extend the BKT phase relative to the impurity-dominated phase.However, this simply reflects the fact that the disorder initially has a stronger effect on the condensate density than on the phase stiffness.

D. Weak Disorder Regime
In the previous section, the primary role of disorder was to mediate interactions between replicated fields via the N = ∞ z and w field propagator Ĝ−1 Naively, it might appear as though the system is blind to disorder in the regimes where ω 0 = 0.In this section, we will show that this is not the case and that fluctuations which are sub-leading in 1/N generate qualitatively the same inter-replica couplings as in the strong disorder regime.To demonstrate this, it will suffice to consider the case where σ 0 = 0.
In the strong disorder regime, the propagator had the double-Lorentzian form characteristic of disordered systems due to the replica-mixing kernel Π (2) σ (p, n).If we can show that this kernel is actually non-zero even when ω 0 = 0, then the same conclusions as before will hold.We can write the kernel extremely generally in the form where a, b, c, d = 1, 2 denote either z or w, respectively, and i, k 1 , k 2 , ℓ 1 , ℓ 2 are replica indices, with implied summation over repeated indices.G ab ij (p) is the exact z, w propagator, and Γ ab i,jk (p, q, k) is, similarly, the exact three-point vertex between σ i and the z and w fields with replica indices j and k; within the large-N method, both these quantities have a perturbative expansion in powers of 1/N .In the N = ∞ limit, that is, at tree level, the propagator is simply the bare propagator Eq. (3.35), and the three-point vertex is independent of momentum and completely replica-diagonal.At this level, the necessary replica off-diagonal terms only appear when ω 0 > 0. However, within the 1/N expansion, both the propagator and the vertex functions receive corrections due to fluctuations of all the collective fields.For example, when σ 0 = 0, the propagator can acquire off-diagonal components which couple the z i and w i fields together due to fluctuations of the σ i .Similarly, when ω 0 = 0, the three-point vertex can acquire contributions which are off-diagonal in the replica indices due to fluctuations of the ζ ij and ω ij .At least at one-loop level, the propagator cannot acquire replica off-diagonal components when ω 0 = 0 (see Fig. 3).A quantitatively correct calculation of the 1/N corrections to physical observables, such as the order parameter correlation function, requires the inclusion of self energy and vertex corrections from every fluctuating field.This is not a particularly instructive calculation, so we will content ourselves with examining the replica-mixing contribution to the three-point vertex, as shown in Fig. 3. To leading order, this involves the order-N contribution to the propagator for the ζ ij and ω ij fields.Expanding Eq. (3.9) to quadratic order in these fields yields where the matrix of kernels in momentum space is (recall that we are setting σ 0 = 0), and there is no dependence on n when ω 0 = 0.This yields a contribution to the three-point function Next, observe that the leading 1/N contribution in an expansion of Eq. (3.36) comes from the two ways of inserting a single copy of the replica-mixing contribution above, Π As noted in the previous section, when deriving the double-Lorentzian propagator we are primarily interested in Π σ,ij (p = 0), and since Γ 12 i,jj (0, −q, q) → 0 as |q| → ∞, the kernel is largely determined by the value of the vertex function at q = 0, which we find to be This expression was obtained by using a Padé approximant of order (0,2) for the (ζ, ω) propagator in Eq. (3.40).This result has two important features: i) It vanishes for ∆ = 1, that is, η 2 = 0. Therefore, some symmetry breaking disorder is necessary to generate replicamixing interactions.ii) It diverges at the phase boundary η 2 c,0 = 4πm 2 (the dependence of m 2 on model parameters will be shifted to order 1/N even if the phase boundary equation appears unchanged in this form).This divergence originates from the infrared divergence in Eq. (3.40) due to gapless ζ ij and ω ij fluctuations.However, a similar IR divergence must occur in the calculation of the replica-diagonal kernel Π (1) σ (p).Therefore, the effective disorder strength will not be infinite.FIG. 4. Schematic phase diagram of the theory with quenched disorder, including the effects of fluctuations.For K 0 < K c there is a featureless symmetric phase (σ 0 = 0) with double-Lorentzian correlations Eq. (3.23).For K c < K 0 < K dis the ground state is a BKT-type vortex plasma (σ 0 ̸ = 0); disorder is RG irrelevant; this region is narrow when N is large.For K 0 > K dis disorder is RG relevant and the (unknown) ground state is dominated by impurities.The light gray dashed line represents the mean-field disorder-driven transition, which does not affect the nature of the ground state.
Finally, the result of this section shows that the order parameter correlation function will still have a double-Lorentzian form in the weak disorder regime ω 0 = 0, albeit with modified parameters m ′ 2 σ and η ′ 2 σ .It is also clear that this analysis carries over to the σ 0 ̸ = 0 regime.Therefore, whether ω 0 is zero or not does not fundamentally alter the nature of the ground state, but simply determines whether the effect of disorder is suppressed by an additional factor of 1/N .All of the previous analysis is summarized in the phase diagram shown in Fig. 4. In closing this section we re-emphasize that the presence of the double-Lorentzian term in the disorder-averaged two-point correlator of Eq. (3.43) implies a stronger infrared singularity in the averaged susceptibility.In turn, this behavior implies the absence of long-range order below four dimensions even in the weak disorder regime, as expected from the Imry-Ma argument [2].

IV. DISCUSSION
In this paper we have introduced a new model, solvable in the large-N limit, to understand the interplay of thermal fluctuations and quenched random disorder in two-dimensional systems with a phase transition in the clean limit.While we were motivated by the existence of charge-ordered states in the cuprate high-T c superconductors, the formalism we have developed is completely general and can, in principle, be applied to any system described by a U (1) order parameter coupled to random field disorder.Previous studies [16] using the large-N technique did not capture the physics of the Berezinskii-Kosterlitz-Thouless phase transition because the order parameter was encoded on a manifold with dimension that scaled with N , which, for any N > 2, does not have a phase transition in 2D in the clean limit.In contrast, we held the U (1) manifold fixed as a subgroup of the larger symmetry group U (N ) × U (N ).Naively, this produced an inconsistency with the Mermin-Wagner theorem [3] in the absence of disorder, and the Imry-Ma theorem [2] in the presence of random field disorder.Throughout this paper, we have shown that including fluctuations about the N = ∞ ground state resolves any apparent contradictions.
Our first main result demonstrated how these fluctuations produce a BKT transition split from the mean-field transition in the clean theory at order 1/N ; the vortex plasma phase is narrow when N is large.Next, we derived the N = ∞ phase diagram of the model as a function of the disorder strength and the coupling between the two CP N components.This revealed a novel disorder-driven transition between a weakly-and a strongly-disordered regime.We then derived an explicit expression for the order parameter correlation function in the strongly-disordered symmetric phase of the theory, finding agreement with the wellknown double-Lorentzian distribution [42,48].In the strongly-disordered phase with a U (1) condensate, we derived an infrared effective theory which had the form of a random field XY model.By mapping our expression to known results [6,8], we were able to derive the phase boundary between a vortex-dominated BKT-like phase and the contentious [12][13][14] impurity-dominated phase.
Finally, we proved that fluctuations are sufficient for ensuring agreement with the Imry-Ma condition in the weakly-disordered regime by showing that the order parameter correlation function again has the double-Lorentzian form.We note that non-analytic contributions that are non-perturbative in the 1/N expansion are not included in our analysis.These include rare configurations of the disorder [49], which may also play a role, for example, by rounding the disorder-driven transition into a broad crossover.
In this paper, we have considered a purely classical theory.This suitably describes static correlations in materials when quantum effects are weak, such as when electronelectron interactions are strong enough to form an insulating CDW state.A complete theory of dynamic correlations in randomly pinned ICDWs necessarily requires the inclusion of quantum fluctuations.Luckily, the approach used in this work lends itself well to such a generalization.Investigating the N = ∞ ground state properties at T = 0 will require little more than changing the dimension from d = 2 to d = 3 in all of the calculations presented, due to the usual quantum-classical correspondence.However, being quenched, disorder also introduces significant temporal non-locality into the system, which will lead to richer behavior of dynamic order parameter correlations.It would also be interesting to study conducting systems, in which case the CDW order parameter will experience Landau damping [50].This would provide a unified theoretical picture of the role of disorder in a wide range of physical systems, and help answer many of the questions raised by modern experiments.In this appendix, we sketch the calculation of the response of the two-component CP N model to a U (1) background gauge field.We start from the effective action with a background field A µ minimally coupled to the (previously global) relative U (1) symmetry The leading order behavior of the sector corresponding to the relative U (1) symmetry involves only θ and A µ .All other couplings are either forbidden by symmetry or sub-leading in 1/N .Expanding the tr ln to quadratic order in these fields yields the following two-point kernels: where Π µν aa (p) is (two times) the well-known electrodynamic response of the CP N model [40,41], which in d = 2 is Since A µ couples to matter (θ), the kernel Π µν AA will have both transverse Π µν T and longitudinal Π µν L responses.One finds that where the kernel Λ(p 2 ) can be evaluated exactly in d = 2, where γ θ is the phase stiffness The transverse vorticity response is Note that since γ θ ∝ ρ 2 0 in the vicinity of the critical point K c , the transverse response does not vanish.Instead, it simply becomes equal to Π µν aa (p) for K 0 ≤ K c .Therefore, in the long-wavelength limit, we have where e 2 = ρ 2 0 /(2m 2 γ θ ) is the effective coupling constant of the probe field and In the above discussion we ignored the fact that the phase field is actually defined mod 2π and that the full computation of the partition function is dominated by the contributions of vortices and anti-vortices.Such contributions can be computed by regarding the local flux of the gauge field A µ as representing vortices and anti-vortices (for a recent discussion see Ref. [51]).It is well understood that the leading contributions to the partition function come from dilute configurations of vortices and anti-vortices, and that in terms of the Cauchy-Riemann dual ϑ of the phase field, which satisfies ∂ µ ϑ = ε µν ∂ ν θ, the effective action is mapped to the sine-Gordon theory [52][53][54].In the present case, where β 2 0 = (4π) 2 γ θ , v 0 = 2e −Ecore /a 2 , E core is the core energy of a vortex, a ∼ µ −1 is an ultraviolet cutoff, and ϑ has been rescaled by 2 √ γ θ to bring the action into canonical form.
Appendix B: Review of O(N ) Model with a Quenched Random Field In this appendix, we review the large-N analysis of O(N ) models in a quenched random field for the purpose of aiding comparison with the main results in this paper.The O(N ) nonlinear sigma model (NLSM) is described by the action where h(x) is a fixed external source.Suppose now that h(x) = h(x) is a random field drawn locally from a Gaussian distribution such that where the overline denotes an average over the disorder configurations, and the variance η 2 represents the strength of the disorder.Using the replica trick formalism, we consider, for where n j and λ j (the Lagrange multiplier imposing the unit vector constraint), for j = 1, . . ., n are the replicated fields corresponding to each factor of Z[h].Because the theory remains quadratic in the replicated O(N ) fields n j , they can be integrated out exactly: where Tr( • ) denotes the functional operator trace as well as the trace over replica indices, and Î is the n × n identity matrix, M is the matrix with a 1 in every entry and diag(λ) = diag(λ 1 , . . ., λ n ) is the diagonal matrix of Lagrange multipliers.The remaining functional integrals over each λ j can be performed using steepest descent, which becomes exact in the limit N → ∞.To make this limit precise, we define g = g 0 /N and η 2 = N η 2 0 , keeping g 0 and η 0 fixed as N → ∞.We then look for a replica-symmetric saddle point where where tr( • ) denotes a trace over only the replica indices.Since M is an idempotent matrix ( M 2 = n M ), finding the inverse matrix in the integrand is straightforward: and hence, taking the replica limit n → 0, While explicitly solving this equation for m 2 requires renormalization of the coupling constant g 0 , it will suffice for our purposes to simply make the following remarks: (i) For any amount of disorder (η 2 0 > 0) in d ≤ 4 there is an infrared singularity in the above integral equation unless m 2 > 0. This implies the absence of any broken symmetry phase, providing an exact, non-perturbative realization of the Imry-Ma argument [2].(ii) The disorder-averaged propagator of the n field has the double-Lorentzian form Appendix C: Review of CP N Model with a Quenched Random Field In the main body of this work, we investigate a two-component generalization of the CP N model.Here, we present a summary of the physics of the simpler model in a quenched random field for pedagogical purposes.
In the CP N model, quenched disorder can only couple to gauge invariant combinations of the scalar field.The natural coupling originates in the Hopf map from O(3) unit vectors to elements of CP 1 , h • n → h a z * α τ a αβ z β , where τ a are the generators of SU (2).Therefore, we consider the partition function where τ a are the generators of SU (N ), h a is a real (N 2 − 1)-component random field transforming under the adjoint (vector) representation of SU (N ) and drawn from the locally Gaussian distribution η 2 is the variance (strength) of the disorder and overlines denote averaging with respect to disorder configurations.For quenched disorder averages, we use the replica trick formalism, where we have used the identity τ a αβ τ a γδ = N δ αδ δ βγ −δ αβ δ γδ [implied summation over repeated indices; working in the convention tr τ a τ b = N δ ab ], and dropped the unimportant constant terms proportional to |z i ||z j | = 1.The quartic interaction between the replicated z j fields can then be decoupled using a Hermitian Hubbard-Stratonovich field ω ij = ω * ji ; since z is a unit vector, we also have ω ii = 0. Therefore, The replicated theory has an enlarged U (1) n gauge symmetry which allows for independent gauge transformations on each replicated field z j (x) −→ e iθ j (x) z j (x), a µ j (x) −→ a µ j (x) − ∂ µ θ j (x).(C.7) As a consequence, the disorder Hubbard-Stratonovich field inherits the transformation ω jk (x) −→ e i(θ j (x)−θ k (x)) ω jk (x), (C.8) transforming as a tensor under U (1) n .Since a replica-symmetric ground state should be invariant under the permutation group S n , we might have assumed this would uniquely constrain ω ij = ω 0 , independent of i, j.However, the manifold of gauge-equivalent ground states is U (1) n /S n ; i.e., ω ij = ω 0 is simply a particular representative of the equivalence class of configurations obtained from ω 0 by gauge transformations.
Given the above considerations, we are free to expand around the replica-symmetric saddle point ω ij (x) = ω 0 .Also taking the replica-symmetric ansatz λ j (x) = m 2 , and fixing a µ j (x) = 0, the replicated effective potential in the N = ∞ limit is 2π) d ln det (q 2 + m 2 + ω 0 ) Î − ω 0 M + n(n − 1) where Î is the n × n identity matrix and M is the matrix with a 1 in every entry.Using the same coupling constant renormalization as in Eq. (2.7) to cure the UV divergence, we recover the physical (renormalized) effective potential by taking the replica limit In d = 2 this evaluates to The values of m 2 and ω 0 are then obtained from the saddle point equations, though a subtlety occurs for ω 0 , as the corresponding equation has two solutions.Note that since ω ij has n(n − 1) complex degrees of freedom, in the replica limit n → 0, the negative number of effective degrees of freedom require us to maximize the effective potential [47].This yields a crossover scale η 2 c = 4πµ 2 e −4π/g R such that We emphasize that this crossover does not break any "physical" symmetries; i.e., those of the non-replicated theory.It is also likely that this apparent mean-field transition would be rounded by [possibly non-perturbative O(e −N )] corrections to the saddle point equations.
Note that this behavior is not an artifact of the absence of a phase transition in the clean d = 2 theory; consider the case of d = 3, where Eq. (C.10) evaluates to It is simple to check that the "clean" symmetry breaking state m 2 = ω 0 = 0 is also never energetically favorable.
For the purposes of this paper, it suffices to note that the 1/N expansion reveals that the Higgs mechanism takes place for η 0 > η c , with n − 1 of the replicated gauge fields a µ j becoming massive; as expected from the form of the disorder, only the replica-diagonal U (1) subgroup is unaffected, with the corresponding gauge field remaining gapless.
Next, the case σ 0 ̸ = 0. Following the same notation as above, we write the effective action in the form θ (x − y; n) Î − Π (2π) 2 G ii zw (q)G ii zw (p + q) − G ii zz (q)G ii zz (p + q) + 2ρ 0 d 2 q (2π) 2 G ii zw (q) + (n − 1)ω 2 0 ∆ 2 d 2 q (2π) 2 G ii zw (q)G ii zw (p + q) − G ii zz (q)G ii zz (p + q) + (n − 2)G ii zw (q)G ij zw (p + q) − (n − 2)G ii zz (q)G ij zz (p + q) + (n − 1)G ij zz (q)G ij zz (p + q) − (n − 1)G ij zw (q)G ij zw (p + q) + (n − 1)ω 0 d 2 q (2π) 2 G ij zz (q) + 4(n − 1)ρ 0 ω 0 ∆ d 2 q (2π) 2 G ii zw (q)G ij zz (p + q) − G ij zw (q)G ii zz (p + q) , (1) order parameter with the large-N technique by solving the model exactly in the N = ∞ limit (II B) and then showing how the BKT transition emerges at order 1/N (II C).In Sec.III we begin by coupling our model to quenched disorder and then derive the phase diagram of the theory at N = ∞ (III B).Finally, we resolve the apparent inconsistencies of the mean field limit with the Imry-Ma condition by incorporating the effects of fluctuations (III C & III D).Section IV presents our conclusions.Technical details and pedagogical reviews of known material are presented in a number of appendices.Appendix A supports the results of Sec.II C, Appendix B reviews the effects of quenched disorder in O(N ) models, Appendix C does the same for the conventional CP N model, and Appendix D supports tht results of Sec.III C.II. TWO-COMPONENT CP N MODEL FOR A U (1) ORDER PARAMETERA.The Model and its SymmetriesTo apply the non-perturbative large-N technique to a theory with a global U (1) symmetry, we wish to construct a model with an internal symmetry group G with dimension scaling as N and a U (1) subgroup which can subsequently be spontaneously broken.A simple approach is a theory of two N -component complex scalar fields z and w, each of which transforms under the fundamental representation of U (N ); the larger symmetry group is then G = U (N ) × U (N ) and each of the U (N ) products contributes a U (1) subgroup.The particular model we have chosen to study is a two flavor generalization of the CP N model, and is described by the action

. 4 )
The coupling to the disorder shown in Eq. (3.3) is manifestly invariant under the local symmetry of diagonal U (1) transformations and transforms non-trivially under the global symmetry of the relative U (1) transformations [see Eqs.(2.2)].Thus, for each realization of the disorder, the random fields explicitly break the relative global symmetry, but it remains unbroken in the ensemble of the distribution of Eq. (3. .10b) Therefore, the η 1 disorder is neutral under all U (1) transformations, while the η 2 disorder explicitly breaks the replicated relative global U (1) n symmetry down to its replica-diagonal U (1) subgroup; a non-zero expectation value of ζ ij and ω ij will also spontaneously break the entire U (1) n × U (1) n symmetry group (see Appendix C for a discussion of the symmetries of the conventional CP N model).This means that the natural replica-symmetric and Z 2 exchange-symmetric saddle point ζ ij (x) = ω ij (x) = ω 0 (1 − δ ij ) is actually a privileged choice and expanding around this specific configuration loses information about the global phase diagram.However, this assumption allows us to obtain a closed-form expression for the disorder-averaged effective potential.We will first derive this effective potential and then .22b)Note that the pole mass of the order parameter is proportional to the square-bracketed term in Eq. (3.22a); substituting the saddle-point values for m 2 and ω 0 recovers the critical value K c (η tot ) above which the σ 0 = 0 state becomes unstable.To compare this propagator with the disorder averaged propagator of an O(N ) model (see Appendix B), we can put the propagator into canonical form by rescaling σ by the coefficient of p 2 .Then, using the small momentum expansion of the kernels,

FIG. 3 .
FIG. 3. One-loop self energy and three-point vertex contributions from fluctuations of ζ ij and ω ij .

ACKNOWLEDGMENTS
MCO thanks J. Gliozzi for useful discussions.This work was supported in part by the US National Science Foundation through the grants DMR 1725401 and DMR 2225920 at the University of Illinois.
Appendix A: Response of Two-Component CP N Model without Disorder to a Background U (1) Gauge Field

4 )
The field ω ij has a simple physical interpretation by comparison with the O(N ) nonlinear sigma model with quenched disorder, Eq. (B.3).Evidently the amplitude of ω ij plays the role of an effective disorder strength for the SU (N ) scalars z.The crucial difference is that ω ij must be allowed to fluctuate to ensure gauge invariance is respected.After rescaling ω ij → ω ij /g, defining g = g 0 /N and η = η 0 /g 0 , and integrating out the z j , one obtains an effective actionS eff /N = Tr ln −diag(D 2 µ [a]) + diag(λ) − ω + tr d d x • )denotes a matrix which is diagonal in replica indices, ω is the matrix with elements ω ij , and Tr includes a trace over functional configurations and replica indices.Before proceeding, we briefly comment on the symmetries of the replicated theory.The original non-replicated theory has a local U (1) gauge symmetry with the transformation rules z(x) −→ e iθ(x) z(x), a µ (x) −→ a µ (x) − ∂ µ θ(x).(C.6)