Meron-cluster simulation of the quantum antiferromagnetic Heisenberg model in a magnetic field in one- and two-dimensions

Motivated by the numerical simulation of systems which display quantum phase transitions, we present a novel application of the meron-cluster algorithm to simulate the quantum antiferromagnetic Heisenberg model coupled to an external uniform magnetic field both in one and in two dimensions. In the infinite volume limit and at zero temperature we found numerical evidence that supports a quantum phase transition very close to the critical values $B_{c}=2$ and $B_{c}=4$ for the system in one and two dimensions, respectively. For the one dimensional system, we have compared the numerical data obtained with analytical predictions for the magnetization density as a function of the external field obtained by scaling-behaviour analysis and Bethe Ansatz techniques. Since there is no analytical solution for the two dimensional case, we have compared our results with the magnetization density obtained by scaling relations for small lattice sizes and with the approximated thermodynamical limit at zero temperature guessed by scaling relations. Moreover, we have compared the numerical data with other numerical simulations performed by using different algorithms in one and two dimensions, like the directed loop method. The numerical data obtained are in perfect agreement with all these previous results, which confirms that the meron-algorithm is reliable for quantum Monte Carlo simulations and applicable both in one and two dimensions. Finally, we have computed the integrated autocorrelation time to measure the efficiency of the meron algorithm in one dimension.


I. INTRODUCTION
When quantum Monte Carlo simulations are perform in frustrated magnetic systems, they suffer from the severe sign problem [1], [2]. Well known representative of this class of systems are antiferromagnetic systems in a non-bipartite lattice, fermions in higher dimensions than one, ferromagnetic or antiferromagnetic systems in a presence of an external magnetic field, among others. In all of these models the Boltzmann factors associated with some possible configurations of the system in the path integral representation can be negative, and therefore cannot be interpreted as probabilities [3]. This feature represents a difficult task to solve when defining the stochastic process. In general the sign problem is of the non-polynomial class [4] and should be solved case by case. Even if one could define positive transition probabilities and include the sign into the observable such that accessibility and detailed balance are fulfilled, fluctuations of the sign in the sampled configurations would generate strong cancellations of the physical observables, which would make it impossible, to measure them accurately using these numerical methods [5].
Different algorithms have been proposed to address this problem, and among them we mention the worm-and the meron algorithms [6], [7]. Both algorithms have the important feature of including clusters in the update procedure, which is essential to beat or at least reduce critical slowing down [8], [9].
In this paper, we perform a numerical simulation of the antiferromagnetic Heisenberg model (AFH) in both one and two dimensions in the presence of an external magnetic field. From the numerical point of view, the dynamical competition between the spin-spin coupling and the coupling to the magnetic field leads to the existence of the sign problem, which has slowed down the development of efficient algorithms. We have implemented the meron cluster algorithm to avoid the sign problem and combined it with the concept of an improved estimator to measure the magnetization, the order parameter in this model. The numerical results agree remarkable well with the corresponding values obtained in article of Chandrasekharan et al [7] in two dimensions and also with the one obtained by using scaling behavior laws and extrapolation of them in the one and two dimensional models [10].
There is another interesting physical feature involved in some physical systems whose numerical simulation is hampered by the notorious sign problem: there is a dynamical competition between the interactions which leads to a non-unique ground state at low temperature as in the quantum antiferromagnetic Heisenberg model in the presence of an external magnetic field. This property leads to a quantum phase transition (QPT) at zero temperature in such systems [11]. In order to find and describe QPTs different observables have been proposed like the entanglement of formation, the concurrence among others [12]. We study the behaviour of the magnetization in the one dimensional AFH spin chain and found numerical evidence of a QPT at the critical magnetic field value B c = 2, which is consistent with the second order critical point conjectured in [13].

