Numerical tests of the large charge expansion

We perform Monte-Carlo measurements of two and three point functions of charged operators in the critical O(2) model in 3 dimensions. Our results are compatible with the predictions of the large charge superfluid effective field theory. To obtain reliable measurements for large values of the charge, we improved the Worm algorithm and devised a measurement scheme which mitigates the uncertainties due to lattice and finite size effects.


Introduction
Conformal field theories (CFTs) play a key role in physics.CFTs are fixed points of the Renormalization Group flow and encode the universal properties of critical points in second-order phase transitions.Additionally, CFTs shed light on some of the mysteries of quantum gravity through the AdS/CFT correspondence.
Local observables in CFTs can be described in terms of a set of dimensionless numbers (the CFT data), i.e. the conformal dimensions and OPE coefficients of the primary fields of the theory [1,2].When the CFT is strongly coupled, the spectrum of low-dimension operators often lacks an organizational principle; consequently one has to resort to numerical methods, such as the bootstrap or Monte Carlo simulations, to compute the corresponding CFT data.
It has been recently realized that sectors with large quantum numbers are often amenable to a perturbative description.Notably, CFTs become weakly coupled in the large spin sector [3], in the sense that the CFT data can be computed in an expansion in inverse powers of the spin via the analytic bootstrap [4,5].Relevant to our work is the large charge expansion [6,7], which applies to operators charged under the internal symmetries of the theory.Both of these expansions, besides being useful per se, allow identifying interesting patterns in the CFT spectrum, with operators naturally organized in Regge trajectories as a function of their quantum numbers. 1et us review the physical picture underlying the large charge expansion.Consider for concreteness a CFT invariant under an internal U (1) symmetry group in three spacetime dimensions.By the state-operator correspondence, an operator with a large U (1) charge is associated with a finite density state for the theory quantized on the cylinder.In [6] it was argued that, in the simplest case, the corresponding state is found in a superfluid phase.In this case, one can describe large charge states via the universal effective field theory (EFT) description for the hydrodynamic Goldstone mode of the superfluid [7]. 2 This allows for the systematic calculation of correlation functions of charged operators, where the derivative expansion coincides with an expansion in inverse powers of the charge.
The superfluid EFT is believed to describe the large charge sector of a vast class of theories.Nonetheless, sometimes other phases are possible, e.g.Fermi spheres in fermionic theories [14,15] or extremal Reissner-Nordstorm black holes in holographic models.Furthermore, in contrast with the large spin expansion, even in specific theories there is no rigorous bootstrap proof of the validity of the superfluid description. 3It is therefore important to check its predictions in theories where explicit calculations are possible.The main purpose of this paper is to provide evidence for the validity of the superfluid EFT in a specific strongly coupled CFT, namely the O(2) model in three dimensions, via Monte-Carlo calculations. 4

