Thermodynamic phase diagram of the competition between superconductivity and charge order in cuprates

We argue that there is a special doping point in the phase diagram of cuprates, such that the condensation of holes into a charge-ordered and into a superconducting phase are degenerate in energy but with an energy barrier in between. We present a Monte Carlo simulation of this problem without and with quenched disorder. While in the clean case charge order and superconductivity are separated by a first-order line which is nearly independent of temperature, in the presence of quenched disorder, charge order is fragmented into domains separated by superconducting filaments reminiscent of the supersolid behavior in $^4$He. The resulting phase diagram is in good agreement with the experiments.


Introduction
There is an overwhelming experimental evidence [2][3][4][5][6][7][8][9] that competition between charge order (CO) and superconductivity (SC) occurs in high-critical-temperature superconducting cuprates.It has been argued [8][9][10][11][12][13] that, under certain circumstances, the superconducting order parameter with U(1) symmetry and a commensurate charge-density-wave (CDW) parameter, with Z 2 (Ising) symmetry, can be encoded in a single-order parameter with higher symmetry.Evidence for this emergent symmetry stems from studies where a non-thermal parameter which couples non-linearly to one of the orders [4,5,[7][8][9]14,15] controls the balance between SC and CO.For example, a uniform magnetic field H disfavours SC with respect to CO.At zero magnetic field, upon reducing the temperature, the correlation length of CO starts to grow at a temperature T QC larger than the superconducting critical temperature T c , as if the system was approaching a charge-ordered state [16].With further lowering the temperature, however, p O (3) Figure 1: a) Phase diagram of YBa 2 Cu 3 O 6+x in the temperature vs. hole doping plane (the figure is adapted from Ref. [1]).The magenta region is the CO induced by a magnetic field.The phase diagram can be seen as a two-dimensional projection of a three-dimensional phase diagram with the magnetic field axis running perpendicular to the plane of the figure.p O (3) indicates the "O(3)" doping at which CO and SC are nearly degenerate.the growth of the CO correlation-length stops near T c , where no CO but rather SC develops.Extrapolating the divergence of the CO correlation length to the superconducting region shows that the temperature T CO of the putative charge-ordered state coincides with the actual threedimensional ordering temperature T 3D CO that is reached once SC is suppressed by a sufficiently strong magnetic field.One fact that points to an approximate O(3) symmetry between the two orders is that the critical temperature for the field-induced CO (and therefore the putative charge-ordering temperature) obeys T 3D CO ≈ T c near a hole doping content, p O(3) ≃ 0.12, hereafter "the O(3) point", in the underdoped region of the phase diagram.As schematically shown in Fig. 1, this has been seen in YBa 2 Cu 3 O 6+x with various probes as nuclear magnetic resonance [1], sound velocity, and Hall effect measurements (see Ref. [17], and references therein).The degeneracy of the ordering temperature at the O(3) point strongly suggests that the tendency towards CO and SC are the same instability manifesting in different channels.
It is useful to visualize this phase diagram with an extra control-parameter axis, such as the magnetic field, perpendicular to the T vs. doping plane.Fig. 2 shows various experimental cuprate phase diagrams (a-c) and compare them with 4 He (d).In panels a and b, the control parameter is the magnetic field, which tunes the energy of the superconducting phase.We will refer to this situation as a SC-driven transition.In contrast, in panel c the control parameter is the isoelectronic doping, favouring stripes [18,19], which will be referred to as a CO-driven transition.The superconducting and charge-ordered phases meet the disordered phase at a so-called "bicritical" point [20], analogous to the bicritical point in 4 He (d).An alternative terminology is that of "triple" point, usually adopted when the transition lines are first-order, as the point where the liquid, gas, and solid phases of a substance meet.
In all the phase diagrams of Fig. 2 a-c, a superconducting foot develops underneath the charge-ordered state at low temperatures [wavy shading in (a-c)], which in Refs.[8,9] was associated with the occurrence of filamentary SC (FSC).This is attributed to a tertius gaudens (rejoicing third) effect.The quenched disorder breaks the CO into domains hosting different variants of CO, related to each other by discrete translations.At the interface between two variants of CO, both get frustrated and SC is stabilized.The same principle is believed [23] to explain the appearance of supersolid phases in 4 He (wavy shading in panel d).Similar effects, dubbed "fragile" SC, have been found in a Landau-Ginzburg-Wilson theory of competing orders [24].
A striking characteristic of the phase diagrams in Fig. 2 is the nearly vertical separation between CO/solid and SC/superfluid phases, once FSC is disregarded.In cuprates, several probes near p O (3) show that, in a wide temperature range, the critical magnetic field H * c , at which a long-range charge-ordered phase stabilizes, is nearly independent of the temperature.This is seen in sound-velocity data [4], resonant inelastic X-ray scattering [15] and nuclear magnetic resonance [1,7] experiments on YBa 2 Cu 3 O 6+x , as shown in Fig. 2a.The transition to the superconducting phase has been determined by the anomaly in the density of states probed by specific heat measurements, which coincides with T c determined from other methods (see Ref. [7], and references therein).A similar phase diagram can be deduced from magnetotransport experiments [8,9,25,26] in La-based cuprates.Here, the lines are not sharp anomalies, probably due to stronger disorder, but the general topology is the same (see Fig. 2b).
Applying a magnetic field is not the only way to tune the balance between CO and SC.An alternative path is the structural enhancement of CO introduced by Tranquada and collaborators [27].In this case, isovalent doping induces a structural distortion which couples with CO.When plotted against the isovalent doping concentration [19,[28][29][30], the phase diagrams bear a striking resemblance with the magnetic-field-controlled phase diagrams for similar hole content (see Fig. 2c).We notice, on passing, that the phase diagram of hole-doped La-based cuprates is characterized by a considerable (and unavoidable) coexistence of different structural phases, i.e., the low-temperature tetragonal (LTT) and the low-temperature orthorombic The onset of superconducting correlations (labelled 2D-SC) was detected by an anomaly in the specific heat (adapted from Ref. [7]).Three-dimensional CO (3D-CO) was detected with X-ray diffraction [3] and sound velocity [4].b) La 2−x Sr x CuO 4 with x = 0.08.The onset temperatures for SC (T ONS ), FSC (T ONSF ), and CO (T CO ) were extracted from magnetoresistance data (adapted from Ref. [9]).The darker blue region corresponds to zero sheet resistance R □ = 0.The red line R Q shows the locus of the quantum of resistance.c) La 1.85− y Nd y Sr 0.15 CuO 4 .Here, isoelectronic doping ( y) favours the LTT, which stabilizes CO in the form of stripes [18].Open circles show T c while solid circles show the structural transition T d which is known to be close to the CO temperature (adapted from Ref. [19]).d) Phase diagram of 4 He (adapted from Ref. [21]).In a)-c) the wavy-shaded regions denote the coexistence of CO and SC which in Ref. [8,9] was attributed to filamentary SC.In d) the wavy shading is the 4 He-analogous supersolid phase observed in Ref. [22].
(LTO) phases [31].Here, the assigned filamentary phase is more prominent, which can be understood as the effect of higher disorder (cf.with Fig. 4a in Ref. [9]).
Notice that in panels a and b of Fig. 2, T c evolves rapidly with the control parameter, while the CO temperature is approximately constant.The situation reverses in panel c.This is due to the different way the control parameter couples to the two competing phases.In panels a and b, the magnetic field destabilizes the superconducting phase having little influence on the CO.Instead, in panel c the structural distortion stabilizes the CO phase and has little effect on the superconducting T c .In the case of 4 He, pressure stabilizes the solid phase, which indeed shows up as a larger slope of the critical temperature vs. pressure line.
Summarizing, the cuprate phase diagrams are in many respects similar to the 4 He phase diagram shown in Fig. 2d, namely: i) similar critical temperatures for SC/superfluid and CO/solid phase, ii) vertical transition lines and iii) FSC/supersolid phase induced by disorder and grain boundary effects.
For 4 He, the control parameter is the pressure P. At low temperatures, and at a critical pressure nearly independent of the temperature, the superfluid phase transforms into the solid phase.The analogy between a continuum system like 4 He and lattice systems as cuprates is not new.Indeed, at least from a theoretical point of view, there is a long tradition [32,33] of modelling 4 He on discrete lattices, very similar in spirit to the model we shall be using for cuprates in this work.Also, analogies between the phase diagram of cuprates and 4 He had been emphasized before [34].
The vertical transition line T SF↔S (P) in panel d between the superfluid (subscript SF) and the solid phase (subscript S), more often drawn as horizontal with exchanged axes, was one of the first arguments put forward by London to advocate for some form of condensation in the early days of research on superfluids.Indeed, according to the Clausius-Clapeyron equation, relating the changes in entropy ∆S and in volume ∆V at a first-order phase transition, the divergent derivative dT SF↔S /dP = ∆V /∆S implies that ∆S = 0, thus identifying the superfluid as a practically zero-entropy state, like a solid phase [35].Analogously, the divergent dT SC↔CO /dH = ∆M /∆S (with ∆M , the change in magnetization) in cuprates implies that at the critical field, the same quasiparticles flip from a momentum-condensed state to a real-space condensed state, that represents two equally ordered (low-entropy) states.
A minimal model to investigate the instability that can occur in the particle-hole or particleparticle channel, i.e., pre-formed electron pairs that are paired in real space (precursors of CO) or in momentum space (Cooper pairs, precursors of SC) is the two-dimensional attractive Hubbard model [36].This model enjoys the property that, exactly at half-filling (one electron per unit cell), CO and SC are degenerate.Moving the electron density away from half-filling usually tilts the balance in favour of SC, unless other interactions (e.g., a nearest-neighbor repulsion) are present.
The attractive Hubbard model can be mapped onto the repulsive Hubbard model [36], which has a spin-density-wave ground state.In this representation, spin-density-wave order along z describes a charge-ordered state, while order in the x y plane describes SC.Order along any other direction maps into a uniform "supersolid" order, where SC and CO coexist.Since the free energy of the repulsive model is invariant with respect to O(3) rotations of the order parameter, one concludes that the charge-ordered, superconducting and supersolid phases are degenerate.In two spatial dimensions, this suppresses the ordering temperature to T = 0, due to the Mermin-Wagner theorem [37].
The critical temperature of underdoped cuprates is believed to be determined by phase fluctuations.Evidence in this direction comes from the famous Uemura plot [34] where T c appears as proportional to the superfluid stiffness.Photoemission data provide further evidence, showing that in the same doping region the superconducting gap does not close at T c but instead gets filled [38].This prompts us to examine the limit of a large Hubbard coupling where amplitude fluctuations are completely neglected and only phase fluctuations are taken into account.
This limit is quite insightful to study the competition between CO and SC, because the model can be mapped onto a Heisenberg model [36] of interacting pseudospins (analogous to Anderson's pseudospins).Here, the pseudospin projection along z, up or down, encodes a double occupied or an empty site on the lattice, respectively, while the in-plane component encodes superconducting correlations.Adding other terms to the model, one can move away from the O(3) point, meaning that either the X Y (SC) or Ising (CO) contribution prevails.The superconducting transition in two dimensions belongs to the Berezinskii-Kosterlitz-Thouless (BKT) universality class, and is characterized by the appearance of a finite stiffness and the binding of vortices and antivortices.
While the Hubbard model is genuinely quantum, it is well known [39] that, provided the ground state is ordered above a microscopic scale (the Josephson correlation length) the system can be described by a classical theory.In the Heisenberg limit, the Josephson correlation length reads ξ J ≈ ħ hc/J s , where J s and c are the stiffness and the zero-temperature spin-wave velocity respectively.Using estimates appropriate for the Heisenberg model [39], yields ξ J = 2πa/C s ≈ 10.9 a with C s a constant, and a the lattice spacing.Therefore, we can use a classical effective lattice spin model to study the competition between CO and SC.Each pseudospin in the lattice model represents a cluster of elementary unit cells with a linear dimension of order ξ J , behaving as a classical variable.
While the above scenario is very appealing to formulating a statistical mechanical description of the competition between CO and SC, a full O(3) symmetry is clearly a drawback of the model.Indeed, for the repulsive Hubbard model, this is a consequence of rotational invariance, but there is no such fundamental symmetry in a generic attractive model.One expects that CO and SC can be tuned to an approximate O(3) symmetry point by a non-ordering field (e.g., p for cuprates), but there is no reason why the barrier between these states should vanish at the O(3) point.In other words, the O(3) symmetry is only approximate, in the sense that the charge-ordered and superconducting phases are still degenerate but are separated by barriers.
Based on the above considerations, we construct an effective XXZ model that accounts for the experimentally observed phenomenology.We study an effective classical spin model on a square lattice, with nearest neighbour exchange interaction and three relevant parameters: an exchange anisotropy, to tilt the balance between the easy-axis (charge) and the easy-plane (superconducting) order; a potential barrier, to remove the unphysical high degeneracy of the O(3) symmetric point; and a random field to mimic disorder.
In the clean system (without disorder), we find that the presence of the barrier allows for a finite-temperature phase transition, otherwise forbidden at the O(3) point.Once disorder is taken into account, CO is fragmented into different domains, resulting in a polycrystalline charge-ordered phase, and FSC sets in as a parasitic phase at the domain boundaries [8,9].
Our analysis is carried out by means of Monte Carlo (MC) simulations, which allow us to study not only the ground state, as in Ref. [9], but the thermodynamic phase diagram itself and the behaviour in temperature of the various physical quantities.
While the above discussion was based on the attractive (negative−U) Hubbard model, which has an s-wave superconducting state, the statistical mechanical model we study can be considered as a simplified version of the CO-SC competition problem in a broader context.Indeed, the X Y model is by itself relevant both for s-wave and d-wave superconductors (see for example Ref. [40]).Furthermore, the possibility of encoding CO and SC order parameters in a single one with higher symmetry has been advocated also in models of cuprates supporting d-wave superconductivity [11][12][13].On the other hand, physical systems may require more complicated order parameters.For example, the d-wave model of Ref. [11] has a larger symmetry group than the negative−U Hubbard model.Furthermore, commensurate CO in cuprates often shows a periodicity of four lattice sites and is unidirectional [41].That makes eight different variants of CO domains related by symmetry, while our model has only two CO variants.Likewise, 4 He would have infinite variants of positional order as the crystal breaks a continuum symmetry.Yet, the phase diagrams shown in Fig. 2 are all very similar, which suggests that in a phenomenological approach, and as a first approximation, these complications can be ignored.
The scheme of this work is the following.In Sec. 2, we discuss the model and methods of investigation.In Sec. 3, we discuss the properties of the model in the absence of disorder, highlighting the role of the potential barrier and using the result of the model in the absence of a barrier (see also Refs.[42][43][44]) as a benchmark.In Sec. 4, we include the effect of disorder and show that this is crucial to promote FSC.Finally, in Sec. 5 we discuss the resulting phase diagram for a three-dimensional system which is in very good agreement with experiments.Our concluding remarks are found in Sec. 6. Appendix A contains details about the physical properties of the clean system, as well as a detailed phase diagram in the case of a larger barrier.Appendix B contains some details about the physical properties of the dirty system, in particular, the behaviour of the superfluid stiffness.

Model and methods
Above the Josephson scale, we can model our system with a classical order parameter.Therefore, as in Refs.[8][9][10], we consider a coarse-grained model of classical pseudospin vectors S R , on the sites R on a square lattice, each representing a region of area ξ 2 J of the quantum system.The new (coarse-grained) lattice spacing is set as a ′ = 1 and the linear size of the lattice is L (i.e., the lattice hosts N = L 2 sites), with periodic boundary conditions.The states with positive or negative pseudomagnetization along the z axis represent two variants of the charge-ordered state, related in the original quantum microscopic model by a translation symmetry, while the in-plane pseudomagnetization describes the superconducting state [8][9][10].In order to lighten the notation, we will henceforth refer to the pseudomagnetization simply as magnetization, not to be confused with the physical magnetization mentioned in connection with the Clausius-Clapeyron argument above.In the following, we set |S R | = 1, and fix the reference frame so that the three Cartesian components of the vector S R are S x R = sin ϕ R cos θ R , S y R = sin ϕ R sin θ R , and S z R = cos ϕ R , in terms of the polar and azimuthal angles θ R and ϕ R .θ R can be identified with the phase of the superconducting order parameter.
The competition between CO and SC is captured by the classical XXZ model (anisotropic Heisenberg model) with an effective barrier potential term and a random field mimicking disorder, where the symbol 〈R, R ′ 〉 specifies that the sum runs over nearest-neighbouring sites.We fix the interaction strength J = 1 unless otherwise specified.Furthermore, we use the anisotropy parameter α ≥ 0 to tune the ground state from being superconducting to charge-ordered.This corresponds to keep constant the ground-state energy of the superconducting state and tune the one of the charge-ordered state, i.e., a CO-driven transition.As we shall see, a simple rescaling of the energy units allows describing a SC-driven transition.
The second term in Eq 1 is the barrier potential, whose height is adjusted by the parameter B. Its role, as anticipated in Sec. 1, is to eliminate the unphysical degeneracy of the charge-ordered and superconducting state with all possible intermediate supersolid phases for α = 1.Note that the introduction of such a term makes our model a phenomenological model instead of an exact mapping of the Hubbard model discussed in the previous section.In the following, we will still call the α = 1 case the "O(3)" or "isotropic" point, keeping in mind that such terminology refers only to the first term of the Hamiltonian in Eq. 1.
The last term in Eq. ( 1) is a random field that mimics impurities coupled to the charge density in a real system.We take h R as independent random variables with a flat probability distribution between −1 and +1.The strength of disorder is controlled by the parameter W .As we shall show, this term is crucial to promote the polycrystalline behaviour of CO for α ≫ 1, as well as the occurrence of FSC in a certain range of anisotropy α ≳ 1, in the form of topologically protected domain walls between regions hosting two different realizations of CO.
In the case where B = W = 0, Eq. ( 1) is the bare XXZ model which has been widely studied in the literature [42][43][44], and whose phase diagram will be used as a benchmark case when discussing the effect of the energy barrier.In the bare model, the anisotropy α allows switching from the BKT universality class, for α < 1, where the ground state is superconducting, to the Ising universality class, for α > 1, where the ground state is charge-ordered.Finally, in the isotropic limit α → 1 +,− , the critical temperature goes to zero logarithmically [45], and at α = 1 no finite-temperature phase transition is possible, according to the Mermin-Wagner theorem [37].The presence of a finite energy barrier separating the three equivalent ground states ϕ R = 0, π/2, π (i.e., S z R = −1, 0, 1) makes the model no longer invariant with respect to O(3) rotations of the order parameter, and the Mermin-Wagner theorem does not apply.Indeed, we find that, for B > 0, the ordering temperature remains finite for all α.Furthermore, since the effect of the barrier persists at finite temperatures, metastability regions appear in the resulting T vs. α phase diagram.
In order to study the physical quantities related to the effective Hamiltonian, Eq. ( 1), as functions of the temperature, we performed large-scale MC simulations, with systems of linear size L ranging from L = 16 up to L = 256.We used the Metropolis and simulated annealing algorithms to optimize the thermalization process: at the highest temperature reached in our calculations the system evolves from an initial configuration of random pseudospins until it reaches its equilibrium state, then the temperature is slightly decreased and a new thermalization starts from the final configuration of the previous step.This process is iterated until the lowest temperature of interest is reached.At each Metropolis step, the whole lattice is updated according to the Metropolis prescription [46], either sequentially updating all the L × L pseudospins or by L × L random choices of the pseudospin to be updated.Thermal averages of any observable O are obtained as the average over N MC (at least 10 3 ) measures, taken τ MC Metropolis steps apart from one another, τ MC being of the order of, or larger than, the autocorrelation time (for the clean system, typically we take τ MC = 30 − 100 Metropolis steps).To account for the thermalization time, we discard the initial N out (at least 10 5 ) Metropolis steps.In the presence of the random field, we also average over N dis realizations of the disorder (henceforth, this average is marked by an overline).

Physical observables
In order to assess a global BKT superconducting transition, we compute the in-plane superfluid stiffness J s , associated with the superconducting phase rigidity and defined as the second derivative of the free energy with respect to a twist of the phase angle δθ , e.g., along the x direction: where Z is the partition function and x is the unit vector in the x direction.
To perform the extrapolation to the thermodynamic limit of the BKT critical point we use the BKT scaling of the superfluid stiffness [47], where T BKT is the BKT critical temperature and we take L 0 as a fitting parameter.Whenever needed, we complement our analysis of the superconducting state, with a closer inspection into the occurrence of vortices in the pattern of the local superconducting order parameter (the in-plane magnetization).A vortex (anti-vortex) is identified whenever a variation of 2π (−2π) of the superconducting phase θ R is found in a closed path around a single plaquette with side equal to the lattice spacing.Defining the superconducting phase difference at site R in the direction of the ν unit vector (ν = x , ŷ) as the notation [•] +π −π meaning that we take the value modulus 2π so that Θ ν(R) ∈ (−π, π], the circulation of the superconducting phase around a plaquette whose center is located at where n R = ±1 is the integer vorticity in the phase angle θ R going around the plaquette [48].Summing over all positive (negative) vorticities per unit length we obtain the density of vortices (antivortices) ρ V > 0 (ρ AV < 0), defining the total vorticity as Concerning the charge-ordered state, we define T CO using the crossing point (as a function of temperature) of the kurtosis of the pseudospin distribution function, i.e., the Binder cumulant Its value at the critical temperature is indeed less sensitive to finite-size effects, as compared to the CO order parameter 〈m z 〉, and unbiased by fitting functions and a priori scaling hypotheses.
In the thermodynamic limit N → ∞, one expects U N → 0 in the high-temperature limit, while U N → 2/3 in the ordered phase for T → 0 [49-51].
3 Clean system

Metastability and spinodal lines
Let us begin to discuss the MC numerical results starting with the clean case in the presence of a finite energy barrier, i.e., B > 0. As can be expected, one prominent effect of the barrier is to introduce metastability in the system.Already at zero temperature, there exists a range of values of α * CO (B) < α < α * SC (B) where both the superconducting and the charge-ordered phase are local minima of the energy [52].Indeed, while α * = 1 marks the first-order phase transition between the two ordered states (SC and CO), the spinodal points α * CO,SC (B) mark the limit of stability of the less stable phase: CO remains metastable down to α * CO < α * , and SC up to α * SC > α * .At zero temperature, the spinodal points α * CO,SC (B) can be calculated analytically as a function of the barrier height B. Indeed, increasing α at T = 0, we can assume that all pseudospins are parallel, i.e., ϕ R = ϕ and θ R −θ R ′ = 0. Up to constant terms in the angle ϕ, the total energy per site from Eq. ( 1) takes the simplified form The resulting energy landscapes for B = 0.2 and B = 2, with varying α, are given in Figs.3a and 3b, respectively.The free energy has a single superconducting global minimum at ϕ = π/2, up to some α = α * CO (B), after which two new local (metastable) minima at ϕ = 0, π appear, corresponding to the two possible CO realizations.Crossing the first-order critical point α * , the situation gets reversed: the new global minima are found at ϕ = 0, π and the ϕ = π/2 configuration becomes a local (metastable) minimum, disappearing at some α * SC (B), after which only the two equivalent charge-ordered states survive.
By substituting ϕ = π/2 and ϕ = 0, π in the second derivative of Eq. ( 5), we obtain the two spinodal points at zero temperature, For B = 0 the model becomes fully symmetric and the zero-temperature spinodal points merge into the zero-temperature transition point.
The study of metastable states at finite temperatures is particularly challenging, both in real experiments and in numerical simulations.If the system is prepared in a metastable state, there is a finite probability that a bubble of the more stable phase nucleates in a finite time and then grows.Therefore, sooner or later the system will transit to the more stable phase.As a consequence, spinodal lines, marking the limit of stability of the metastable phase, require time-domain considerations to be defined.A metastable phase is well-defined only if it persists at least for its equilibration time, otherwise, it can not be considered a thermodynamic phase.The line where the equilibration and the nucleation times become equal defines the dynamical spinodal line [53].While these considerations allow for a rigorous definition of spinodals, here we circumvent the problem of estimating the nucleation times and take a more pragmatic approach, which is enough for our purpose.We then proceeded as follows.For a small enough size, during a simulation, the system may get trapped in one metastable state for a long time (compared to the equilibration time of the state) and sporadically change to a different state, where it gets trapped again.Evolving the system long enough, we can construct a reliable effective probability density function P eff (m z ) for the order parameter associated with CO, that is encoded in By definition, the free energy of the system is a) b) Figure 3: Energy landscapes at T = 0 [Eq.( 5)], for various α and: a) B = 0.2; b) B = 2. Notice that in this case even for α = 0 the charge-ordered phase remains as a metastable minimum, and there are two possible realizations (two equivalent minima).
Notice that when drawing the histogram in m z one automatically takes into account the correct measure for the probability distribution.It is easy to check that, for B = 0, the probability distribution is flat at α = 1, P(m) = 1 2 , consistent with a flat free energy.
The numerical identification of the spinodal lines and the equilibrium coexistence line of SC and CO can be obtained by studying the form of F (m z ) as a function of the parameters (see App. A.3 for more details).For a fixed value of α, the first-order transition temperature T SC↔CO , is defined as the temperature at which the absolute minimum of the free energy changes from m z = ±1 to m z = 0.A spinodal temperature T sp is defined by the local minimum of a metastable phase becoming an inflection point.Thus, once we extracted all the effective free energies F (m z ) at each temperature, we can infer the full phase diagram.
Of course, it must be borne in mind that the condition m z = 0 cannot distinguish a superconducting state from a charge-disordered state.To construct the phase diagram one has to complement the previous study with the computation of the superfluid stiffness.