II. THE MERON-CLUSTER ALGORITHM
We consider a system of 1/2 quantum spins defined on both a lattice of one dimension as well as on a square lattice, with site labels i, j and equipped with periodic boundary conditions. The spin operators fulfil the SU(2) algebra of commutations: whereŜ a i is the a component of the spin operator defined on site with label i. The Hamilton operator describes the interaction between nearest neighbours on the lattice, where J > 0 is the coupling constant of the antiferromagnetic interaction and B represents an external constant magnetic field.
We will compute the partition function Z of the system written as a path integral. The path integral representation is obtained by decomposing the Hamiltonian into (2d+1)Ĥ i pieces, where d = 1,2 is the system dimension and by using the Trotter-Suzuki formula [14], [15]: The important point to note here is that all the individual terms present inĤ i commute with each other, but two differentĤ i andĤ j do not commute. This allows introducing (2d + 1)N identities between each exponential factor in order to derive the path integral. Here ǫ = β/N represents the elemental lattice spacing in the extra dimension of length N.
Once these formal manipulations have been performed, one obtains a (d+1)-dimensional lattice of classical Isingtype of spins. Using the eigenstates |↑ and |↓ of theσ 3 operator, the transfer matrices turn out to be: The non-zero matrix elements ofT 1 (up to an irrelevant overall pre-factor e ǫJ/4 ) are: while the corresponding non-trivial matrix elements ofT B are: where in these equations, the bra vectors ↑| and ↓| are the dual eigenstates ofσ 3 defined on the time slice t, and the ket vectors |↑ and |↓ are the eigenstate defined on the time slice t + 1 of the d + 1 dimensional lattice.
At this point it is worth mentioning that we have chosen the quantization direction perpendicular to the one of the magnetic field, such that one can still update whole spin clusters with probability 1/2. This choice is essential because a flip probability depending on the value of the magnetic field would rarely lead to possible flips of magnetized clusters for strong or weak fields, which would destroy the efficiency of any algorithm.
The non-zero elements of the transfer matrices define the allowed spin-configurations of the auxiliary lattice of d + 1 dimensions. From Eq. (2.7) and (2.8) it follows that a sign problem arises for J > 0 and B < 0 respectively. The minus sign appearing in F can be easily avoided by making the choice B > 0 which does not affect the physical contents described by the Hamiltonian. Nevertheless, for an anti-ferromagnetic interaction, J is necessarily positive, which leads to a minus sign in the matrix element C and prevents interpreting this term as a Boltzmann factor (or probability) in the simulation process. From the physical point of view this sign problem is generated by the competence between the anti-ferromagnetic interaction and the interaction with the external magnetic field B, which leads to a frustrated ground state.
In order to get rid of the sign problem it is useful first to rewrite the matrix elements as a product of the sign of the elemental configurations (plaquettes) times the Boltzmann's factor associated: where We now briefly explain how the meron-cluster algorithm proposed originally in [5] gets rid of the sign problem. The main idea is to compute observables in the original system by computing averages using the modified system with action S[s] which does not suffer from the sign problem: It follows that the original expectation value O is calculated as the ratio of two exponentially small modified expectation values OSign and Sign . Due to the strong cancellations coming from the sign fluctuations both terms appearing on this ratio are very small compared to their statistical errors and therefore it is in practice impossible to measure them accurately. The solution to this apparent problem involves a two step procedure: in a first step the algorithm matches any contribution with one sign with another paired contribution with the opposite sign, which will not give a net contribution to the observable. Only a few unmatched relevant contributions remain. In a second step, the algorithm discard these paired contributions by using a Metropolis decision, and includes only those contributions to the observables that are relevant (see discussion of the improved estimators in the next section). This step is called reweighting, and aims to suppress multi meron configurations, which should not be sampled in the stochastic process as they do not contribute to the physical observables. In the modified system without the minus signs it is possible to define a cluster algorithm by introducing break-ups associated to the active plaquettes and build a Markov chain which fulfils detailed balance and the micro-causality condition in a similar way as it is performed with the loop cluster algorithm [16], [17]. In fact, in the modified system the weight are always positive and therefore can be interpreted as Boltzmann's factors, which leads to the break-ups displayed in figure 1. The probabilities are defined by: P A || = P C = = P F : = 1, P B || = exp (−ǫJ/2)/ cosh (ǫJ/2) = 1 − P B = and P E : = 1 − P E | = tanh (ǫ|B|/2). By joining together the resulting break-ups on each plaquette different strings are produced, which we will call clusters. Due to the " : " break-up defined for the E y F configurations, the clusters will not correspond necessarily to closed loops. Nevertheless, each spin s of the d + 1 dimensional lattice will belong to only one cluster and each cluster will be updated independently with probability 1/2.

