Critical $1$- and $2$-point spin correlations for the $O(2)$ model in $3d$ bounded domains

We study the critical properties of the $3d$ $O(2)$ universality class in bounded domains through Monte Carlo simulations of the clock model. We use an improved version of the latter, chosen to minimize finite-size corrections at criticality, with 8 orientations of the spins and in the presence of vacancies. The domain chosen for the simulations is the slab configuration with fixed spins at the boundaries. We obtain the universal critical magnetization profile and two-point correlations, which favorably compare with the predictions of the critical geometry approach based on the Yamabe equation. The main result is that the correlations, once the dimensionful contributions are factored out with the critical magnetization profile, are shown to only depend on the distance between the points computed using a metric found solving the fractional Yamabe equation. The quantitative comparison with the corresponding results for the Ising model at criticality is shown and discussed. Moreover, from the magnetization profiles the critical exponent $\eta$ is extracted and found to be in reasonable agreement with up-to-date results.


Introduction
The universality class of 3d models with U (1) [or O (2)] symmetry finds broad application in different areas of physics, thanks to its ability to describe various phenomena ranging from high-energy models and lattice gauge theories [1,2] to superfluids and superconductors [3], while at the same time being amenable to be studied by both analytical and numerical techniques.
The most studied lattice counterpart to the O(2) field theory is the XY spin model, which can describe the chiral transition of lattice QCD in the presence of two staggered quark flavors [4,5]. The XY Hamiltonian also provides a model for liquid helium on a lattice [6,7], which allows it to describe the superfluid transition of 4 He across the λ line. This phase transition enjoys extremely accurate experimental studies [8], thanks to a weakly singular compressibility of the fluid and the possibility of performing experiments in space [9], to reduce gravitational effects that would broaden the transition. However, the experimentally determined exponent of the specific heat α (and therefore, the exponent ν as well) differs from the most accurate numerical results, obtained from Monte Carlo simulations [10], by a whopping eight standard deviations. Successive conformal bootstrap results [11] confirm the Monte Carlo values, meaning that this discrepancy is still an open problem. Another interesting phenomenon displayed by the XY model is the appearance of vortices and vortex lines: in 2d, vortices drive the celebrated Berezinskii-Kosterlitz-Thouless transition [12,13], while vortex loops appear in the 3d XY model [14], providing a topological mechanism for its order-disorder transition an additional tool to tackle the problem of the specific heat of helium [15] and duality in (2+1)-d systems at finite temperature [16]. The 3d XY model is such a rich topic that a lot has yet to be discovered: in particular, we are going to focus on the properties of the critical XY model in a slab geometry.
There are many reasons to focus one's attention to systems with boundaries [17,18]: first of all, it is necessary to compare theoretical predictions with experimental results or interpret numerical data, which come from finite systems. Techniques such as finite-size scaling can be used to extract critical exponents and other important quantities. The field is also quite rich: a given bulk universality class corresponds to several surface universality classes. In some cases, boundaries are not just a nuisance, but are rather the source of the phenomenon being considered. The most well-known case is the thermodynamic Casimir effect: it has been studied thoroughly to obtain universal scaling function and Casimir amplitudes [19]. The interest in the Casimir effect has shifted the attention from semiinfinite systems [20][21][22][23] to the geometry of a slab [24][25][26][27][28].
The XY model in a slab, in particular, can model the Casimir effect in Helium near the superfluid transition. Different kinds of boundary conditions can be used to describe the various surface universality classes [29]. In an ordinary transition, the bulk orders in the presence of a disordered surface: this corresponds to Dirichlet boundary conditions and is used to describe a pure 4 He fluid, whose wave function vanishes at the boundary [30]. Enhancing the surface coupling up to a critical value causes the boundary and the bulk to order at the same temperature; this special transition is characterized by Neumann boundary conditions [31]. Finally, if the boundary coupling is enhanced further, or a magnetic field is added at the boundaries, the surface will order at higher temperatures than the bulk. The transition of the bulk in the presence of an ordered surface is called extraordinary transition, and is associated with a diverging order parameter at the boundaries. In the presence of two ordered boundaries, we must then distinguish the case in which the respective spins are aligned (labeled ++) from the case where they are antiparallel (+−). The latter case describes experiments with mixtures of two liquids, attracted to opposite surfaces [32][33][34]. The former case has also been studied experimentally [35]. We are going to focus on the ++ case for the XY model. It should be added that, while measuring bulk quantities in an extraordinary transition, it does not matter how the surface order was achieved, whether it is by stronger couplings or by magnetic fields localized on the boundaries [18]. In our case, surface order can only be obtained through surface (magnetic) fields, since the surfaces are 2d and the XY model has a continuous symmetry: stronger couplings alone could not produce nonvanishing magnetization. We refer to [36] for a recent detailed analysis of the boundary phase diagram of the 3d O(n) model for varying n.
Our aim is both to obtain magnetization profiles and data for two-point correlations, and to apply the theory introduced in [37] to the XY universality class, to confirm the validity of the fractional Yamabe approach. To this end, we will describe the model that was used for numerical simulation, followed by a summary of [37]. Then we will see how the theoretical predictions provided by the fractional Yamabe equation match our numerical results.