Phase diagram for B = 0.2
The resulting phase diagram is reported in Fig. 4. For comparison, we also report the phase diagram for B = 0 (light-blue line).Cyan circles and purple triangles refer, respectively, to the critical temperatures T BKT and T CO calculated using the scaling laws of J s [Eq.(2)] and U N [Eq.( 4)].The points in green along the line T = 0 are the analytical results: the two squares at α * CO = 0.6 and α * SC = 1.4 are the spinodal points, calculated as described in Sec. 3 [see Eq. ( 6)], and the green square at α * = 1 is the value at which the free energy has three equivalent minima at m z = 0, ±1 (first-order phase transition at T = 0).
The presence of the barrier with B = 0.2 shifts the superconducting transition line to higher temperatures, and to a higher value of the anisotropy parameter α at T > 0, α B = 1.04.For slightly larger values of α, we recover the CO (Ising) transition line, which is shifted downwards with respect to the case B = 0. Note that the nature of the transition, and hence the shift of α B , can be also inferred by looking at the susceptibility of the in-plane (SC) and of the out-of-plane (CO) spin magnetizations.See App.A.1 and A.2 for more details.
The non-trivial consequence of the shift of α B towards higher values translates into a positive slope of the first-order line, being α B > α * , indicating that entropy slightly favours SC over CO.In agreement with our previous discussion, CO is indeed the stable phase at low temperatures for α ≳ 1, then, with the increase of the temperature, the system switches to SC eventually ending in a disordered state for higher temperatures (see App. A.3).We predict that in a very clean system close to the O(3) point, the very interesting phenomenon of SC stabilized by temperature could be seen.

