1. Introduction
Since the advent of the theoretical description of classical and quantum phase transitions (QPTs), long-range interactions between degrees of freedom challenged the established concepts and propelled the development of new ideas in the field [
1,
2,
3,
4,
5]. It is remarkable that only a few years after the introduction of the renormalisation group (RG) theory by K. G. Wilson in 1971 as a tool to study phase transitions and as an explanation for universality classes [
6,
7,
8,
9,
10,
11], it was used to investigate ordering phase transitions with long-range interactions. These studies found that the criticality depends on the decay strength of the interaction [
1,
2,
3]. It then took two decades to develop numerical Monte Carlo (MC) tools capable to simulate basic magnetic long-range models with thermal phase transitions following the behaviour predicted by the RG theory [
12,
13]. The results of these simulations sparked a renewed interest in finite-size scaling above the upper critical dimension [
12,
14,
15,
16,
17,
18,
19] since "hyperscaling is violated" [
13] for long-range interactions that decay slowly enough. In this regime, the treatment of dangerous irrelevant variables (DIVs) in the scaling forms is required to extract critical exponents from finite systems.
Meanwhile, a similar historic development took place regarding the study of QPTs under the influence of long-range interactions. By virtue of pioneering RG studies [
20,
21], the numerical investigation of long-range interacing magnetic systems has been triggered [
22,
23,
24,
25,
26,
27,
28,
29,
30,
31,
32,
33,
34,
35,
36,
37,
38]. In particular, Monte Carlo based techniques became a popular tool to gain quantitative insight into these long-range interacting quantum magnets [
22,
25,
26,
29,
30,
31,
32,
33,
34,
35,
36,
37,
38,
39,
40]. On the one hand this includes high-order series expansion techniques, where classical Monte Carlo integration is applied for the graph embedding scheme, allowing to extract energies and observables in the thermodynamic limit [
25,
29]. On the other hand, there is stochastic series expansion quantum Monte Carlo [
39], which enables calculations on large finite systems. To determine physical properties of the infinite system, finite-size scaling is performed with the results of these computations. Inspired by the recent developments for classical phase transitions [
15,
16,
17,
18,
19,
41], a theory for finite-size scaling above the upper critical dimension for QPTs was introduced [
32,
34].
When investigating algebraically decaying long-range interactions
with the distance
r and the dimension
d of the system, there are two distinct regimes: One for,
(strong long-range interaction) and another one for
(weak long-range interaction) [
5,
42,
43,
44,
45]. In the case of strong long-range interactions, common definitions of internal energy and entropy in the thermodynamic limit are not applicable and standard thermodynamics breaks down [
5,
42,
43,
44,
45]. We will not focus on this regime in this review. For details specific to strong long-range interactions, we refer to other review articles such as Refs. [
5,
42,
43,
44,
45]. For the sake of this work, we restrict the discussion to weak long-range interaction or competing antiferromagnetic strong long-range interactions, for which an extensive ground-state energy can be defined without rescaling of the coupling constant [
5].
The interest in quantum magnets with long-range interactions is further fueled by the relevance of these models in state-of-the-art quantum-optical platforms [
5,
46,
47,
48,
49,
50,
51,
52,
53,
54,
55,
56,
57,
58,
59,
60,
61,
62,
63,
64,
65,
66,
67,
68,
69,
70,
71,
72,
73,
74,
75,
76,
77,
78,
79,
80,
81,
82,
83,
84,
85,
86,
87]. To realise long-range interacting quantum lattice models with a tunable algebraic decay exponent, one can use trapped ions which are coupled off-resonantly to motional degrees of freedom [
5,
81,
82,
83,
84,
85,
88]. Another possibility is to couple trapped neutral atoms to photonic modes of a cavity [
5,
86,
87]. Alternatively, one can realise long-range interactions decaying with a fixed algebraic decay exponent of six or three using Rydberg atom quantum simulators [
46,
47,
48,
49,
50,
51,
52,
53,
54,
55] or ultracold dipolar quantum atomic or molecular gases in optical lattices [
56,
57,
58,
59,
60,
61,
62,
63,
64,
65,
66,
67,
68,
69,
70,
71,
72,
73]. Note, in many of the above listed cases it is possible to map the long-range interacting atomic degrees of freedom onto quantum spin models [
5,
52,
89]. Therefore, they can be exploited as analogue quantum simulators for long-range interacting quantum magnets and the relevance of the theoretical concepts transcends the boundary between the fields.
From the perspective of condensed matter physics, there are multiple materials with relevant long-range interactions [
90,
91,
92,
93,
94,
95,
96,
97,
98,
99,
100,
101,
102,
103,
104,
105,
106]. The compound LiHoF
4 in an external field realises an Ising magnet in a transverse magnetic field [
102,
103,
104,
105]. A recent experiment with the two-dimensional Heisenberg ferromagnet Fe
3GeTe
2 demonstrates that phase transitions and continuous symmetry breaking can be implemented by circumventing the Hohenberg-Mermin-Wagner theorem with long-range interactions [
106]. This material is in the recently discovered material class of 2d magnetic van der Waals systems [
107,
108]. Further, dipolar interactions play a crucial role in the spin-ice state in frustrated magnetic pyrochlore materials Ho
2Ti
2O
7 and Dy
2Ti
2O
7 [
90,
91,
92,
93,
94,
95,
96,
97,
98,
99,
100,
101].
In this review, we are interested in physical systems described by quantum spin models, where the magnetic degrees of freedom are located on the sites of a lattice. We concentrate on the following three paradigmatic types of magnetic interactions between lattice sites: First, Ising interactions where the magnetic interaction is oriented only in the direction of one quantisation axis. Second, XY interactions with a -symmetric magnetic interaction invariant under planar rotations and, third, Heisenberg interactions with a -symmetric magnetic interaction invariant under rotations in 3d spin space. In the microscopic models of interest a competition between magnetic ordering and trivial product states, external fields, or quasi long-range order leads to QPTs.
In this context, the primary research pursuit revolves around how the properties of the QPT depend on the long-range interaction. The upper critical dimension of a QPT in magnetic models with non-competing algebraically decaying long-range interactions is known to depend on the decay exponent of the interaction for a small enough exponent, and decreases as the decay exponent decreases [
20,
21]. If the dimension of a system is equal or exceeds the upper critical dimension, the QPT displays mean-field critical behaviour. At the same time standard finite-size scaling as well as standard hyperscaling relations are no longer applicable. Therefore, these systems are primary workhorse models to study finite-size scaling above the upper critical dimension. In this case the numerical simulation of these systems is crucial in order to gauge novel theoretical developments. Further, systems with competing long-range interactions do not tend to depend on the long-range nature of the interaction [
23,
24,
25,
26,
29,
30,
32]. In several cases long-range interactions then lead to the emergence of ground states and QPT which are not present in the corresponding short-range interacting models [
27,
30,
54,
55,
109,
110,
111,
112].
In this review, we are mainly interested in the description and discussion of two Monte Carlo based numerical techniques, which were successfully used to study the low-energy physics of long-range interacting quantum magnets, in particular with respect to the quantitative investigation of QPTs [
22,
25,
29,
30,
31,
32,
34,
35,
36,
37,
38,
40]. We chose to review this topic due to our personal involvement with the application and development of these methods [
25,
29,
30,
31,
32,
34,
35]. On one hand, we explain in detail how classical Monte Carlo integration can enhance the capabilities of linked-cluster expansions (LCEs) with the pCUT+MC approach (a combination of the perturbative unitary transforms approach (pCUT) and MC embedding). On the other hand, we describe how stochastic series expansion (SSE) quantum Monte Carlo (QMC) is used to directly sample the thermodynamic properties of suitable long-range quantum magnets on finite systems.
This review is structured as follows. In
Section 2 we review the basic concept of a QPT in a condensed way, focusing on the details relevant for this review. We define the quantum critical exponents and the relations between them in
Section 2.1. Here, we also have the first encounter with the generalised hyperscaling relation which is also valid above the upper critical dimension where conventional hyperscaling breaks down. As the SSE QMC discussed in this review is a finite-system simulation, we discuss the conventional finite-size scaling below the upper critical dimension in
Section 2.2 and the peculiarities of finite-size scaling above the upper critical dimension in
Section 2.3. In
Section 3, we summarise the basic concepts of Markov chain Monte Carlo integration: Monte Carlo sampling, Markovian random walks, stationary distributions, the detailed balance condition, and the Metropolis-Hastings algorithm. We continue by introducing the series-expansion Monte Carlo embedding method pCUT+MC in
Section 4. We start with the basic concepts of a graph expansion in
Section 4.1 and introduce the perturbative method of our choice, the perturbative continuous unitary transformation method, in
Section 4.2. We introduce the theoretical concepts for setting up a linked-cluster expansion as a full graph decomposition in
Section 4.3 and subsequently discuss how to practically calculate perturbative contributions in
Section 4.4 and
Section 4.5. We prepare the discussion of the white-graph decomposition in
Section 4.6 with an interlude on the relevant graph theory in
Section 4.6.1 and
Section 4.6.2 and the important concept of white graphs in
Section 4.6.3. Further, in
Section 4.7 we discuss the embedding problem for the white-graph contributions. Starting from the nearest-neighbour embedding problem in
Section 4.7.1, we generalise it to the long-range case in
Section 4.7.2 and then introduce a classical Monte Carlo algorithm to calculate the resulting high-dimensional sums in
Section 4.7.3. This is followed by some technical aspects on series extrapolations in
Section 4.8 and a summary of the entire workflow in
Section 4.9. In the next section the topic changes towards the review of the SSE, which is an approach to simulate thermodynamic properties of suitable quantum many-body systems on finite systems at a finite temperature. First, we discuss the general concepts of the method in
Section 5. We review the algorithm to simulate arbitrary transverse-field Ising models introduced by A. Sandvik [
39] in
Section 5.1. We then review an algorithm used to simulate non-frustrated Heisenberg models in
Section 5.2. After the introduction to the algorithms, we summarise techniques on how to measure common observables in the SSE QMC scheme in
Section 5.3. Since the SSE QMC method is a finite temperature method, we discuss how to rigorously use this scheme to perform simulations at effective zero temperature in
Section 5.4. We conclude this section with a brief summary of path integral Monte Carlo techniques used for systems with long-range interactions (see
Section 5.5). To maintain the balance between algorithmic aspects and their physical relevance, we summarise several theoretical and numerical results for quantum phase transitions in basics long-range interacting quantum spin models, for which the discussed Monte Carlo based techniques provided significant results. First, we discuss long-range interacting transverse-field Ising models in
Section 6. For ferromagnetic interactions this model displays three regimes of universality: A long-range mean-field regime for slowly decaying long-range interactions, an intermediate long-range non-trivial regime, and a regime of short-range universality for strong decaying long-range interactions. We discuss the theoretical origins of this behaviour in
Section 6.1.1 and numerical results for quantum critical exponents in
Section 6.1.2. Since this model is a prime example to study scaling above the upper critical dimension in the long-range mean-field regime, we emphasise these aspects in
Section 6.1.3. Further, we discuss the antiferromagnetic long-range transverse-field Ising model on bipartite lattices in
Section 6.2 and on non-bipartite lattices in
Section 6.3. The next obvious step is to change the symmetry of the magnetic interactions. Therefore, we turn to long-range interacting XY models in
Section 7 and Heisenberg models in
Section 8. We discuss the long-range interacting transverse-field XY chain in
Section 7 starting with the
-symmetric isotropic case in
Section 7.1, followed by the anisotropic case for ferromagnetic (see
Section 7.2) and antiferromagnetic (see
Section 7.3) interactions which display similar behaviour to the long-range transverse-field Ising model on the chain discussed in
Section 6. We conclude the discussion of results with unfrustrated long-range Heisenberg models in
Section 8. We focus on the staggered antiferromagnetic long-range Heisenberg square lattice bilayer model in
Section 8.1 followed by long-range Heisenberg ladders in
Section 8.2 and the long-range Heisenberg chain in
Section 8.3. We conclude in
Section 9 with a brief summary and with some comments on next possible steps in the field.
2. Quantum Phase Transitions
This review is part of the special issue with the topic "Violations of Hyperscaling in Phase Transitions and Critical Phenomena". In this work, we summarise investigations of low-dimensional quantum magnets with long-range interactions targeting in particular quantum phase transitions (QPTs) above the upper critical dimension, where the naive hyperscaling relation is no longer applicable. In this section, we recapitulate the relevant aspects of QPTs needed to discuss the results of the Monte Carlo based numerical approaches. First, we give a general introduction to QPTs. After that, we discuss in detail the definition of critical exponents and relations among them in
Section 2.1 as well as the scaling below (see
Section 2.2) and above (see
Section 2.3) the upper critical dimension.
Any non-analytic point of the ground-state energy of an infinite quantum system as a function of a tuning parameter
is identified with a QPT [
113]. This tuning parameter can, for instance, be a magnetic field or pressure but not the temperature. Quantum phase transitions are a concept of zero temperature as there are no thermal fluctuations and all excited states are supressed infinitely strong such that the system remains in its ground state. There are two scenarios how a non-analytic point in the ground-state energy can emerge [
113]: First, an actual (sharp) level crossing between the ground-state energy and another energy level. Second, the non-analytic point can be considered as a limiting case of an avoided level crossing.
In this review, we are interested in second-order QPTs, which fall into the second scenario. At a second-order QPT, the relevant elementary excitations condense into a novel ground state while the characteristic length and time scales diverge. Apart from topological QPTs involving long-range entangled topological phase, such continuous transitions are described by the concept of spontaneous symmetry breaking. On one side of the QPT the ground state obeys a symmetry of the Hamiltonian, while on the other side this symmetry is broken in the ground state and a ground-state degeneracy arises.
Following the idea of the quantum-to-classical mapping [
113,
114],
d-dimensional quantum systems can be mapped in the vicinity of a second-order QPT to models of statistical mechanics with a classical (thermal) second-order phase transition in
dimensions. In many cases, the models obtained from a quantum-to-classical mapping are rather artificial [
113]. However, such mappings often allow to categorise QPTs in terms of universality classes and associated critical exponents by the non-analytic behaviour of the classical counterparts [
10,
11,
113,
115,
116]. The mapping further illustrates that the renormalisation group (RG) theory is also applicable to describe QPTs [
10,
11,
113,
115,
116].
In the RG theory, each QPT belongs to a non-trivial
1 fixed point of the RG transformation [
10,
11]. Critical exponents are connected to the RG flow in the immediate vicinity of these non-trivial fixed points [
10,
11,
113]. The concept of universality classes arises from the fact that different microscopic Hamiltonians can have a quantum critical point that is attracted by the same non-trivial fixed point under successive application of the RG transformation [
10,
11]. Due to this, the QPTs in these models have the same critical exponents.
Another remarkable result of the RG theory is the scaling of observables in the vicinity of phase transitions. Historically, the theory of scaling at phase transitions was heuristically introduced before the RG approach [
117,
118,
119,
120,
121,
122,
123]. The latter provided the theoretical foundation for the scaling hypothesis [
6,
7]. The main statement of the scaling theory is that the non-analytic contributions to the free energy and correlation functions are mathematically described by generalised homogeneous functions (GHF) [
123]. A function with
n variables
is called a GHF, if there exist
with at least one being non-zero and
such that for
The exponents
are the scaling powers of the variables and
is the scaling power of the function
f itself. An in-depth summary of the mathematical properties of GHFs can be found in
Appendix B. The most important properties of GHFs are that its derivatives, Legendre transforms, and Fourier transforms are also GHFs. As we will outline in
Section 2.2, the theory of finite-size scaling is formulated in terms of GHFs and relates the non-analytic behaviour at QPTs in the thermodynamic limit with the scaling of observables for different finite system sizes. In this, the variables
are related to physical parameters like the temperature
T, control parameter
, symmetry-breaking field
H and also irrelevant, more abstract parameters that parameterise microscopic details of the model like the lattice spacing. Later in this section, we will define irrelevant variables in the context of RG and GHFs.
Another aspect relevant for this work is that quantum fluctuations are the driving factor with QPTs [
113]. In general, fluctuations are more important in low dimensions [
116]. The universality class of QPTs for a certain symmetry breaking depends on the dimensionality of the system.
An important aspect regarding this review is the so-called upper critical dimension
. The upper critical dimension is defined as a dimensional boundary such that for systems with dimension
, the critical exponents are those obtained from mean-field considerations. The upper critical dimension is of particular importance for QPTs in systems with non-competing long-range interactions. For sufficiently small decay exponents of an algebraically decaying long-range interaction
, the upper critical dimension starts to decrease as a function of the decay exponent
[
20,
32,
34]. In the limiting case of a completely flat decay (
) of the long-range interaction resulting in an all-to-all coupling, the model is intrinsically of mean-field type and mean-field considerations become exact. For a certain value of the decay exponent, the upper critical dimension becomes equal to the fixed spatial dimension, and for decay exponents below this value, the dimension of the system is above the upper critical dimension of the transition [
20,
32,
34]. This makes long-range interacting systems an ideal test bed for studying phase transitions above the upper critical dimension in low-dimensional systems. In particular, long-range interactions can make the upper critical dimension accessible in real-world experiments as the upper critical dimension of short-range models is usually not below three.
Although phase transitions above the upper critical dimension display mean-field criticality, they are still a matter worth studying, since naive scaling theory describing the behaviour of finite systems close to a phase transition (see
Section 2.2) is no longer applicable [
16,
19,
124]. Moreover, the naive versions of some relations between critical exponents, as discussed in
Section 2.1, do not hold anymore [
15,
16]. The reason for this issue are the dangerous irrelevant variables (DIVs) in the RG framework [
125,
126,
127]. During the application of the RG transformation, the original Hamiltonian is successively mapped to other Hamiltonians which can have infinitely many couplings. All these couplings, in principle, enter the GHFs. In practice, all but a finite number of these couplings are set to zero since their scaling powers are negative which means they flow to zero under renormalisation. These couplings are therefore called irrelevant. This approach of setting irrelevant couplings to zero can be used to derive the finite-size scaling behaviour as described in
Section 2.2. However, above the upper critical dimension, this approach breaks down because it is only possible to set irrelevant variables to zero if the GHF does not have a singularity in this limit [
125]. Above the upper critical dimension such singularities in irrelevant parameters exist, which makes them DIVs [
126]. We explain the effect of DIVs on scaling in
Section 2.3.
2.1. Critical Exponents in the Thermodynamic Limit
As outlined above, a second-order QPT comes with a singularity in the free energy density. In fact, also other observables experience singular behaviour at the critical point in the form of power law singularities. For instance, the order parameter
m as a function of the control parameter
behaves as
in the ordered phase. Without loss of generality, the system is taken to be in the ordered phase for
and the notation
means that
is approaching
from below, i. e. it is approaching in the ordered phase. In the disordered phase
, the order parameter by definition vanishes such that
. The observables with their respective power-law singular behaviour, that is characterised by the critical exponents
and
z, are summarised in
Table 1 together with how they are commonly defined in terms of the free energy density
f, the symmetry-breaking field
H, that couples to the order parameter, and the reduced control parameter
.
One usually defines reduced parameters like
r that vanish at the critical point not only to shorten the notation but also to express the power-law singularities independent of the microscopic details of the specific model one is looking at. While the value of
depends on these details, the power-law singularities are empirically known to not depend on the microscopic details but only on more general properties like the dimensionality, the symmetry that is being broken and, with particular emphasis due to the focus of this review, on the range of the interaction. It is therefore common to classify continuous phase transitions in terms of universality classes. These universality classes share the same set of critical exponents. In terms of RG, this behaviour is understood as distinct critical points of microscopically different models flowing to the same renormalisation group fixed point, which determines the criticality of the system [
6,
7,
10]. Prominent examples for universality classes of magnets are the 2d- and 3d-Ising (
symmetry), 3d XY (
symmetry), 3d Heisenberg (
symmetry) universality classes [
113]. It is important to mention that the dimension in the classifications are referring to classical and not quantum systems and they should not be confused with each other. In fact, the universality class of a short-range interacting non-frustrated quantum Ising model of dimension
d lies in the universality of the classical
-dimensional Ising model.
There are only a few dimensions for which a separate universality class is defined for the different models. For lower dimension, the fluctuations are too strong in order for a spontaneous symmetry breaking to occur. In case of the classical Ising model, there is no phase transition for 1d, while for the classical XY and Heisenberg models with continuous symmetries there is not even a phase transitions for 2d due to the Hohenberg-Mermin-Wagner (HWM) theorem [
129,
130,
131]. This dimensional boundary is referred to as lower critical dimension
dlc. The lower critical dimension is the highest dimension for which no transition occurs, i. e.
dlc for the Ising model and
dlc for the XY and Heisenberg model. For higher dimensions
, the critical exponents of the mentioned models do not depend on the dimensionality anymore and they take on the mean-field critical exponents in all dimensions. The underlying reason is that with increasing dimensions the local fluctuations become smaller due to the higher connectivity of the system [
132]. This has been also exploited in diagrammatic and series expansions in
[
133,
134,
135]. This dimensional boundary, at which the criticality becomes the mean-field one, is called upper critical dimension
duc. Usually, the upper critical dimension is too large to realise a system above its upper critical dimension in the real world. However, long-rang interactions can increase the connectivity of a system in a similar sense as the dimensionality. Sufficiently long-ranged interaction can therefore lower the upper critical to a value that is accessible in experiments.
Finally, it is worth mentioning that the critical exponents are not independent from each other but obey certain relations [
128], namely
The first relation in Eq. (3) is the so-called hyperscaling relation whose classical analogue (without
z) was introduced by Widom [
10,
136]. The Essam-Fisher relation in Eq. (4) [
137,
138] is reminiscent of a similar inequality that was proven rigorously by Rushbrooke using thermodynamic stability arguments. Eq. (5) is called Widom relation. The last relation in Eq. (6) is the Fisher scaling relation which can be derived using the fluctuation-dissipation theorem [
10,
128,
138]. Those relations were originally obtained from scaling assumptions of observables close to the critical point which were only later derived rigorously when the RG formalism was introduced to critical phenomena [
10,
128]. Due to these relations, it is sufficient to calculate three, albeit not arbitrary, exponents to obtain the full set of critical exponents.
The hyperscaling relation Eq. (3) is the only relation containing the dimension of the system and is therefore often said to break down above the upper critical dimension where one expects the same mean-field critical exponents independent of the dimension
d [
10]. It therefore deserves special focus in this review since the long-range models discussed will be above the upper critical dimension in certain parameter regimes. Personally, we would not agree that the hyperscaling relation breaks down above the upper critical dimension, but we would rather call Eq. (3) a special case of a more general hyperscaling relation
with the pseudo-critical exponent Ϙ ("koppa") [
34]
Below the upper critical dimension, the general hyperscaling relation therefore relaxes to Eq. (3). Above the upper critical dimension the relation becomes
which is independent of the dimension of the system. For the derivation of this generalised version of the hyperscaling relation for QPTs, see
Section 2.3 or Ref. [
34]. The derivation of the classical counterpart can be found in Ref. [
15] and is reviewed in Ref. [
41].
2.2. Finite-Size Scaling below the Upper Critical Dimension
Even though the singular behaviour of observables at the critical point is only present in the thermodynamic limit, it is possible to study the criticality of an infinite system by investigating their finite counterparts. In finite systems, the power-law singularities of the infinite system are rounded and shifted with respect to the critical point, e.g. the susceptibility with its characteristic divergence at the critical point
is deformed to a broadened peak of finite height. The peak’s position
is shifted with respect to the critical point
. A possible definition of a pseudo-critical point of a finite system is the peak position
. As the characteristic length scale of fluctuations
diverges at the critical point, the finite system approaching the critical point will at some point begin to "feel" its finite extent and the observables start to deviate from the ones in the thermodynamic limit. As
diverges with the exponent
like
at the critical point, the extent of rounding in a finite system is related to the value of
. Similarly, the peak magnitude of finite-size observables at the pseudo-critical point will depend on how strong the singularity in the thermodynamic limit is, which means it depends on the respective critical exponents
,
,
, and
. The shifting, rounding, and the varying peak magnitude are shown for the susceptibility of the long-range transverse-field Ising model in
Figure 1. This dependence of observables in finite systems on the criticality of the infinite system is the basis of finite-size scaling.
In a more mathematical sense, the relation between critical exponents and finite-size observables has its origin in the renormalisation group (RG) flow close to the corresponding RG fixed point that determines the criticality [
139]. Close to this fixed point, one can linearise the RG flow so that the free energy density and the characteristic length
become generalised homogeneous functions (GHFs) in their parameters [
123,
127,
128,
140,
141]. For a thorough discussion of the mathematical properties of GHFs we refer to Ref. [
123] and
Appendix B. This means, the free energy density
f and characteristic length scale
as functions of the couplings
and the inverse system length
obey the relations
with the respective scaling dimensions
,
, and
governing the linearised RG flow with spatial rescaling factor
around the RG fixed point, at which all couplings vanish by definition. All of those couplings are relevant except for
u which denotes the leading irrelevant coupling [
10,
113]. Relevant couplings are related to real-world parameters that can be used to tune our system away from the critical point like the temperature
T, a symmetry-breaking field
H, or simply the control parameter
r. The irrelevant couplings
u do not per se vanish at the critical point like the relevant ones do. However, they flow to zero under the RG transformation and are commonly set to zero in the scaling laws
by assuming analyticity in these parameters. The generalised homogeneity of thermodynamic observables can be derived from the one of the free energy density
f. For example, the generalised homogeneity of the magnetisation
can be derived by taking the derivative of
f with respect to the symmetry-breaking field
H.
By investigating the singularity of
in
r via Eq. (13), one can show that the scaling power
of the control parameter
r is related to the critical exponent
by
[
113]. For this, one fixes the value
of the first argument to
in the right-hand side of Eq. (13) by setting
such that
Analogously, further relations between the scaling powers and other critical exponents can be derived by looking at the singular behaviour of the respective observables in the corresponding parameters. Overall, one further gets
From these equations, one can already tell that the critical exponents are not independent from each other. In fact, the scaling relations
,
and
(see Eqs. (3) (5)) can be derived from Eq. (16) and
. By expressing the RG scaling powers
in terms of critical exponents, the homogeneity law for an observable
with a bulk divergence
is given by
In order to investigate the dependence on the linear system size
L, the last argument in the homogeneity law is fixed to
by inserting
. This readily gives the finite-size scaling form
with
being the universal scaling function of the observables
. The scaling function
itself does not depend on
L anymore, but in order to compare different linear system sizes one has to rescale its arguments. To extract the critical exponents from finite systems, the observable
is measured for different system sizes
L and parameters
close to the critical point
. The
L-dependence according to Eq. (19) is then fitted with the critical exponents
,
,
and
z as well as the critical point
, which is hidden in the definition of
r, as free parameters. It is advisable to fix two of the three parameters
to their critical values in order to minimise the amount of free parameters in the fit. For example, with
and only
one can extract the two critical exponents
and
alongside
. For further details on a fitting procedure, we refer to Ref. [
32]. If one knows the critical point, one can also set
and look at the
L-dependent scaling
directly at the critical point to extract the exponent ratio
. There are many more possible approaches to extract critical exponents from the FSS law in Eq. (19) [
142,
143,
144]. For relatively small system sizes, it might be required to take corrections to scaling into account [
143,
144].
2.3. Finite-Size Scaling above the Upper Critical Dimension
In the derivation of finite-size scaling below the upper critical dimension, it was assumed that the free energy density is an analytic function in the leading irrelevant coupling u and therefore one can set it to zero. However, this is not the case above the upper critical dimension anymore and the free energy density f is singular at . Due to this singular behaviour u is referred to as a dangerous irrelevant variable (DIV).
As a consequence, one has to take the scaling of
u close to the RG fixed point into account. This is done by absorbing the scaling of
f in
u for small
u into the scaling of the other variables [
127]
up to a global power
of
u. This leads to a modification of the scaling powers in the homogeneity law for the free energy density [
127]
with the modified scaling powers [
34,
127]
In the classical case [
127],
yL* was commonly set to 1 by choice. This is justified because the scaling powers of a GHF are only determined up to a common non-zero factor [
123]. However, for the quantum case [
34], this was kept general as it has no impact on the FSS.
As the predictions from Gaussian field-theory and mean-field differed for the critical exponents
,
, and
but not for the "correlation" critical exponents
,
z,
, and
[
145], the correlation sector was thought not to be affected by DIVs at first [
127,
142,
145]. Later Q-FSS, another approach to FSS above the upper critical dimension, pioneered by Ralph Kenna and his colleagues, was developed for classical [
15,
19,
124] as well as for quantum systems [
34] which explicitly allowed the correlation sector to also be affected by DIV. In analogy to the free energy density, the homogeneity law of the characteristic length scale is then also modified to
with
in order to reproduce the correct bulk singularity
. A new pseudo-critical exponent Ϙ ("koppa")
is introduced. This exponent describes the scaling of the characteristic length scale with the linear system size. This non-linear scaling of
with
L is one of the key differences to the previous treatments above the upper critical dimension in Ref. [
127].
Analogous to the case below the upper critical dimension, the modified scaling powers
can be related to the critical exponents:
By using the mean-field critical exponents for the
quantum rotor model, one gets restrictions for the ratios of modified scaling powers
Furthermore, one can link the bulk scaling powers
, and
to the scaling power
yL* of the inverse linear system size [
34]
by looking at the scaling of the susceptibility in a finite system [
34,
127]. This relation is crucial for deriving a FSS form above the upper critical dimension as the modified scaling power
yL* or rather its relation to the other scaling powers determines the scaling with the linear system size
L. For details on the derivation, we refer to Ref. [
34]. We want to stress again that the scaling powers of GHFs are only determined up to a common non-zero factor [
123]. Therefore, it is evident that one can only determine the ratios of the modified scaling powers but not their absolute value. The absolute values are subject to choice. Different choices were discussed in Ref. [
34], but these choices rather correspond to taking on different perspectives and have no impact on FSS or the physics.
From Eq. (32) together with Eqs. (26) and (27), a generalised hyperscaling relation
can be derived. This also determines the pseudo-critical exponent
Finally, we can express the modified scaling powers in the FSS law for an observable
with power-law singularity
For
below the upper critical dimension, Eq. (36) relaxes to the standard FSS law Eq. (19). The scaling in temperature has not yet been studied for finite quantum systems above the upper critical dimension. However, in Ref. [
34] it was conjectured that
based on Eq. (32), which is also in agreement with
z being the dynamical critical exponent that determines the space-time anisotropy
as we will shortly see. This means that the finite-size gap scales as
with the system size [
34]. Of particular interest is the scaling of the characteristic length scale above the upper critical dimension, for which the modified scaling law Eq. (36) also holds with
. Hence, the characteristic length scale in dependence of the control parameter
r scales like
with the scaling function
. Directly at the critical point
, this leads to
. Comparing this with the scaling of the inverse finite-size gap
verifies that
z still determines the space-time anisotropy. Prior to Q-FSS [
15,
17], the characteristic length scale was thought to be bound by the linear system size
L [
127]. However, this was shown not to be the case by measurements of the characteristic length scale for the classical five-dimensional Ising model [
14] and for long-range transverse-field Ising models [
34].
For the latter, the data collapse of the correlation length according to Eq. (37) is shown in
Figure 2 as an example.
3. Monte Carlo Integration
In this section we provide a brief introduction to Monte Carlo integration (MCI). We focus on the aspects of Markov chain MCI as the basis to formulate the white-graph Monte Carlo embedding scheme of the pCUT+MC method in
Section 4 and the stochastic series expansion (SSE) quantum Monte Carlo (QMC) algorithm in
Section 5 in a self-contained fashion. MCI is the foundation for countless numerical applications which require the integration over high-dimensional integration spaces. As this review has a focus on "Monte Carlo based techniques for quantum magnets with long-range interactions", we forward readers with a deeper interest in the fundamental aspects of MCI and Markov chains to Refs. [
146,
147,
148].
MCI summarises numerical techniques to find estimators for integrals of functions
over an integration space
using random numbers. The underlying idea behind MCI is to estimate the integral, or the sum in the case of discrete variables, of the function
f over the configuration space by an expectation value
with
sampled according to a probability density function
(PDF) and the function
reweighted by the PDF. A famous direct application of this idea is the calculation of number "pi" which is in great detail discussed in Ref. [
148].
In this review, MCI is used for the embedding of white graphs on a lattice to evaluate high orders of a perturbative series expansion or to calculate thermodynamic observables using SSE. In both cases, non-normalised relative weights
within a configuration space
arise which are used for the sampling of the PDF
P
being oftentimes not directly accessible. In the context of statistical physics,
is often chosen to be the relative Boltzmann weight
of each configuration
. While this relative Boltzmann weight is accessible as long as
is known, the full partition function to normalise the weights is in general not.
In order to efficiently sample the configuration space
according to the relative weights, the methods in this review use a Markovian random walk to generate
. Let
be the random state of a random walk at a discrete step
n. The state
at the next step is randomly determined according to the conditional probabilities
(transition probabilities). This transition probabilities are normalised by
Markovian random walks obey the Markov property, which means the random walk is memory free and the transition probability for multiple steps factorises into a product over all time steps
with
the start configuration. We require the Markovian random walk to fulfil the following conditions: First, the random walk should have a certain PDF
defined by the weights
in Eq. (
39) as a stationary distribution. By definition,
is a stationary distribution of the Markov chain if it satisfies the global balance condition
Second, we require the random walk to be irreducible, which means that the transition graph must be connected and every configuration
can be reached from any configuration
in a finite number of steps. This property is necessary for the uniqueness of the stationary distribution [
147]. Lastly, we require the random walk to be aperiodic (see Ref. [
147] for a rigorous definition). Together with the irreducibility condition, this ensures convergence to the stationary distribution [
147].
There are several possibilities to design a Markov chain with a desired stationary distribution [
146,
147,
148,
149,
150,
151]. Commonly, the Markov chain is constructed to be reversible. This means that it satisfies the detailed balance condition
which is a stronger condition for stationarity of
P than the global balance condition in Eq. (42). One popular choice for the transition probabilities
that satisfies the detailed balance condition is given by the Metropolis-Hastings algorithm. Most applications of MCI reviewed in this work are based on the Metropolis-Hastings algorithm [
149,
150]. In this approach, the transition probabilities
are decomposed into propositions
and acceptance probabilities
as follows
The probabilities to propose a move
can be any random walk satisfying the irreducibility and aperiodicity condition. By inserting the decomposition of the transition probabilities Eq. (
44) into the detailed balance condition Eq. (
43), one obtains for the acceptance probabilities
where in the last step it was used that the unknown normalisation factors (see Eq. (
39)) of the PDF cancel. The condition in Eq. (
45) is fulfilled by the Metropolis-Hastings acceptance probabilities [
150]
For the special case, that the proposition probabilities are symmetric
, Eq. (46) reduces to the Metropolis acceptance probabilities
As an example, we regard a classical thermodynamic system with Boltzmann weights
given by the energies of configurations and the inverse temperature to give an intuitive interpretation of the Metropolis acceptance probabilities in Eq. (
47). The proposition to move from a configuration
to a configuration
with a smaller energy
is always accepted independent of the temperature. On the other hand, the proposition to move to a configuration
with a larger energy than
is only accepted with a probability depending on the ratio of the Boltzmann weights. If the temperature is higher, it is more likely to move to states with a larger energy. This reflects the physics of the system in the algorithm, focusing on the low-energy states at low-temperatures and going to the maximum entropy state at large temperatures.
4. Series-Expansion Monte Carlo Embedding
In this section we provide a self-contained and comprehensive overview of linked-cluster expansions for long-range interacting systems [
25,
29,
30,
31,
34,
35] using white graphs [
152] in combination with Monte Carlo integration for the graph embedding. First, we introduce linked-cluster expansions (LCEs) and discuss perturbative continuous unitary transformations (pCUT) [
153,
154] as a suitable high-order series expansion method. We then establish an adequate formalism for setting up LCEs and discuss the calculation of suitable physical quantities in practice. With the help of white graphs we can employ LCEs for models with long-range interactions and use the resulting contributions in a the Monte Carlo algorithm to deal with the embedding problem posed by long-range interactions. This approach we dub pCUT+MC.
4.1. Motivation and Basic Concepts
The goal of all our efforts is to calculate physical observables in the thermodynamic limit, i.e., for an infinite lattice
, using high-order series expansions. The starting point of every perturbative problem is
where the Hamiltonian describing the full physical problem
can be split up into an unperturbed part
that is readily diagonalisable and a perturbation
associated with a perturbation parameter
that is small compared to the energy scales of
. We aim to obtain a power series up to a maximal reachable order
as an approximation of a desired physical quantity
where the coefficients
are to be determined by series expansion. We want to use the information contained in the power series to infer properties of the approximated function
[
155]. The cost of determining the coefficients is associated with an exponential growth in complexity with increasing order [
155]. Hence, calculations are performed with the help of a computer programme. Obviously, the computer cannot deal with an infinitely large lattice. Instead, we must look at finite cut-outs consisting of a finite set of lattice sites that are connected by bonds (or links) symbolising the interactions of the Hamiltonian on the lattice. We term these cut-outs
clusters. If two clusters
A and
B do not share a common site or conterminously do not have a link that connects any site of A and B with each other (
), then the cluster
is called a
disconnected cluster. Otherwise, if no such partition into disconnected clusters
A and
B exists (
), the cluster
C is called
connected. We can define quantum operators
(e.g., a physical observable) on these clusters just as on the infinite lattice.
There are essentially two ways of performing high-order series expansions. The first one is the naive approach of taking a single finite cluster
[
153,
154,
156,
157,
158] and designing it such that the contribution of
coincides with the contributions on the infinite lattice
up to the considered order in the perturbation parameter. The cluster needs to be chosen large enough such that the perturbative calculations up to the considered order are not affected by the boundaries of the cluster. Another way of performing calculations is to construct the operator contribution on a cluster – coinciding with the infinite lattice contributions up to a given order – by decomposing it into all possible contributions on smaller clusters [
155,
159,
160,
161,
162,
163,
164,
165,
166,
167,
168,
169,
170]. Now the contributions on many but smaller clusters must be determined and added up to obtain the contribution on the infinite lattice
In contrast to the previous approach we willingly accept boundary effects for the many subclusters. Such a cluster decomposition is known to be computationally more efficient because it suffices to calculate the contributions on the Hilbert space of the smaller clusters reducing the overhead of non-contributing processes and it also suffices to do the calculations only on a few clusters as many give identical contributions due to symmetries. This can significantly reduce the overhead.
However, there are subtleties about the validity of performing calculations on finite clusters, e.g., when setting up a linked-cluster expansion (linked-cluster means only connected clusters contribute) the operator
must satisfy a certain property, namely cluster additivity. The quantity
is called
cluster additive if and only if the contribution on disconnected clusters
solely comes from the contributions on its individual connected clusters
A and
B. This means we can simply add the contributions of
A and
B from the smaller connected clusters to obtain the one for the disconnected cluster, i.e.,
as illustrated schematically in
Figure 3.
We can also understand cluster additivity in the language of perturbation theory where acting on bonds in every order forms a cluster of perturbatively active bonds. If such a cluster is connected, we call these processes linked. So cluster additivity simultaneously means that only linked processes will contribute. Cluster additivity is at the heart of the linked-cluster theorem, which states that only linked processes will contribute to the overall contribution in the thermodynamic limit. To set up a linked-cluster expansion we want to exploit cluster additivity and the linked-cluster theorem so that we can "simply" add up the contributions from individual connected clusters to obtain the desired power series in the thermodynamic limit.
An example of a cluster additive quantity is the ground-state energy
. Imagine we want to calculate the ground-state energy of a non-degenerate ground-state subspace, then cluster additivity is naturally fulfilled
and we can calculate the ground-state energy on
from its individual parts
A and
B. However, cluster additivity is not satisfied in general. We can construct a counterexample by considering the first excited state with energy
. For example, the first excitation above the ferromagnetic ground state of the transverse-field Ising model in the low-field limit is a single spin flip dressed with quantum fluctuations induced by the transverse field [
113]. We usually refer to such excitation as
quasiparticles qp. Here, we cannot add the contributions on clusters
A and
B to obtain the excitation energy on cluster
How to set up a linked-cluster expansion for intensive properties is not obvious and it seemed out of reach after the introduction of linked-cluster expansions in the 1980s [
159,
160,
161]. Only several years later, it was noticed by Gelfand [
162] that additivity can be restored for excited states when properly subtracting the ground-state energy. This approach was later generalised to multiparticle excitations [
154,
156,
164,
165] and observables [
154,
163].
In the following, we first introduce a perturbation theory method that maps the original problem in Eq. (
48) to an effective one. We will show that the derived effective Hamiltonian and observables satisfy cluster additivity. In the subsequent section we make use of the property and show how we can set up a linked-cluster expansion for energies of excited states and observables by properly subtracting contributions from lower energy states.
4.2. Perturbation Method—Perturbative Continuous Unitary Transformations
The first step towards setting up a linked-cluster expansion is to find a perturbation method that satisfies cluster additivity, which is generically not given and a non-trivial task [
171]. Here, we use perturbative continuous unitary transformations (pCUT) [
153,
154] that transform the original Hamiltonian perturbatively order by order to a quasiparticle-conserving Hamiltonian, reducing the original many-body problem to an effective few-body problem. We start discussing how to solve the flow equation to obtain the pCUT method and show afterwards how the Hamiltonian decomposes into additive parts that can be used for a linked-cluster expansion.
We strive to solve the usual problem of perturbation theory of Eq. (
48). The unperturbed part
can be easily diagonalised exactly with a spectrum that has to be equidistant and bounded from below. Additionally, the perturbation
must be a sum of operators
containing all processes changing the energy by
n quanta and – if properly rescaled – corresponds to the same number of quasiparticles
n. The goal of the pCUT method is to find an optimal basis in which the many-body problem of the original Hamiltonian reduces to an effective few-body problem. For that, we introduce a unitary transformation depending on a continuous flow parameter
ℓ and define
In the limiting case
we require
to recover the original Hamiltonian and for
we require
so that the unitary transformation maps the original to the desired effective Hamiltonian. We can rewrite the unitary transformation as
where
is the anti-hermitian generator generating the unitary transformation and
the ordering operator for the flow parameter. Taking the derivatives of Eq. (
55) and Eq. (
56) in
ℓ, we eventually arrive at the flow equation
Flow equations have been studied for quite some time in mathematics and physics with a variety of applications [
172,
173,
174,
175,
176,
177,
178,
179,
180]. It was Knetter and Uhrig [
153] who proposed a perturbative ansatz for the generator of continuous unitary transformations along the lines of Mielke [
180], introducing the quasiparticle generator (also known as the "MKU generator") for the pCUT method
where the indices
refer to blocks of the Hamiltonian labelling the quasiparticle number. Diagonal blocks
contain all processes conserving the number of quasiparticles
i, while offdiagonal blocks
contain all processes changing the quasiparticle number from
i to
j. The reasoning behind the ansatz can be explained by looking at
, where processes
in
are assigned the opposite sign of the inverse processes
and therefore the idea is to "rotate away" offdiagonal blocks by the unitary transformation during the flow of
ℓ, while processes that do not change the quasiparticle number are not transformed away due to
but get renormalised during the flow. Consequently, in the limit
we obtain an effective Hamiltonian
that is block diagonal in
n.
This idea is depicted in
Figure 4. Next, we make a perturbative ansatz for the Hamiltonian during the flow
with the notation
and
being undetermined real functions. We introduce
distinct expansion parameters instead of just a single
to keep the notation as general as possible because below in
Section 4.6.3 about white-graphs we will need multiple expansion parameters to encode additional information. Inserting Eq. (
59) and Eq. (
58) into the flow equation (
57), we can solve the equation perturbatively order by order as we get a recursive set of differential equations for
.
To recover the original Hamiltonian
, we have to demand the correct initial conditions
for
and
for
. We can solve the differential equations (c.f. Ref. [
153]) exactly for
, yielding
with
being exact rational coefficients and the restriction
making the products
quasiparticle conserving [
153]. Hence, the commutator of the effective Hamiltonian with the unperturbed diagonal part of the original Hamiltonian vanishes (
). Note, so far the effective Hamiltonian (
64) is model independent. It only depends on the overall structure of Eq. (
54). The generic form of the Hamiltonian comes at the cost of an additional normal ordering usually by applying the Hamiltonian to a cluster. Of course, it could also be done explicitly by using the hard-core bosonic commutation relations but the former approach can be handled much easier by a computer programme. Yet, we achieved our goal of obtaining a block-diagonal Hamiltonian
where
is the effective irreducible Hamiltonian of
n quasiparticle processes (see also
Section 4.3). Let us emphasise again, that this block diagonal structure allows us to solve the
n quasiparticle blocks individually which significantly reduces the complexity of the original many-body problem to an effective one that is block-diagonal in the quasiparticle number
n.
If we want to calculate an effective observable we can make an ansatz along the same lines [
154]. We insert the perturbative ansatz
with undetermined functions
. The operator product is defined as
Inserting exactly the same generator (
58) and the ansatz for the observable in Eq. (
66) instead of the Hamiltonian into the flow equation (
57), we arrive at
with
by solving the resulting set of differential equations for
[
154]. Note that the last sum does not contain a restriction
and therefore – in contrast to the effective Hamiltonian – effective observables are not (necessarily) quasiparticle conserving.
We have just derived the effective form of the Hamiltonian and observables in the pCUT method that have a very generic form depending only on the structure of the perturbation of Eq. (
54). As already stated, the model dependence of our approach comes into play when performing a linked-cluster expansion by applying the effective Hamiltonian or observable to finite clusters. But how do we know if the effective quantities are cluster additive?
We follow the argumentation of Refs. [
181,
182] by looking at the original Hamiltonian (
48) that trivially satisfies cluster additivity as long as all bonds represent a non-vanishing term in
between sites ("bond equals interaction"). Thus, the Hamiltonian on a disconnected cluster
is cluster additive because
and
are non-interacting. Here, we denote the restriction of the Hamiltonian
to a cluster
C as
. We can further insert this property into the flow equation
Here, we used the property of
that it commutes on disconnected clusters and the fact that the Hamiltonian is continuously transformed during the flow starting from
. Therefore, the derivative and the commutator can be split up acting on each cluster individually and preserving cluster additivity during the flow. Consequently, the effective Hamiltonian
in the limit
is cluster additive as well. The same proof holds for effective observables. Another more physical argument is that the effective pCUT Hamiltonian can be written as a sum of nested commutators of
T-operators [
152,
167]. For instance, considering the perturbation
, the effective Hamiltonian looks like
Splitting up
into local operators acting on bonds
l, the nested commutators vanish for processes that are not linked. Hence the linked-cluster theorem is fulfilled and the effective Hamiltonian is cluster additive. To emphasise the linked-cluster property the generic effective Hamiltonian is often written as
where the sum over
C runs over all possible connected clusters with maximal
k bonds (
) [
152]. The notation
will be clarified in
Section 4.6.2 where graphs are formally introduced. In this context it is simply the set of bonds of a connected cluster
C and
the number of bonds in this set. The condition
arising from the linked-cluster theorem ensures that the cluster consisting of active links and sites during a process must match with the bonds and sites of the connected cluster
C. For observables the generalised condition
holds, where the index
x can either refer to a site (local observable) or a link (non-local observable) and we have
Although we showed that the effective Hamiltonian and observables are cluster additive and therefore fulfil the linked-cluster theorem, to set up a linked-cluster expansion, there are important subtleties remaining when we restrict the effective Hamiltonian and observables to the quasiparticle basis that we need to address before we can discuss how to do the calculations in practice.
4.3. Unraveling Cluster Additivity
In this subsection we need to clarify how we can use the cluster additive property of the effective pCUT Hamiltonian and observables to set up a linked-cluster expansion not only for the ground-state energy but also for energies of excited states. Many aspects of this section are based on the original work of Ref. [
154], in which a general formalism was developed how to derive suitable quantities for the calculation of multiparticle excitations and observables. We further develop this formalism by inferring the concept of cluster additivity to the quasiparticle basis, introducing the notion of particle additivity. The term "additivity" in this context was recently introduced by Ref. [
171].
We start by recalling that the effective Hamiltonian is block diagonal and we can write the Hamiltonian operator as a sum of irreducible operators of
n quasiparticle processes
We can express the
n quasiparticle processes in second quantisation in terms of local (hard-core) bosonic operators
creating and
annihilating a quasiparticle at site
i. When considering quantum magnets like we do in this review a hard-core repulsion comes into play allowing only a single quasiparticle at a given site [
154]. For instance, in the ferromagnetic ground-state of the 2D Ising model an elementary excitation is given by a single spin flip that can be interpreted as a quasiparticle excitation [
113]. Obviously, at most one excitation on the same site is allowed. Different particle flavours
can also be accounted for by incorporating an additional index
of the operator
. To keep the notation simple, we will drop this additional index in the following. The irreducible operators in second quantisation and normal-ordered form then read
Written in normal order the meaning of these processes is directly clear when acting on states in the quasiparticle basis. The prefactors
are to be determined by applying the effective pCUT Hamiltonian (
64) to a appropriately designed cluster
C. Let us consider the quasiparticle basis on a connected cluster
C
where the number
n specifies the number of particle excitations and the indices
denote the positions of the
n (local) excitations. The effective pCUT Hamiltonian is quasiparticle conserving, so let us restrict it to
N particle states, like when evaluating its matrix elements in this basis. If we evaluate by acting on a state with fewer particles than particles involved in the process then the irreducible operator annihilates more particles than there are and the contribution is zero, that is
for
. When we determine the action of
for
this allows us to determine all prefactors
defining the action of
on the entire unrestricted Hilbert space. Second quantisation presents a natural generalisation of the Hamiltonian restricted to a finite number of particles to an arbitrary number of particles [
154]. We can construct
for
from
since the latter completely defines the action
on the entire Hilbert space.
Although everything seems fine so far and we cannot wait to set up a linked-cluster expansions, let us tell you that it is not. We finished the motivation in
Section 4.1 with the statement that we cannot simply add up contributions for energies of excited states (c.f Eq. (
53)). The reason is that although we showed that
is cluster additive, the irreducible operators restricted to the
N particle basis
are in fact not. To grasp a better understanding about the abstract concept of cluster additivity and why setting up a linked cluster expansion for higher particle channels usually fails, let us consider the following basis on a disconnected cluster
where
represents all possible
n-particle states living on a cluster
C. While there is only one way to decompose the zero-particle states
on the disconnected cluster
, one-particle states
decompose into two sets of states with a particle on cluster
A (
) and a particle on
B (
). For two-particle excitations
there are three possibilities to distribute the particles. In general,
N-particle states have the form
with
and there are
possibilities to decompose the states.
When restricting Eq. (
51) to these
N-quasiparticle states, a cluster-additive Hamiltonian must decompose as
where we introduce the notation
restricting the Hamiltonian to all
n-quasiparticle states on a cluster
C. The direct sum in Eq. (
79) is introduced to emphasise that for a cluster additive Hamiltonian there must not be any particle processes between the two disconnected clusters
A and
B. The Hilbert space on the disconnected cluster
can be seen as the natural extension of the Hilbert spaces on cluster
A and
B and we can define the operators on the clusters
A and
B in terms of the Hilbert space on
as a tensor product
The issue with Eq. (
79) is that when we restrict the particle basis to
N on the disconnected cluster
there are contributions from lower particles channels coming from the
possibilities to distribute the
N particles on the two clusters. For example, if we look at the one-particle space
we see that besides the one-particle contributions we get additional zero-particle contributions. The left part of the direct sum stems from acting on
and the right side from acting on
. There are always two possibilities to distribute the one-particle excitation on the two clusters where the other cluster is always unoccupied which gives an additional zero-particle contribution. This fact is schematically illustrated in
Figure 5a.
This is not the desired behaviour for a linked-cluster expansion. We would like the general notion of Eq. (
51) to directly translate to the particle-restricted basis which is illustrated in
Figure 5b. In words, our goal is to find a notion of (cluster) additivity in the restricted particle basis such that we can simply add up the
N-particle contributions of individual clusters without caring about lower particle channels. We define
particle additivity as
where we demand the other contributions from Eq. (
79) to vanish. The crucial thing to notice is that irreducible operators (
76) in the particle basis
are in fact particle additive [
154,
171]. In the following we will show that this is indeed the case. First, we remember that for the ground-state energy we can trivially add up the contributions. Starting from the definition for cluster additivity (
79), we have
Second, from restricting the decomposition of Eq. (
75) to the
particle channel, we can express the irreducible one-particle operator as
We recall that by calculating
, we can automatically derive
on the entire Hilbert space which subsequently defines us
. Therefore, by inserting the definition for cluster additivity (
79), we obtain
where we used the definition of Eq. (
82) in the last line. Hence, we have proven Eq. (
83) for
. The above prove can be readily extended to
.
We achieved our goal of finding a notion of cluster additivity in the particle basis which we termed particle additivity. We can determine the desired particle additive quantities by using the subtraction scheme
that comes from Eq. (
75) by restricting it to a
N particle basis. This is an inductive scheme starting from
calculating the irreducible additive quantity
. This result can be used to calculate the subsequent irreducible additive quantity
for
. Then for
we use the results from
and
to calculate
and so on. Again, it is important that
completely defines the operator
and therefore any
.
When considering effective observables, the particle number is no longer conserved and more types of processes are allowed. We need to generalise Eq. (
75) for effective observables by introducing an additional sum over
d that is the change in the quasiparticle number. An effective observable thus decomposes as
where
are irreducible contributions [
154]. When writing them in second quantisation, we have
We can directly see that the
d quasiparticles are created because there are
d additional creation operators. When
d is negative
d quasiparticles are annihilated. We can infer a notion for cluster additivity and particle additivity along the same lines
To determine the additive parts we can use an analog subtraction scheme as described in Eq. (
87), that can be denoted as
If we want to calculate
then we have to inductively apply Eq. (
92).
There are several things to be noted at this point. First, not all perturbation theory methods satisfy cluster additivity (
79) and in this case we cannot write operators as a direct sum anymore. There will be quasiparticle processes between one and the other cluster changing the number of particles on each cluster [
155,
162]. This is sketched in
Figure 6 by the presence of an additional term.
When falsely performing a linked-cluster expansion, it can be noticed immediately that the approach breaks down. A symptom of non-cluster additivity is the presence of contributions of lower orders than expected from the number of edges of the graph. When calculating reduced contributions in a linked-cluster expansion, we subtract only contributions of connected subgraphs which leaves non-zero contributions of disconnected clusters when the perturbation theory method is not cluster additive. However, there are notable exceptions when a linked-cluster expansion for energies of excited states is still correct even though the perturbation theory method is not cluster additive [
163,
164,
165,
183,
184,
185,
186]. This is only possible when the considered excitation does not couple with a lower particle channel, i.e., lower-lying states are described by a distinct set of quantum numbers lying in another symmetry sector [
155,
165,
187]. For instance, consider the elementary quasiparticle excitation in a high-field expansion for the TFIM, then the structure of the perturbation is
and therefore the first excited state does not couple directly with the ground state (there is no
which is due to symmetry). If one wants to draw a comparison, we can think of this to be similar to the case when calculating the excitations in Density Matrix Renormalisation Group. It is no problem to target an excited state if it is in a different symmetry sector than the ground state but if it is in the same symmetry sector described by same set of quantum numbers then the Hamiltonian needs to be modified to project out the ground state [
188].
Recently, a minimal transformation to an effective Hamiltonian was discovered that preserves cluster additivity. This method called "projective cluster additive transformation" [
171] can be used analogously and is even more efficient for the calculation of high-order perturbative series. In this review, however, we stick to the well established pCUT method.
4.4. Calculating Cluster Contributions
At this point, we may ask how to evaluate physical quantities on finite clusters in practice. To evaluate these quantities we must evaluate them in the quasiparticle basis. In general, when setting up a cluster expansion may it be non-linked or linked, it is important to subtract contributions from all possible subclusters to prevent over counting. Mathematically for a quantity
(
or
), this can be written as
where the sum runs over all real subclusters
in cluster
C and we call the resulting quantity
reduced. Starting from the smallest possible cluster (e.g., a single bond between two sites) this formula can be inductively applied to determine the reduced quantity on increasingly big clusters. An essential observation to make is that reduced operators vanish on disconnected clusters if the operator
is additive by construction since we subtract all contributions from individual subclusters. As the linked-cluster theorem applies, we can set up a cluster expansion
of connected clusters but we need to consider reduced quantities
to prevent over counting. For a light notation we will drop the bar in the sections below as we will only consider reduced contributions on graphs anyway.
Now we are ready to look at the problem from a more practical point of view. From the previous subsection we know how cluster additivity translates into the particle basis and how to construct particle additive parts, namely the irreducible quasiparticle contributions. We decompose the effective Hamiltonian into its irreducible contributions by explicitly calculating:
Again, consider the effective Hamiltonian in second quantisation made up of hard-core bosonic operators annihilating (creating) quasiparticles and the quasiparticle counting operator occurring in the unperturbed Hamiltonian . We also consider a connected cluster C and we denote n quasiparticle states on this cluster as with the quasiparticles on the sites to . Note, for multiple quasiparticle flavours or multiple sites within a lattice unit cell this notation can be generalised to by introducing additional indices . To lighten the notation in the following we stick to the former case. Let us consider the three lowest particle channels:
- n = 0
We can directly calculate the ground-state energy
on a cluster
C as it is already additive
as can be seen from Eq. (
95).
-
n = 1
To calculate the irreducible amplitudes
associated with the hopping process
in
, we need to subtract the zero-particle channel as can be seen from Eq. (96). However, we only need to subtract the ground-state energy if the hopping process is local,
since the ground-state energy only contributes for diagonal processes. Thus, we calculate
In the two-particle case we have to distinguish between three processes: pair hoppings (
with four distinct indices), correlated hoppings (
) and density-density interactions (
). The free quasiparticle hopping is already irreducible and nothing has to be done, but for the correlated hopping contribution we have to subtract the free one-particle hopping. In case of the two-particle density-density interactions we need to subtract the local one-particle hoppings as well as the ground-state energy as this process is diagonal (c. f. Eq. (97)). Therefore we calculate
An analog procedure can be applied for effective observables. Here, we need to determine the irreducible contributions for a fixed
d. The subtraction scheme is given by
For
we recover exactly the same subtraction procedure as before. It is straightforward to generalise this procedure for
. Let us specifically consider the example we will encounter in the next section, when calculating correlations for the spectral weight. The effective observable is applied to the unperturbed ground state
. Hence, there are only contributions out of the ground state (
) and the effective observables decomposes into
Since only processes contribute, nothing needs to be subtracted and the effective observable is irreducible and already particle additive.
Although for these types of calculations the lowest orders can be analytically determined by hand, the calculations usually become cumbersome quickly and must be evaluated using a computer programme to push to higher perturbative orders (in most cases a maximal order ranging from 8 to 20 is achievable). Such a programme reads the information of the cluster (site and bond information), the bra and ket states, the coefficient lists
or
from Eqs. (
64) and (
68) up to the desired order, the structure of the Hamiltonian
, the local
-operators from the perturbation
and if necessary the observable. The input states as well as the
-operators should be efficiently implemented in the eigenbasis of
bitwise encoding the information as known for instance from exact diagonalisation. If possible, the calculation should be done with rational coefficients for the exact representation of the perturbative series up to a desired order. The routine of the programme is then to iterate through the operator sequences from the coefficient list
or
and to consecutively apply the
-operators by systematically iterating over all bonds of the cluster and calculating the action of the operator, saving the intermediate states for the action of the next operator in the sequence. As intermediate states are superpositions of basis states, they are saved in associative containers (maps in C++ or dictionaries in python)
where
is the bit representation of a basis state. The key of the associative container is the basis state
and the associated value the prefactor
. The bitwise representation of the basis states
as well as the
-operators allow for a fast access and modification during the calculation.
So far, we have introduced the pCUT method for calculating the perturbative contributions on clusters. We demonstrated that the resulting effective quasiparticle-conserving Hamiltonian is cluster additive and showed how to extract particle-additive irreducible contributions. We introduced a subtraction scheme for the effective Hamiltonian and observables that can be easily applied to calculate additive quantities to set up a linked-cluster expansion. Finally, we briefly discussed an efficient implementation of the pCUT method.
4.5. Energy Spectrum and Observables
Having established the basic theoretical framework for the pCUT method we want to give a short overview over the physical quantities that are most frequently calculated with this approach. To this end we assume that we consider a single suitably designed cluster for the calculations instead of setting up a LCE as full graph decomposition. We will see how we can calculate the desired quantities without thinking about the abstract concepts necessary for a linked-cluster expansions and with the insights from this section it will be easier to recognise what we are aiming at. Here, we first consider the energy spectrum of the Hamiltonian. We derive both the control-parameter susceptibility as the second derivative of the ground-state energy as well as the elementary excitation gap from the effective one-quasiparticle Hamiltonian. Second, we consider observables. In the pCUT approach often spectral weights are calculated that are derived from the dynamic structure factor which is of great importance for inelastic neutron scattering experiments.
4.5.1. Ground-State Energy and Elementary Excitation Gap
Following the above described recursive scheme, we start with the zero-quasiparticle channel assuming a non-degenerate unperturbed ground state which is the situation in all applications discussed in this review. The ground-state energy can be directly calculated from the cluster as in Eq. (
100). We consider a suitable designed cluster that is large enough to accommodate all fluctuations of a given maximal order
and has periodic boundary conditions to correctly account for translational invariance. We calculate
and obtain a high-order perturbative series
in the expansion parameter
, which is valid in the thermodynamic limit up to the given maximal order
. We can extract the ground-state energy per site by dividing
by the number of sites of the cluster
. By taking the second derivative, we obtain the control-parameter susceptibility
We are usually interested in the quantum-critical properties of the model and therefore analyse the behaviour about the quantum-critical point
. The control-parameter susceptibility shows the diverging power-law behaviour
with the associated critical exponent
as we know from
Table 1.
Turning to the one-quasiparticle channel, we calculate the hopping amplitudes following Eq. (
101). Here, we use open boundary conditions, again with a cluster large enough to accommodate all fluctuations contributing to the hopping process. Note, in our notation we denote
for hopping amplitudes on a graph or cluster level and the character
a for processes in the thermodynamic limit. As for the ground-state energy we can directly infer the contribution in the thermodynamic limit if contributing fluctuations do not feel finite-size effects and thus we use
in the following. Further, we can generalise our notation to multiple particle types or larger unit cells containing multiple sites as mentioned earlier by introducing additional indices
,
. We denote a hopping from unit cell
to
with
and within the unit cells from
to
. We calculate
by fixing
, due to translational symmetry. The effective one-quasiparticle Hamiltonian in second quantisation then is
Applying the Fourier transform for a discrete lattice
we can diagonalise the resulting Hamiltonian in momentum space
introducing the operators
that diagonalise the matrix
. The eigenenergies
are the associated bands of the one-quasiparticle dispersion. In case of a trivial unit cell or a single particle flavour, the dispersion matrix becomes a scalar and we can directly express the single-banded dispersion as
where the sum over
is restricted to symmetry representatives and we assumed real hopping amplitudes
. We determine the hopping amplitudes
as a perturbative series by performing the calculations on the properly designed cluster. Note, even for a single cluster we need to perform a subtraction scheme because in order to get the irreducible contribution
explicitly we need to subtract the ground-state energy from
. When we determine the elementary excitation gap at the minimum of the lowest band
we can as well extract the gap directly as such a series. The gap closing shows a power-law behaviour
about the critical point with the critical exponent
.
4.5.2. Spectral Properties
Neutron scattering is a powerful method resolving spatial and dynamical structures in condensed matter physics since thermal neutrons have a de Broglie wavelength of similar length scale as interatomic distances and their energy is of the same order of magnitude as typical excitations [
189]. By measuring the change in momentum and kinetic energy in inelastic neutron scattering experiments determining the dynamic response, not only static properties like magnetic order can be resolved but also dynamic properties like spin-spin correlations [
190]. The dynamic response
can be determined as it is proportional to the cross section
of inelastic neutron scattering [
189,
190]. We follow the derivations in Refs. [
155,
191,
192] and start with the definition of the
dynamic response
which is the space and time Fourier transform of the spin correlation function
with
and
referring to the thermal expectation value. In the limit of vanishing temperature
, the expectation value simplifies to
with
being the ground state. Then we call
the
dynamic structure factor. We introduce a complete set of energy eigenstates
where
denotes a set of quantum numbers, for instance
, where
n is the number of quasiparticles and
the lattice momentum. Writing the dynamic structure factor in terms of these energy eigenstates, this yields the dynamic structure factor in the spectral form as a sum
where
is called
exclusive structure factor or
spectral weight associated with the quantum numbers
. We insert
into the correlation function and switch to the Heisenberg picture, where
. The spectral weight then reads
which after integrating over time t yields
If we consider the case
or some observable with
, that could be a linear combination of spin operators, then the above expression further simplifies to
By the above equations we can identify the spectral weights
2 as
These quantities are usually visualised as a heat map where the dispersion is plotted against the momentum and the value , that is the intensity of the scattering signal associated with , is colour-coded.
In the pCUT approach we want to reformulate the observable in terms of an effective one
. Here, we want to restrict ourselves to the one-quasiparticle spectral weight. In Ref. [
191] you can also find a formulation for the two-quasiparticle case. For the 1qp spectral weight,
are the quantum numbers defining one-quasiparticle states and we denote
in the following. Since the pCUT method is a perturbative approach we want to reformulate the problem in the language of
states
where we introduce the momentum states
with additional index
denoting a quantum number like a flavour of the excitation or denoting a site in a unit cell. The momentum states are defined via Fourier transform
Inserting these identities, we obtain
where we defined
. For a problem with translational invariance, we can fix the site
of the observable and introduce
, which yields
where we used
. Note, the form of the effective observable is
, exactly as denoted above in Eq. (
108). As long as the processes
and
are equivalent by means of the model symmetries we can use
and further simplify
to
So, in the pCUT method, we determine
as a perturbative series with which we can determine
. Note, this observable is already an irreducible contribution according to Eq. (
92) and no subtraction scheme needs to be performed. Because we consider an observable
the initial particle number is
and nothing needs to be subtracted. To extract the quantum-critical behaviour, we need to evaluate the expression at the critical momentum
as
yielding a perturbative series for this quantity as well. It shows a diverging critical behaviour that goes as
with the associated critical exponent
. After determining the three quantities
,
, and
described above, we can use the extrapolation techniques described in
Section 4.8 to extract estimates for the critical point
and the associated critical exponents. However, we continue with the description of a linked-cluster expansion as a full graph decomposition for long-range interacting systems. The next step on the way is to formally introduce graphs, their generation and the associated concept of white graphs.
4.6. White Graph Decomposition
In this section we first give a brief introduction to graph theory which is the basis to understand how linked-cluster expansions as a graph decomposition work. Then, we discuss how to generate all topologically distinct graphs up to a certain number of edges corresponding to the maximal number of active edges at a given perturbative order. We conclude this section by explaining the concept of white graphs where edge colours are ignored as a topological attribute in the classification of graphs. White graphs are essential for tackling long-range interactions and therefore every graph decomposition in this review is in fact a white graph decomposition.
4.6.1. Graph Theory
So far, we already defined clusters as a cut-out of the infinite lattice with a finite set of sites and a set of bonds connecting those sites. More generally, only considering the topology of clusters without restricting them to the geometry of lattices, we can define a
graph as a tuple
consisting of a (finite) set of vertices
and a (finite) set of edges
. An edge
consists of a pair of vertices
and these vertices
are called
adjacent. The
degree of a vertex is the number of edges connecting it to other vertices of the graph. In the following, we only consider
undirected,
simple, and
connected graphs which means there are neither directed edges, multiple edges between two vertices, nor loops (no edge that is joining a vertex to itself) and there always exists a path of edges connecting any two vertices of a graph [
155,
193,
194]. As an example we depict a graph
with
and
in
Figure 7.
A
subgraph of a graph
– we write
– is defined as a subset of
and
[
193,
194,
195]. We call
a proper subgraph if
and
are proper subsets of
and
.
To set up a full graph decomposition it is essential to define how to distinguish different graphs. If there exists an isomorphism between two graphs we call them
topologically equivalent otherwise they are
topologically distinct. A
graph isomorphism is a bijective map
between the vertex sets of two graphs, such that
So, an isomorphism preserves adjacency and non-adjacency, i.e.,
and
are adjacent if and only if the vertices
and
are adjacent for any vertices
[
193,
194,
196]. A special case are
graph automorphisms , which are maps of a graph on itself, so we have
In other words a graph automorphism is a permutation of the vertex set preserving adjacency [
194,
196]. The number of graph automorphisms
of a given graph gives the number of its symmetries and we call it the
symmetry number of a graph [
155,
195]. Examples of a graph isomorphism and automorphism are depicted in
Figure 7. Further, if
the mapping (
134) is injective instead of a bijective such that
and we call it
subgraph isomorphism or
monomorphism [
195,
197,
198].
An example of a monomorphism is depicted in
Figure 8a. Monomorphisms will later become of utmost importance as they give the number of embeddings of the subgraph onto a graph which is the infinite lattice.
To account for hopping processes during the embedding of graphs we need to assign additional attributes to graphs like colouring their vertices [
195]. Then, a coloured graph
is a tuple
, where
is a set of pairs
with
and
a the colour attribute. In
Figure 8b an example is depicted with
. We can extend the above definitions for isomorphisms and automorphisms and say that they must preserve the vertex colour, i.e., only vertices of the same colour can be mapped onto each other. Of course, this reduces the cardinality
and therefore the symmetry number
associated to the coloured graph. We can later exploit the colour information for matching hopping vertices of the graph with the actual hopping of the quasiparticle on the lattice. Note, also colouring edges of graphs with attributes
is useful to distinguish different types of interactions of a Hamiltonian on the graph level and very similar properties as for coloured vertices hold. As stated above two graphs are topologically distinct if there does not exist a graph isomorphism between the two graphs. Thus, for coloured graphs the colour attribute serves as another topological attribute. The importance of this will become apparent later in
Section 4.6.3. For a more elaborate overview of graph theory in the context of linked-cluster expansions, we recommend Ref. [
195].
4.6.2. Graph Generation
For the graph generation of simple connected graphs we need to define an ordering between all graph isomorphs such that it is possible to pick out a unique representative of all graph isomorphs that is called
canonical representative [
199]. One challenge lays in efficiently generating graphs as the number of connected graphs grows exponentially with the number of edges and the idea behind every algorithm must be to restrict the number of possible graph isomorph candidates when generating new edge sets by permutation. Most popular, there is McKay’s algorithm [
199,
200] that exploits the degree information of the vertices and sets up a search tree for vertices with the same degree and uses the ordering to check if the canonical representative of the current graph already exists. We recommend using open-source libraries for the graph generation and calculation of additional symmetries. For instance, there is "nauty" [
201] and "networkX" [
202] or the "Boost Graph Library" [
203].
There are various conventions for how to write the graph information to a file. One could simply save the site and edge set as lists or save its
adjacency matrix, where the rows and columns refer to the sites and a non-zero entry in the corresponding matrix element marks the existence of an edge between the sites [
155,
193,
194]. Here, we suggest to simply use a
bond list. Each entry in the list contains the edge information of the graph, denoted as
, where
with
adjacent to
e and
is just a number associated to the edge
e. The number
can be interpreted as a specific expansion parameter corresponding to this bond. In the simplest case there is just a single expansion parameter and number
for all edges. Assigning multiple expansion parameters becomes especially important in the next
Section 4.6.3. Usually, the symmetry number
of a graph is calculated on the fly when generating a set of graphs and should be saved as well. In our workflow, we save the generated graphs into bond lists and create a list containing all connected graphs (e.g., the filenames) along with their symmetry number
. These types of lists suffice for the calculation of the ground-state energy. When calculating 1qp irreducible processes for the dispersion or for the spectral weight observable, we can think of these processes breaking the graph symmetry. Therefore, after the graph generation, we consider all possible processes on a graph and assign colour attributes to the start and end vertices. Due to the symmetry of the processes we assign the same colour for hoppings and distinct colours for processes of the spectral weight.
3 We calculate the symmetry number
associated with the coloured graph. In the end, we create a list, where each entry contains the graph, its symmetry number
, a representative process (start and end vertex) and the associated symmetry number of the coloured graph
, counting the number of processes that give the same contributions as the representative process due to symmetries.
After generating the graphs and the lists, we can employ perturbation theory calculations on the graph level viewing graphs as abstract clusters with vertices as lattice sites connected by bonds. A programme as described above reads in the graph and process information and repeatedly performs the pCUT calculations on every graph. The resulting graph contributions must be added up in a way such that the information of the lattice geometry is restored by weighting the graph contributions with embedding factors. That is, how many ways a graph can be fitted onto a lattice apart from translational invariance (c.f. Ref. [
155]). For conventional embedding of contributions from models with just nearest-neighbour interactions, the number of graphs can be reduced as the lattice geometry puts a restriction on the graphs. For example, graphs containing cycles of odd length cannot be embedded on the square lattice. In contrast, for long-range interactions no such restriction exists and therefore the lattice can be seen as a single fully connected graph where every site interacts with each other. Hence, we have to generate every simply connected graph up to a given number of edges as all graphs contribute. Before we turn to the embedding problem in more detail, we will first deal with the challenges long-range interactions pose and address the problem by using white graphs.
4.6.3. White-Graphs for Long-Range Interactions
In many cases in perturbation theory there may be more types of expansion parameters, e.g., due to the geometry of the lattice as can be seen in the
n-leg Heisenberg ladder [
152]
with an expansion parameter
associated with the rungs and another one
associated with the legs. The interaction
is given by XY-interactions and the series expansion is done about the Ising limit
. Setting up a graph decomposition for this model, the canonical approach would be to associate each expansion parameter with a distinct edge colour, blue for
and purple for
.
In
Figure 9a all graphs with two edges and two colours (one for
and one for
) are depicted on the left. It is necessary to incorporate the edge colour information as the graphs can only be embedded correctly on the infinite lattice when the edge colour matches the bond type of the lattice. Another common type of perturbation problem is where the perturbation splits into different types of interactions although associated with the same perturbation parameter. For instance, the problem can look like
which is essentially the form of the Hamiltonian when performing a high-field series expansion for Kitaev’s honeycomb model [
205] where each site has one
x-,
y-, and
z-type Ising interaction bond to one of its neighbours on the honeycomb lattice, respectively [
204]. Here, we associate the three different interaction types with three types of edge colours blue, green, and purple for
x-,
y-, and
z-bonds. See
Figure 9b (left) for an illustration of all three graphs with two edges and three colours (there are only three possibilities because colours must alternate due to the constraints posed by Kitaev’s honeycomb model). In both cases the edge colour is an additional topological attribute of the graph leading to exponentially more graphs with the number of colours which becomes relevant when pushing to high orders in perturbation.
In case of long-range interactions such an approach becomes unfeasible. A Hamiltonian with long-range interactions is of the form
where
is the distance between interacting sites,
d the dimension of the lattice, and
the long-range decay exponent and
. Applying the conventional approach from above we would associate each of the infinitely many perturbations
with its own edge colour. The only obvious way to resolve this problem would be to truncate the sum over
only considering very small distances. Instead, the use of
white graphs [
152] can be a solution to problems of this kind [
25,
29]. The idea is to change our view how to tackle these Hamiltonians. We ignore the edge colours of the graph – significantly reducing the number of graphs for a given order – and instead encode the colour information in the expansion parameters on the graph level in a more abstract way. This is done by associating each edge
e of the graph
with a different "abstract" perturbation parameter
such that we can track how often
-operators acted on each edge
e yielding a multivariable polynomial
with the sum over all monomial contributions. The individual monomial contributions consist of a prefactor
and its monomial dependency
.
comprises a product of expansion parameters
associated with edge
e and their respective integer powers
tracking how often each bond was active during the calculation. Let us emphasise, we simply wrote the white graph contribution for
,
, or
explicitly as a polynomial
. It is only later during the embedding of these abstract generalised contributions when the proper link colour is reintroduced by substituting the expansion parameters by the actual expansion parameters for each realisation on the actual lattice. This is the origin of the name "white graphs" because during the calculation on the graph the colour of the links is unspecified (but encoded in the multivariable polynomial) and only during the embedding the colour of the edges is reintroduced. In
Figure 9 we illustrate the white-graph concept in the middle. On the right, we depict how to recover the polynomial of the conventional graph contribution from the abstract white graph contribution for the models in Eq. (
137) and Eq. (
138) by correctly substituting the abstract expansion parameters. The main difference between these two models (
137) and Eq. (
138) is that for different interaction types multiple parameters are associated with each edge. To account for three types of interaction flavours we have to consider three expansion parameters
with
per edge
e.
For long-range interactions we can straight-forwardly apply the exact same white-graph scheme. We substitute the abstract expansion parameter with the correct algebraically decaying interaction strength depending on the distance between the interacting sites on the lattice [
25,
29]. For this, we have to use the substitution
with
and
. It is also possible to incorporate multiple interaction types but the substitution for the different interaction types must then be done before the embedding as it was done in Refs. [
31,
35,
206].
So far, we explained how to resolve the problem of infinitely many perturbation parameters by introducing white graphs. We managed to reduce the number of expansion parameters from infinity to the number of edges, i.e., the order of the perturbation. Yet, the polynomial in Eq. (
140) still grows exponentially with the number of expansion parameters. Following the description in Ref. [
152], we can further mitigate this issue by an efficient representation of white-graph contributions. The abstract concept is to track the relevant information in a quantity
as a generalised monomial with the property
and use
as an abstract parameter that encodes the tracked information. In Eq. (
140) the monomial quantity
is just
, tracking how often each edge was active during the calculation by associating each edge with its own abstract expansion parameter. We can generalise the expression of Eq. (
109) for states comprising additional information
Instead of a simple superposition of states
with associated prefactor
as in Eq. (
109), we have an additional superposition over all monomials
comprising the information of the all distinct processes encoded in
leading to the state
. We can make use of this factorisation property by using nested containers. The key of the outer container is the basis state
and the value contains an inner container with the monomial
as the key and the prefactor
as the value. Thus, the action of a
-operator on a state
can be calculated independently of the action on the monomial
. For a flat container we would have to calculate the action on the same state
multiple times. To further improve the efficiency of the calculation we can directly encode into the programme which edge of a graph was active in a bitwise representation. Using the information of the number of edges of a graph and tracking the applied perturbative order during the calculation we can neglect subcluster contributions on the fly and reduce the computational overhead even further. Therefore, we can directly calculate the reduced contribution on the graph without the need of explicitly subtracting subcluster contributions.
Wrapping things up, we explained how the use of white graphs can be applied to models with long-range interactions resolving the issue of infinitely many graphs or expansion parameters at any finite perturbative order. Instead of associating each perturbation as another topological attribute in the classification of graphs we associate an abstract expansion parameter with each edge of a graph and only during the embedding on the lattice we substitute these expansion parameters with the actual bond information of long-range interacting sites. An efficient representation of white graphs can further help to reduce the computational overhead.
4.7. Monte Carlo Embedding of White Graphs
We discussed ways to set up a linked-cluster expansion, either by designing a single appropriately large cluster hosting all relevant fluctuations up to a given perturbative order or by setting up a linked-cluster expansion as a full graph decomposition, where the latter is the more efficient way to perform high-order series expansions. Of course, a quantity
must be cluster additive in the first place. To decompose
on the lattice
into many smaller subcluster contributions, we can add up the contributions
where
are reduced contributions on a cluster
C to prevent overcounting. It is not necessary to calculate the contribution of every possible cluster since many clusters have the same contribution. It suffices to calculate the contribution of a representative cluster, only containing the relevant topological information, that is a graph defined by its vertex set
and its edge set
. We can see a graph
as representing an equivalence class
, whose elements are all clusters
C realising all possible embeddings of a graph on the lattice [
195]. Figuratively, all elements of the equivalence class are related by translational and rotational symmetry or are different geometric realisations on the lattice. We can now split the sum over all clusters into one sum over all possible graphs on the lattice and another one over all elements in the equivalence class
where it suffices to calculate
once for all
[
195]. We are left counting the number of elements
C in the equivalence class such that we can write
where the embedding factor
is simply a number counting the number of embeddings [
195]. The important point in our line of thought is that we can calculate the quantity
only once on the graph level and multiply the contribution with a weight that is the embedding factor
and sum up the resulting contribution for all graphs to obtain the desired quantity
in the thermodynamic limit. We are essentially left with determining the embedding factors after calculating the graph contributions
. Depending on the topology of the graph, the number of possible embeddings is different and therefore also the embedding factor. When calculating one quasiparticle processes for the 1qp dispersion or the 1qp spectral weight it is important that we account for the graph vertices’ colour attributes for the definition of the equivalence class. If conventional graph contributions are considered we directly obtain the correct contribution
. If there are multiple physical parameters or different interaction flavours the graph edges have to be colour matched with the bonds of the lattice. If white graphs are used we do not have to match any colours and can be ignored but the white graph contributions have to be evaluated appropriated substituting the abstract expansion parameters with the correct physical parameter for each embedding. Regarding continuative reading on the embedding, we want to point out the standard literature for linked-cluster expansions [
152,
155,
187] as well as Ref. [
195]. Nonetheless, we show in the following how the embedding procedure works for models with nearest-neighbour interactions specifically on the example of the ground-state energy and one-quasiparticle processes yielding the 1qp dispersion of the Hamiltonian. Likewise, we derive the 1qp spectral weight. With these findings we can turn to the embedding problem of long-range interacting models and eventually describe the Monte-Carlo algorithm as a solution to this problem.
4.7.1. Conventional Nearest-Neighbour Embedding
Ground-state energy: We start with the simplest case of calculating the ground-state energy on the infinite lattice
. We know that
is cluster additive and that the ground-state energy is an extensive quantity so we can directly calculate
In words this means we have to multiply the ground-state energy contributions on the graphs with the correct embedding factor and add up the resulting weighted contributions for every graph. For the definition of the embedding factor we follow the formalism introduced be Ref. [
195] and write
The embedding factor
is the number of subgraph isomorophism (monomorphism) of the graph
on the lattice
divided by number of graph automorphisms, i.e., the symmetry number of the graph. It is necessary to divide by the symmetry number
because the number of monomorphisms is in general not equal to the number of subgraphs because there may be multiple monomorphisms that map to the identical subgraph [
195]. To properly account for this, we have to divide by the number of automorphisms.
In
Figure 10 you see an example for the embedding problem. We can recognise the fact that there are multiple monomorphisms by the presence of arrows illustrating the ambiguity of mapping onto a specific subgraph embedding. Further, we always have to consider reduced graph contributions subtracting all subgraph contributions as
Here,
is the embedding factor for subgraphs
with respect to the considered graph
. Note again, with an efficient white graph implementation reduced contributions are calculated on the fly without the need of explicit subtraction. We state the subtraction scheme for completeness. Going back to the embedding problem for the infinite lattice
, the embedding factor will be extensive as the ground-state energy is extensive as well. Usually, we calculate the ground-state energy per site
, where
N is the number of sites and fix an arbitrary edge of the graph on the lattice. We then have
where the normalised embedding factor is
denoted with a lower case
w, where
q is the coordination number of the lattice (the number of neighbours of any site) [
155,
195]. One can think about fixing sites as colouring two adjacent vertices on the graph and two adjacent sites on the lattice and considering the monomorphism with respect to the additional colour attributes (for instance
) of the coloured graph
and lattice
. In
Figure 10 we depict the number of embeddings for a given graph on the square lattice and show how the embedding factor is calculated for this example. Note, it would be equally valid to just fix a single site. Then one would not have to account for the coordination number
q and the number of monomorphism would be larger by a factor of
q making it computationally more expensive. In the end we obtain the ground-state energy (per site) as a high-order perturbative series in an expansion parameter
up to a computationally feasible maximal order
as in Eq. (
110).
1qp dispersion: We now turn to the one-quasiparticle channel and calculate the hopping amplitudes. For a quasiparticle hopping
on the lattice with additional hopping within the unit cell from
to
or a quasiparticle changing its flavour from
to
, we need to calculate
where we fix the initial vertex
and end vertex
of the quasiparticle process on the graph to the initial site
and end site
in real space. Due to translational symmetry we can fix site
. This can be formally achieved by assigning colours to the initial and end sites on the graph and on the lattice. The second sum goes over all representative hopping processes
as described above in
Section 4.6.2. The embedding factor here is similar to Eq. (
150) and reads
where the colour attribute comes from fixing the sites to the hopping. We account for the reduced symmetry of coloured graphs stemming from the representative hopping by multiplying with the symmetry number of the coloured graph
.
The embedding factor is then calculated with respect to the coloured graph
and lattice
. Again, we need to be careful as we have to consider reduced and additive contributions, i.e., we have to determine
for each contribution on a graph
. After having determined the ground-state energy per site
and the hopping amplitudes
on the lattice
, we derive the effective one-quasiparticle Hamiltonian (
113). As we have seen in
Section 4.5, we can readily derive the one-quasiparticle gap as a series in the perturbation parameter
as in Eq. (
117).
1qp spectral weight: Lastly, we can do the same for one-quasiparticle spectral weights. We calculate the process amplitudes
Note that the process
does not satisfy
in general as a quasiparticle is created at
and subsequently hops to
. On the other hand 1qp hopping processes of the Hamiltonian fulfil
as long as the hopping amplitudes are real. While for hopping processes of the Hamiltonian we can use the same colour (for start and end vertex) as a topological attribute, here we must use two different colours leading to a smaller symmetry number of the coloured graphs. On the graph level only subgraph contributions must be subtracted
With the contributions
the spectral weight can be determined with Eq. (
131) and evaluating this quantity for example at the critical momentum
, we again obtain a series (
132) in the perturbation parameter
.
For the conventional embedding problem we can also consider several generalisation, which we want to briefly mention. First, we could consider quasiparticle processes between different particle types. In a 1qp process a particle flavour could change from
to
and the graph contribution would be denoted by including the additional indices
and
. The rest of the formalism is identical. Second, we may want to consider different interaction types like in Eq. (
138) or more expansion parameter like in Eq. (
137). In such a case, when considering coloured graphs, we also have to consider their coloured edges that must be matched with the appropriate bonds on the lattice. Then, the coloured graphs
are given by the tuple
where
are the edge colour attributes and
the vertex colour attributes. Now, every embedding factor
w needs to be determined using coloured graphs with respect to
and
. When using white graphs the additional edge colour can be ignored and the embedding is as before but we need to appropriately substitute the abstract expansion parameters with the actual physical expansions parameters.
By now, we have everything together to calculate the ground-state energy, the 1qp dispersion, and the 1qp spectral weight in the thermodynamic limit from a linked-cluster expansion set up as a full-graph decomposition. Therefore, we need to calculate the embedding factor, multiply them with the associated graph contribution and add up the resulting weighted contributions. The embedding factor can be determined using available graph libraries like the Boost Graph library [
203] as only automorphisms and monomorphisms (with colours) need to be determined. For long-range interactions however, we cannot just simply calculate the embedding factor because for every graph there are infinitely many possible embeddings even when accounting for translational invariance.
4.7.2. Embedding for Models with Long-Range Interactions
In this section we only consider a single physical perturbation parameter because there cannot be an additional functional dependency on a parameter when we later perform numerical Monte Carlo summation. If we have more than one physical perturbation parameter we have to sample the additional parameters and perform the Monte Carlo summation for each parameter ratio like it was done in [
31,
35,
206] for anisotropic XY interactions, distinct ladder interactions, and for XXZ interactions. We also need to use white graphs to account for long-range interactions. In the following, we restrict to a single quasiparticle flavour due to simplicity and a trivial unit cell as we have not generalised the algorithm to larger unit cells yet. The starting point to describe the embedding problem for long-range interactions are the embedding formulas given in the previous section in Eqs. (
149), (
151), and (
154) for the ground-state energy per site, 1qp hopping processes and the 1qp spectral weight respectively, which we have to rewrite and adapt for long-range interactions.
Ground-state energy: Starting with the embedding of the ground-state energy per site (
149), we can write
We have replaced
with
in the sum to emphasise that the graphs are not restricted to the nearest-neighbour lattice geometry anymore due to the long-range interactions making the lattice a fully connected graph. From the first to second line, we decomposed the ground-state energy contribution from graphs into contributions of individual orders
. Note, the minimum order of the ground-state energy is
. From the second to third line we introduced a second sum over sites
and restricted the sum over all graphs to a sum over graphs with a fixed number of vertices
. Next, we can reformulate the embedding factor in Eq. (
147) as
where we replace the expression for the number of monomorphisms (divided by
N) with a sum over all possible configurations
. A configuration is nothing else than the current embedding of graphs. But, as we will see later, we can calculate the contributions of multiple graphs simultaneously as the Monte Carlo sum only depends on the number of sites (we use "site" also as a synonym for "vertex position"). When using the word configuration we think about it as the current set of vertex positions on the lattice. The sum over all configurations comprises individual sums for each vertex over all lattice sites excluding configurations where vertices overlap as shown in the next subsection in Eq. (180).
4 At first sight we have not gained anything from rewriting the embedding factor. However, the white-graph contributions
still need to be replaced with the correct colour, i.e., the general expansion parameters need to be substituted with the algebraically decaying long-range interaction depending on the current configuration. In reality, the contribution
depends on the current configuration
c. Thus, replacing the expression with an explicit sum is necessary as the contribution for each configuration is different and
cannot just be a number. The substitution must look like
where the index
m of the sum runs over all monomials of the contribution
and the product is over all edges
of the graph
with the adjacent vertices
. The power law in the product arises from substituting the expansion parameters
on the edges
e with the appropriate algebraically decaying long-range interaction of the current embedding (c.f. Eq. (
140)). The adjacent vertices
are embedded on the lattice sites
and
with the distance
and the multiplicity
comes from the power
of the associated expansion parameter
. This way we can reduce the many expansion parameters from the white-graph contribution to a single physical perturbation parameter
and by reordering the expression (
156) of the ground-state energy, we have
We can define
where
is the integrand function and
denotes the associated sum over all possible embeddings on the lattice that will be evaluated using a classical Monte Carlo algorithm. Since Monte Carlo runs are usually done for a batch of different seeds we introduce an additional sum over seeds averaging the Monte Carlo runs, which we denote as
Eventually, this yields the expression
To express the ground-state energy
as a perturbative series, we write
where we have to sum up the contributions of multiple Monte Carlo runs to obtain the series prefactors
for
. The zeroth prefactor
is simply given by the ground-state energy of
.
1qp dispersion: We now turn to extracting the one-quasiparticle dispersion. In Eq. (
116) we have seen that the dispersion can be analytically determined in terms of the hopping amplitudes of Eq. (
151) up to some perturbative order. For nearest-neighbour interactions it is an analytic function in
, however, for long-range interactions there are infinitely many hoppings possible at any order so we can neither explicitly determine the hopping amplitudes
nor is it possible to have the dispersion as an analytical function in
. This would introduce an functional dependence in the integrand that we cannot sample. Instead, we will have to evaluate the dispersion for certain values
to get an explicit series in the perturbation parameter
. Evaluating Eq. (
116) at
and inserting Eq. (
151), we can write
where we rewrote the sum over
by introducing the function
Again, we can split
into contributions of individual orders
and split the graph set into subsets of fixed number of sites
, yielding
Note, here the second sum goes until while for the ground-state energy it runs only until . The maximal number of sites at a given order is because graphs with edges can maximally have sites. For the ground-state energy fluctuations every quasiparticle that is created has to be annihilated again, so at order a process can only touch maximal edges which restricts the sum to sites.
Now, we argue that we can drop the sum over
by thinking differently about the embedding problem for the dispersion. The information of the start and end vertex of the hopping process is encoded into vertex colours and when finding the subgraph monomorphisms for the embedding on the infinite lattice
the colours of the vertices must match with the coloured sites on the lattice, i.e., the hopping vertices are fixed to the hopping sites on the lattice. Since the long-range interactions allow any hopping – i.e., of any distance – at any order, it is not useful to think in this picture. Instead, we should think about the embedding problem analogous to the one for the ground-state energy, where no such hopping constraint exists and the embedding factor
is simply proportional to a sum over all configurations. This is valid as we let the sum over all lattice sites and account for constraints on
by multiplying with the symmetry number of the coloured graph
. The relevant hopping information of the vertices, that was previously fixed by coloured vertices, is anyway encoded into the cosine terms. Hence, we can make the substitution
where we account for the reduced symmetry of the graph due to the hopping by multiplying with the symmetry number of the coloured graph
. As before, we need to substitute the general white-graph contribution with the actual algebraically decaying long-range interactions of the current embedding
Inserting Eq. (
168) into Eq. (
167) we end up with the expression
For a lighter notation we again define the integrand function and the Monte Carlo sum
Introducing an average over a batch of seeds for the MC sum
, we obtain
The perturbative series of the dispersion evaluated at
can then be expressed as
The sum prefactors for can be determined by summing up the individual contributions from the Monte Carlo runs.
1qp spectral weight: Last, the evaluation for the spectral weight observable is analogous to the 1qp dispersion. The integrand and Monte Carlo sum are defined as
which we use to calculate
We determine by again calculating the series prefactors for and then determine the 1qp spectral weight with .
It should be noted that we have to perform
Monte Carlo runs for a series of order with . This means the number of runs grows quadratically with the maximal order .
So far, we have derived the necessary formalism for how to express the embedding problem of models with long-range interactions but we have not talked about how to evaluate the Monte Carlo sums for the integrand functions f. In the next section we investigate how to evaluate such sums by introducing a suitable Monte Carlo algorithm.
4.7.3. Monte Carlo Algorithm for the Long-Range Embedding Problem
We are left with evaluating the Monte Carlo sum
which runs over all configurations
of graphs. The embeddings on the lattice depend only on the number of vertices of a graph
and there is no constraint by the edge set as in the nearest-neighbour case because every site of the lattice interacts with any other site making it a fully connected graph. The only restriction for models with long-range interactions is that the vertices of the graph are not allowed to overlap as they do not self interact and an overlap would imply infinite interaction strength resulting from the term
. In conclusion, we can write the sum over all possible configurations as number-of-sites-many sums over all possible lattice positions
where the primed sums
for any vertex at position
is the short notation for excluding overlaps with any other vertex position. For example, for contributions with three sites on a one-dimensional chain, this sum would look like
Due to the overlap constraint these are nested high-dimensional summations over the integrand functions which are in general hard to solve. The dimensionality of the summation is given by because higher dimensions introduce additional sums for each component. If we wanted to evaluate the Monte Carlo sum in two dimensions for contributions with 8 sites, which already occurs in eighth order perturbation theory, we would have to determine the integral value of 16 nested sums over the integrand function . This makes it clear, that the evaluation of such sums becomes challenging very quickly. In the following we use the short notation interchangeable for both the ground-state energy integrand and the 1qp process integrands .
The first approach to tackle the problem of evaluating these sums using conventional numerical integration techniques was pioneered by S. Fey and K. P. Schmidt already in 2016 [
25]. They managed to successfully determine the phase diagram and critical exponents from the closing of the gap of the long-range transverse-field Ising model on a one-dimensional chain with ferro- and antiferromagnetic interactions. While they were successful with their approach over a large range of decay exponents in one dimension, the extraction of the critical properties for small decay exponents was challenging. The two-dimensional problem was out of reach with this approach as the number of nested sums doubles and the integral converges significantly slower. Here, Monte Carlo integration came into play as it is known to be a powerful integration technique for high-dimensional problems where conventional integration fails. The reason behind the slow convergence of such high-dimensional sums is often that the configuration space where the integrand mainly contributes to the integral is significantly smaller than the entire configuration space. In 2019 Fey et al. [
29] introduced a Markov-chain Monte Carlo algorithm to efficiently sample the relevant configuration space. They were able to determine the quantum-critical properties of the long-range TFIM on two-dimensional lattices to even higher precision than previously for the one-dimensional chain extending the accessible range of decay exponents without having to forfeit perturbative orders in higher dimensions.
In the following, we describe the Markov-chain Monte Carlo algorithm introduced by Ref. [
29] to evaluate the high-dimensional nested sums. To sample the relevant configuration space efficiently, we use importance sampling with respect to some convenient probability weight
with respect to a configuration
c and the associated partition function
. We can insert an identity to Eq. (
179) or Eq. (180) and rewrite it as
where
can be interpreted as the probability of being in configuration
c. The integrand now reads as the contribution
(we dropped the order
and momentum
as indices to lighten the notation) from configuration
c multiplied with its probability, which allows us to write the sum as the expectation value
of
with respect to the weight
. We later call this sum "target sum". Since the partition function
Z is not known a priori, we introduce also a "reference sum"
over a reference function
. We require this sum to be analytically solvable to avoid introducing an additional source of error. We denote its analytical expression as
. The reference function
should behave similar to the integrand function of interest
. This means that reference sum and target sum should have considerable contributions in the same area of configuration space and their asymptotic behaviour should be similar as well. Although we could make in principle an arbitrary choice for the reference function the latter properties guarantee to lead to good convergence. In one dimension we choose the reference integrand as
with
and the reference integrand exponent
that is a free simulation parameter. We can solve the reference sum as follows
where we accounted for
possibilities to randomly embed the vertices by ordering the indices of the sums and
is the Riemann
function. One major difference between the reference and the target sum is that in the target sum many different graph contributions contribute. In fact the reference sum above is exactly the contribution of order
of a chain graph with
vertices and the contribution from the associated target sum is the same up to a linear factor. In higher dimension we cannot choose a contribution proportional to the one of a chain graph anymore since it cannot be solved analytically. Instead we make simplifications to the reference sum and require that the reference sum is still good enough to capture the same properties as the target sum. We choose to decouple the dimensions in the reference integrand
and explicitly allow overlaps in the reference sum, such that it can be solved analytically as follows
Although the exponent
can be chosen freely, we want to achieve similar asymptotic behaviour as the target integrand, therefore we choose
in one dimension as the reference sum exactly behaves like the target integrand of a chain graph. In two dimension we made some simplifications to the reference sum and we have to adopt the parameter to
for
,
for
, and
for
. This is by no means a rigorous choice but it empirically proofed to produce good convergence [
29,
207]. Solving Eq. (
183) for
Z and inserting it into Eq. (
182), we obtain
We use the analytic expression of Eq. (
185) in 1d or Eq. (
187) in 2d for the reference sum
. We got rid of the partition function
Z and now can use this expression in our Monte Carlo run to determine the sum
using the analytic expression of the reference sum
while tracking the running averages in the numerator and denominator expressions.
We are left with just one missing ingredient that is the choice of the probability weight
. For our choice to be a probability weight, it must fulfil
and we want the weight to be the largest if both the reference and the target integrand contribute the most. An obvious choice may be the quadratic mean
that is always
and rotationally invariant in
and
. However, we also want both quantities to contribute equally to the probability weight on average over all configurations. As the contributions of target and reference sum may differ significantly, we introduce a factor for rescaling
that can be estimated in an in-advance calibration run. We then use an adjusted probability weight
for the actual Monte Carlo run. The weight needs to be evaluated at every Monte Carlo step to track the running averages of the numerator and denominator in Eq. (
188).
Now, we have everything together to describe the actual Monte Carlo algorithm. We employ a Markov-Chain Monte Carlo algorithm where we need to sample the configuration space
according to the probability weight
in Eq. (
191). Each configuration with a non-zero weight must be in principle accessible by the Markov chain. On the one hand, we propose a high acceptance rate of the Monte Carlo steps to sample the large configuration space efficiently not staying in the same configuration for too long. On the other hand, we want to sample closely confined configurations with rather small distances between vertices more often than configurations that are farther apart (configurations with large distances between vertices) such that we capture the asymptotic behaviour of the model. The interaction strength decays algebraically with the distance between vertices leading to smaller contributions for configurations in which sites are far apart. What we call a confined configuration therefore depends on the decay exponent of the long-range interaction
. In the algorithm we have the free exponent parameter
(for the reference sum) and
(for probability distributions) that can be changed to tweak this behaviour but are usually chosen similar or equal to
. In two dimensions we had to adapt the values of
to obtain a similar asymptotic behaviour for reference and target sum due the approximations we made. An optimal choice of these parameters ensures fast convergence of the Monte Carlo sum. The current embedding configuration should be represented as an array container where the entries are the positions of the graph vertices. In one dimension entries are simple integers while in higher dimensions the position needs to be represented as a mathematical vector. For small decay exponents
, very large distances can occur between vertices from time to time which need to be squared when calculating the absolute distance. This can lead to an integer overflow and therefore the use of 128 bit integer may be considered. Further, we define functions in the programme for the target integrand
and for the reference integrand
, where the current configuration is passed as a parameter and the function returns the contribution from the integrand evaluated for the current configuration.
Turning back to the sampling scheme, the idea of the Markov chain is to interpret the vertex positions on the lattices as random walkers. We randomly select a graph vertex and then draw a new position from a probability distribution such that the move fulfils the detailed balance condition. In each Monte Carlo step we perform the following two moves:
-
shift move:
-
This Monte Carlo move is implemented to introduce confined random fluctuations to the current configuration independent of the strength of the algebraically decaying long-range interactions. It is especially important for larger decay exponents
when the configurations are much more likely to be confined. First, we randomly select a vertex
drawn from a discrete uniform distribution with
. Second, for the fluctuation we draw a shift value
from a discrete uniform distribution
. In one dimension we have to draw a single time and in higher dimensions we draw repeatedly for each component. Subsequently, we add the shift to the position of the selected vertex and propose the position
We might have proposed a position that is already occupied by another vertex so we have to check for overlaps. In one dimension we reset the proposed position to the original one if there is an overlap while in higher dimensions we explicitly allow overlaps. As we remember from above, this distinction is also present in the reference sums in one dimension in Eq. (
185) compared to higher dimensions in Eq. (
187). If an overlap occurs in dimensions higher than one then the target summand is explicitly set to zero such that these configurations can not contribute (otherwise the sum would become infinity). Then, we calculate the Metropolis acceptance probability
by determining the probability weights
of the current and the proposed configuration. The result of the target and reference function calls should be saved into variables to prevent redundant and expensive function calls at each Monte Carlo step. Note, the transition weights
cancel out as we draw only from uniform distributions. Last, the minimum function is implemented by drawing a random number
and we accept the proposed move if
and update the current configuration. An example of such a shift move is depicted in
Figure 11a.
-
rift move:
-
In contrast to the previous move that should introduce fluctuations to the configuration independent of the current one and independent of the long-range interaction strength, "rift moves" are introduced to better capture the correct asymptotic behaviour induced by the algebraically decaying interactions. The moves are able to propose very large distances between vertices but are also able to do the opposite closing the "rift" between vertices when the configuration is split in essentially two clusters. At first, we select a site
from the vertex set with discrete uniform probability
, explicitly excluding the last site. In one dimension we can order the vertex set such that the first vertex is the one with the smallest positional value and the last the one with the largest value, so we order by
where
are vertex indices and
,
the associated sites on the lattice. The same ordering was also done when we solved the reference sum in Eq. (
185). In higher dimensions a similar ordering comes at a much higher computational cost so we stick to the vertex numbering given by the array indices, i.e., the order is
. Here, it is also important that the vertex labelling of the reference sum coincides with the labelling of the chain graph. To capture the physical asymptotics of the system, we draw random values from a
-function distribution. In one dimension we draw from
yielding a power-law distribution with
with the free exponent parameter
. We choose
for obvious reasons. The distance to the next vertex is given by
and
is the proposed distance drawn from the
distribution. Since we ordered by the position and only selected sites in
it is sufficient to draw positive values only. We shift all indices
according to
In higher dimensions we have no such ordering and therefore extend such a distribution to negative values – we refer to it as "double-sided"
-function distribution – and draw random values from
for each component. Note, the additional one is introduced to prevent divergence when sites overlap. After drawing the new distance
we shift all vertices componentwise with
according to
The underlying idea is that if there is a large distance between two vertices
and
we can close the "rift" of the entire configuration instead of introducing a new one between
and
. The transition weights for this move are given by
With these we can calculate the Metropolis-Hastings acceptance probability in one dimension