(N + 1)-state clock model
The spins in the XY model take value in [0, 2π). Continuous variables are, however, inconvenient for numerical simulations, as generating them and computing Monte Carlo weights takes longer. Universality comes to our aid, providing a way to bypass this difficulty: we can obtain the same critical exponents and correlation functions while considering a different model in the same universality class as the XY model. A good candidate is the N -state clock model: it has the same Hamiltonian as the XY model, but the spins can take only N possible values, the vertices of a regular N -sided polygon: It has been shown [38] that the 3d clock model can be seen as an XY model with an additional crystal field, which reduces the symmetry of the system from O(2) to Z N : at the critical point, this field becomes irrelevant as long as N > 4, meaning that it does not alter the universal properties we are seeking. In general, however, models with Z N symmetry may have different properties than models with O(2) symmetry, as seen in [39].
In any numerical simulation, the free energy per site and other quantities contain a leading term which is the thermodynamic limit, plus various finite-size corrections. Following [10], we take measures to reduce these corrections: the Hamiltonian is further modified by introducing vacancies. Placing a spin in any given site lowers the energy of the system by some amount D; by tuning this parameter, one can find a point D * where the leading corrections to scaling vanish. The final Hamiltonian is then To fully define the model, it is also necessary to specify the probability measure of each spin: the (0, 0) state is chosen to have the same weight as all the other states combined, so the partition function is then (2.4)

Monte Carlo simulation
The geometry of the system is a slab of sizes (L + 1) × 6L × 6L (in the x, y, and z directions respectively), with periodic boundary conditions in the parallel (y and z) directions and fixed boundary conditions in the transverse direction. The values used for L are 32, 48, 64, 96, 128. The parallel directions have been chosen to be six times larger than the transverse one as they provide an excellent representation for the system with infinite transverse directions for the observables under scrutiny in this work. The model was simulated at the bulk critical temperature, using the values of β c = 0.56379622(10) and D = 1.02 obtained in [10]. We also used the same three kinds of Monte Carlo moves: two types of single-spin flip moves, and a cluster update.
The quickest spin flip move consists in proposing to empty the chosen site if it contains a spin, or fill it if empty, with a spin uniformly chosen out of the N values. This move, however, does not guarantee ergodicity, which is why we also use the second, more standard, type of move. In this second case, once a site is picked, regardless of the value of the spin we generate a new spin according to (2.3): en empty site half the time, a random spin, uniformly chosen out of the N possible ones, the other half. In either case, the move is then accepted or refused with the usual criterion of the Metropolis algorithm.
One of the well-known problems of simulating systems at the critical point is the critical slowing down effect: since the autocorrelation time diverges, even after many single-spin flips we still have samples correlated to the previous ones, and therefore unfit for new measurements. This problem can be ameliorated through moves that update large (correlated) chunks of the lattice at once, using the Wolff algorithm [40]. First, we select a reflection axis by picking a vector r = cos 2πm N , sin 2πm N uniformly in m ∈ {0, . . . , N − 1}. Then, we pick a random spin, mark it as part of the cluster and flip it: s i → s i − 2( s i · r) r. Next, we try and add to the cluster every neighbor of s i , according to the bond probability (2.5) The sites that have been added to the cluster are flipped, and the procedure continues in the same way. Once a spin has been added to the cluster and flipped, it cannot be selected or flipped again. To save time, the bond probabilities are computed once, at the beginning of the simulation. A difference worth highlighting, with respect to the case with periodic boundary conditions, is that the cluster may reach a spin on the boundary. To overcome this difficulty we use the recent extension of the Wolff algorithm to the inclusion of arbitrary fields recently proposed in [41]. When the cluster touches the boundary all spins on both boundary planes must be treated as a single spin S B , which gets flipped according to the same rule. All boundary spins are added to the cluster, which can then keep growing from each boundary spin.The value of the new boundary spin is then stored to be used in future measurements. For instance, whenever the magnetization of a yz plane, at a given distance x from a boundary, is computed, this must be taken as the angle between the spins on the plane and the current boundary spin (2.6) The data from the simulation are used to test the hypothesis introduced in [37], which is summarized below.

