Finite density 2d O(3) sigma model: dualization and numerical simulations

The action of the 2d O(3) non-linear sigma model on the lattice in a bath of particles, when expressed in terms of standard O(3) degrees of freedom, is complex. A reformulation of the model in terms of new variables that makes the action real is presented. This reshaping enables us to utilize Monte Carlo simulations based on usual importance sampling. Several observables, including the correlation function and the mass gap, are measured.


Introduction
The physics of dense matter is crucial for understanding the drastic changes underwent by the Universe during the first minutes after the Big Bang. The extreme conditions that prevailed in that period can be partially reproduced in experiments with heavy-ion collisions where the effects of dense matter are probed.
The presence of a density of matter affects the dynamical laws that govern the physical systems. The study of these effects by numerical simulations is usually hindered by the so-called sign problem. This problem consists in the fact that upon adding a coupling with an external chemical potential, the functional that weighs field configurations may not furnish a positive number, thus ruining any attempt to use that functional as a probability for importance sampling in a Monte Carlo procedure.
Several methods have been proposed to overcome those difficulties in Monte Carlo simulations of Quantum Chromodynamics (QCD), see for instance Ref. [1]. The quest for better strategies has raised the interest in 2d toy models that are afflicted by similar problems [2][3][4][5][6][7][8][9]. A technique that has proven to be particularly appealing to evade the sign problem is dualization: the model is recast in terms of new (dual) variables in such a way that the new action turns out to be real.
No general recipe exists for transforming ordinary into dual variables. Every model has to be studied on its own and the strategy to get a dual representation may differ significantly for different models. Moreover, not even a unique dual representation exists for a given model. In fact, among the few different possible dual representations available for a given model, some may be more advantageous than others for simulation purposes or some might even present a so heavy slowing down that it makes the dual version, albeit real, useless.
Another type of drawback that often appears in dual formulations is that certain observables cannot be disentangled from the probability weight in such a way that their expectation values have to be extracted as the ratio of expectation values with different Hamiltonians, which in general gives rise to extremely inefficient numerical calculations. Such difficulties typically arise in the evaluation of non-local observables like correlation functions.
In the present paper we apply the dualization idea to the 2d O(3) non-linear sigma model in presence of a chemical potential. The standard action of the model on a d-dimensional lattice Λ ∈ Z d without a chemical potential is given by and the partition function reads In 2d this model possesses a particle triplet with a spontaneously generated mass gap m [10]. The value of this mass has been verified numerically [11]. The model is asymptotically free and presents a rich topology [12]. All those properties make this model a close relative of QCD. The physical interest of the 2d O(3) non-linear sigma model goes well beyond the above list of properties. In the first place the non-linear sigma model reproduces reliably several qualitative traits of ferromagnetic materials. Another reason for casting relevance to the model is that it is involved in the development of the resurgence program [13,14].
Our purpose is to construct a real action for the 2d O(3) non-linear sigma model at finite density expressed in terms of dual variables and in such a way that it allows us to determine vacuum expectation values with usual Monte Carlo methods.
To describe the theory at finite density one introduces a chemical potential µ in the original action as an external source for the total third component of the angular momentum [10]. Following [15] and using the spherical parameterization the action (1) can be written as where the more general situation with anisotropic chemical potentials µ ν , ν = 1, 2 is considered. The conventional physical situation is recovered if we put µ 1 = µ, µ 2 = 0. Three different routes to introduce dual variables are possible, and it is not obvious a priori which one is preferable and under what circumstances may it be so. The first route relies on the fact that the chemical potential is introduced in the This action is of an XY type with a fluctuating coupling and in the presence of µ ν . As shown in [16] in the case of the Villain formulation of the XY model and for β ν (x) = β for all x, the conventional dual transformations performed with a nonvanishing chemical potential lead to a positive definite Boltzmann weight. Moreover, if the coupling β ν (x) is positive for all x, then it is straightforward to prove that exactly the same transformations lead to a positive dual weight for all O(N) models. This conclusion has been explicitly demonstrated in [17] for the 2d O(3) model, and the proof can be readily extended to all O(N) models. The second route consists in Taylor expanding the Boltzmann weight and integrating over the original degrees of freedom. The dual variables appear as flux variables subject to certain constraints, and the dual weight can be proven to be positive for O(N) models [18,19]. The resulting dual theory can be simulated by a worm algorithm, and a number of thermodynamic quantities have been computed along this route [18][19][20].
In the third and last route one constructs the dual theory by expanding the Boltzmann weight in hyperspherical harmonics on O(N) and integrating out the original variables. This program has been accomplished for the 2d O(3) model in [17]. However, the full positivity of the resulting dual weight remains to be proven. It is important to stress that, at least in the context of these two-dimensional models, the dual formulation appears as the only reliable tool to investigate the properties of the model. Indeed, the results of Ref. [20] show that alternative approaches to the sign problem -the reweighting and the complex Langevin-have certain drawbacks and lead to incorrect results in some regions of the β-µ plane.
One of the main purposes of the present work is to develop a dual formulation applicable not only for computing thermodynamic quantities, but also for extracting long-distance quantities like the correlation functions. With such a formulation in hand, one could be able to reliably address the question of the hypothetical Berezinskii-Kosterlitz-Thouless (BKT) transition in O(N) models at non-vanishing chemical potential. Our approach is essentially the first route described above. We reformulate the O(N) model in terms of the link formulation for the Abelian (sub)group. In this formulation our results can be straightforwardly extended to all O(N) models in any dimension. Here, for the sake of simplicity, we consider only the two-dimensional O(3) sigma model.
The paper is organized as follows. In the next Section we construct the dual formulation. Firstly we introduce the link representation, then we obtain two alternative dual Boltzmann weights, both of which are positive. To study the real effectiveness of the above Boltzmann weights in Monte Carlo simulations, we test one of them by calculating numerically several observables. We are particularly careful to individuate any slowing down during the simulations. In Section 3 we outline the procedure employed during the simulations and list the observables that have been studied. Specifically we calculate the particle density to find the expected threshold at µ = m, verify that the mass gap extracted from a correlation function is insensitive to µ, and evaluate the energy as a function of β and µ. The results are presented in Section 4 and some conclusive remarks in Section 5.