Dirty system
We now discuss the role of disorder.The localizing effect of impurities is not expected to significantly alter the BKT transition found at α < 1 [54].For α > 1, a study at zero temperature has shown that the effect of disorder is to break CO into a polycrystalline state [8][9][10].This can be seen in Fig. 5, where we show low-temperature snapshots (T = 0.001) of the MC simulations for increasing α.The colour code maps the CO order parameter ranging from S z = +1 (blue), through S z = 0 (no CO, yellow), to S z = −1 (red).We remind the reader that S z = ±1 encodes two variants of the CO, e.g., with maxima of the charge density located at two different lattice positions, connected by translational symmetry of the lattice.It was argued in Refs.[8][9][10] that at the boundary of such domains CO gets frustrated and FSC emerges.Indeed, as CO fluctuations are enhanced by increasing α, the superconducting condensate gradually loses its two-dimensional nature by forming thinner and thinner filamentary structures.As a consequence, as the superconducting cluster gets narrower, a smearing of the BKT signatures is expected.
To probe the BKT transition, we monitor the superfluid stiffness J s and its scaling according to Eq. 2. However, this will not be sufficient in our discussion because of the gradual broadening of the BKT jump of J s , along with the gradual violation of the BKT scaling relation, Eq. ( 2).Moreover, a substantial fraction of in-plane pseudospins will survive also in the charge-ordered region of the phase diagram, as it is already visible from the snapshots in Fig. 5. Thus, in the dirty system, we will also study the disorder-averaged superconducting correlation function of two pseudospins separated by r , defined as Note that the average over many disorder realizations (indicated by the overline) restores spatial isotropy at large distances.Indeed, as it is well known, one of the hallmarks of the BKT topological phase transition is encoded in the peculiar behaviour of the correlation function: in the thermodynamic limit, where J is the stiffness at T = 0 and a V is the characteristic size of a vortex core [55].The infinite correlation length of superconducting fluctuations ξ x y at T ≤ T BK T cannot be probed by numerical simulations of a finite system.Instead, in a MC simulation, one has ξ x y ∝ L. It should also be noted that, as a consequence of the presence of out-of-plane fluctuations, Eq. ( 7) acquires an extra factor sin ϕ R sin ϕ R+r with respect to the standard BKT case, in which provides instead a good estimate of the short-range SC still present in the system.It is worth noting that this quantity again does not contain the information about the coherence of the condensate, which is encoded in the θ variable.We checked, however, that for all the values of α we investigated, the spins belonging to the same superconducting cluster are indeed coherent at least at low temperatures.Moreover, the presence of disorder results in the charged ordered state losing its long-range order, so that the resulting magnetization and the corresponding Binder cumulant in Eq. ( 4) cannot be used to define T CO .For α > 1, the system gradually evolves toward a random-field Ising model and the ground state appears as a rather inhomogeneous landscape, characterized by large charge-ordered puddles.Therefore, we will use the CO correlation function defined as to characterize the behaviour of the charge-ordered state.C zz (r) is expected to decay exponentially as ∼ e −|r|/ξ z , with ξ z a fitting parameter.The random-field Ising model in two dimensions has no finite critical temperature [56] and is characterized by a finite low-temperature correlation length that grows exponentially with reducing the strength of the random field [57].The presence of the barrier term in Eq. ( 1) further suppresses transverse pseudospin fluctuations at low temperatures and enhances the clustering of up/down charge-ordered regions, even at α ≳ 1, thereby favouring the polycrystalline behaviour up to a finite temperature T CO .
We estimate T CO assuming the CO correlation length to behave as for T approaching T CO from above (without getting too close to it), with the critical index ν = 1 of the clean Ising model [58].The idea is that, starting from high temperatures, one can follow the critical behaviour of the clean Ising model, down to a temperature at which the system crosses over to the non-critical behaviour of the random field Ising model and the correlation length saturates to a finite value that determines the typical size of the clusters.