Background and summary
Let us begin discussing some of the predictions of the large charge expansion, with a particular focus on the features which are specific to the superfluid EFT.The scaling dimension D(Q) of the operator O Q with lowest dimension at fixed charge Q is given by [6,7] (1.1) The coefficients c 3/2 and c 1/2 in (1.1) are model-dependent Wilson coefficients, while c 0 ≃ −0.0937.The leading behaviour at large Q follows from dimensional analysis.The existence of an expansion in 1/Q is less trivial, but it is not specific to the superfluid EFT only; for instance, the result for a Reissner-Nordstrom black hole admits a similar expansion (see e.g.[18]).The Q 0 contribution c 0 ≃ −0.0937 is associated with the Casimir energy of the Goldstone and it is thus a specific prediction of the superfluid phase [19].
The superfluid EFT also predicts other observables in terms of the same Wilson coefficients.For instance, it predicts that the primary operator with the next-to-lowest dimension with charge Q has spin 2 and scaling dimension D(Q) + √ 3 [6].Importantly for this work, the EFT also predicts OPE coefficients of the operators O Q [7,[20][21][22], defined as (1. 2) The structure of the prediction depends on whether one considers three large charge operators or two large charge operators and one with a small charge.We refer to these predictions, respectively, as OPE in Regime I and Regime II.The results are • Regime I: the prediction for the OPE of three large charge operators reads [21] λ 3) where Q 1 , Q 2 ≫ 1, and f (y) is given in terms of the solution of a non-linear PDE.We shall only need its value at y = 1, f (1) ≃ 0.996.Notice that this prediction depends only on c3 /2 to leading order in the derivative expansion.
• Regime II: the prediction for the OPE of two large charge operators and a small charge operator is [7,20,22] where ) is a novel Wilson coefficient, associated with the operator matching for O Q 2 in terms of the superfluid Goldstone [7].As in (1.1), the leading scaling with Q 1 can be inferred from dimensional analysis [16].The subleading correction depends on c3 /2 and it is a specific prediction of the superfluid EFT.The numerical coefficient multiplying c3 /2 does not admit a simple analytic expression.It was computed in [20] from the (small) shift of the superfluid saddle-point (equivalent to a tadpole diagram) due to the charge Q 2 sourced by the operator insertion. 5e now discuss former tests of the large charge expansion in CFTs.The validity of the EFT has been unambiguously demonstrated in several perturbative theories, see e.g.[13,[23][24][25]; particularly relevant for us are the results for large charge operators in the O(N ) model in the ε-expansion [26,27] and at large N [28,29] It is of course harder to study directly strongly coupled CFTs.[30] initiated the Monte-Carlo study of the large charge sector of the O(2) model computing the scaling dimension D(Q) for Q = 1, 2, . . ., 12. To perform the calculations the authors applied the worm algorithm [31] to the worldline formulation of the classical O(2) sigma-model [32] (see appendix A for details).Remarkably, the numerical results agree with the theoretical prediction (1.1) up to Q = O(1), providing a determination of the coefficients c3 /2 and c 1 2 from the fit (with c 0 assumed as input).Similar calculations have been performed in the O(N ) models for N = 4 [33,34] and N = 6, 8 [35].
The results of [30] provide strong evidence for the existence of a 1/Q expansion for D(Q).This is a very nontrivial result, but, as commented earlier, it is not necessarily specific to the superfluid EFT (even if it is admittedly hard to think of alternative descriptions for the large charge sector of the O(2) model).The main goal of this work is to test specifically the superfluid EFT by studying OPE coefficients of charged operators.Below we give a brief summary of our results.
First, in sec. 2 we compute the scaling dimension D(Q) for charges up to Q = 19, thus extending the pre-existing results for Q ≤ 12 [30].The results are plotted in fig. 2. Our measurements are compatible with those of [30] and provide improved estimates for the values of the Wilson coefficients c3 /2 and c 1 2 in the O(2) model, cfr.eq.(2.4).Unfortunately, we are not able to reach the precision needed to obtain a reliable estimate for the coefficient c 0 in eq.(1.1); our results are nonetheless compatible with the theoretically predicted value.
In order to test the superfluid EFT, in sec 3 we study the OPE coefficients in eq. ( 1.3) and (1.4).Notice that extracting three-point functions from Monte-Carlo simulations is significantly more involved than computing two-point functions.In Regime I we computed the OPE coefficient for Q 1 = Q 2 = 1, 2, 3, 4, while in Regime II we obtained results for Q 1 = 1, 2, 3, 4, 5 with Q 2 = 1 (fixed).The results are shown in fig.8. Despite the relative smallness of the charges we find good agreement between the numerical results and the EFT predictions, in both regimes.In particular, from the extrapolation of the result in Regime I to larger values of the charges we extract the coefficient c3 /2 , finding remarkable agreement with the value extracted from the measurement of the scaling dimension.The comparison is shown in fig.9. From the results in Regime II we measure the value of the coefficient C(1) in eq.(1.4), see fig. 10.The estimate for c3 /2 extracted from Regime II is encouragingly compatible with the one obtained from D(Q) and the OPE coefficient in Regime I, but uncertainties are too large for our analysis to be conclusive; see fig.11.
Overall our results provide encouraging evidence for the validity of the superfluid EFT in the large charge sector of the O(2) model, but additional data would be helpful to unambiguously confirm the EFT description.In sec. 4 we further speculate on the implications of our findings and comment on possible future directions.
To compute the correlation functions numerically we used the worm algorithm.We introduced two technical improvements with respect to the strategy of [30].First, we introduced the continuous time update step, which reduces the computational time; details are given in appendix A. Additionally, we devised an improved procedure to take the continuum limit.To this aim, we carefully analysed lattice effects, combining numerical experiments and conformal perturbation theory; some details are given in appendix C.

Conformal dimension of lightest charged scalar operator
Measuring 2pt functions of operators with large conformal dimensions is challenging.On the one hand, the 2pt function decays quickly when the distance between the operators increases.On the other hand, measurements at short distances are contaminated by large lattice effects.3), for α = 2 and L ∈ [16,32].Notice that the region where there are significant deviations from a constant value, on the left of the vertical line, gets smaller as L increases.This is expected from lattice effects.We include the previous Monte-Carlo result with error bars at 1σ, taken from Tab.I in [30], as well as the bootstrap result [36].
To make progress, [30] introduced a method that does not require sampling directly 2pt functions.They measure the difference between the conformal dimensions of operators with consecutive charges, ∆(Q) ≡ D(Q + 1) − D(Q), which scales as Q 1/2 instead of Q 3/2 .This is achieved by rewriting the 2pt function This is useful because, due to the worldline formulation of the O(2) model, R q (x), can be sampled directly by computing the expectation value of operators with Q = 1 in a background charge distribution R q (x) = e iθ(0) e −iθ(x) where e iθ (e −iθ ) represents a charge 1 (−1) operator in the nonlinear sigma model and the subscript indicates that the expectation value is computed in the presence of charge (q − 1) at the origin 0 and charge −(q − 1) at position x.Check App.A.2 for details.
Having reconstructed the 2pt function, the next step is to extract ∆(Q).The naive approach is to compute R q (x) for different values of x, take the log and fit the slope.Unfortunately, this cannot be done systematically due to both finite-size effects and lattice effects.Lattice effects are due to the discrete nature of the lattice.This introduces another distance scale, the lattice spacing a, such that in the region where x/a ∼ O(1), the discrete Worm steps), as well as the results of the extrapolation to L = ∞, obtained using α = 2 in eq.(2.3).We test the leading behaviour of eq.(1.1) by demonstrating the linearity of D(Q)/Q 1/2 .We reproduce the previous results of [30].As explained at the end of the section, we performed best-fits of the coefficients c3 /2 and c 1 2 in eq.(1.1), restricting to charges Q ≥ Q min for different choices of Q min .The blue dashed curve is obtained using the coefficients extracted from the average of all the fits with Q min = 4, . . ., 10.The orange dashed curve is the best fit obtained using only charges Q ≥ 8. Remarkably, the two lines are almost indistinguishable and almost overlap with all the data points.In the table below we list the data used in the fits presented in this section.Check app.B for a detailed discussion.
nature of the lattice spoils the CFT predictions.We set a = 1 in the following unless specified otherwise.In summary, there is an intermediate region, where x/a ≫ 1 and x/L ≪ 1/2, such that the continuum infinite size CFT predictions hold.As we show next, we are able to drop the second restriction through a choice of observable that eliminates finite-size effects.
Systematic errors coming from finite size effects are eliminated by computing the ratio between 2pt functions measured for different lattice sizes but at the same relative position 3) The right-hand side holds as long as lattice effects are negligible.Notice that these relations are independent of the position x.Both ratios are insensitive to finite size effects since these are parameterized by the relative position x/L.Deviations from a constant value (as a function of x) are a proxy for lattice effects.
In fig. 1 we present the measurements of ∆(1), obtained through eq. ( 2.3).For small values of x there are deviations as expected.These are the lattice effects.They disappear when x ≈ 7, identified by the vertical lines.In particular, we can also reliably measure the correlation function for distances x < L/2, which we will use to optimize our measurements.In the constant region, to the right of the vertical lines, we have an independent estimate of ∆(1) for each position.From these, we estimate the error bars on the final result.
The generalization of this analysis to different values of Q provides reliable measurements of ∆(Q), which are not susceptible to finite-size effects and do not require fits.Similarly, we also obtain accurate estimates of the errors introduced by lattice effects.The measurements of ∆(Q) are independent of the lattice size.We explicitly checked this for larger charges.The results for L = 32 and L = 64 match within the statistical uncertainties, but L = 32 has a smaller error 6 and we focused on this system size.
We measured D(Q) up to Q = 19, see fig. 2. As remarked in the introduction, this represents a considerable improvement with respect to the existing results, which stopped at Q = 12 [30].To obtain results for such high values of the charge we used a continuoustime update step (see App. A).Additionally, we sampled the ratio in eq. ( 2.3) for relatively small values of the distance x7 , while [30] performed all measurements for x ∼ L/2.Indeed, as explained earlier, lattice effects are negligible already for x ≳ 7, with the most precise measurements obtained for 7 ≲ x ≪ L/2. 8et us now discuss the comparison with the theoretical prediction eq.(1.1).In doing so, we face some important questions: What is the theoretical error of the large charge expansion?Is this error also under control for small values of Q?
The large charge expansion is believed to be an asymptotic expansion [37].The series in eq.(1.1) thus includes both perturbative terms, suppressed by inverse powers of Q, 9 as well as non-perturbative corrections, which are exponentially suppressed at large charge ∼ e −# √ Q . 10Most importantly, as typical with asymptotic expansions, the series is not expected to converge to the exact result upon including infinitely many terms; rather, the large N analysis of [37] suggests that the large charge expansion of the scaling dimension D(Q) admits an optimal truncation after n ∼ √ Q terms.
While the subtleties associated with the asymptotic nature of the series are unimportant for very large charges, they make it challenging to estimate the accuracy of the   expansion for our data.In particular, we do not expect the theoretical error to be a simple function of Q when Q ∼ O(1).
To (partially) account for these effects, we performed fits of c3 /2 and c1 /2 in eq.(1.1) using charges in the range [Q min , 19] for different values of Q min .We indeed expect that the theoretical error decreases with Q.The results are shown in fig. 3. The results are independent of Q min for Q min ≳ 3, suggesting that the truncated asymptotic expansion is trustworthy beyond this value of the charge.By averaging over the results for Q min ∈ [4,8] we obtain c3 /2 = 0.339(1) c1 /2 = 0.25(1). (2.4) These values are compatible with the previous estimates of [30] c3 /2 = 0.337(3) and c1 /2 = 0.27(4).
In fig. 4 we plot the difference between the numerical data and the best-fit curve obtained using the averaged parameters.In the same plot, we also show the deviation for the best-fit curve obtained setting Q min = 8.For small charges, there are systematic deviations, but they become smaller than the numerical uncertainties for Q ≳ 4.This justifies a posteriori the choice of fitting in the range Q min ∈ [4,8].It is remarkable that also for small charges the relative deviations are rather small.For instance, the best-fit curve obtained for Q min = 8 extrapolated to Q = 1 agrees with the measurement within a 4% relative error.
Notice that we did not try to fit the value of c 0 , which we held fixed at its theoretical value c 0 ≃ −0.0937.Indeed, as we argue below, the data are compatible with this value, but the increasing numerical uncertainties with the charge make it impossible to obtain a reliable estimate for c 0 or other subleading coefficients.
To justify the compatibility of c 0 = −0.0937with the Monte Carlo data, it is convenient to define the following quantity where The quantity A(Q) is so defined to be independent of c3 /2 and c 1 2 when evaluated using eq.(1.1).Its 1/Q expansion reads 11 A similar sum rule was introduced in [6].
The results for A(Q) are presented in fig. 5. On the one hand, measurements for large charges do not achieve a sufficient level of precision to extract c 0 .On the other hand, the results for small charges, for which the precision is high, are subject to unknown theoretical errors. 12Our results are nonetheless compatible with the theoretical value.