Dual representation
Action (4) becomes complex if any of the µ ν 's is non-zero. However, an action with a positive Boltzmann weight can be constructed by using, in the spirit of [18,19], a dual representation of the model with the action (4). Nevertheless, our strategy is somewhat different from [18,19] and relies on the use of the so-called link formulation [21,22]. In fact, this approach is similar to the one used in [16,17]. One of its advantages is that it can be readily extended to any O(N) model, in any number of space dimensions.

Link formulation for Abelian subgroup
To build a dual representation, we start from the following partition function with the action (4), The integration can be done in a number of ways. Here we use the fact that the dependence of (4) on the angles φ(x) is only through their differences (U(1) variables). This allows to make a change of variables and rewrite the partition function in terms of the link angles φ(l) The procedure generates local and global constraints known as Bianchi identities on the link variables [21,22]. The local identity constrains the allowed configurations of φ(l) on every plaquette p of the lattice and can be embedded into the partition function in the form of a periodic delta-function as Global identities constrain two holonomies winding through the lattice in periodic directions. They have the form Then, it is easy to prove the following equality in any number of dimensions d: One should keep in mind that when d > 2 not all local Bianchi identities are independent.
Applying this approach to the 2d O(3) model one gets, after integration over link variables, where L ν are the values of the lattice size in the two directions, I r(l) is the modified Bessel function of first kind, and The plaquettes p 1 and p 2 have the link l = (x; ν) in common. The two-point correlation function between the origin (denoted by a nought 0) and a point R in the parameterization (3) reads with where the expectation values · · · are evaluated with (6). In the link formulation Γ 2 (R) is given by where η(l) is defined just below and C R is any lattice path connecting the points 0 and R. Introducing a set of sources ζ = {h 1 (x), h 2 (x), η(l)} and integrating out link variables one gets where h 1 (x) = 1 for x = 0, R and h 1 (x) = 0 otherwise and We have introduced here notations ζ = (0, h 2 (x), η(l)) and ζ ′ = (0, h 2 (x), −η(l)), where h 2 (x) = 1 for x = 0, R and h 2 (x) = 0 otherwise; η(l) = 1 if l = (x; ν) ∈ C R , η(l) = −1 if l = (x − e ν ; ν) ∈ C R and η(l) = 0, otherwise. The partition function utilized in (15) and (16) is given by with B η (l) = I r(l)+η(l) (β sin α(x) sin α(x + e ν )) .
As it stands, the expression (17) is much more general and allows to compute correlations of any kind as