III. IMPROVED ESTIMATORS
The sign of a configuration Sign[s] was defined as the product of the signs of the elemental configurations or plaquettes which it is made of. This distinction can be translated to the clusters itself by defining a meron cluster as a cluster whose flip would change Sign[s]. This definition leads to the distinction of the configuration space between the zero-meron section, which consists of only non-meron clusters, and the meron sector, which consists of at least one meron cluster. Using the meron concept allows gaining an exponential factor in statistics and moreover it allows to compute improved estimators by measuring global properties of the cluster as it was shown in [18]. In particular the magnetization is an interesting physical observable. Due to symmetry properties of the physical system described by the Hamiltonian in Eq. (2.2), the only non-trivial component of the magnetization is the one parallel to the external magnetic field. The average of this component can be very efficiently measured by computing the temporal winding number of the closed-loops which result from joining the open string clusters [18]: where W j is the temporal winding number associated to the closed loop j. If the loop is not made by open string clusters, then necessarily it does not contribute to the observable W = 0. N corresponds to the number of meron strings appearing in the considered configuration [s]. The kronecker's delta constraints the configuration space to the zero-meron sector, only from which the magnetization gets non-vanishing contributions.
Since we are interested in the magnetization, it is in principle unnecessary to generate any configuration that contains meron-clusters. This is a crucial point of the meron-cluster algorithm as it allows to gain an exponential factor in statistics by restricting the simulation to only the zero-meron sector, which is exponentially small compared to the whole configuration space. In spite of this enormous statistical advantage it has been turned out advantageous to allow some meron configurations which do not contribute to the magnetization but speed up the stochastic dynamics by reducing the autocorrelation time. This step was performed by introducing an additional Boltzmann-type factor q < 1 for each meron such that if a new configuration with N ′ merons is proposed starting from a configuration containing N merons the new configuration is accepted with conditional probability p = min[1, q N ′ −N ]. As it was shown in [18], adding this extra Boltzmann factor in order to suppress multi-meron clusters for the re-weighting only modifies the algorithm but the physics remains the same.