OPE coefficients
Let us now discuss how to measure the OPE coefficients (1.2).We consider in particular the OPE coefficient in regime I (cfr.eq.(1.3)) for 11 Here we restored two-subleading orders in the expansion of the scaling dimension D(Q): (2.7) The Q −1 term in this expression represents the contribution to the Casimir energy from higher derivative corrections to the Goldstone dispersion relation, see [20] for details. 12 Perhaps relatedly, the extrapolation of the large N analysis of [37] suggests that non-perturbative terms in the series (1.1) are of the same order of c0 for Q ∼ O(1).where f (1) ≃ 1 and we included the first two subleading terms, which are multiplied by two unknown coefficients α1 /2 and α 0 , for future reference. 13For the OPE in regime II, we take where we also included an extra-subleading term, which depends upon a new Wilson coefficient β −1 . 14For the sake of concreteness, in the following we discuss how to measure the OPE coefficient (3.1).A similar discussion applies to the OPE in regime II.
In Monte Carlo simulations operators are normalized differently than in the CFT literature.In Monte Carlo, 2pt functions at coincident points are normalized to 1, while in the CFT literature 2pt functions are normalized to 1/|x| 2D(Q) asymptotically.Therefore, to extract the OPE coefficient as defined in eq.(1.2), we measure a suitable ratio between the 3pt function and 2pt functions.The ratio of interest is where on the right-hand side we expressed it in terms of the OPE coefficient (1.2), using the continuum infinite size CFT prediction.
In order to measure 3pt functions with the Worm algorithm, we need to rewrite them as some combination of 2pt functions in the presence of background charges, as in eq.(2.2). 1513 The Q 1/2 term depends upon a subleading Wilson coefficient of the EFT which does not contribute to the scaling dimension D(Q); therefore we cannot compute its value from the estimates obtained in sec.2. The Q 0 term, α0, is instead independent of Wilson coefficients, analogously to the c0 term in eq.(1.1).In principle, its value could be computed from the one-loop fluctuation determinant around the saddle-point of [21].In practice, this calculation is technically challenging and we treat α0 as an unknown parameter.We, therefore, write the OPE coefficient as: The task of measuring OPE coefficients reduces to the measurement of 2pt functions in the presence of background charges; these can be efficiently sampled in terms of ratios (as in eq.(2.1)) using the strategy outlined in the previous section.Further details are given in App.A.2.We now discuss lattice and finite-size effects.First, it is useful to determine the region where lattice effects are negligible.To this aim, we consider the ratio of three-point functions at different lattice sizes but at the same relative position: where we defined the 3pt function on the lattice as the value for the exponent γ Q in eq.(3.5) follows from scale invariance in the CFT.Analogously to eq. (2.3), the ratio (3.5) is independent of x when lattice effects are negligible.
In fig.6 we plot the difference between the exponent γ Q extracted from numerical measurements for different values of x and the CFT prediction in eq.(3.5).We show explicitly the results for Q = 1, 2, 3 (the plot for Q = 4 is analogous) and L = 32.The plot clearly shows that the region where lattice effects are negligible decreases with the charge.This makes it challenging to perform measurements for large values of the charge.Differently than with 2pt functions, it is not possible to eliminate completely systematic errors from finite size effects. 16Therefore, to extract the OPE coefficient we fix x and study the L dependence of the ratio of 3-and 2-pt functions in eq.(3.3).The idea is that, as long as x is outside the region where lattice effects are relevant, the extrapolation to L → ∞ is unaffected by lattice corrections.Moreover, in the limit x/L → 0, at fixed x, the lattice 3pt function should be well described by the continuum infinite size prediction.Thus, this provides a direct measurement of the OPE coefficients.We show the results for Q = 1 in fig. 7.While for x ∈ [1, 2, 4, 6] there are deviations, all results for x > 6 converge to the same value, within uncertainties.The error bars for the final result are estimated from the dispersion of the intercept.
For larger values of Q one obtains similar plots, but the data has larger uncertainties and the extrapolation for L → ∞ does not converge as nicely.This is due to the increased significance of lattice effects, as shown in fig.6.To improve the precision and accuracy, we used conformal perturbation theory to parameterize both lattice corrections and finite size effects in eq.(3.3).We obtain an expression as a series in powers of a/x, a/L and L/x, see App.C for details.We improve on the naive linear extrapolation by fitting the coefficients of these powers.Notice that in this approach we perform a unique multidimensional-fit with all the data (i.e. for all values of x and L), rather than performing separate linear extrapolations for each value of x.
The numerical values for the OPE coefficients in regimes I and II are shown in fig.8.We show both the results of the linear extrapolation, as well as those obtained from accounting lattice and finite size corrections, as explained in the previous paragraph.Both methods yield the same results for the OPE coefficients Regime II.Also for the OPE coefficients in regime I the two measurements are compatible, but the results accounting In regime I, the OPE coefficients extracted from the linear extrapolation are in cyan, and OPE coefficients extracted using lattice corrections are in dark blue.In regime II, both methods yield the same results.The OPE coefficient for Q = 1 is the same in both regimes.Notice that the OPE coefficient in regime I grows much faster with Q than the one in regime II; this is in qualitative agreement with the EFT predictions.The data in this plot is in the table below.for lattice corrections, in dark blue, are consistently larger than those obtained from the linear extrapolation, in cyan.This suggests that our results for λ Q,Q,−2Q may be underestimated.Unfortunately, to further investigate this issue we would need to perform more precise simulations at larger distances, which are currently beyond our reach. 17Our results are tabulated in the table in fig.8.The range of charges that can be sampled in regime I is limited, as the statistical errors grow quickly with the charge.In regime II, it would be possible to go further, but this would not improve the precision with which we can measure the coefficient c3 /2 , as we will explain in the following paragraphs.
Let us now analyze our data.To this aim, for the OPE coefficient in regime I, we perform two separate fits: one for the first three coefficients in eq.(1.2) and all data points, and one for the first two Wilson coefficients only and the results for Q ≥ 2. The results are shown in Tab. 1.
Clearly, our analysis is limited by the small number of data points.Nonetheless, the estimates in table 1 suggest that the leading Wilson coefficient should lie in the range c3 /2 ≈ 0.4 ± 0.1, which is compatible with the value c3 /2 ≈ 0.34 obtained from the measurements of scaling dimensions in the previous section.To appreciate this point better, in fig. 9 1: Result for the fits of the OPE coefficient in regime I with different free parameters and for different values of the charges in eq.3.1.We show the fits with both the data obtained from the linear extrapolation and with the inclusion of lattice corrections; notice the latter are more precise.show our results for log(λ

Conformal dimension Linear extrapolation Lattice corrections
The value of c3 /2 can be extracted from the asymptotic behaviour of this quantity.The value obtained with the linear extrapolation is smaller than the one with lattice corrections.The range of charges is small hence it is difficult to extrapolate to infinite charge.Nevertheless, the results are compatible with asymptotically approaching c 3/2 ≈ 0.34.We now discuss the results for the OPE coefficient in regime II.We begin by testing the leading behaviour of the OPE coefficient.In fig.(10), we show λ 1,Q,−Q−1 /Q D(1)/2 .If the EFT prediction eq.(3.2) holds, this ratio should approach a constant for large Q.The plot shows that this is indeed the case.
Having checked the leading behaviour, we focus on studying the sub-leading correction and measuring c3 /2 .Unfortunately, fits are not very precise due to large correlations between C(1) and c3 /2 in the range of charges available.For reference, the fits including β −1 1.14(6) -0.2(3) 0.18(7) 1.02(2) 0.47 (5) --Table 2: Result for the fit of the OPE coefficient in Regime II, with different parameters, in eq.(3.2).  and excluding the subleading coefficient β −1 (cfr.eq.(3.2)) are given in table 2. Clearly the fit including the coefficient β −1 makes the uncertainties too large for the results to be meaningful.The second fit is compatible within 2σ with the estimate c3 /2 ≈ 0.34.

Conformal Dimension Monte Carlo
To study directly the first sub-leading contribution in eq.(3.2), we consider the following ratio (3.8) The EFT predicts that eq.(3.8) asymptotes to c3 /2 for large Q.Our numerical results are shown in fig.11.They are compatible with the value of c3 /2 obtained in sec.2, represented by the black dashed line.However, similarly to the analysis c 0 in the previous section, the uncertainties increase rapidly with Q, making it impossible to obtain a reliable estimate.These uncertainties also represent an obstacle towards improving our results with measurements of the OPE coefficients at larger values of Q.

Conclusions and outlook
In this work, we used Monte-Carlo calculations to test the validity of the superfluid EFT for describing the large charge sector of the O(2) model.Our results were already summarized in the introduction.Here we instead discuss the implications of our findings and potential future directions.
The most surprising aspect of the result for D(Q) is the effectiveness of the large charge predictions also for Q = O(1).For instance, the extrapolation of the best-fit curve obtained from the data with Q ≥ 8 reproduces the measured value of D(1) with 2% accuracy.Overall, our results confirm that this phenomenon persists also for OPE coefficients.A partial justification for the accuracy of the large charge expansion for the scaling dimension was given in [37] via resurgence analysis of the large N result of [28].It might be interesting to perform similar analyses for OPE coefficients.
While our results are encouraging, the scarcity of data points does not allow us to draw unambiguous conclusions about the validity of the superfluid EFT.It is therefore important to obtain more data.Unfortunately, obtaining OPE coefficients for higher values of the charges is beyond the reach of our current Monte-Carlo algorithm.For instance, we estimate that obtaining the OPE coefficient in Regime I for Q 1 = Q 2 = 5, with an uncertainty of 10%, would require 10 CPU years.
A simpler target might be the Monte-Carlo calculation of the scaling dimension of the lightest charged operator with spin J = 2. 18 Notice that the cubic symmetry group of the lattice naturally allows to represent operators up to spin J < 4. 19 As commented in the introduction the superfluid EFT predicts the scaling dimension of the spin 2 operator to be D(Q) + √ 3.However, it is unclear whether one should expect this result to converge for small values of Q, as it happens for scalar operators.Indeed it is expected that the large charge sector of the O(2) model admits a rich phase diagram as a function of the ratio J/Q [39,40].We hope to report about progress in this direction in the future.
It would also be interesting to explore alternative methods to compute the spectrum of charged operators in the O(2) model.An intriguing possibility is provided by the fuzzysphere regularization of [17], which allows directly computing the spectrum of the theory on the cylinder. 20We were also informed of interesting numerical bootstrap results for charged operators in the O(3) model [41].
The accuracy of the large charge expansion is reminiscent of the success of the Regge relation for the QCD and Yang-Mills spectrum [42], and more recently of the remarkable results of the large spin bootstrap in the O(N ) models [43,44].In both of these examples, the results are (partially) explained by the analyticity properties of the spectrum [8,45,46]

A Monte Carlo
This section describes the Monte Carlo method and measurement strategies employed.It describes the world line formulation, the Worm algorithm and the procedures required to express the correlation functions as averages that can be efficiently estimated with Monte Carlo.We describe an improved Worm update, the continuous time update, that guarantees the Worm tail always moves.
The lattice Hamiltonian of the O(2) model is where the field θ x is defined at the cubic lattice nodes.Simulations can be performed in this representation.However, it is more efficient to use a world-line representation [32], where the node variables are mapped into edge variables using where I k (β) is the modified Bessel function of the first kind and β is the inverse temperature.We work at the critical temperature of the three dimensional O(2) model, β = 0.4541652 [47].Since each k is associated with a pair (θ x , θ x+aρ ), these live on the edge of the lattice connecting x to x + aρ.After this rewriting, the path integral over θ can be performed explicitly and the partition function becomes a sum over all possible values k for all the edges where the sum is over all possible configurations of edge variables and the product is over nodes, n, and the edges connected to it.The world-line formulation brings two significant improvements.The first is the possibility of using the Worm algorithm [31], which has one of the smallest dynamical critical exponents [48].The second is that correlation functions can be reinterpreted as the partition function in the presence of some background charge where the subscript i (Q i ) x i indicates that the correlation function is computed with the partition function ) is interpreted as charge conservation at each node, the extra term i (Q i ) x i can be interpreted as a source/sink of charge, or in other words, a background charge distribution.