Dual Boltzmann weight 1
The dual Boltzmann weight can be read off either from (10) or from (17). It has the form Clearly, it is strictly positive in the integration domain over α(x) and hence it can be used for the numerical simulations. In this case we have

Dual Boltzmann weight 2
A different dual Boltzmann weight can been constructed by performing a complete integration over the remaining original degrees of freedom. We begin with the formula π 0 dα F (cos α, sin α) = s=±1 π/2 0 dα F (s cos α, sin α) .
This introduces Ising-like partition/correlation function with fluctuating coupling Equation (17) becomes where (taking into account that h 1 (x) = 1 for x = 0, R and zero otherwise) The dual transformation can be performed in a standard way by introducing the link variables z(l) = s(x)s(x + e ν ) (this time the global constraints on holonomies can be omitted) to obtain where z(p) = l∈p z(l). Summation over z(l) leads to the following representation on the dual lattice (now C R is a path connecting the points 0 and R and consisting of links dual to the original links): By inserting (28) into (25), we get the expression Here B η (l) = I r(l)+η(l) (β sin α(p) sin α(p ′ )) , where Finally, one can integrate out the α(p) angles. This can be done by Taylor expanding the factor e ±βJν (x) and by using either the series representation for the Bessel function or the multiplication theorem for the Bessel function to decouple the sin α(p) factors from their argument. The first approach is somewhat simpler and leads to the following result: where B(x, y) is the beta-function and a(p) and b(p) are given by The partition function can be obtained from the last expression if we put h 1 (p) = h 2 (p) = η(l) = 0 and extend the second product in the last line to all links of the lattice The Boltzmann weights of all three representations (10), (29) and (34) are positive. Moreover, the dual weights are free of constraints. It means, in particular, that one can compute not only local quantities (those which can be presented as derivatives of the partition function), but also long-distance quantities like the two-point correlation function, because it can be presented as an ordinary expectation value.

Simulation details
We have two real Boltzmann weights, (10) and (34). However, being real is not the only condition they must satisfy to be handy. In conjunction with an adequate simulation algorithm, they should also avoid slowing down. For this reason we have tested one of the two weights. The tests have been done with (10) because, as it includes a few dynamical variables, we expect that it is more economic and thereby requires less computational effort. We have simulated the action (10) on a square lattice Λ ∈ Z 2 of lateral extents L 1 , L 2 with periodic boundary conditions. The variables are α(x) ∈ [0, π], r(p) ∈ Z and q 1 , q 2 ∈ Z where x indicates sites and p plaquettes. A single Monte Carlo step consists in updating separately variables α(x), r(p) and q 1 , q 2 . All updatings were performed by the Metropolis algorithm [23]. The angle variables α(x) were refreshed by proposing a brand new value of cos α(x) and applying the usual Metropolis test on it. Plaquette variables r(p) and global variables q 1 , q 2 were updated after randomly choosing a new value that differs from the old one by at most ±∆ units (we took ∆ = 5 in r(p) and q 1 , q 2 ).
Next we introduce the observables measured in this work. In order to enhance decorrelation, single measurements were evaluated on configurations separated by 10 Metropolis updatings. (i) Since the number N of particles is we deduce that the particle density operator is or, equivalently, Both observables were measured for testing purposes as both (37) and (38) should provide the same result. The test was successfully passed.
(iii) The correlation was measured by a wall-wall correlation function. It was constructed according to the definition (22). Since the path followed to join the walls is irrelevant (on average), we chose it as the simplest one guaranteeing computing efficiency: along the x 1 -axis. For clarity, in Fig. 1 the 16 paths composing a wall-wall correlation at distance 2 from site (2, 0) to site (4, 0) on a 4 × 4 lattice are shown. The total wall-wall correlation function at distance 2 consists in summing the contributions from all the above paths and averaging the result over all possible ways the walls can be placed at distance 2: erecting the walls at sites (1, 0) and (3, 0), at (2, 0) and (4, 0) -this is the one shown in Fig. 1-, at (3, 0) and (1,0), and at (4, 0) and (2, 0), having used periodicity on the last two. In general, that average contains L terms on a L × L lattice.