Phase diagram for B = 0.2, W = 5
We present our results for barrier height B = 0.2 and disorder strength W = 5, to explore the effect of disorder in a situation when the first-order transition between the two phases would be nearly vertical in the clean case (see Sec. 3.2).
The phase diagram T vs. α is reported in Fig. 6, where the T BKT points (cyan dots) are calculated using the BKT scaling law of J s , while T CO (purple) is computed from the fit of 1/ξ z in the temperature range where it exhibits a linear behaviour.For α < 1, the superconducting state is not much affected by the disorder, except for a small suppression of the superfluid stiffness (see App. B).Indeed, according to the Harris criterion [54], the presence of spatially uncorrelated disorder does not alter the universality class of the topological phase transition, nor induces smearing of the superfluid-stiffness jump [59,60].
Instead, non-trivial features are expected for values α ≳ 1 where spatially correlated disorder emerges from the interplay between the competition and the presence of impurities.Our two striking results are indeed found on the α ≥ 1 side of the phase diagram: i) the observation of FSC for α ≳ 1; and ii) the formation of a polycrystalline charge-ordered phase when α ≫ 1.
In Fig. 7 we report the rescaled superfluid stiffness J s (L, T ) at various L (the grey line is the critical line 2T /π), following the scaling law in Eq. (2).We find that the BKT scaling law works quite well up to α = 1.1 (see Fig. 6) although we observe deviations in the non-universal features of the phase transition.The vertical line and the grey shaded area correspond to T BKT ± σ T BKT (the error being calculated as to include the smearing of the jump and the uncertainty on the fitting parameter L 0 ).
It is not surprising that for α = 1 (panel a) we still observe a pretty clear jump of the stiffness, although the clustering of small charge-ordered regions in the system can already be observed (see the corresponding snapshot in Fig. 5).In fact, as we showed in Sec.3.2, in the clean system the potential barrier stabilizes the superconducting state up to α B = 1.04 .
Going towards α = 1.05 (panel b) we can still observe a well-defined crossing of J s with the critical line, whereas in the clean system this value of α already corresponded to a chargeordered global minimum of the free energy.
Finally, for α = 1.1, we still find a finite superfluid stiffness, and yet the usual BKT scaling relation, Eq. ( 2), has noticeable deviations as one can realize by scrutiny of Fig. 7c.Rescaling the curves according to Eq. ( 2) leads to a spread of crossing points with the critical line 2T /π.Notwithstanding that, compared with the unscaled curves shown in Fig. 19c of App.B, the different curves here exhibit convergence to a small region of temperatures showing approximate scaling.Within the BKT scenario, we obtain T BKT = 0.31 ± 0.08.This is consistent with the (negative) minimum value of the derivative of the superfluid stiffness with respect to temperature (see Refs. [42,61] and Fig. 19d in App.B).To take into account uncertainties in the definition of the BKT critical temperature, we considered a conservative confidence interval highlighted in grey in Fig. 6c.The difficulties in applying scaling relations in this case can be  linked to the emergence of new length scales, presumably related to the geometrical structure of the system, along with the typical sizes of vortex-antivortex pairs in the BKT theory.
It is very interesting that in this regime the stiffness decreases with decreasing temperature.This can be seen as a remnant of the entropy-induced SC observed in the clean system, which disfavours SC at low temperatures.Such an anomaly may be measured in samples close to the O(3) point and is another prediction of this work.
Note that the downward curvature in the low-temperature limit is not the consequence of finite-size effects, since a 256×256 lattice with periodic boundary conditions already provides a reasonable size to observe reliable thermodynamic quantities.
Let us now discuss the CO correlation lengths ξ z , used to define the polycrystalline chargeordered phase for α ≥ 1.15 (see purple dots in Fig. 6).The inverse correlation lengths 1/ξ z as a function of temperature are displayed in Fig. 8 for sizes L = 64, 128, 256.The dashed lines correspond to the linear fits, and the different colours refer to the system size L, as shown in the legend.As one can see, the linear decrease in temperature of 1/ξ z deviates towards a constant plateau when the temperature is lowered below a certain value.The CO critical temperature T CO (of the clean Ising model universality class) is defined as the intercept of the linear fit of 1/ξ z performed at large enough temperature.Note that all the curves show no sign of scaling and the low-temperature saturation value of 1/ξ z , for fixed values of the energy barrier B and the disorder strength W , only depends on the anisotropy parameter α.That signals the presence of an intrinsic length scale related to the clusters.A small downward deviation from linearity at high temperature, observed in the studied temperature interval when α ≥ 1.5, signals that the system is exiting the critical regime of the clean Ising model with further increasing T .We would have observed the same deviation for α < 1.5, at higher temperatures.We conclude this section by observing that the curves in Fig. 8 resemble the behaviour of the full width at half maximum of the CDW peak [proportional to (ξ z ) −2 ] probed in YBa 2 Cu 3 O 7−δ and Nd 1+x Ba 2−x Cu 3 O 7−δ , by means of resonant inelastic X-ray scattering in Ref. [16].There, the extrapolated T CO coincides with the temperature at which CO would occur once SC is suppressed by a magnetic field [1,14], while the saturation at low temperature, in the absence of a magnetic field, signals that CO competes with SC, and SC is more stable.Here, instead, the CO temperature obtained by this criterion near the O(3) point is much lower than the asymptotic value at large α.We will come back to this important point in Sec. 5.