IV. NUMERICAL RESULTS
Now we present numerical results obtained by using the meron cluster algorithm briefly explained in the last two sections. We studied the anti-ferromagnetic Heisenberg model defined both in one as well as two spatial dimensions. The spin-spin interaction constant was set equal to J = 1, while the coupling B for the external magnetic field was varied in the interesting physical region 0 < B < 6. As a first application we studied the one dimensional Heisenberg model in the low-temperature region. In Fig. 2 the magnetization density is shown as a function of the external magnetic field for four different lattice sizes L = 60, 80, 160 and 240 and for three values of the inverse of the temperature β, 15, 20 and 30. As a comparison, we include the results obtained in [10], which are based on exact diagonalization on small lattice sizes, guessed scaling relations and on the Bethe ansatz to extract the thermodynamic limit of the magnetization at zero temperature. Both results agree very accurately in the whole range of the external magnetic field 0 < B < 2 which shows that indeed the finite size effects are not really important in one dimension. As the external magnetic field goes to zero the magnetization goes also to zero. in agreement with the Mermin-Wagner theorem. Close to the threshold value B = 2 the magnetization saturates to its maximum value 1/2 which corresponds to a perfect alignment of the spins. These results show that the meron-cluster algorithm represents a powerful tool to simulate frustrated magnetic systems close to zero temperature and encourages to use it to study the structure of quantum phase transitions. respectively; the diamond data corresponds to L = 80 for inverse temperature β = 20, the inverted triangle data corresponds to L = 160 for inverse temperature β = 30. and the emplty circle data corresponds to L = 240 at inverse temperature β = 30. All these simulations were performed using ǫ = 0.05. The dot line corresponds to Fabricius et al. data reported in [10], which was obtained based on scaling behaviour and on the Bethe Ansatz for the magnetization in the thermodynamic limit L → ∞ and T = 0.
As it is well known the antiferromagnetic Heisenberg (AFH) spin chain at zero temperature in the presence of an external magnetic field has a quantum phase transition (QPT) at the critical value B c = 2. At this critical value the correlation length diverges and some physical quantities like the first derivative of the entanglement entropy and the susceptibility are singular which are general features of QPT [11] . Above this critical value the ground state becomes ferromagnetic in the thermodynamic limit and below it the physics is governed by magnons and the aniferromagnetic order is preferred by the system.
In [10] it was shown that the magnetization curve develops a square root singularity as the magnetic field approaches its critical value B c = 2. Quantum entanglement as well as the concurrence have been widely use to find QPT in quantum systems [12]. In [13] the same physical system is studied by using exact diagonalization combined with the infinite time-evolving block decimation technique to compute the entanglement, the ground state energy and the magnetization as a function of the external magnetic field. It is shown that the quantum entanglement is sensitive to the subtle changing ground state and that it can be used to describe the magnetization and the QPT. In particular, they show numerical evidence supporting that the magnetization displays a square root singularity as B → 2, and its derivative becomes discontinuous at B c = 2 beyond it the system becomes ferromagnetic. It is further conjectured that this point corresponds to second order QPT.
In Fig. 3 the first derivative of the magnetization density is shown for the AFH spin chain for L = 80, 160 and 240 sites and inverse temperature β = 20 and β = 30. The behaviour of these curves allows to support the existence of a discontinuity very close to B c = 2 as L goes to infinity and at zero temperature, in agreement with [13]. In two dimensions we have studied the antiferromagnetic Heisenberg model in an external magnetic field defined on spin ladders and on square lattices. In Fig. 4 the magnetization density is displayed as a function of the external magnetic field B for two different square lattices of lattice sizes L = 12 and L = 16, for an inverse temperature β = 15. As a comparison, the numerical results obtained in [10] for the magnetization in the thermodynamic limit and at zero temperature are also displayed. In this reference a judicious law for the scaling behaviour for the ground state energy based on numerical results obtained by exact diagonalization on small lattices is proposed, and by differentiation, the magnetization in the thermodynamic limit and at zero temperature is attained. At low external magnetic field both results agree with each other with high accuracy and as B increases slightly deviations appear, which can be conjectured due to the lack of an exact analytical expression for the scaling behaviour of the ground state stressed in that paper. For larger values of the magnetic field B > 4 the magnetization value saturates to its maximum value 1/2, corresponding to a complete spin alignment. As the magnetic field goes to zero B → 0 the magnetization also goes to zero, again in agreement with the Mermin-Wagner theorem. Compared to the d = 1 case, the finite size effects are larger for large B-values, which is well known. We also computed the magnetization of the antiferromagnetic Heisenberg model on small lattices and high temperature values. In Fig. 5 the magnetization density as a function of the external magnetic field is plotted for a square lattice with size L = 4 at high temperature values T = 0.5, 1.0, 1.5 and 2.0. As a comparison, the sub-plot shows the results from exact numerical diagonalization obtained in [10] for the same parameter values as the ones used in our simulations. Within a remarkable accuracy both curves agree with each other, which shows how accurate and trustfully the meron-cluster algorithm works in the whole parameter region considered to perform numerical simulations. Antiferromagnetic spin ladders are interesting physical systems as they interpolate between single d = 1 spin chains and d = 2 quantum antiferromagnets. Their low-energy dynamics are accurately described by (1+1)-d quantum field theories, which can be exactly solved analytically using the Bethe Ansatz method. This correspondence was explicitly used in [7] to compare numerical results obtained by the meron-cluster algorithm with analytical predictions for the magnetization of spin ladder systems. We used these results also to compare with the results of our numerical simulations.
In Fig. 6 the magnetization as a function of the external magnetic field is displayed for a spin ladder of size 4 × 20 at inverse temperature β = 15. As a comparison, numerical results obtained for the same physical system by Chandrasekharan et al [7] is plotted. A remarkable agreement between both results is obtained in the whole range of the parameter values, which amounts to a correct implementation of the meron-cluster algorithm.  [7] for the same physical system at the same β value.
The magnetization density as a function of the external magnetic field displayed in logarithmic scale is shown in Fig. 7 for an antiferromagnetic spin ladder (L = 20, L ′ = 4) at low temperature β = 15. As already mentioned, spin ladders are spatially quasi one dimensional systems whose low-energy physics is governed by d = (1 + 1) dimensional quantum field theories. According to this connection the uniform magnetic field B corresponds to a chemical potential µ = B/c in the (1 + 1) − d quantum field theory, where c is the spin-wave velocity. Using the Bethe Ansatz allows obtaining an exact solution for the magnetization in this quantum field theories [19], which we will compare with numerical results along the lines proposed in [7]. If one constraints the simulations to an even number of coupled spin 1/2 chains, the corresponding quantum field theory describes an asymptotically free theory with a non-perturbatively generated mass-gap m, and the (1 + 1) − d O(3) effective action for the spin ladder contains no topological terms.
It is known that close to the threshold region µ ≈ m numerical results are very sensitive to both finite-size as well as finite-temperature effects. Nevertheless, in the region µ no too large the system is accurately described by a dilute fermion gas with magnetization density given by: For chemical potential slightly above the threshold value the magnetization in the thermodynamic limit can be written as: All these analytic expressions will be used to compare with our following numerical results. The magnetization density given by Eq. (4.1) is shown in Fig. 7 as a function of the external magnetic field (dot-dashed line) in logarithmic scale. Eq. (4.2) is represented as solid line and corresponds to its asymptotic value very close to the threshold region µ ≈ m. We have used the numerical values m = 0.141 and c = 1.657 known from references [20], [21] respectively. Finally, the circle-segmented line represents the numerical values obtained by the meron-cluster simulation. It can be concluded that the numerical results agree remarkably well close to the threshold region with the behaviour predicted by the analytical expressions of the effective field theory. The deviation of the numerical results from the behaviour depicted by Eq. (4.1) is due to the fact that this equation holds only for sufficiently large number L ′ >> 1, which is not fulfilled in the spin-ladder we have simulated. Our numerical data agrees perfectly with the results obtained in [7] within an error less than 2 percent.

