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 $\sim r^{-(d+\sigma )}$ with the distance r and the dimension d of the system, there are two distinct regimes: One for, $\sigma \le 0$ (strong long-range interaction) and another one for $\sigma >0$ (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₄ 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₃GeTe₂ 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₂Ti₂O₇ and Dy₂Ti₂O₇ [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 $U\left(1\right)$-symmetric magnetic interaction invariant under planar rotations and, third, Heisenberg interactions with a $SU\left(2\right)$-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 $U\left(1\right)$-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.

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 $\lambda$ 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 $d+1$ 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-trivial1 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 $f(x_1,x_2,...,x_n)$ is called a GHF, if there exist $a_1,a_2,...,a_n\in \mathbb{R}$ with at least one being non-zero and $a_f\in \mathbb{R}$ such that for $b\in {\mathbb{R}}_{+}$

The exponents $a_1,a_2,...,a_n$ are the scaling powers of the variables and $a_f$ 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 $x_i$ are related to physical parameters like the temperature T, control parameter $\lambda$, 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 $d_{\mathrm{uc}}$. The upper critical dimension is defined as a dimensional boundary such that for systems with dimension $d\ge d_{\mathrm{uc}}$, 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 $1/r^{d+\sigma }$, the upper critical dimension starts to decrease as a function of the decay exponent $\sigma$ [20,32,34]. In the limiting case of a completely flat decay ($\sigma =-d$) 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.

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 $f:\mathcal{C}\to \mathbb{R}$ over an integration space $\mathcal{C}$ 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 ${\omega }_i\in \mathcal{C}$ sampled according to a probability density function $P:\mathcal{C}\to {\mathbb{R}}_{\ge 0}$ (PDF) and the function $\tilde{f}\left(\omega \right)=\frac{f\left(\omega \right)}{P\left(\omega \right)}$ 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 $\pi \left(\omega \right)$ within a configuration space $\mathcal{C}$ arise which are used for the sampling of the PDF P being oftentimes not directly accessible. In the context of statistical physics, $\pi \left(\omega \right)$ is often chosen to be the relative Boltzmann weight $e^{-\beta E\left(\omega \right)}$ of each configuration $\omega$. While this relative Boltzmann weight is accessible as long as $E\left(\omega \right)$ is known, the full partition function to normalise the weights is in general not.

In order to efficiently sample the configuration space $\mathcal{C}$ according to the relative weights, the methods in this review use a Markovian random walk to generate $\{{\omega }_1,...,{\omega }_m\}$. Let ${\omega }_n$ be the random state of a random walk at a discrete step n. The state ${\omega }_{n+1}$ at the next step is randomly determined according to the conditional probabilities $T(\omega \to {\omega }^{\prime })$ (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 ${\omega }^{\left(0\right)}$ the start configuration. We require the Markovian random walk to fulfil the following conditions: First, the random walk should have a certain PDF $P\left(\omega \right)$ defined by the weights $\pi \left(\omega \right)$ in Eq. (39) as a stationary distribution. By definition, $P\left(\omega \right)$ 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 $\omega \in \mathcal{C}$ can be reached from any configuration ${\omega }^{\prime }\in \mathcal{C}$ 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 $T(\omega \to {\omega }^{\prime })$ 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 $T(\omega \to {\omega }^{\prime })$ are decomposed into propositions $\tilde{T}(\omega \to {\omega }^{\prime })$ and acceptance probabilities $p_{\mathrm{acc}}(\omega \to {\omega }^{\prime })$ as follows

The probabilities to propose a move $\tilde{T}(\omega \to {\omega }^{\prime })$ 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 $\tilde{T}(\omega \to {\omega }^{\prime })=\tilde{T}({\omega }^{\prime }\to \omega )$, Eq. (46) reduces to the Metropolis acceptance probabilities

As an example, we regard a classical thermodynamic system with Boltzmann weights $e^{-\beta E\left(\omega \right)}$ 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 $\omega$ to a configuration ${\omega }^{\prime }$ with a smaller energy $E\left({\omega }^{\prime }\right)<E\left(\omega \right)$ is always accepted independent of the temperature. On the other hand, the proposition to move to a configuration ${\omega }^{\prime }$ with a larger energy than $\omega$ 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.

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.