Short-range and filamentary superconductivity
We discuss now the survival of a filamentary superconducting cluster and the presence of shortrange SC in the polycrystalline charge-ordered side of the phase diagram.Whereas in the case α = 1.1 it was still possible to define a BKT transition, albeit with a certain degree of uncertainty, when the anisotropy parameter is increased up to α ≥ 1.15 no T BKT can be defined from the crossing point.In particular, at α = 1.15 one can see from Fig. 9a that J s vanishes already at L = 128.The point α = 1.15 deserves, however, more attention, since it displays both a finite critical temperature T CO and a short-range coherence of the superconducting cluster, as indicated by the finite value of J s at L = 16, 32, 64.The substantial fraction of superconducting pseudospins can also be observed by comparing the snapshots in Fig. 5 (yellow component), in particular, those corresponding to α = 1.1, 1.15, 1.2.Therefore, for α ≥ 1.15, although the nearly one-dimensional nature of the superconducting cluster does not allow for the binding of vortex-antivortex pairs [62], the finite residual superconducting component can still exhibit some short-ranged stiffness.It is worth noting, once again, that in our coarse-grained model the spacing of the pseudospin lattice a ′ corresponds to the Josephson scale, i.e., a ′ ≈ ξ J ≈ 11a (see Sec. 1), meaning that L = 16 corresponds to about (16•11) 2 ≈ 31000 atoms.Such a large coherent region should have a strong impact in transport properties.While we have not computed the resistivity, it is clear that it will be quite small, as large patches of coherent regions will short-circuit the sample.We speculate that a broad transition should be observed with a large drop of the resistivity to a small but finite value.We thus include in our phase diagram the temperature at which we find short-range SC (blue symbols) in the FSC region of our phase diagram of Fig. 6.Those points indicate the crossing of J s at L = 16 with the universal critical line 2T /π.A finite, even if exponentially small, stiffness is found up to values α = 1.2.This behaviour of SC with T c ≈ 0 is reminiscent of transport experiments in cuprates [9,25].
In order to get a more quantitative idea of short-range SC, in Fig. 9b we show the average in-plane component 〈sin ϕ〉 as a function of temperature, for different values of α.Note that in the standard BKT model this should be identically equal to one at T = 0.However, the presence of out-of-plane (corresponding to CO) fluctuations renormalizes it to a lower value, that decreases with increasing α.We indicate with gray crosses and with hexagons the two critical temperatures, respectively T BKT (α < 1.15) and T CO (α ≥ 1.15), which we discussed above and presented in the phase diagram (Fig. 6).For α < 1, where SC is well described within the BKT scenario and no spatially correlated disorder emerges, 〈sin ϕ〉 increases quite monotonically with lowering the temperature.As α = 1 (in orange) we still observe the monotonic increase of 〈sin ϕ〉 with decreasing T , and from superfluid stiffness computations we know that the system still exhibits quite clear BKT signatures (see Sec. 4.1 and App.B).
For 1 < α < 1.15, at high T we observe first a slow decrease of 〈sin ϕ〉 with increasing the temperature followed by an inflection point at T infl ≳ T BKT .This range of the control parameter α lies inside the region of the phase diagram that we labelled with FSC in Fig. 6.We stress again that up to α < 1.1, it is still possible to define the BKT temperature from the jump of the superfluid stiffness, which is smeared out but still clearly visible.For α = 1.1 (light green), instead, the BKT scaling law starts showing deviations, and we observe a downturn of 〈sin ϕ〉 at T down < T BKT .This may be related again to the entropically favoured SC of the clean case, which might also be the cause of the downturn of J s at very low temperatures.
The curve for α = 1.15 (dark green) highlights again an interesting crossover scenario, which presents a filamentary pattern, clearly visible in the snapshots, but no long-range stiffness (see Fig. 9a).In this case, the decrease at high temperature follows a behaviour similar to the one found for α = 1.1, but with no inflection point, down to T = 0.2 (indicated with a red dot).By further lowering the temperature, 〈sin ϕ〉 becomes steeper.The absence of an inflection point in 〈sin ϕ〉 might be a proxy that the entropically favoured superconducting state is now suppressed by the large CO fluctuations, although a small J s survives at finite L, becoming exponentially small with increasing the size.In fact, by comparing 〈sin ϕ〉 for α = 1.15 with the corresponding J s , one can observe that short-range SC is still present (L < 128).Note that the curve for L = 16 of Fig. 7 has a maximum for T = 0.3, then decreases by further lowering T , crossing the critical line 2T /π at T = 0.2 (red dot).We point out that at the lowest temperature, T = 0.001, a substantial in-plane residue survives, 〈sin ϕ〉 = 0.35, exhibiting phase coherence.Even increasing the anisotropy up to α = 1.2, the superconducting fraction is still about 20%.We thus include 〈sin ϕ〉 at T = 0.001 and L = 256 in our phase diagram in Fig. 6 (right axis, in green) to stress out the presence of a macroscopic superconducting residue, that can show signatures in transport experiments even if it lacks long-range coherence.
In order to investigate the role of superconducting phase fluctuations in this crossover filamentary state, we analyze the correlation length ξ x y .In Fig. 10 we present ξ x y for α = 1, 1.05, 1.1, 1.15 and L = 32, 64, 128, 256, as found by fitting the correlation function in Eq. ( 7).For α < 1.15, we find that ξ x y ∼ L, thus following the expected BKT scenario, thereby justifying the BKT analysis discussed above.Note that the black-dashed lines mark T BKT , extractred from the crossing of J s (Fig. 7).For α = 1.15, although J s vanishes at L = 128, we can still observe BKT-like features of ξ x y : the saturation value at low temperatures is in fact increasing with L, with some slowing down for L = 256.Again, we mark the temperature T = 0.2 with a dashed red line: this temperature corresponds to the maximum of J s found at L = 16 and to the change of slope in the decrease of 〈sin ϕ〉.The behaviour of ξ x y (T ) at this temperature seems to suggest the occurrence of an "avoided" superconducting state, reminiscent of the findings in transport experiments [9,25].