A.1 Worm algorithm
The Worm algorithm consists of two steps.An update step generates new configurations, Alg. 1.Then, a measurement step extracts the desired correlation function, see Alg. 2.
1. Pick a random site x h (head site).Define x t = x h (tail site).
3. Let k be the current flowing through the bond connecting x t to x t + σ ρ.If:

If:
• update is accepted: • update is not accepted: x t = x t .
5. If x t = x h the update ends.Else go to step 2.
Algorithm 2: Measurement step: e iθ(x) e −iθ(y) .1. Start two counters c x = 0 and c y = 0, that count the number of times the head and the tail coincide and the number of times the tail is at y, respectively.2. Define the head and tail sites as x: x t = x h = x.
3. Perform step 2, 3 and 4 of the update step.4. Whenever x t = y, increment the counter c y = c y + 1. 5. When x t = x, increment the counter c x = c x + 1. Otherwise go back to step 3. 6. Repeat the previous steps N times.The expectation value is given by cy /cx.
In the measurement step, the head of the Worm transports charge 1, and it generates configurations with charge insertions at its head and tail.Then the ratio between the number of times the head and tail are at the positions at which the correlation function is being measured and the number of times the head and the tail are at the same position is an estimate of the correlation function.This algorithm can be improved by modifying step 3 and removing step 4. Instead of uniformly choosing a direction and then choosing to accept it, we can immediately move in a given direction with a probability that makes it equivalent to choosing a direction and then accepting or not that direction.We denote this algorithm by continuous time update, contrary to the Metropolis-like update of the original algorithm.This is achieved by choosing to move through a given edge with probability P (σ ρ), given by where P (σ ρ | n) is the probability of being at position n, proposing to move in the direction σ ρ and accepting the proposal (as described in steps 2 and 3 of the update step).Since k P (σ ρ) = 1, the tail always moves.If there is a counter associated with the position n, then its value should be incremented by 1 /( σ ρ P (σ ρ|n)), the expected time, in the original algorithm, the tail is at n before moving.In essence, we are replacing a stochastic step in the algorithm with its exact solution, which reduces the statistical errors.This update can also be understood as the heat bath step, from which the detailed balance follows.The comparison with the previous algorithm, shown in fig.12, demonstrates that the continuous time update represents a significant improvement over the standard algorithm.