Fractional Yamabe equation
The main idea is that a bounded system at the critical point alters its original (euclidean) metric in order to be as uniform as possible. The first requirement for this new metric is that the distance between any point in the system and the boundary be infinite, so that parts of the system that were close to the boundary while using a flat metric now appear to be on the same footing as ones lying deeper in the bulk. Next, we need to preserve local properties: the system must then locally appear euclidean, so we only allow conformal changes of the flat metric: Distances between a point x and x measured with the metric g will be written as D g (x, x ). Requiring homogeneity in the system is implemented by making the Ricci scalar curvature constant (the reader is reminded that . . , d} and repeated indices are summed over): Choosing κ > 0 would mean having a space with no boundaries, such as a d-dimensional sphere S d , so κ must be negative. Without loss of generality, we can set κ = −1. Negative curvature spaces (in their prototypical realization H d ) indeed possess boundaries that are infinitely apart matching our desiderata for uniformization. This requirement turns into a constraint on γ(x), called the Yamabe equation [42] ( There are a few cases where the Yamabe equation has simple solutions: e.g., for a ball of radius ρ in any dimension, one finds γ = ρ 2 −|x| 2 2ρ [37,43]. Solutions for the case of a slab are found in [44]. This equation, however, is only expected to be valid at the mean field level, since it does not take into account the anomalous dimension of the fields. Indeed connections between the mean-field case and the Yamabe equation have been explored and successfully tested via numerical simulations in [44]. To give a hint, notice that the the exponent d−2 2 and d+2 2 correspond to the scaling dimensions of the order parameter field of a free theory in d dimensions and its conjugate field respectively: the Yamabe equation is therefore equivalent to the saddle-point equation of an O(n) model at the upper critical dimension. When anomalous dimensions are present, as in our case, the Yamabe equation must be modified into the fractional Yamabe equation (the dependence on the anomalous dimension of γ has been highlighted by introducing a subscript) The introduced modification to Yamabe equation is in principle easy and stems from the assumption that γ (∆ φ ) (x) constitutes a local gauge for measuring lengths. The subtleties involved in the definition of the suitable fractional power of the laplacian (− ) s call upon the development of an ambient construction adapted for bounded domains as done in Appendix B of [37]. The fractional Yamabe equation (3.4) can also be given some geometric meaning: it can be interpreted as the problem of finding a metric making R (s) , the so-called fractional Q-curvature (of order s = d 2 − ∆ φ ), constant. In this framework the ordinary scalar curvature R is proportional to R (1) , so for clarity one may refer to (3.3) as the integer Yamabe equation. To get some intuition on why they could appear in a critical anomalous system consider that a sphere of radius ρ has R (s) ∝ ρ −2s . This makes R (s) at least a reasonable candidate to rule the space dependence of a quantity with arbitrary real scaling dimension. For more details on the fractional Yamabe problem please consult [45]. In Figure 1 we collect solutions in the slab geometry of the integer Yamabe problem and of the fractional Yamabe for a range of anomalous dimensions relevant for this work as derived in [37].
Once the Yamabe equation is solved and the function γ(x) is found, it can be used to express correlation functions [37]: to this end, one notices that one-point functions and the scale factor transform similarly under a dilation of the system Ω → b Ω:  for the XY universality class [11].
Therefore, once γ(x) is known, all one-point functions are determined up to a multiplicative constant α: As an example, for a half-space in any dimension, with x 1 > 0 and {x 2 , . . . , x d } ∈ R d−1 , the solution to the integer Yamabe equation is simply γ(x) = x 1 , so The hyperbolic half-space is peculiar in the sense that the same γ(x) also solves the fractional Yamabe equation for any ∆ φ , meaning that (3.7) is also valid below the upper critical dimension.
Higher-order correlations are not completely determined; however, this approach predicts that they will contain a dimensionful prefactor for every field, and will also depend on the distances between the points. Based on our conjecture, the bounded critical system no longer exhibits a euclidean metric: for this reasons, distances between points are to be computed with the metric g we determined (3.1): Both (3.7) and, for the case of the 2d Ising model, (3.8), reproduce known results from Conformal Field Theory, as shown in section 4.3. In higher dimensions, they also reproduce what is known for specific domains that are interiors of balls, possibly of infinite radius, reducing to semi-infinite systems (see [46,47] for applications in the spirit of this work). For any other domain the outlined approach provides new, testable predictions that have already been numerically verified in the Ising model [37]. The hypotheses for one-point and two-point functions will now be tested against the result of the Monte Carlo simulations for the system under scrutiny belonging to the XY universality class.  To test (3.6), we employ the clock model in a d = 3 slab defined in Section 2. We measure the average of the order parameter field, i.e., the magnetization, on various planes parallel to the boundaries, obtaining a function of the transverse distance. In order to improve readability in this section and in the following we slightly change the notation since we will be dealing just with models in three (two) dimensions: the coordinates of a point p will be denoted by {x, y, z} ({x, y}). If there is a preferred direction it will be chosen to be the x coordinate. The magnetization profile will thus be m(x) with x transverse to the slab. This profile can then be fitted by (4.1) The fit parameters are a multiplicative constant α, the scaling dimension of the field ∆ φ and the extrapolation length λ. This last parameter is introduced to match numerical profiles for spin systems with continuous profiles from field theory: on the lattice, the magnetization does not diverge. A divergence would appear if the lattice profile were continued λ sites beyond either boundary [17,29]. We notice that the extrapolation length remains roughly constant as the system size is increased, meaning that larger sizes capture a greater fraction of the continuum profile.
In Figure 2, we plot the magnetization profiles for different system sizes, appropriately rescaled by multiplying each profile by L ∆ φ . The x coordinate has also been rescaled, to take into account the extrapolation length of each size. First of all, we notice a clear collapse of the different profiles onto the same curve, proving that we are at the critical point. For each system size, the fit using (4.1) gives us a value of the scaling dimension ∆ φ , as shown in Table 1. As our final estimate, we take the average of the different values, obtaining ∆ φ = 0.5206 (15). (4.2) This is compatible, to one standard deviation, with the most precise estimate ∆ CB φ = 0.519088 (22), obtained through conformal bootstrap [11]. This result should be intended as a check of the theory, since obtaining critical exponents is far from the only use of the critical geometry approach. Details on the data analysis are found in Appendix A.  Table 1: Scaling dimensions obtained from the fit for each system size. As final value, we take their average: ∆ φ = 0.5206 (15).

Correlation functions
Testing (3.8) is slightly more subtle because we cannot determine the function F. Once the correlation data are divided by the product of one-point functions, i.e., they will depend only on the distance between the points computed with our metric g. This means that the points will collapse onto a single line if plotted as a function of this distance, while they will appear scattered if plotted as a function of any other (nonequivalent) distance. The data we collected is averaged over the parallel directions: φ(x, y, z)φ(x , y + ∆y, z + ∆z) . (4.4) To limit the volume of data, we measured correlations only in a 15 × 15 × 15 grid: x, x = (i + 1)L/16, ∆y, ∆z = iL/16, where i = 0, . . . , 14. Once we take x ≥ x , we consider the points C(x, x , ∆y, ∆z) and C(x, x , ∆z, ∆y) to be the same data point, and exclude coinciding points, we get 7672 independent correlators. The fractional Yamabe distance between the points was computed via Surface Evolver [48] using as metric g = δ/γ 2 (∆ CB φ ) . The correlation ratio as a function of the distance D g (x, y) are shown in the left pane of Figure 3: we see a clear collapse onto a single line. The value chosen for the scaling dimension ∆ φ affecting γ (∆ φ ) (x) and in turn the distances D is ∆ CB φ ; however taking our best estimate (4.2) makes actually no visible difference. To see what would happen for a different, wrong from the outset, metric, we also plotted the same points as a function of 2d Ising (CFT) Figure 3: : on the right, plotted as a function of the euclidean distance: since that is not the "correct" metric, the points are scattered. On the right, the same points plotted as a function of the distance D g (x, x ), computed with the metric which solves the fractional Yamabe equation for ∆ φ = ∆ CB φ , the most accurate estimate. The points collapse onto one line. In the inset a zoom of the region where the ratio displays more pronounced variations is shown.
the euclidean distance between them, in the right pane of Figure 3: a collapse is absent and the point are scattered.
To quantify the goodness of the collapse, we fit the points in 3 with the function f (ξ) = 1 + a 1 e −b 1 ξ + a 2 e −b 2 ξ , and then compute the mean square displacement between the points and the curve (n d.o.f are the independent degrees of freedom): We obtain σ = 0.0077 for L = 32, σ = 0.0036 for L = 64 and σ = 0.0029 for L = 128: the fact that displacements become smaller for larger sizes lends quantitative support to the collapse. Notice that the chosen fitting function f correctly captures the long (mutual) distance limit of the correlation ratio, that is r → 1 as D → ∞. Concerning the short distance properties, the operator product expansion structure should be retrieved; this, however, is not the focus of the present investigation, since short distance behavior is not crucially affected by the curved space properties we are interested in, as well as being challenging to retrieve from lattice approaches.

Comparison with the Ising model
As we mentioned, the function F(D g ) of the distance between the points is not known in d = 2 and our approach cannot determine it without using the input from the numerical simulations. Figure 3 shows the existence of F through the collapse of the data. One sees that we can obtain a precise numerical estimate of it.
To put in context the result obtained for the F of the 3d XY model, we can now set it side by side with the same function for the Ising model. In Figure 4, we compare the correlation ratio r we obtained for the XY class with the corresponding profiles for the 3d Ising model (obtained in [37]). We see qualitatively similar curves in the two cases, with greater correlation ratio for the XY model with respect to the Ising one. The quality of the collapse in 3d both for the XY and the Ising models is similarly convincing. Notice that the metric for the Ising model is different from the one for the XY model, since they have different scaling dimensions, producing distinct fractional Yamabe equations and yielding in turn different collapse distances. As expected, both correlation ratios tend to 1 for large D.
To analyze the role of the dimension d, we consider the 2d Ising case, which is an important example of the validity of our approach when applied to d = 2. Obtaining the function F for the 2d Ising model is straightforward via boundary Conformal Field Theory (CFT). In particular, the two-point function for the Ising model on a half-plane x > 0, with diverging magnetization at the boundary, is known to be [49]: where α is a constant. This clearly contains the product of one point functions φ({x, y}) = α x −1/8 , and τ is indeed a function of the hyperbolic distance between the two points, since so τ = coth 2 (D/2). This means that (4.6) is in the form (3.8). 1 One can then conformally map the half-plane into a strip, and obtain a result which is still in the form (3.8): additionally, while one-point functions and the distances between the points change with the mapping, the function F remains the same. Figure 4 is the main result of this work, as it shows that the critical geometry approach is valid for different models. It appears that the 2d case has weaker correlations (i.e., smaller r) than 3d, and -comparing models living in the same space dimension -the XY model has stronger (normalized) correlations than the Ising one.
Simulations for the 1-point magnetization of the 4d Ising model have been recently performed [44]. Preliminary results for the 2-point spin correlations show that, at the considered sizes, collapse takes place similarly to Figure 4 and the correlation ratio r appears to be different from 1, although tending to 1 for large distances, as it should. A careful analysis of the logarithmic corrections at the upper critical dimension and more numerical data are needed to draw further conclusions.
In future works, it would be interesting to compare the results in Figure 4 with the corresponding results for 3d O(n) models with n ≥ 3 and for models at the upper critical dimension. Finally, a crucial point is to confirm that, as is the case in 2d and as the Yamabe approach seems to suggest, the curve of Figure 4 in 3d does not depend on the bounded domain one chooses, be it a slab, a half-space, or any other.

Conclusions
We have shown that the introduction of a curved metric in a bounded domain can reliably describe correlation functions of statistical fields for the 3d XY universality class. Numerical simulations of a 3d improved clock model on a slab geometry have provided magnetization profiles consistent with the critical geometry approach [37,44], as well as a clear collapse of the 2-point functions. Our main results are summarized in Figure 4, where the structure of correlators has been determined and the dependence on the curved metric distance is found for the 3d XY model and compared with the Ising model in 3d and 2d. A value of ∆ φ has been determined, and is consistent with the most recent estimates.
An open question, suggested by the results of the present paper, is whether there is a way to determine the function of the curved distance which appears in (3.8), to determine 2-point correlations completely. One could also investigate whether correlations in different O(n) models (such as the Ising, XY and Heisenberg models) could by described by a boundary-independent function dependent on n. Moreover, since in this paper we argued that 1-and 2-point spin correlation functions in a bounded domain are described by the fractional Yamabe equation, it would be interesting to investigate possible generalizations of the Yamabe equation in connection to models with vectorial order parameters and more general symmetry groups.
The well-known 3d models (such as the Ising and XY models, as well as percolation) have a rather small anomalous dimension, meaning that ∆ φ is close to the mean field value 1/2, which is described by the integer Yamabe equation [37,44]. This suggests another line of research: it may be possible to obtain an approximate solution of the fractional Yamabe equation as a perturbation of the corresponding integer equation. This would provide a 1 Remember that the hyperbolic metric is gij = δ ij x 2 and it solves the Yamabe problem for the half-plane, and actually solves the fractional Yamabe problem since it makes the half-plane fully homogeneous so all sensible curvatures are constant. This is discussed in detail in [37]. perturbative solution to a fractional differential equation without the need to compute fractional derivatives. Finally, the study of critical correlation functions and exponents through the critical geometry approach based on the fractional Yamabe equation is not limited to unitary models. A famous example of a statistical model lacking unitarity is percolation, which is therefore a subject worth considering for future studies.
Acknowledgements: The authors thank M. Hasenbusch for useful discussion and correspondence. GG is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC 2181/1 -390900948 (the Heidelberg STRUCTURES Excellence Cluster). GG also acknowledges QSTAR for hospitality during completion of this work. Stimulating discussions with several participants of Bootstat 2021, held in May in Institut Pascal, Université Paris-Saclay and online, are also acknowledged.

A Simulation and data analysis details
To ensure the validity of the simulations described in the main text, the following tests have been performed: • The critical temperature obtained in [10] has been checked by computing the Binder cumulant [50] in a system without fixed boundaries.
• The system has also been simulated at β c + σ βc and at β c − σ βc (with β c = 0.56379622 and σ βc = 10 −7 ): the resulting magnetization profiles are compatible with the ones obtained at β c , i.e. they are within the respective statistical errors. This means that, for the studied sizes, the uncertainty on the critical temperature does not affect our results. However, further increasing the size L, the uncertainty on the location of the critical point could affect the results, and supplementary studies would be necessary.
• The profiles have been obtained for various ratios between the length in the parallel directions and in the transverse direction. The final value we chose is 6, since increasing it further does not change the profiles beyond statistical uncertainties.
• Additional small-size simulations were performed using only single-spin flips, without cluster update, as a cross-check: once again, we verified that the results were compatible with the ones obtained from the employed optimized algorithm, which includes cluster updates.
For any system size L, two-to-four-day simulations were performed on up to 160 Intel(R) Xeon(R) CPUs E5-2680 v2 cores running at 2.80GHz. After the magnetization data have been obtained, an additional step is needed before performing the fit. The points closest to the boundary are most affected by finite-size effects. Therefore, despite having smaller errors than the central points, a few of them have to be discarded. In order to determine how many to discard in an unbiased way, as well as to avoid a sharp distinction between discarded and included points, we introduce a window function w(x). The weight of each point in the fit is given by the square of the ratio between this function and the error of that point. This function starts off from 0 at the boundary, ramps linearly to 1 around a movable point t, and maintains the value 1 up to the center of the slab. To determine the location of the point t, we start from t = −1 (the boundary point) and gradually move towards t = 0. For each value of t we compute the χ 2 of our data, and the corresponding p-value. We stop once the p-value reaches the reference value of p = 0.95.
Once the values ∆ φ (L) have been obtained, as seen in Table 1, the final value is calculated as their average; the error on ∆ φ is the standard deviation σ of the set {∆ φ (L)}. Since their errors are too small for them to be compatible with one another, σ was not divided by the square root of the number of data points.