V. CONCLUSIONS
In this article we have simulated numerically the antiferromagnetic quantum Heisenberg model in a uniform external field in one and two spacial dimensions by using the meron-cluster algorithm proposed in [7]. Our results agree remarkably well with previous results obtained in [10] by using suitable scaling behaviour formulae for the energy eigenvalues and performing numerically a thermodynamic limit. We have also simulated the quantum spin ladders in a magnetic field studied in [7], and found a perfect agreement with their numerical results with numerical differences less than a few percent. The present study allows to confirm that the meron-cluster concept applied to quantum spin systems in the presence of an arbitrary external magnetic field represents an efficient solution to the sign problem, which is present in both fermionic and bosonic models. Its efficiency encourages to apply it to the study of systems having quantum phase transitions since in these systems there are a dynamical competition between the different interactions, which produces a non-unique ground state at very low temperature. This dynamical competition leads to a severe sign problem when trying to perform numerical simulations. As a first application, we have used the meron-cluster algorithm to simulate the one dimensional antiferromagnetic Heisenberg model in an external magnetic field close to the critical magnetic field B c = 2 and found numerical evidence for a discontinuous first derivative of the magnetization at this value in the infinite volume limit and at zero temperature, which is consistent with a second-order quantum phase transition conjectured in [13]. Nevertheless, and in order to address this question new and more accurate simulations should be performed, which we plan to carry out in future investigations.