A.2 Ratio between correlation functions
First, let us show how operator insertions can be interpreted as background charges.Consider the correlation function e iQθ(x) e −iQθ(y) .Then, by expanding the definition of the expectation value and using the identity the path integral over θ can be performed analytically, yielding e iQθ(x) e −iQθ(y) = ⟨i,j⟩ ⟨i,j⟩ where D i ≡ ρ∈ 1, 2, 3 k i,i+aρ − k i−aρ,i and ⟨i, j⟩ is the sum over first neighbors.
We are now ready to study ratios between correlation functions.Consider the ratio appearing in eq.(2.1) To implement this on the lattice, it suffices to generate an initial configuration of link variables satisfying at every point and then perform the standard worm update.This can be generalized to higher point functions.The general rule is that charge insertions in the denominator are removed from the numerator and added to the background.Thus, the 3pt function in eq.(3.4) can be rewritten as

B Finite size scaling analysis
In this appendix, we study the dependence of the numerical measurements of the conformal dimension on the size of the system.The data presented in the table below fig. 2 are obtained using the procedure described here.In fig.13, we show the numerical results for the differences ∆(Q) between conformal dimensions for L = 32, 48 and 64 and their extrapolation to L = ∞, and present the comparison with the available results in the literature.For Q = 1, 2, 3, 4, finite-size effects are relevant and bigger than statistical uncertainties.In particular, the measurement for Q = 1 and L = 32 is incompatible with the bootstrap result [49] and with previous Monte Carlo results [50,51], while our extrapolated value is compatible.
Figure 13: Extrapolations of the difference between consecutive conformal dimensions to L = ∞.In blue we present our data for L = 32, 48, 64.In orange we present the linear extrapolation to L = ∞ and in red the results in [30].In black, we present the bootstrap results [49].In green, we present the results from [50], for Q = 1, and [51] for Q = 2, 3, 4.Not all of the results exist for all the charges.
As charges become larger, systematic errors become less relevant and, for Q = 7, the results for L = 32 are compatible with the extrapolated results, within the uncertainties of the latter.
In light of the above analysis, we made the following choice for the data presented in tab.3: for charges Q ≤ 7, we use the extrapolated data and uncertainties; for Q > 7 we use the data for L = 32 with doubled uncertainties.We made this choice because, for L = 64 and charge Q ≳ 10, there are systematic errors due to the small statistics.Such systematic errors prevent us from obtaining reliable extrapolations to L = ∞.From fig. 13 we observe that the statistical errors are roughly comparable with the systematic ones for Q ≥ 7, therefore justifying the choice of doubling the statistical uncertainties to estimate the overall uncertainty.