Phase diagram for a three-dimensional system
Starting from our detailed analysis of the two-dimensional model, we can make a comparison between the theoretical phase diagram Fig. 6 and the experimental phase diagrams of the CO-SC competition of Fig. 2. The most striking difference is that the CO temperature is strongly suppressed near the O(3) point, both when compared with experiments and with the clean case.Also, FSC does not develop an evident foot in the CO region, although local SC regions are present.To a large extent, both deficiencies can be ascribed to the low dimensionality of the model.To show this, we compute a phase diagram of a three-dimensional system assuming an interlayer coupling both for the CO and the SC phases.As it is known from studies of the layered Heisenberg model, the transition temperatures T 3D SC and T 3D CO can be estimated from the superconducting and CO correlation length [63] solving the following equations, meaning that the interlayer energy associated with a correlated region of area (ξ x y,z ) 2 is of the order of the critical temperature.
Figure 11: (a) CO-driven and (b) SC-driven phase diagrams.Critical temperatures are normalized with respect to J and αJ, respectively.The 3D critical temperatures are calculated from the correlation length according to Eq. ( 8) with J x y ⊥ = 0.01 and J z ⊥ = 0.1.Other parameters are as in Fig. 6.We specify that the FSC point in α = 1.5 is extrapolated from the vanishing of the in-plane superconducting component (see purple curve in Fig. 9b).For comparison, the grey lines in panel a) refer to the twosystem.
The superconducting and charge-ordered interlayer couplings, J x y ⊥ and J z ⊥ respectively, are not known.In view of the quasi two-dimensional nature of cuprates we take J z ⊥ = 0.1 and J x y ⊥ = 0.01.The much smaller value of the superconducting coupling with respect to the charge-ordered one is by the fact that CO is coupled by the long-range Coulomb interaction while SC is coupled by Josephson tunnelling through the insulating layers, and one expects a large difference between these two scales.For the rest, these parameters are rather arbitrary, but the qualitative form of the phase diagram is not expected to be sensitive to the precise value of the couplings.
Fig. 11a shows the resulting phase diagram.The dotted grey lines refer to the original two-dimensional system for comparison.Panel a follows our previous convention of measuring energy and temperature in units of the superconducting scale J so that the scale of the CO state changes with α.This corresponds to the CO-driven transition mentioned in the introduction.In panel b, by simply rescaling our energy units, we derive the phase diagram for the SC-driven transition.Here, the CO energy scale is, by definition, constant.
We see that, indeed, the phase diagram of the three-dimensional system bears a strong resemblance with the experimental phase diagrams for the SC-driven case (Fig. 2a,b) and COdriven case (Fig. 2c).Now the bicritical point is at a temperature of the order of the ordering temperature and the FSC foot extends more into the CO region.
Comparing phase diagrams of Figs.2 and 11 allows one to estimate J/k B ≈ 100 K as a reasonable parameter for cuprates.The barrier height can not be very large otherwise the slope of the SC-CO transition would be more pronounced than in the experiment (cf.Fig. 4 and Fig. 16 in App.A).Disorder also tends to change the slope of the transition [9] and depends on the sample quality.The value taken here is appropriate for not-too-disordered samples.

Conclusions
We used Monte-Carlo simulations to solve a statistical mechanical model of a two-dimensional system presenting competition between SC and CO, both in the absence and in the presence of quenched disorder.We computed thermodynamic quantities, correlation functions, and thermodynamic phase diagrams.
In a clean system, the competition mechanism generates metastability regions in the phase diagram, bounded by two spinodal lines and encompassing the first-order phase transition line.As the temperature increases, the region of metastability shrinks to a single point, which coincides (within numerical accuracy) with the bicritical point, where the charge-ordered, superconducting and disordered phases meet.The first-order line separating the charge-ordered phase and the superconducting phase is rather steep for low values of the barrier height (B = 0.2), indicating that the two phases have similar entropy, as one can check using the Clausius-Clapeyron relation.We can thus make a comparison with the case of 4 He [32,33], where the almost vertical line separating the solid and superfluid phases led to the hypothesis that superfluidity was a low-entropy phase, as a crystal, fuelling explanations based on condensation in momentum space rather than in real space.
A closer inspection shows that the first-order line is not exactly vertical, and a re-entrance appears, thus showing that near α ≳ 1 one can make a transition from the charge-ordered phase to the superconducting phase by increasing the temperature.This means that the superconducting phase has actually slightly higher entropy than the charge-ordered phase.A posteriori, this result is reasonable as the charge-ordered state has two gapped transverse modes while the superconducting state has one gapped mode and one Goldstone mode.Thus, just considering low-lying excitations near T = 0, it is reasonable that the superconductor can have larger thermal fluctuations and entropy.Interestingly, the re-entrant behaviour of the superconducting phase is also reminiscent of the phase diagram of 4 He, in which a range of pressures is found where the solid 4 He, if heated, transits to its superfluid state before becoming a simple liquid.In 4 He, however, this happens in the high-temperature part of the phase diagram, while here we observe it at low temperatures.In fact, in the low-temperature region, the slope of our phase diagram and the one of 4 He have opposite sign.We speculate that this qualitative difference is due to the fact that in our case the charge-ordered state has no Goldstone modes while in the case of 4 He the crystal has sound (Goldstone) modes.
The metastability regions and the first-order transition line disappear when quenched disorder is considered, giving rise instead to a phase-separated region where FSC appears.Indeed, moving from the superconducting to the charge-ordered regime we find the gradual disappearance of the two-dimensional superconducting phase towards a polycrystalline chargeordered phase with the tuning of the anisotropy parameter α.As the BKT signatures disappear, one-dimensional-like superconducting patterns still survive inside the polycrystalline chargeordered phase.
As mentioned in the introduction, in our model, up and down pseudospins encode only two possible realizations of CO corresponding to a checkerboard pattern in a bipartite lattice.This is a simplification of cuprates where, for example, non-magnetic charge stripes with a fourfold periodicity have four CO variants for each orientation, yielding 16 possible "colours" of CO patterns.Furthermore, cuprates have spin order at low doping whereas magnetic fluctuations are completely neglected in our model.Still, our simplified two-colour CO model captures many subtleties of the phase diagram near the O(3) point.
In Refs.[8][9][10] it was already proposed that disorder may have a peculiar effect in the coexistence region discussed above, turning the metastable superconducting state into a stable state, where FSC is topologically protected at the boundaries between different chargeordered domains, in agreement with the tentative phase diagram proposed for cuprates in Ref. [8].Such a phase diagram was purely based on the peculiarities of the resistance curves as a function of the temperature, with varying magnetic field and doping, and showed that SC can develop at low temperatures even when at high temperature the system is well inside the charge-ordered region of the phase diagram.In this work, we provided a solid background to the above scenario, showing that within our minimal model for the competition between CO and SC, Eq. ( 1), a random magnetic field has exactly the effect of promoting the fragmentation of the charge-ordered state into domains exhibiting the two different realizations of CO.In the domain wall, they frustrate each other resulting in the stabilization of the superconducting state.Once FSC is suppressed (by increasing the temperature and/or the non-thermal parameter α), only polycrystalline CO remains.The polycrystalline CO is characterized by large puddles with different realizations of CO, and it resembles the complex landscape of charge-ordered domains experimentally observed in cuprates [64].
The filamentary superconducting foot we find in the two-dimensional phase diagram (Fig. 6) is relatively small compared with experiments (Fig. 2).Also, the CO temperature is strongly suppressed close to the O(3) point in the presence of quenched disorder.Both features are cured considering an interlayer coupling, yielding a phase diagram nicely resembling experiments, both in the CO-and SC-driven case.
Very near the O(3) point in the dirty case, we find that the superfluid stiffness has a non-monotonic behaviour as a function of T with a maximum at intermediate temperatures (Fig. 7c).It would be very interesting to observe this effect as it would be a signature of entropically favoured SC.
The presence of some one-dimensional-like superconducting patterns persisting on the CO side of the phase diagram can have a striking effect on the macroscopic observables, such as specific heat [7] or spin susceptibility [65], and particularly on transport measurements [8,25,26].