Results
The Monte Carlo simulations of (10) were performed on square lattices L 1 = L 2 ≡ L. We used L = 20, 30 and 40. We employed such small sizes for testing purposes only. Moreover, our algorithm slowed down on larger lattices (L = 40 or larger). This problem manifests itself in a poor variability of the q ν variables. Therefore we carefully verified that such a variability and the corresponding fluctuations were present and regarded this verification as an indispensable condition to accept the results of every simulation run. We typically collected some 10 7 measurements for each run. Our strategy to construct dual Boltzmann weights allows us to employ them to calculate every observable, including correlation functions. Such functions permit to extract the mass gap. This gap can also be determined (i) as the value of the chemical potential µ at which non-zero particle numbers start off and (ii) as the value of µ at which the energy density detaches from its value at zero particle density. Therefore our procedure enables us to cross-check the results for the mass gap.
Assuming the following theoretical form for the wall-wall correlation function, Here m eff (R) is given in units of the inverse lattice spacing. We expect that m eff (R) exhibits a plateau for large enough R. The height of that plateau is the value of the mass gap resulting from our simulations. We call m MC this height. In Fig. 2 we summarize our findings for two values of β and µ. This plot follows after measuring the correlation function on 2 · 10 7 configurations obtained in a 20 × 20 lattice. Similar results have been derived for other β and µ (not shown). In all cases the expected independence on the chemical potential is apparent.
In Fig. 3 the energy density (40) for zero particle density is displayed as a function of beta. Again L = 20. Data exhibit a satisfactory agreement with the weak coupling and with the strong coupling expansion [26] E = y + 2y 3 + 12 5 These agreements constitute further positive tests of our procedure. Figure 4 shows the energy density E as a function of β and µ on a 20 × 20 lattice. As expected, the value of µ at which E detaches from its value at µ = 0 is µ = m MC . Even though we have used rather small lattice sizes, the coincidence between thresholds and effective mass m MC is manifest because the energy density operator has negligible size effects. Figure 5 shows the dependence of the particle density on µ for several β. Each point is the average of 10 6 measurements. In principle also the thresholds in this plot should coincide with the plateaux in Fig. 2. This coincidence is not so manifest at L = 20 because, contrary to the energy density, n has strong size effects [20]. We have verified this assertion by repeating the study for L = 30 and L = 40 at β = 1.2. Figure 6 blatantly exhibits that size dependence. While the threshold at L = 20 is at about µ ≈ 0.2, at L = 40 it has shifted to about µ ≈ 0.3. The point along the line for β = 1.2 in Figure 4 at which data detach from the value of the energy at zero density agrees reasonably well with the threshold in Figure 6 for L = 40. Figure 6 also displays the consequences of the above-mentioned slowing down. Indeed, although the three lines have been obtained with the same statistics, namely 2 · 10 7 measurements, all points on the line for L = 40 and some of the points along the line for L = 30 present much larger error bars. They are due to long correlations among data. Such errors can be reduced only at the cost of increasing exorbitantly the statistics.

Discussion
We have derived two different positive Boltzmann weights for simulating the 2d O(3) non-linear sigma model at non-zero chemical potential. These weights permit the evaluation of any observable average. In particular we have extracted the correlation functions and from them the mass gap. These calculations have been confronted with the determination of the mass gap from the behaviour of the particle and energy densities as functions of µ. The comparison was successful. Also the energy density matches the weak and strong coupling expansions (where they should) for µ = 0. We were forced to do all these tests on small lattice sizes (L=20) because the algorithm slows down in the refreshment of variables q ν as the size grows beyond approximately L = 40. Tempered Monte Carlo and other improvements did not manage to speed up the simulation at those lattice sizes.       described here to the 2d O(3) non-linear sigma model with a topological θ-term. In this case the sign problem is even more severe. In the past we obtained a positive weight for the model at non-zero θ by mapping it to a certain dual version of the SU(2) principal chiral model [27,28]. However, the worm algorithm applied for this dual model performed rather poorly during the Monte Carlo tests [28]. The approach described in this paper can be generalized to the model which includes the θ-term. Whether or not one can construct the positive Boltzmann weight in this case or, at least, significantly reduce the sign problem remains to be verified.