C Lattice corrections
The analysis we perform in sec.3 relies on the hypothesis that both lattice and finite size effects are small.However, we observe significant lattice effects for large charges, see fig. 6, and thus linear extrapolation is no longer enough.Given that we are already using the largest lattice size and distances that we can reasonably simulate, we now study lattice and finite size effects.We start by identifying the projection of a lattice operator into continuum operators.Next, we identify neutral scalar irrelevant operators that can be added to the action.These need to be irrelevant otherwise the system flows away from this fixed point.Their goal is to encode the information about the existence of a lattice.Finally, our space has the topology of a torus, such that the 2pt functions are not fully fixed by symmetry.Thus, we only have access to correlation functions when operators are close, |x − y| ≪ L, and the OPE quickly converges.Since this is related to the UV behaviour of theory, it is not sensitive to the global topology of the system.
For the sake of simplicity, we will only present explicit computations for the 2pt functions with Q = 1.Generalizations to other charges or higher-order correlation functions are straightforward.We use the standard CFT notation.To make a connection with the rest of the paper, the reader should keep in mind that O lat (x) = e iθ(x) , ∆ ϕ = D(1).
The lightest charged scalar lattice operator can overlap with all operators that have the same charge and are scalar under the cubic subgroup (i.e.operators whose spin s ≡ 0 mod 4).We will just consider the two lightest operators ϕ and □ϕ, check Tab.The correlation functions on the right-hand side are computed on the unperturbed CFT.
In flat space, the 2pt and 3pt functions are known.In the torus, they are not and the only tool available is the OPE.This means we can only study the short-range behaviour of these correlation functions.In the OPE of ϕ † × ϕ we will only include the lightest neutral scalar operator 22ϕ † (x) × ϕ(0) ∼ |x| −2∆ ϕ I + |x| ∆s λ ϕ † ϕs s(0) .which depends on vacuum expectation value (VEV) of s and s ′ and in the integrated 2ptfunction on the torus of s and s ′ .These are unknown in general, hence, instead of focusing on computing them explicitly, we extract their dependence on the dimensionful parameter L. Let us go case by case: • ⟨I⟩ CFT = 1, by definition of the identity operator.
• ⟨s (0)⟩ CFT = α 1 L ∆s , where α 1 is a dimensionless parameter.L appears raised to the power of the conformal dimension of s since it is the only dimensionful parameter available 23 .
• T 3 d 3 z ⟨s ′ (z)⟩ CFT = ⟨s ′ (0)⟩ CFT T 3 d 3 z = α 2 L ∆ s ′ −3 , where we used translation invariance to bring the expectation value of s ′ out of the integral.By acting with □ on this, we obtain the second term in eq.(C.2).
Corrections to the OPE coefficients are obtained using the same ideas, but the derivations are significantly more cumbersome.As such, we only show the end result.Thus, the lattice estimation of the OPE coefficient appearing on the right-hand side of eq.(3.3), here denoted as λ This relation depends on the charges appearing on the left-hand side of eq.(3.3) through OPE coefficients of the type λ Q,−Q,s and multiplicative factors of ∆ Q (these will never appear on the exponents of x or a).The expansion (C.8) is independent of the charges of the operators, up to the undetermined coefficients.