A.1 Effect of the barrier potential B
To analyse the tendency to order in the in-plane (SC) or out-of-plane (CO) direction, we compute the mean-square magnetization, where m ν = 1 N R S ν R is the magnetization per unit surface area calculated at each MC step.Note that the mean-square magnetization is directly related to the charge (ν = z) and superconducting (ν = x, y) susceptibilities [44], In the absence of long-range order 〈m ν 〉 = 0, therefore, for a BKT system in the thermodynamic limit, χ ν = N χ ν /T .Of course, in numerical calculations, the system is always finite and never reaches the real thermodynamic limit, preventing the vanishing of 〈m x, y 〉.In the following, we will use the quantity χ x y ≡ 1 2 ( χ x + χ y ), to monitor the superconducting correlations.While χ x y can be seen as a proxy of the superconducting susceptibility, in the charge-ordered phase, instead, the order parameter is nonzero, so that χ z and χ z are not simply proportional.In this case, to monitor the response of the CO correlations, we use both the susceptibility χ z , which is the true response of the system to an external field, and χ z .
From the superconducting χ x y and charge-ordered χ z functions at different values of the barrier one can have a first insight about the stabilization on the BKT state operated by the potential barrier.In Fig. 12 we show χ x y and χ z functions, defined as in Eq. (A.1) for different values of B > 0. When the barrier is present (B > 0), we observe a sizable χ x y .The temperature at which the superconducting response significantly rises is an increasing function of B (see Fig. 12a).At lower temperatures, the mean-square magnetization tends to a finite value, indicating the stabilization of superconducting correlations.Indeed, contextually, χ z  presents a peak and is driven to zero at low temperatures.This behaviour is characteristic of the bare XXZ model with α < 1, i.e., in the superfluid region of the phase diagram.In the presence of the barrier, we find that the same results also persist for a small range of α > 1.
Thus, the effect of the barrier is to shift the bicritical point (α B , T B ) to α B > 1.However, as discussed in the main text, at T = 0 the superconducting and the charge-ordered phases are degenerate at α = 1.This implies that the first order line T CO↔SC , which by definition starts at (α = α * ≡ 1, T = 0) and ends at the bicritical point, must have a positive slope.This indicates that for a small range α ≳ 1 and lowering the temperature, one has the sequence of phases: disorder → SC → CO.Thereby, two spinodal lines starting from (α B , T B ) and terminating at points α * CO,SC (B) and T = 0 appear.
A.2 B = 0.2: Mean-square magnetizations and susceptibilities In this section, we discuss the case B = 0.2.In Fig. 13 we plot the superconducting χ x y (panel a) and charge-ordered χ z (panel b) mean-square magnetization, as well as the susceptibility χ z (panel c), for different values of the anisotropy parameter, in the range 0.1 < α < 2.
The superconducting mean-square magnetization χ x y grows monotonously by lowering the temperature for values of the anisotropy parameter as large as α = 1.04, i.e., above the isotropic Heisenberg limit.For the same range of anisotropy, χ z shows a maximum and then drops nearly to zero at lower temperatures.Clearly, thus, in this region the superconducting phase prevails at low temperature, so that, at some temperature T BKT , the system transitions from a high-temperature disordered state to a superconducting state.
For α > 1.04, the situation gets reversed with the charge correlations growing monotonically and the superconducting ones getting suppressed.This behaviour is coherent with the results found in the bare XXZ model (B = 0), where χ z decreases with T for α ≪ 1 while for  α ≲ 1 it displays a peak at T ≃ T BKT [44], as a precursor of the Ising transition that is found for α > 1.
Observing the yellow curve corresponding to α = 1.04, in Figs.13a-c, one can understand the importance of considering all three quantities χ x y , χ z and χ z .As a matter of fact, the CO susceptibility χ z in Fig. 13c presents a peak at α = 1.04, although smeared with respect to the peaks for α ≥ 1.05; concurrently, χ x y , Fig. 13a, shows that a superconducting state is present at α = 1.04.The doubt about whether the system has a superconducting or charge-ordered ground state is solved by looking at χ z , Fig. 13b, in which the α = 1.04 curve grows with T following the typical behaviours of the charge-ordered states α ≥ 1.05, but then decreases below a temperature T = 0.65, at which χ z ∼ 37.The BKT scaling of J s for α = 1.04 allowed us to extract T BKT = 0.575 ≲ T .
As soon as α ≥ 1.05, the main response of the system is in the out-of-plane direction (CO), as it is clear looking at χ z , Fig. 13b, and at the susceptibility χ z , Fig. 13c.

A.3 B = 0.2: First order transition and spinodals
The first-order transition line (T SC↔CO , marked by the yellow squares) and the spinodal lines (T sp red diamonds) in the phase diagram are obtained by constructing the effective distribution function P eff (m z ), as discussed in Sec. 3.
Since the flip of the whole phase is a very rare event, we need to take a very small system in order to have enough flips to consider the system at equilibrium within a reasonable simulation time.As a proof of principle and in order to have an approximate map of the spinodal lines, we take L = 4 and construct P eff (m z ) using histograms of m z measured at each MC step.The system is evolved for N MC = 5 × 10 5 after N out = 5 × 10 5 , with τ MC = 50.
We report in Fig. 14a the minima of the free energy F min (T ) as a function of the temperature, for the case α = 1.04,where the superconducting state is marked in red and the charge-ordered state in green.The crossing point between the two curves is the first-order critical temperature T SC↔CO .We see that, in agreement with our previous discussion, CO is the stable phase at low temperatures, then, with rising the temperature, the system switches to SC and then reaches the disordered state.
In panels b-e of Fig. 14, we report the histograms at the temperatures T = 0.25, 0.35, 0.60, 0.65, where the distribution of m z is in turquoise (left axis) and the corresponding free energy F (m z ) is in magenta (right axis).At T = 0.25, F (m z ) displays three minima, Fig. 14b, the global ones being at m z = ±1 (corresponding to CO).By increasing the temperature, at T = 0.35, the three minima become equivalent, Fig. 14c, while at T ≥ 0.35 the global minimum is at m z = 0 (corresponding to SC).To define the spinodal temperature, at each α, we fitted the data F (m z ) in the region around m z = 0.5, and looked for the temperature at which the free-energy curvature changes from downward to upward.It can be seen, by comparing panels d and e in Fig. 14 how the two minima at m z = ±1 disappear when the temperature is increased from T = 0.60 (panel d) to T = 0.65 (panel e), where the curvature near m z = ±1 appears to be flat.

A.4 B = 2: Phase diagram
In order to emphasize the shifting of α B towards higher values, we consider a barrier parameter B = 2.The search for metastable states and first-order lines within the protocol described in Sec.3.1 becomes harder and harder with increasing B. In the case of large B, we follow a different protocol to numerically estimate the spinodal points.Instead of constructing the free energy landscape, we prepare the system in the two metastable states, i.e., the superconducting metastable state above α * and the charge-ordered one below α * .Starting from a very low temperature (T = 0.001), we compute the superfluid stiffness and the CO square Figure 15: a) Superfluid stiffness rescaled according to Eq. ( 2) (in the label of the vertical axis we defined f (L, L 0 ) = 1 + [2 ln(L/L 0 )] −1 for brevity) for α = 0.5 and b) mean-square charge-ordered magnetization 〈m 2 z 〉 for α = 2 at various L (color code as indicated in the legend).Circles correspond to the usual cooling down protocol; lines stand for the results obtained when the system is heated up starting from its metastable state, being a) all spins parallel and oriented along the z-axis (CO metastable phase), b) all spins parallel on the x y plane (superconducting metastable phase).Note that the temperature at which the system jumps from the local to the global minimum state is not strongly dependent on the system size L. The error bars are calculated using the bootstrap resampling method with 100 datasets and blocks of size 100.magnetization respectively in the charge-ordered and superconducting metastable states.By heating the system via a simulated annealing procedure, we thus define the spinodal points as the temperature at which J s and 〈m 2 z 〉 jump from zero to their finite value, checking that this temperature is not strongly dependent on the system size L. The absence of a significant size dependence can also be viewed as a confirmation of the spinodal points extracted with the free energy protocol.As an example, we report in Fig. 15 the superfluid stiffness, rescaled according to Eq. ( 2), for α = 0.5 (panel a) and the mean-square magnetization 〈m 2 z 〉 signalling CO, for α = 2 (panel b).The dots shown in the plot are computed cooling down the system from a random configuration at a given size L while thick lines stand for the heating up process from the metastable state.As one can see, the jump from the metastable state at low temperature to the ground state is not strongly dependent on the system size L.
The phase diagram in Fig. 16 shows that the BKT line survives for values α ≫ 1, up to α B = 1.325.The two spinodal points at T = 0, calculated according to Eq. ( 6 As in App.A.2, we rely on mean-square magnetizations and CO susceptibility at various α to have a quick and comprehensive view of the re-entrance of the superconducting phase.In Fig. 13 we plot the superconducting ( χ x y ) and charge-ordered ( χ z ) mean-square magnetization (panels d and e, respectively), and the susceptibility χ z (panel f), at different values of the anisotropy parameter, in the range 0.1 < α < 2, for B = 2.As one can see, χ x y is significant for values of anisotropy α ≤ 1.325, well above the isotropic Heisenberg limit.The situation gets reversed as soon as α ≥ 1.35, where the main response of the system is in the out-of-plane direction (corresponding to CO).As in Sec.3.2, the susceptibility χ z shows precursor peaks of the charge-ordered state found at α ≤ 1.35, down to α = 1.3,where a very broad peak can be observed.
The case α = 1.35 highlights again the possibility of having a superconducting state stabilized by entropic effects.Indeed, upon cooling, the superconducting mean-square order parameter, light-blue curve in Fig. 13f, follows the BKT behaviour lowering T , similarly to the curve at α = 1.325 (green curve), down to T = 0.8.By further lowering T to T = 0.775, χ x y drops down by a factor ∼500 and, correspondingly, χ z is increased by a factor ∼ 800.It is worth noting how the peak in χ z for this value of anisotropy is still very smeared as for α < 1.35, thus leaving no doubt about the nature of the ground state.
To fully describe the properties of the anomalous transition found at α = 1.35, we looked at the evolution in temperature of the total density of vortices, given in Eq. ( 3).The Tdependence of ρ V, tot , is shown in Fig. 17a.In the BKT scenario, ρ V,tot is supposed to be exponentially suppressed as the temperature is lowered towards T BKT , as a consequence of the binding of vortex-antivortex pairs.This happens up to α = 1.325 (shown as a benchmark with light colours in Fig. 17a), where the suppression of ρ V, tot coincides with the appearance of a finite J s (light colours in Fig. 17b).The α = 1.35 curve to follow this trend in the hightemperature regime.Crossing the temperature T ≈ 0.8, a sudden proliferation of free vortices is observed.This indicates an anomalous transition from an almost BKT-like superconducting state at high temperatures, turning into a charge-ordered state below T ≈ 0.8, in agreement with the trend found in χ x y and χ z .Note that such anomalous behaviour is also detected by finite-size effects in the superfluid stiffness plotted in Fig. 17b.At high temperatures, the paramagnetic phase seems to be on the verge of undergoing a BKT transition, as it is visible from the tails of J s , while instead at T ≃ 0.76 (vertical dashed line) the system develops CO and J s drops to zero.

B Dirty system: Superfluid stiffness
For the sake of completeness, we show the superfluid stiffness and their BKT critical jump for all the values of α considered to construct our phase diagram (Fig. 6).In Fig. 18 it is possible to observe the validity of Harris criterion when addressing the BKT transition for α < 1, where the disorder leaves J s almost unaffected.The only appreciable effect is seen in the suppression of both the saturating value of J s for T → 0 (see panels a, b and c), which is lowered to 0.75 for α = 0.9 while the critical temperature is only very slightly decreased.This can be appreciated looking at panels d, e, and f, where we show the relative crossing points with the universal critical line 2T /π, indicating with a vertical line the corresponding T BKT .A first consequence of the random field is indeed visible in the smearing of this crossing at α = 0.9, highlighted in grey.
In panels a, b and c of Fig. 19 we present the superfluid stiffness in the filamentary region of our phase diagram, namely α = 1, 1.05, 1.1.The suppression of J s caused by the emergence of the correlated disorder is much more visible here.In particular, we highlight the fact that, while for α = 1, 1.05 the scaling law still produces reliable results (see panels a and b of Figs.18 and 7), this does not seem to be the case for α = 1.1 (see panels a and b of Figs.19  and 7).However, the extrapolated T BKT is consistent with the minimum found for its derivative ∂ J s /∂ T [42,61].This is shown in panel d where the vertical line is at T BKT and the grey area highlights the estimated error.T BKT with its error, extracted using the BKT scaling law (see Fig. 7).

Figure 2 :
Figure 2: Phase diagrams showing competition between SC/superfluidity and CO/solid phases.a) YBa 2 Cu 3 O y with y = 6.67.The onset of superconducting correlations (labelled 2D-SC) was detected by an anomaly in the specific heat (adapted from Ref. [7]).Three-dimensional CO (3D-CO) was detected with X-ray diffraction[3] and sound velocity[4].b) La 2−x Sr x CuO 4 with x = 0.08.The onset temperatures for SC (T ONS ), FSC (T ONSF ), and CO (T CO ) were extracted from magnetoresistance data (adapted from Ref.[9]).The darker blue region corresponds to zero sheet resistance R □ = 0.The red line R Q shows the locus of the quantum of resistance.c) La 1.85− y Nd y Sr 0.15 CuO 4 .Here, isoelectronic doping ( y) favours the LTT, which stabilizes CO in the form of stripes[18].Open circles show T c while solid circles show the structural transition T d which is known to be close to the CO temperature (adapted from Ref.[19]).d) Phase diagram of 4 He (adapted from Ref.[21]).In a)-c) the wavy-shaded regions denote the coexistence of CO and SC which in Ref.[8,9] was attributed to filamentary SC.In d) the wavy shading is the 4 He-analogous supersolid phase observed in Ref.[22].