Figure 1 :
Figure 1: ∆(1), extracted with eq.(2.3), for α = 2 and L ∈[16,32].Notice that the region where there are significant deviations from a constant value, on the left of the vertical line, gets smaller as L increases.This is expected from lattice effects.We include the previous Monte-Carlo result with error bars at 1σ, taken from Tab.I in[30], as well as the bootstrap result[36].

Figure 2 :
Figure2: Here we plot both the measurements of D(Q)/Q 1/2 obtained for L = 32(10 7  Worm steps), as well as the results of the extrapolation to L = ∞, obtained using α = 2 in eq.(2.3).We test the leading behaviour of eq.(1.1) by demonstrating the linearity of D(Q)/Q 1/2 .We reproduce the previous results of[30].As explained at the end of the section, we performed best-fits of the coefficients c3 /2 and c 1

2 Figure 3 :
Figure 3: Best-fit values of c3 /2 (left) and c1 /2 (right) as a function of the minimum charge included in the fit.For larger values of Q min the error bars are larger than the plotted range.The coloured regions represent the 1σ interval quoted on (2.4).

Figure 4 :
Figure 4: Difference between the data and the best-fit curves with the averaged parameters in eq.(2.4) (red) and with Q min = 8 (blue).For large values of Q the uncertainties exceed the range displayed in the plot.