BFigure 4 :
Figure 4: Phase diagram in the T vs. α plane for the competition between SC and CO, modelled with an XXZ model with a barrier height B = 0.2.The light blue lines refer to the bare XXZ model, for comparison.Cyan dots are the T BKT points calculated with the BKT scaling law [Eq.(2)]; purple triangles refer to the Ising (CO) transition, for which the T CO points are found by using the Binder cumulant U N [Eq.(4)]; the first order transition line T SC↔CO (yellow squares) and the spinodal points T sp (red diamonds) are computed from the effective free energies F (m z ) and the locations of its minima F min (T ) (more details in App.A.3); green points at T = 0 are calculated analytically [Eq.(5)].

B = 0. 2 ,Figure 6 :
Figure6: T vs. α phase diagram for B = 0.2 and W = 5 (light blue, purple and blue symbols and lines, left axis) and average in-plane component 〈sin ϕ〉 at T = 0.001 (green circles, right axis); the complete temperature dependence of 〈sin ϕ〉 can be found in Fig.9b.T BKT points (cyan) are computed using Eq.2; T CO points (purple) are computed from the linear fitting of 1/ξ z ; short-range SC points (blue) refers to the temperature at which J s for L = 16 crosses the critical line 2T /π.The errorbars are calculated from the standard deviation of independent disorder configurations.The grey symbols show the phase diagram of the clean system (W = 0), for comparison.

Figure 9 :
Figure 9: a) Superfluid stiffness J s at α = 1.15,B = 0.2, and W = 5, for different system sizes.MC parameters are N MC = 2 × 10 3 , τ MC = 100, N out = 2 × 10 5 ; N dis are 25 (L = 16), 20 (L = 32, 64) and 15 (L = 128, 256).The errorbars are calculated from the standard deviation of N dis independent disorder configurations.b) Average in-plane component 〈sin ϕ〉 as a function of temperature for various anisotropies α at L = 256.Crosses are T BKT points and hexagons are the T CO points.The red dot highlights the short-range SC in both panels.

15 Figure 10 :
Figure 10: superconducting correlation length ξ x y as a function of temperature for L = 32, 64, 128, 256 at various anisotropy parameter α = 1, 1.05 1.1, 1.15.The vertical dashed black lines signal the critical temperatures T BKT , the dashed red line in the panel α = 1.15 at T = 0.2 is the same temperature indicated (also in red) in Fig. 9; the grey dashed line at T = 0.3 correspond to the maximum of J s for L = 16.

Figure 12 :
Figure12: a) Superconducting χ x y and b) charge-ordered χ z mean-square magnetizations vs. temperature T for the isotropic case α = 1 at various heights of the potential barrier B. While for the bare Heisenberg model (B = 0) no transition is possible, the presence of a barrier allows for a BKT transition at α = 1.We used L = 128, N MC = 2 × 10 3 , N out = 5 × 10 4 and τ MC = 100.The error bars are calculated using the bootstrap resampling method with 100 datasets and blocks of size 100[66].

Figure 14 :
Figure 14: Effective free energies and probability distributions of m z for a system of linear size L = 4 with α = 1.04, b = 0.1.a) Local minimum of the free energy as a function of temperature F min (T ).The red curve corresponds to the minimum F (m z ≈ 0) while the green line corresponds to the minimum F (m z ≈ 1).The crossing temperature between the two lines marks the first order transition T SC↔CO .b-e) Effective probability density function P(m z ) (turquoise) and free energy F (m z ) (magenta) for: b) T = 0.25, c) T = T SC↔CO = 0.35.d) T = 0.60, and e) T = T sp = 0.65.The free energies F (m z ) at each temperature were constructed from the distribution of m z within N MC = 5 × 10 5 , τ MC = 5 × 10 5 and τ MC = 50.

Figure 16 :
Figure 16: Phase diagram in the T vs. α plane for the competition between SC and CO, modelled with an XXZ model with a barrier height B = 2.The light blue lines refer to the bare XXZ model, for comparison.Cyan dots are the T BKT points calculated with the BKT scaling law [Eq.(2)]; purple triangles refer to the Ising (CO) transition, for which the T CO points are found by use of the Binder cumulant U N [Eq.(4)]; red diamonds indicate the spinodal points; α * is calculated analytically; the yellow dashed line is a guide to the eye that sketches the expected first-order phase transition.

Figure 17
Figure17: a) The total density of vortices and antivortices ρ V,tot [Eq.(3)] as a function of the temperature, for α = 1.35,B = 2, shows a re-entrant phase as ρ V, tot seems to decay exponentially lowering T down to a temperature T c where vortices suddenly proliferates.T c , marked with the vertical dashed line, was deduced from the Binder cumulant U N .b) The same trend is also confirmed by the finite-size effects in J s .Error bars are calculated using the bootstrap resampling method with 100 datasets and blocks of size 100.The case α = 1.325, showing typical BKT feature, is plotted in lighter colours as a benchmark.

Figure 19 :
Figure 19: Superfluid stiffness for a) α = 1, b) α = 1.05, c) α = 1.1.Error bars are calculated from the standard deviation of independent disorder configurations.Black lines are the universal critical line 2T /π.d) First derivative with respect to the temperature ∂ J s /∂ T for α = 1.1.The vertical line and the grey shaded area indicateT BKT with its error, extracted using the BKT scaling law (see Fig.7).