Figure 5 :
Figure 5: Results for A(Q), defined in eq.2.5.Notice the large uncertainties despite the precision of the measurements of ∆(Q).

Figure 6 :
Figure 6: Difference between the value of γ Q measured from Monte-Carlo and the theoretical prediction 2D(Q) + D(2Q) (see eq. (3.5)).The result is obtained from the ratio of correlation functions sampled at L = 32 and L = 64, at the same relative position.We used the values D(Q) and D(2Q) measured in the previous section for the theoretical prediction.

Figure 7 :
Figure 7: Linear extrapolation of the OPE coefficient λ Q,Q,−2Q to L → ∞ at fixed x and Q = 1.

Figure 8 :
Figure8: Numerical results.In regime I, the OPE coefficients extracted from the linear extrapolation are in cyan, and OPE coefficients extracted using lattice corrections are in dark blue.In regime II, both methods yield the same results.The OPE coefficient for Q = 1 is the same in both regimes.Notice that the OPE coefficient in regime I grows much faster with Q than the one in regime II; this is in qualitative agreement with the EFT predictions.The data in this plot is in the table below.Q 1 2 3 4 5 Figure8: Numerical results.In regime I, the OPE coefficients extracted from the linear extrapolation are in cyan, and OPE coefficients extracted using lattice corrections are in dark blue.In regime II, both methods yield the same results.The OPE coefficient for Q = 1 is the same in both regimes.Notice that the OPE coefficient in regime I grows much faster with Q than the one in regime II; this is in qualitative agreement with the EFT predictions.The data in this plot is in the table below.Q 1 2 3 45Regime I -linear extrapolation 1.252(6) 1.73(2) 2.46(3) 3.51(7) --Regime I -lattice corrections 1.254(2) 1.77(3) 2.51(6) 3.7(3) --Regime II 1.250(9) 1.40(2) 1.54(4) 1.64(5) 1.70(6)

Figure 9 :
Figure 9:  Comparison between the numerical results for eq.(3.7) and the predicted value c3 /2 = 0.34(1).The black line is the value of c3 /2 obtained in sec.2, and its width is the uncertainty.The coloured full lines are the fits and the coloured regions around their uncertainty.The fits converge to the black line for larger values of Q.

Figure 10 :
Figure 10: Coefficient of leading behaviour OPE coefficient in regime II.
. It remains an important open question whether the CFT data of the O(2) model enjoy similar analyticity properties as a function of the charge.GRID FEUP.They also thank Centro de Física do Porto funded by Portuguese Foundation for Science and Technology (FCT) within the Strategic Funding UIDB/04650/2020.JP is supported by the Simons Foundation grant 488649 (Simons Collaboration on the Nonperturbative Bootstrap) and the Swiss National Science Foundation through the project 200020 197160 and through the National Centre of Competence in Research SwissMAP.

Figure 12 :
Figure 12: Variance of R Q (16) as a function of Q, defined in eq.(2.1).The comparison is made at a fixed CPU time.The histograms show the distribution of R Q (16), as obtained from the measurement step described in alg. 2.