# Bosonized 1D quantum systems through enhanced Event-Chain Monte Carlo

Oscar Bouverot-Dupuis, Alberto Rosso, Manon Michel

► To cite this version:

Oscar Bouverot-Dupuis, Alberto Rosso, Manon Michel. Bosonized 1D quantum systems through enhanced Event-Chain Monte Carlo. Physical Review B, 2025, 112 (3), pp.035148. 10.1103/k43n-yz82. hal-04994186v2

HAL Id: hal-04994186

https://hal.science/hal-04994186v2

Submitted on 17 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Bosonized one-dimensional quantum systems through enhanced event-chain
Monte Carlo

Oscar Bouverot-Dupuis,¹,² Alberto Rosso,¹ and Manon Michel³

¹Université Paris Saclay, CNRS, LPTMS, 91405, Orsay, France

²IPhT, CNRS, CEA, Université Paris Saclay, 91191 Gif-sur-Yvette, France

³Laboratoire de Mathématiques Blaise Pascal UMR 6620,
CNRS, Université Clermont-Auvergne, Aubière, France.

(Dated: August 13, 2025)

We design an enhanced event-chain Monte Carlo algorithm to study 1D quantum dissipative systems, using their bosonized representation. Expressing the bosonized Hamiltonian as a path integral over a scalar field enables the application of Monte Carlo algorithms developed for classical systems. Specifically, we focus on a dissipative XXZ spin chain, exhibiting critical slowing down, minima degeneracy and long-range interactions. Addressing all three bottlenecks, we design an algorithm that combines local persistent event-chain Monte Carlo moves with global cluster moves, in a O(1)-complexity implementation. Through systematic performance analysis, we show that such an algorithm outperforms traditional Metropolis algorithms by more than a magnitude factor and is competitive with current state-of-the-art Quantum Monte Carlo algorithms. We then use this approach to determine the dissipative spin chain's phase diagram, thereby reinforcing prior analytical predictions.

# I. INTRODUCTION

The world of one dimensional (1D) quantum systems has, since nearly a century ago, intrigued by its phenomenology which sets it apart of higher dimensional quantum systems. The peculiarity of 1D quantum systems lies in the negligible role played by particle statistics. Intuitively, this stems from the fact that rotations used to exchange particles and probe their statistics cannot be performed in 1D systems. This idea is at the core of the bosonization mapping [1–6] which provides a unifying description of bosonic and fermionic 1D systems in terms of bosonic collective modes, hence its name. A bosonized system can then be represented in the path integral formalism by a classical scalar field in dimension (1 + 1), which, for fermions or spins, is a considerable simplification compared to the Grassman algebra or spin coherent states usually required to treat such problems [7]. This procedure has been mainly used to perform renormalization group (RG) studies of the bosonized field-theory, leading to a better understanding of various 1D quantum phenomena [8–11].

Among all 1D quantum systems, dissipative quantum systems have provided a privileged framework for studying the physics of open quantum systems. Dissipation is typically introduced by coupling the system to many bosonic degrees of freedom acting as a bath. The latter are usually simple enough that they can be exactly traced out in the path-integral formalism to yield an effective retarded interaction

for the system [12]. Dissipative quantum dynamics were first studied in the context of systems with a single degree of freedom such as a spin coupled to a bosonic bath in the so-called spin-boson model [13–15] or a Josephson junction coupled to a resistive environment [16–18]. More recently, attention has been drawn to the many-body counterparts of these dissipative systems. As previously stated, most of the analytical progress has come from studying the bosonized Hamiltonians [19–23]. However, the potential of leveraging advanced computational techniques, such as classical Monte Carlo (MC) algorithms, for a direct investigation of the (1 + 1)D bosonized field theory remains underappreciated. Most of the numerical effort [24–26] put into the study of 1D dissipative quantum system has indeed been performed directly on the microscopic Hamiltonians using standard numerical schemes for quantum systems such as exact numerical diagonalization [27], the density-matrix renormalization group [28, 29], or Quantum Monte Carlo (QMC) methods [30, 31]. However, it is in most cases impossible to quantitatively compare results obtained from the microscopic and bosonized Hamiltonians as the relation between microscopic and bosonized couplings is usually not known exactly. This begs the question of simulating directly the bosonized field theory. Such an approach was for instance followed to study the resistively shunted Josephson junction [17] where a MC algorithm coupling global cluster moves [32] to local Metropolis moves [33] was proposed to study the underlying (0 + 1)D field theory. A more
recent example can be found in the study of an XXZ spin chain subject to local dissipation. This model, which can be seen as a many-body generalization of the spin-boson, was studied by [22, 34] with an emphasis on analytical techniques, although a brief numerical analysis was conducted using unadjusted Langevin dynamics. This numerical analysis was restricted to characterizing the bulk of the phases and was not powerful enough to probe the phase transition. The bosonized dissipative spin chain indeed provides a prototypical example of a 1D dissipative system hard to simulate: it suffers from severe critical slowing-down, has an infinite number of degenerate classical minima, and a dissipative long-range interaction which is computationally costly to deal with.

In this work we address the simulation of the bosonized dissipative XXZ spin chain by combining state-of-the-art classical MC algorithms and targeting each of the three previously identified bottlenecks. We propose an event-chain Monte Carlo (ECMC) algorithm [35] enhanced with global cluster moves and benefiting from a $O(1)$-complexity implementation [36]. The ECMC is a class of algorithms which are non-reversible as they break the detailed balance and only satisfy the global balance. By generating persistent and correlated moves, ECMC algorithms produce numerical accelerations, in particular near phase transitions, and have therefore proven useful to simulate various particle [35, 37-40] and spin models [41, 42], polymers [43, 44], or more recently hard disk packings [45]. Coupling local moves to global cluster moves is similar to overrelaxation moves and effective in addressing the minima degeneracy bottleneck. It has already been used to study, for example, the $\phi^4$ theory [46] or the sine-Gordon model [47], with the local moves generated by standard Metropolis updates. Here, we show how to adapt this enhancement by cluster moves to the local persistent ECMC moves. Finally, we exploit the factorized structure of the acceptance rate present in both the ECMC and cluster moves to implement them at a reduced $O(1)$-complexity cost thanks to the clock method [36]. Our results show that such an enhanced-ECMC algorithm is at least comparable to state-of-the-art QMC approaches in the simulated system sizes and observables. Therefore, while still being straightforward to code, such an algorithm can be used to precisely determine the phase diagram of the dissipative spin chain.

The paper is organized as follows: Section II introduces the microscopic model for the dissipative XXZ spin chain and its mapping onto a bosonic field theory. We also give a brief summary of known results.

Section III describes the developed algorithm based on an ECMC scheme enhanced by cluster moves. The performance of the algorithm is then tested against other standard algorithms and a clear scaling reduction is shown. Eventually, in Sec. IV, we use our algorithm to obtain a precise phase diagram of the dissipative XXZ spin chain. A brief discussion of the results and concluding remarks are made in Sec. V. Throughout the article, we work with $\hbar = k_B = 1$.

A public version of the code is available [48].

## II. MODEL

### A. Microscopic model

We study an XXZ spin chain coupled to a set of local baths as in Refs. [22, 49, 50]. The XXZ spin chain is a 1D spin-$\frac{1}{2}$ chain of N spins, total length L, lattice spacing $a = L/N$, and whose Hamiltonian is

$$ \hat{H}_S = \sum_{j=1}^{N} J_z \hat{S}_j^z \hat{S}_{j+1}^z - J_{xy} (\hat{S}_j^x \hat{S}_{j+1}^x + \hat{S}_j^y \hat{S}_{j+1}^y). \quad (1) $$

Each site j is coupled through $\hat{S}_j^z$ to a local set of quantum harmonic oscillators $\{\hat{X}_{\gamma,j}\}$ as

$$ \hat{H}_{SB} = \sum_{j=1}^{N} \hat{S}_j^z \sum_{\gamma} \lambda_{\gamma} \hat{X}_{\gamma,j}, \quad (2) $$

and the baths are governed by the Hamiltonian

$$ \hat{H}_B = \sum_{j=1}^{N} \sum_{\gamma} \frac{\hat{P}_{\gamma,j}^2}{2m_{\gamma}} + \frac{1}{2} m_{\gamma} \Omega_{\gamma}^2 \hat{X}_{\gamma,j}^2. \quad (3) $$

The total Hamiltonian of the dissipative system is therefore $\hat{H} = \hat{H}_S + \hat{H}_{SB} + \hat{H}_B$. We assume periodic boundary conditions for the spin chain and restrict our analysis to the zero magnetization sector.

### B. Effective field theory

This model can be rewritten as a bosonic field theory [22] by first bosonizing the spin chain [2, 3, 6], and then integrating out the bath degrees of freedom à la Caldeira-Leggett [12, 13]. Out of completeness, we provide in the following a more heuristic semiclassical derivation but omit most mathematical details.
<figure><img src="image_2.png" /><figcaption>Figure 1. Schematic representation of the model. An XXZ spin chain has its spins coupled to independent and identical collections of harmonic oscillators. The orientation of the j-th spin is described by the angle $2\phi(x_j = ja)$ it makes with the alternating z-axis.</figcaption></figure>

First, anticipating the fact that the dissipative XXZ spin chain has a transition to an antiferromagnetic ordered phase [51], we semi-classically represent spins through a field $\phi(x)$ as

$$ \hat{S}_j^z \sim (-1)^j \cos[2\phi(x_j)], \quad (4) $$

with $x_j = ja$ the positions of the spins. The field $2\phi$ is simply interpreted as the angle between a spin and the alternating z axis, as illustrated in Fig. 1. Using this mapping, the canonical partition function $Z = \text{Tr} e^{-\beta \hat{H}}$ at inverse temperature $\beta = 1/T$ can be expressed in terms of a Euclidean field theory $Z = \int \mathcal{D}\phi e^{-S[\phi]}$. The spin chain Hamiltonian $\hat{H}_S$ contains interactions with coupling $J_{xy}$ that want to order the spins in the XY plane. Since $[\hat{S}^{x/y}, \hat{S}^z] \neq 0$, this amounts to scrambling $\hat{S}^z$ (and thus $\phi$) so the associated action is

$$ S_{LL} = \int \frac{dx \, d\tau}{2\pi K} \left[ u \left(\partial_x \phi(x, \tau)\right)^2 + \frac{1}{u} \left(\partial_\tau \phi(x, \tau)\right)^2 \right], \quad (5) $$

where $x \in [0, L]$ denotes the spatial coordinate and $\tau \in [0, \beta]$ is the imaginary time coordinate. This is the so-called Luttinger liquid (LL) action with parameters $K$ and $u$. Next, the interactions in $\hat{H}_S$ with coupling $J_z$ are derived from Eq. (4) by noting that $\hat{S}_j^z \hat{S}_{j+1}^z \sim -\frac{1}{2} \cos[4\phi(x_j)] + \text{cst.}$ This yields

$$ S_g = -\frac{g}{2\pi^2} \int \frac{dx \, d\tau}{a\tau_c} \cos[4\phi(x, \tau)], \quad (6) $$

with the short-time cutoff $\tau_c = a/u$. Finally, the baths generate an effective retarded spin-spin interaction through a dissipative kernel $K(\tau) \sim \alpha/\tau^{1+s}$ with $\alpha$ the spin-bath coupling strength and $s$ the bath exponent. Following the nomenclature introduced by Caldeira and Leggett, $s = 1$ is an ohmic bath, while $0 < s < 1$ is a subohmic bath and $s > 1$ is a superohmic bath. Since the baths are local in

space and couple to $\hat{S}_j^z$, this final contribution adds up to $S_{LL} + S_g$ to give the total effective action

$$ S = \int \frac{dx \, d\tau}{2\pi K} \left[ u \left(\partial_x \phi(x, \tau)\right)^2 + \frac{1}{u} \left(\partial_\tau \phi(x, \tau)\right)^2 \right] \\ - \frac{g}{2\pi^2} \int \frac{dx \, d\tau}{a\tau_c} \cos[4\phi(x, \tau)] \quad (7) \\ - \frac{\alpha}{2\pi^2} \int \frac{dx \, d\tau}{a\tau_c^{1-s}} \frac{\cos[2\phi(x, \tau)] \cos[2\phi(x, \tau')]}{|\tau - \tau'|^{1+s}}, $$

where the last integral is over $|\tau - \tau'| > \tau_c$ to avoid any short-time divergence. This action is that of the sine-Gordon model with an additional long-range interaction that is commonly encountered in dissipative quantum systems [16, 23, 26, 52]. The effective parameters $u, K, g, \alpha, s$ can be related to the microscopic ones through bosonization (see [22]).

### C. Summary of known results

The field theory (7) has been previously studied by means of a renormalization group (RG) analysis [22] and a mapping to a generalized Coulomb gas [34]. These works showed the existence of two phases, a Luttinger liquid (LL) and an antiferromagnet (AFM), separated by a Berezinsky—Kosterlitz—Thouless (BKT) transition. The LL is a quasi-long-range ordered phase described by the RG fixed point

$$ S_{LL}[\phi] = \int \frac{dx \, d\tau}{2\pi K_R} \left[ u_R (\partial_x \phi)^2 + \frac{1}{u_R} (\partial_\tau \phi)^2 \right], \quad (8) $$

where $u_R, K_R$ are renormalized couplings and $\phi$ is implicitly evaluated at $(x, \tau)$. This phase is expected from RG arguments to exist for $K_R \ge K_R^c = \max(1/2, 1 - s/2)$. For $K_R < K_R^c$, the system becomes an AFM and is well-approximated by the massive action

$$ S_{AFM}[\phi] = \int \frac{dx \, d\tau}{2\pi K_R} \left[ u_R (\partial_x \phi)^2 + \frac{(\partial_\tau \phi)^2}{u_R} + \frac{u_R}{\xi^2} \phi^2 \right], \quad (9) $$

where $\xi$ is the correlation length which sets the asymptotic decay of the two-point function $\langle\phi(x, \tau)\phi(0)\rangle_{AFM} \sim e^{-\sqrt{x^2+(u_R\tau)^2}/\xi}$.

From Eq. (4), it appears natural to define the staggered magnetization averaged over space and imaginary-time, as the operator

$$ m = \cos(2\phi(x, \tau) - 2\phi(x, \tau)), \quad (10) $$
where $\overline{f(x, \tau)} = (L\beta)^{-1} \int dx d\tau f(x, \tau)$. It was shown that $m$ acts as an (infinite-order) order parameter [22]. The subtraction of the average $\overline{\phi(x, \tau)}$ ensures the existence of $\langle m \rangle$, in spite of the ill-defined nature of the LL where $\phi$ is formally defined up to a constant. Using Eqs. (8,9), the finite size scaling of the order parameter $\langle m \rangle$, where $\langle O[\phi] \rangle = Z^{-1} \int D\phi O[\phi] e^{-S[\phi]}$, can be established in both phases for a system of size $L \times \beta$ with $\beta \simeq L$,

$$ \langle m \rangle_{\text{LL}} = \left( \frac{L_0}{L} \right)^{K_R} \underset{L \to \infty}{\longrightarrow} 0, \quad (11) $$

$$ \langle m \rangle_{\text{AFM}} = \left( \frac{L_0}{\xi} \right)^{K_R}, \quad (12) $$

where $L_0$ is a short-distance cutoff that does not depend on $K_R$.

### D. Discretized field-theory

For the field theory (7) to be well-defined, one has to introduce a short-distance cutoff [53] which we do by putting the field theory on a lattice. Since we want to investigate the zero-temperature and thermodynamic ($\beta, L \to \infty$) properties of the system, we consider a finite system of dimensions $L \times \beta$ as large as possible. Moreover, because the phase transition of the model belongs to the BKT universality class which has a dynamical exponent $z = 1$, we take $L = u\beta$ to mitigate finite-size effects. The continuous field $\phi(x, \tau)$ is then replaced by its discretized version $\phi_i$ such that, with $\hat{x}$ and $\hat{\tau}$ the unit vectors in space and time directions,

$$ \phi_{i=n\hat{x}+m\hat{\tau}} = \phi(x = n\alpha, \tau = m\tau_c), \quad i \in [1, N]^2, \quad (13) $$

with $N = L/a = \beta/\tau_c$. The periodic boundary conditions require that $\phi_{i+N\hat{x}} = \phi_i$ and $\phi_{i+N\hat{\tau}} = \phi_i$. With this discretized field, the action (7) becomes

$$ \begin{aligned} S(\phi) = & \sum_i \frac{1}{2\pi K} \left[ (\phi_i - \phi_{i+\hat{x}})^2 + (\phi_i - \phi_{i+\hat{\tau}})^2 \right] \\ & - \sum_i \frac{g}{2\pi^2} \cos(4\phi_i) \\ & - \sum_i \sum_{\substack{k=-\lfloor N/2 \rfloor \\ k \neq 0}}^{\lfloor(N-1)/2\rfloor} \frac{\alpha}{2\pi^2} \frac{\cos(2\phi_i) \cos(2\phi_{i+k\hat{\tau}})}{k^{1+s}} \\ & \equiv \frac{1}{2} \sum_i \left[ \sum_{j \sim i} S_q^{i,j}(\phi) + \sum_{k=-\lfloor N/2 \rfloor}^{\lfloor(N-1)/2\rfloor} S_c^{i,i+k\hat{\tau}}(\phi) \right], \end{aligned} \quad (14) $$

where $S_q^{i,j}(\phi)$ (resp. $S_c^{i,j}(\phi)$) is the pairwise quadratic (resp. cosine) interaction, $j \sim i$ denotes the nearest neighbors $j$ of $i$ and the bounds of the sum over $k$ correspond to the system's periodicity.

The equilibrium probability distribution is then $\pi(\phi) = e^{-S(\phi)/Z}$ with the partition function $Z$ defined as $Z = \int d\phi e^{-S(\phi)}$ with $d\phi = \prod_i d\phi_i$ the measure over fields. Although the partition function is formally divergent, $\phi$ being defined up to a constant, this is not an issue for our sampling procedure (see Appendix A for more details).

### III. ALGORITHM

An efficient algorithm sampling according to the action (14) must overcome three main computational bottlenecks, from both computational complexity and dynamical origins. In order to show how they arise in practice, let us consider a usual classical MCMC scheme sequentially generating single-site updates, e.g. $\phi_i \to \phi_i + \varepsilon$ with $i$ randomly picked among all sites and $\varepsilon$ some random increment (possibly depending on $\phi_i$) such that proposing $\phi_i \to \phi_i + \varepsilon$ and $\phi_i + \varepsilon \to \phi_i$ are equiprobable. At each step, such an update is accepted according to the following Metropolis acceptance ratio,

$$ p_{\text{Met}}(\phi, (i, \varepsilon)) = \exp \left( - \left[ \sum_{j \sim i} \Delta_{\varepsilon} S_q^{i,j}(\phi) + \sum_{k=-\lfloor N/2 \rfloor}^{\lfloor(N-1)/2\rfloor} \Delta_{\varepsilon} S_c^{i,i+k\hat{\tau}}(\phi) \right]_+ \right), \quad (15) $$

with $[a]_+ = \max(0, a)$ the positive part function and the $\Delta_{\varepsilon}$ variations of the quadratic/cosine pairwise interactions $S_q^{i,j}$ are

$$ \Delta_{\varepsilon} S_q^{i,j}(\phi) = \frac{\varepsilon^2}{2\pi K} + \frac{\varepsilon}{\pi K} (\phi_i - \phi_j), \quad (16) $$

$$ \begin{aligned} \Delta_{\varepsilon} S_c^{i,i+k\hat{\tau}}(\phi) = & \delta_{k,0} \frac{g}{\pi^2} \sin(2\varepsilon) \sin(4\phi_i + 2\varepsilon) \\ & + (1 - \delta_{k,0}) \frac{2\alpha}{\pi^2} \sin(\varepsilon) \frac{\sin(2\phi_i + \varepsilon) \cos(2\phi_i + k\hat{\tau})}{|k|^{1+s}}, \end{aligned} \quad (17) $$

with $\delta_{k,l}$ the Kronecker delta. This algorithm exhibits the following bottlenecks:

(1st bottleneck) Independent and local single-site updates most often result in diffusive dynamics that converge slowly. This issue is particularly severe in presence of strong and long-range correlations that are reinforced by the vicinity of the phase transition.
(2nd bottleneck) Multimodality, arising from the symmetry $S(\phi) = S(\phi + n\pi/2)$, $n \in \mathbb{Z}$ further exacerbates the slow convergence of the diffusive process.

(3rd bottleneck) Computing the acceptance probability (15) for the bath-induced long-range interactions $S_c^{i,i+k\hat{\tau}}$ requires $O(N)$ operations.

To tackle these three bottlenecks, we propose a variant of the non-reversible event-chain Monte Carlo (ECMC) algorithm [35] enhanced by a coupling to cluster-type moves. This algorithm significantly decreases critical slowing down, although not completely suppressing it. We furthermore leverage the *strict-extensive* nature of the long-range interactions to implement both types of moves in a $O(1)$ computational complexity [36]. The key idea underlying these three elements is to consider a *factorized* Metropolis acceptance probability [35], namely

$$p_{\text{FMet}}(\phi, (i, \varepsilon)) = \prod_{j \sim i} \exp(-[\Delta_{\varepsilon} S_{\text{q}}^{i,j}(\phi)]_{+}) \times \prod_{k=-\lfloor N/2 \rfloor}^{\lfloor(N-1)/2\rfloor} \exp(-[\Delta_{\varepsilon} S_{\text{c}}^{i,i+k\hat{\tau}}(\phi)]_{+}). \quad (18)$$

The factorized Metropolis filter is a valid acceptance ratio, as it obeys detailed balance. Although it yields on average a smaller acceptance rate than the standard Metropolis one, its product structure makes the interactions independent. Here, this factorization allows to treat separately the different interaction terms, which is key for efficiently implementing non-reversible schemes such as ECMC [35], cluster moves [32] and general computational complexity reduction schemes [36]. In the following, we detail the three points and eventually compare the performance of our enhanced ECMC algorithm with a standard Metropolis algorithm, with and without enhancement by cluster moves, and a standard ECMC algorithm.

## A. ECMC

### 1. Generating a non-reversible and rejection-free dynamics

An ECMC algorithm [35] is a continuous-time, rejection-free and irreversible sampling scheme. While directly identifying with a continuous-time piecewise deterministic Markov process (PDMP) [54], such algorithms were first introduced as the infinitesimal limit $\varepsilon \to 0$ of a persistent Markov

<figure><img src="image_3.png" /><figcaption>Figure 2. The two elementary steps of the ECMC algorithm. Top: a field component $\phi_i$ (in yellow) performs a ballistic shift until an event is triggered by one of its neighbors drawn in white. Bottom: when an event occurs, the neighbor involved in the event becomes the new moving field component.</figcaption></figure>

chain based on successive fixed $(i, \varepsilon)$ updates coupled to a factorized Metropolis acceptance rate (as in Eq. (18)). Infinitesimal updates $(i, \varepsilon)$ are proposed until a rejection eventually happens and the update direction is resampled and given another value $(j, \varepsilon')$, so that the correct equilibrium distribution is left invariant.

Such schemes were first introduced to upgrade the dynamics from diffusive and slow ones, as typically generated by Metropolis schemes, which is exactly our first bottleneck. However, the ECMC and the linked PDMP formalism allow for great versatility, which can be exploited to couple such continuous local moves to global discrete ones, in order to address the second bottleneck linked to the multimodality. This is discussed in Section III B.

For a general system, we recall that an ECMC algorithm is formally built from factorized $(i, \varepsilon)$ Metropolis updates following three distinct steps:

1. the state space $\Omega$ is augmented to $\Omega \times \mathcal{V}$ where $v = (i, e) \in \mathcal{V}$ encodes the current ballistic motion with $e$ the sign of $\varepsilon = e|\varepsilon|$. The target distribution over the augmented space is taken of the form $\phi(\phi)\mu(v)$ in order to retrieve $\pi(\phi)$ easily. The additional variable $v$ is commonly referred to as a *lifting variable* [55].

2. the infinitesimal update regime $\varepsilon \to 0$ limit is considered. The binomial law generating a
rejection from the factorized acceptance rate $p_{\text{FMet}}(\phi, (i, \varepsilon))$ transforms into a Poisson process of rate

$$ \lambda(\phi, (i, e)) = \lim_{\varepsilon \to 0} \frac{1 - p_{\text{FMet}}(\phi, (i, \varepsilon))}{|\varepsilon|}. \quad (19) $$

3. rejections according to the rate $\lambda(\phi, (i, e))$ trigger $v$-updates, called events, that ensure the global balance condition in the augmented state space $\Omega \times \mathcal{V}$.

Since the introduction of the ECMC method, it has been successfully applied to different systems thanks to the flexibility enabled by the freedom of choice of the lifting variable and its update scheme. Adapting this procedure to our system, the lifting variable is $v = (i, e) \in \mathcal{V} = [1, N]^2 \times \{-1, 1\}$ and we chose for $\mu(v)$ to be the uniform distribution over $\mathcal{V}$. In the infinitesimal regime $\varepsilon \to 0$, the factorized acceptance probability (18) gives rise to the Poisson rate,

$$ \begin{aligned} \lambda(\phi, (i, e)) &= \lim_{\varepsilon \to 0} \frac{1 - p_{\text{FMet}}(\phi, (i, \varepsilon))}{|\varepsilon|} \quad (20) \\ &= \sum_{j \sim i} \lambda_{\text{q}}^{i, j, e}(\phi) + \sum_{k = -\lfloor N/2 \rfloor}^{\lfloor (N-1)/2 \rfloor} \lambda_{\text{c}}^{i, i+k\hat{\tau}, e}(\phi), \end{aligned} $$

with $\lambda_{\text{q/c}}^{i,j,e}(\phi) = \lim_{\varepsilon \to 0} \frac{[\Delta_{\varepsilon} S_{\text{q/c}}^{i,j}]_+}{|\varepsilon|} = [e \partial_{\phi_i} S_{\text{q/c}}^{i,j}(\phi)]_+$ such that, from Eqs. (16,17),

$$ \partial_{\phi_i} S_{\text{q}}^{i,j}(\phi) = \frac{1}{\pi K} (\phi_i - \phi_j) \quad (21) $$

$$ \begin{aligned} \partial_{\phi_i} S_c^{i, i+k\hat{\tau}}(\phi) = \delta_{k,0} \frac{2g}{\pi^2} \sin(4\phi_i) \\ + (1 - \delta_{k,0}) \frac{2\alpha \sin(2\phi_i) \cos(2\phi_{i+k\hat{\tau}})}{\pi^2 k^{1+s}}, \end{aligned} \quad (22) $$

Thanks to the factorization (18), this Poisson process is a superposition of Poisson sub-processes of rates $\lambda_{\text{q/c}}^{i,j,e}$ corresponding to the different interaction terms. Therefore, an event triggered by the total Poisson process of rate $\lambda$ corresponds to the event of a single interaction term. As a rejection can be understood as an energy increment excess, a correct redirection $v \to v'$ should ensure a balanced redistribution of this excess. Here, and as is usually done in ECMC schemes [35, 41], we first leverage the pairwise symmetry

$$ [e \partial_{\phi_i} S_{\text{q}}^{i,j}(\phi)]_+ = [-e \partial_{\phi_j} S_{\text{q}}^{i,j}(\phi)]_+, \quad (23) $$

which relates the energy increments obtained by varying $\phi_i$ and $\phi_j$, so as to propose

$$ (i, e) \to (j, e) \text{ for an event triggered by } \lambda_{\text{q}}^{i, j, e}. \quad (24) $$

For the events triggered by $\lambda_c^{i, i+k\hat{\tau}, e}$, a similar pairwise symmetry does not exist. However, we can exploit the probability flow conservation itself [56, 57]. It relies here on the invariance of the uniform distribution for $v = (i, e)$ under the (1D) rotational symmetry $e \leftrightarrow -e$, and yields,

$$ \begin{aligned} & \sum_{e} ([e \partial_{\phi_i} S_c^{i, i+k\hat{\tau}}(\phi)]_+ + [e \partial_{\phi_{i+k\hat{\tau}}} S_c^{i+k\hat{\tau}, i}(\phi)]_+) \\ &= \sum_{e} ([-e \partial_{\phi_i} S_c^{i, i+k\hat{\tau}}(\phi)]_+ + [-e \partial_{\phi_{i+k\hat{\tau}}} S_c^{i+k\hat{\tau}, i}(\phi)]_+), \end{aligned} \quad (25) $$

This simple relation tells us, for instance, that an energy excess $[e \partial_{\phi_i} S_c^{i, i+k\hat{\tau}}(\phi)]_+$ arising from updating $\phi_i$ along $e$ can be redistributed over opposite $(-e)$-updates of $\phi_i$ and updates of $\phi_{i+k\hat{\tau}}$ along $e'$, according to their respective energy decrease $[e \partial_{\phi_i} S_c^{i, i+k\hat{\tau}}(\phi)]_+$ and $[-e' \partial_{\phi_{i+k\hat{\tau}}} S_c^{i, i+k\hat{\tau}}(\phi)]_+$. The stochastic sampling of the updated field, whether $\phi_i$ or $\phi_{i+k\hat{\tau}}$, prevents systematic backtracking, which would otherwise occur if $\phi_i$ were always kept fixed, necessitating repeated $e/-e$ flips that undo previous moves. Therefore, we propose, at events triggered by $\lambda_c^{i, i+k\hat{\tau}, e}$,

$$ (i, e) \to \begin{cases} (i, -e) & \text{with probability } p_c^{ik}(\phi) \\ (i + k\hat{\tau}, e_c^{ik}(\phi)) & \text{otherwise} \end{cases}, \quad (26) $$

with probability

$$ p_c^{ik}(\phi) = \frac{|\partial_{\phi_i} S_c^{i, i+k\hat{\tau}}(\phi)|}{|\partial_{\phi_j} S_c^{i, i+k\hat{\tau}}(\phi)| + |\partial_{\phi_{i+k\hat{\tau}}} S_c^{i, i+k\hat{\tau}}(\phi)|} \quad (27) $$

and,

$$ e_c^{ik}(\phi) = -\text{sign}(\partial_{\phi_j} S_c^{i, i+k\hat{\tau}}(\phi)), \quad (28) $$

the direction which lowers the energy. From Eq. (22), this amounts to selecting $j = i$ with probability

$$ p_c^{ik}(\phi) = \frac{|\tan(2\phi_i)|}{|\tan(2\phi_i)| + |\tan(2\phi_{i+k\hat{\tau}})|}, \quad (29) $$

and taking $j = i + k\hat{\tau}$ otherwise.

For $k = 0$, the scheme (26) comes down to $(i, e) \to (i, -e)$, which agrees with the trivial symmetry $[e \partial_{\phi_i} S_c^{i, i}(\phi)]_+ = [(-e) \partial_{\phi_i} S_c^{i, i}(\phi)]_+$ analogous to Eq. (23).

Such $v$-updates (24,26) satisfy the infinitesimal limit of the global balance, through the satisfaction of a *lifted* global balance [56, 57],

$$ \lambda(\phi, (j, -e')) = \sum_v \lambda(\phi, v) Q((\phi, v), (j, e')), \quad (30) $$
where $Q$ is the Markov kernel summarizing all the possible $v$-updates (24,26),

$$Q((\phi, (i, e)), (j, e')) = \sum_{i' \sim i} \frac{\lambda_{\mathrm{q}}^{i, i', e}(\phi)}{\lambda(\phi, (i, e))} \delta_{i', j} \delta_{e, e'} \\ + \sum_{k=-\lfloor N/2 \rfloor}^{\lfloor (N-1)/2 \rfloor} \frac{\lambda_{\mathrm{c}}^{i, i+k\hat{\tau}, e}(\phi)}{\lambda(\phi, (i, e))} (p_{\mathrm{c}}^{ik}(\phi)\delta_{i, j}\delta_{-e, e'} \\ + (1 - p_{\mathrm{c}}^{ik}(\phi))\delta_{i+k\hat{\tau}, j}\delta_{e^{ik}(\phi), e'}). \quad (31)$$

The validity of the lifted global balance is checked using Eqs. (23,25).

The ECMC thus generates a deterministic dynamics determined by $v$ that is resampled according to (31) following (24,26) at events happening at a rate (20) (see Fig. 2). It has now proven more efficient and robust to directly describe the generated Markov process by a PDMP [54, 57]. The entire PDMP can be characterized through its infinitesimal generator $\mathcal{A}$ [58]. It acts on smooth test functions $f$ as

$$\mathcal{A}f(\phi, v) = \lim_{t \to 0^+} \mathbb{E}_{\varphi_t} \frac{f(\varphi_t(\phi, v)) - f(\phi, v)}{t}, \quad (32)$$

where $\varphi_t$ is the time-evolution operator generated by a PDMP over a time $t$. Owing to the definition of the PDMP, the generator of the considered ECMC scheme is

$$\mathcal{A}f(\phi, v) = \underbrace{\frac{\partial_{\phi_i} f(\phi, v)e}{\text{deterministic drift}}}_{\text{event rate}} + \underbrace{\lambda(\phi, v) \sum_{v' \in \mathcal{V}} \underbrace{Q((\phi, v), v')}_{\text{jump } v \to v'} [f(\phi, v') - f(\phi, v)]}_{\text{event rate}} \quad (33)$$

with $\lambda$ and $Q$ respectively defined in (20) and (31). The derivation of the lifted global balance [56, 57] is straightforward as the invariance of $\pi$ is equivalent to

$$\sum_{v \in \mathcal{V}} \frac{1}{2N^2} \int_{\Omega} d\phi \mathcal{A}f(\phi, v)\pi(\phi) = 0. \quad (34)$$

The lifted global balance (30) is obtained after some integration by parts and its satisfaction is checked in Appendix B. The fact that $e$ and $-e$ balance each other out in (30) is the mark of the non-reversibility of the generated process. More details on PDMP and ECMC can be found in [54, 57]. In particular, we consider here a translational transport by updating a single field $\phi_i$ along $e$, but more general flows could be used [57]. Also, other factorization schemes than the one in (18) and/or choices

of $Q$ are possible. For instance, we have explored schemes which factorize the quadratic interactions in groups of 2, 3 or 4 pairwise interactions, thus associating a single event rate to each group instead of each pairwise interaction. We have also tried to separate the cosine interaction as $2\cos(2\phi_i)\cos(2\phi_j) = \cos(2\phi_i - 2\phi_j) + \cos(2\phi_i + 2\phi_j)$ which, factorizing both terms separately, can be simply dealt with by proposing the energy redirection $(i, e) \to (j, e)$ for the former, and $(i, e) \to (j, -e)$ for the latter. However, these schemes were always slower than the one detailed in the above.

Finally, ECMC schemes often include a refreshment mechanism, where the variable $v$ is resampled at some fixed time or following some homogeneous Poisson process. Both these schemes, among others, do not impact the invariance condition (30) [54]. Such refreshment mechanisms may prove useful to ensure ergodicity and sometimes even improve on the convergence [35, 59–61]. In the following, we opted for a refreshment event every fixed $T_{\text{Ref}}$, with $T_{\text{Ref}}$ much larger than the typical time between events, but much smaller than the total simulation duration.

## 2. Practical implementation

To numerically simulate this PDMP, one has to first compute the event times associated with the event rates, then perform the deterministic shift until the smallest event time, and finally update the lifting variable using the Markov kernel (31). When moving a field component $\phi_i$, one thus has to compute an event time for each of the 4 quadratic neighbors and for each of the $N$ cosine neighbors. A quadratic rate $[e\partial_{\phi_i} S_{\mathrm{q}}^{i,j}(\phi)]_+$ leads to an event time $t_{\mathrm{q}}^{i,j}$ which is a random variable such that

$$P(t_{\mathrm{q}}^{i,j} \ge t) = \exp \left[ - \int_{0}^{t} [e \partial_{\phi_i}(t')] S_{\mathrm{q}}^{i,j}(\phi(t')) \right]_{+} dt', \quad (35)$$

where $\phi_j(t) = \phi_j + \delta_{i,j} e$ is the field evolution imposed by the lifting variable $v = (i, e)$. It can be sampled by drawing $\nu \sim \text{ran}(0, 1)$ (the uniform distribution over $[0, 1]$) and solving

$$\nu = \exp \left[ - \int_{0}^{t_{\mathrm{q}}^{i,j}} [e \partial_{\phi_i}(t)] S_{\mathrm{q}}^{i,j}(\phi(t)) \right]_{+} dt \right]. \quad (36)$$
<caption>**Algorithm 1: Enhanced ECMC**</caption>

Input $\phi, n_{\text{sample}}, T_{\text{Ref}}, T_{\text{Sample}}, T_{\text{Cluster}}$;
Sample = {};
$i \leftarrow \text{choice}[1, \dots, N^2]$; // Lifting variable
e $\leftarrow$ choice$[-1, 1]$;
$t_r \leftarrow T_{\text{Ref}}$;
$t_s \leftarrow T_{\text{Sample}}$;
$t_{\text{Cluster}} \leftarrow T_{\text{Cluster}}$;

while True do

j_q $\leftarrow \underset{j \sim i}{\operatorname{argmin}}(t_q^{i,j} \leftarrow \text{Eq. 37})$;
t_q $\leftarrow \underset{j \sim i}{\operatorname{min}}(t_q^{i,j} \leftarrow \text{Eq. 37})$;
t_c, k_c $\leftarrow$ Alg. 3;
t_min = \min(t_c, t_q, t_{\text{Cluster}}, t_r);
if $t_s < t_{\text{min}}$ then // Outputting a sample
    $\phi_i \leftarrow \phi_i + e t_s$;
    Sample $\leftarrow$ Sample $\cup \{\phi\}$;
    $\phi_i \leftarrow \phi_i + e(t_{\text{min}} - t_s)$;
    $t_s \leftarrow t_s - t_{\text{min}} + T_{\text{Sample}}$;
    $n_{\text{sample}} \leftarrow n_{\text{sample}} - 1$;
    if $n_{\text{sample}} = 0$ then
        | Break
else
    $t_s \leftarrow t_s - t_{\text{min}}$;
    $\phi_i \leftarrow \phi_i + e t_{\text{min}}$;
    if $t_c = t_{\text{min}}$ then // Cosine event
        $t_r \leftarrow t_r - t_{\text{min}}$;
        $t_{\text{Cluster}} \leftarrow t_{\text{Cluster}} - t_{\text{min}}$;
        if $\text{ran}(0, 1) < p_c^{ik}(\phi)$ then
            $e, i \leftarrow -e, i$;
        else
            $e, i \leftarrow e_c^{ikc}(\phi), i + k_c \hat{\tau}$;
    end
else if $t_q = t_{\text{min}}$ then // Quadratic event
    $t_r \leftarrow t_r - t_{\text{min}}$;
    $t_{\text{Cluster}} \leftarrow t_{\text{Cluster}} - t_{\text{min}}$;
    $e, i \leftarrow e, j_q$;
else if $t_{\text{Cluster}} = t_{\text{min}}$ then // Cluster event
    $t_r \leftarrow t_r - t_{\text{min}}$;
    $t_{\text{Cluster}} \leftarrow T_{\text{Cluster}}$;
    $\phi \leftarrow$ Alg. 2;
else // Refreshment event
    $t_r \leftarrow T_{\text{Ref}}$;
    $t_{\text{Cluster}} \leftarrow t_{\text{Cluster}} - t_{\text{min}}$;
    $i \leftarrow$ choice$[1, \dots, N^2]$;
    $e \leftarrow$ choice$[-1, 1]$;

end

Return Sample

Using the definition of $\partial_{\phi_i} S_q^{i,j}$ in Eq. (21) leads to

$$t_q^{i,j} = e(\phi_j - \phi_i) + \sqrt{[e(\phi_i - \phi_j)]_+^2 - 2\pi K \ln \nu}. \quad (37)$$

Similar analytic expressions can be found for the $N$ cosine event times $t_c^{i,j}$. However, computing all

of these $N$ event times is unnecessarily costly and can be drastically reduced by thinning the cosine interactions [36, 62], as detailed in Section III C and Alg. 3.

Eventually, as the generated PDMP process is continuous in time, samples are outputted every fixed time interval $T_{\text{Sample}}$, corresponding here to fixed traveled distance as $e$ is of constant norm. As for the refreshment mechanism, it is possible to output samples according to exponentially distributed random times. A pseudo-code that sums up such ECMC implementation can be found in Alg. 1. It already describes the enhancement by cluster moves, introduced at fixed times $T_{\text{Cluster}}$, that we present in the next Section.

## B. Enhancing ECMC by cluster moves

We now directly address the second dynamical bottleneck created by the symmetry $S(\phi) = S(\phi + n\pi/2)$ with $n \in \mathbb{Z}$. It can be mitigated by efficiently moving the field from a minimum $\phi = n\pi/2$ of the action to another. The idea we propose is then to couple the continuous, rejection-free and non-reversible but local updates of standard ECMC schemes with rejection-free and global but reversible and discrete moves. Here, we leverage the symmetry $S(\phi) = S(n\pi/2 - \phi)$ to design cluster moves based on the reflection

$$R_n \phi_i = n\pi/2 - \phi_i. \quad (38)$$

As the $R_n$s are all involutions that leave the system globally invariant, this enables a correct cluster algorithm [32, 63]. Note that the reflection $R_n$ exploits the shift symmetry $\phi \to \phi + n\pi/2$ and the particle-hole symmetry $\phi \to -\phi$. While the shift symmetry (at least modulo $\pi$) is a consequence of the compact nature of the boson $\phi$ and applies to any bosonized system, the latter particle-hole symmetry may be absent in more general systems. This would render the cluster algorithm inefficient, but we expect that it could be replaced by a worm-type algorithm [64]. In the following, we consider a single-cluster construction as in the Wolff algorithm [65] and illustrated in Fig. 3. Since the reflections $R_n$ do not change the initial value of the field up to some modulo $\pi/2$, such cluster moves are not ergodic on their own and are more akin to overrelaxation moves.

We now recall the single-cluster recursive construction [65]. First, one uniformly picks an initial node $i$ and adds it to the cluster. Any node $j$ interacting with $i$ is then added to the cluster according
<figure><img src="image_4.png" /><figcaption>Figure 3. Basic principle of the cluster algorithm. Top: a cluster (pink dots) is selected and a flipping plane (in translucent pink) is set at one (or in between two) of the nearby potential minima located at $n\pi/2$, $n \in \mathbb{Z}$. Bottom: the entire cluster is flipped in respect to the flipping plane.</figcaption></figure>

to some probability $p_{i,j}(\phi)$. Once every node $j$ interacting with $i$ has been checked, the nodes interacting with the newly added nodes are checked, and so on until the process terminates. The entire cluster is then reflected according to $R_n$.

The probability $p_{i,j}(\phi)$ corresponds to a rejection of some factorized Metropolis filter. Here, we consider a factorization different from (18), and group all interactions between two fields $\phi_i, \phi_j$ together, i.e.

$$p_{i,j}(\phi) = 1 - \exp(-[\Delta S_{\text{pair}}^{R_n}(\phi_i, \phi_j)]_+), \quad (39)$$

where we have introduced the total pairwise action change $\Delta S_{\text{pair}}^{R_n}(\phi_i, \phi_j)$ defined as

$$\Delta S_{\text{pair}}^{R_n}(\phi_i, \phi_j) = \sum_{j' \sim i} \delta_{j,j'} \Delta_{\varepsilon} S_{q}^{i,j'}(\phi) + \sum_{k=-\lfloor N/2 \rfloor}^{\lfloor (N-1)/2 \rfloor} \delta_{j,i+k\hat{\tau}} \Delta_{\varepsilon} S_{c}^{i,i+k\hat{\tau}}(\phi), \quad (40)$$

with $\varepsilon = R_n\phi_i - \phi_i$.

There remains the choice of the reflection operator $R_n$, i.e. the value $n$. A reflection operator $R_n$ is typically useful if it flips fields $\phi_i$ such that $R_n\phi_i$ is not too far from $\phi_i$. Indeed, if the initial field component $\phi_i$ is flipped far away from the rest of the field, the probability $p_{i,j}(\phi)$ will likely be large, thus adding many sites to the cluster. This will tend to flip the entire system, resulting in a trivial, but costly, move.

We therefore condition the choice of the reflection $R_n$ on the initial field component $\phi_i$ by uniformly picking

$$n \in \lfloor 2m - 2, 2m - 1, 2m, 2m + 1, 2m + 2 \rfloor, \quad (41)$$

where $m = \lfloor \frac{2}{\pi}(\phi_i + \frac{\pi}{4}) \rfloor$ is the index of the minimum closest to $\phi_i$. Such a flip $R_n$ shifts $\phi_i$ of maximum $\pm 2$ minima, which we have numerically determined to be optimal (see Appendix C). A proof that this algorithm satisfies the detailed balance condition, thus proving its validity, can be found in Appendix D.

Eventually, we incorporate the cluster moves into the PDMP scheme by performing them at fixed time intervals $T_{\text{Cluster}}$. The PDMP then resumes by first randomly drawing the lifting variable $v$. As the cluster moves leave the correct Gibbs measure $\pi$ invariant, such fixed-time additions to the PDMP dynamics also leave $\pi$ invariant [54]. Furthermore, as for the event-time generation, the long-range interactions are dealt through the complexity-reduction clock method detailed in the next Section. A pseudo-code that sums up the cluster generation process can be found in Alg. 2, along with the complexity-reduction implementation.

### C. $O(1)$-complexity reduction

We now address the third and last bottleneck of the $O(N)$-computational complexity linked to the direct computations of the $N$ event times $t_c^{i,j}$ in the PDMP dynamics and the $N$ interacting fields to check in the cluster moves. Both stem from the cosine interactions in the long-range dissipative term. Fortunately, both the events and cluster moves are based on factorized rates as in Eqs. (39,20) which enable the use of complexity reduction schemes [62, 66–68] unified in the clock method [36]. In Ref. [36], it was shown that the achievable acceleration separates models into three classes defined from the strength of their long-range interactions and/or of the frustration. Our system belongs to the so-called strict-extensivity class [36], i.e.

$$\sum_{\substack{k=-\lfloor N/2 \rfloor \\ k \neq 0}}^{\lfloor (N-1)/2 \rfloor} \max |\Delta_{\varepsilon} S_c^{i,i+k\hat{\tau}}| \\ \leq \sum_{\substack{k=-\lfloor N/2 \rfloor \\ k \neq 0}}^{\lfloor (N-1)/2 \rfloor} \frac{2\alpha}{\pi^2} \frac{1}{|k|^{1+s}} \sim N^0. \quad (42)$$

Therefore, we can achieve the maximal $O(N)$ to $O(1)$ (i.e. $O(N^0)$) complexity reduction. We now
<caption>**Algorithm 2:** Complexity reduction for cluster generation</caption>

**Input** $\phi$;
$i_0 \leftarrow \text{choice}[1, \dots, N^2]$; // initial node
$C \leftarrow \{i_0\}$; // Add $i_0$ to cluster C
$S \leftarrow \{i_0\}$; // Add $i_0$ to the stack S
$m \leftarrow \lfloor \frac{2}{\pi}(\phi_{i_0} + \frac{\pi}{4}) \rfloor$;
$n \leftarrow \text{choice}[2m - 2, 2m - 1, 2m, 2m + 1, 2m + 2]$;
**while** $S \neq \emptyset$ **do**

$i \leftarrow \text{Choice}(S)$;
$S \leftarrow S \setminus \{i\}$;
**for** $j \sim i, j \notin C$ **do** // Step 1: short-range
$p_{i,j}(\phi) \leftarrow \text{Eq. (39)}$;
**if** $\text{ran}(0, 1) < p_{i,j}(\phi)$ **then**
$C \leftarrow C \cup \{j\}$;
$S \leftarrow S \cup \{j\}$;
**end**

**end**
$n \sim f(n, \lambda_{\text{Clu}}^B)$; // Step 2: long-range
**for** $l = 1, 2, \dots, n$ **do**
  Pick
    $k_c \in [-\lfloor N/2 \rfloor, \lfloor (N-1)/2 \rfloor] \setminus \{-1, 0, 1\}$
    according to $\lambda_{k_c}^B / \lambda_{\text{Clu}}^B$ using a Walker
    table;
    $j \leftarrow i + k_c \hat{\tau}$;
  **if** $j \notin C$ and $\text{ran}(0, 1) < p_{i,j}(\phi) / p_{k_c}^B$ **then**
    $C \leftarrow C \cup \{j\}$;
    $S \leftarrow S \cup \{j\}$;
  **end**
**end**

**end**
**for** $i \in C$ **do**
  $\phi_i \leftarrow R_n \phi_i$
**end**
**Return** $\phi$;

<caption>**Algorithm 3:** Complexity reduction for event generation by $\lambda_c^{i,j}$</caption>

**Input** $(\phi, (i, e))$;
$t_c \leftarrow 0$;
**while** True **do**
  $\nu \sim \text{ran}(0, 1)$;
  $t_c \leftarrow t_c - \ln \nu / \lambda_{\text{EC}}^B$; // Time to next bound
  **event**
    Pick $k_c \in [-\lfloor N/2 \rfloor, \lfloor (N-1)/2 \rfloor]$ according to
    $\lambda_{k_c}^B / \lambda_{\text{EC}}^B$ using a Walker table;
  **if** $\text{ran}(0, 1) < \lambda_c^{i, i+k_c \hat{\tau}, e}(\phi_j + e t_c \delta_{i,j}) / \lambda_{k_c}^B$
    **then**
      **Break**; // Generated true event
  **end**

**end**
**Return** $(t_c, k_c)$

describe the implementation for cluster moves and event computations. Further details about the generality of the clock method and its application to a general factorized Monte Carlo procedure, including any factorized Metropolis scheme, can be found in [36].

Applied to cluster moves, the clock method is used to find in $O(1)$-complexity all sites to add to the cluster from a given site $i$. From (39), a long-range neighbor $j$ of site $i$ is added to the cluster with probability $p_{i,j}(\phi)$. For now, $p_{i,j}(\phi)$ explicitly depends on the current field configuration $\phi$ so we need $O(N)$ operations to check all neighbors. We decide to deal with the nearest neighbors $j \sim i$ by explicitly computing $p_{i,j}(\phi)$, and only use a complexity-reduced implementation for distant neighbors $j = i + k\hat{\tau}$, $|k| > 1$. Following the clock method [36, 66–68], such a distant neighbor is added to the cluster if it is first accepted with the bound probability $p_k^B \ge p_{i,i+k\hat{\tau}}$ and then accepted with the resampling probability $\frac{p_{i,i+k\hat{\tau}}}{p_k^B}$. In practice, we choose the tightest possible bound

$$p_k^B = 1 - \exp\left(-\frac{2\alpha}{\pi^2} \frac{1}{|k|^{1+s}}\right). \quad (43)$$

Since the $p_k^B$'s do not depend on $\phi$, checking the bound acceptances of all neighbors of $i$ in $O(1)$ time is now a standard procedure [66–68]. In the following, we follow the ideas of Fukui and Todo [67] to sample the bound probabilities $p_k^B$. The neighbor $i + k\hat{\tau}$ is added with probability $p_k^B$, which is expressed as drawing an integer $n_k$ from the Poisson distribution $f_k(n_k) = \frac{e^{-\lambda_k^B} (\lambda_k^B)^{n_k}}{n_k!}$ of rate $\lambda_k^B = -\ln(1 - p_k^B) = \frac{2\alpha}{\pi^2} \frac{1}{|k|^{1+s}}$ and adding site $i + k\hat{\tau}$ if $n_k > 0$. Owing to the properties of Poisson distributions, the variable $n = \sum_{|k|>1} n_k$ follows a Poisson distribution of rate $\lambda_{\text{Clu}}^B = \sum_{|k|>1} \lambda_k^B$. Each of the $n$ events generated is then assigned using a Walker alias table [69, 70] to the neighbor $i + k\hat{\tau}$ with probability $\lambda_k^B / \lambda_{\text{Clu}}^B$. The complete procedure is thus:

* Draw $n$ events from the Poisson distribution of global rate $\lambda_{\text{Clu}}^B$.

* Assign each event to an interaction at distance $k$ ($|k| > 1$) with probability $\lambda_k^B / \lambda_{\text{Clu}}^B$ using a Walker table.

* Resample each event as actually adding site $i + k\hat{\tau}$ with probability $\frac{p_{i,i+k\hat{\tau}}}{p_k^B}$.

It is clear that the complexity identifies with the average number of events generated. This number
<figure><img src="image_5.png" /><figcaption>Figure 4. Integrated auto-correlation time $\tau_{\text{int}}$ for the staggered magnetization $m = \cos(2\phi_i - 2\phi_j)$ deep in the AFM (left), at the phase transition (middle), and deep in the LL phase (right). The other parameters are $s = 0.5$, $\alpha = 1$ and $g = 1$. The critical value $K = 0.85$ at the transition was extracted from the results of Sec. IV B. The error bars were obtained by averaging over 40 runs with independent gaussian initial conditions $\phi_i \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0, 1)$.</figcaption></figure>

is $\lambda_{\text{Clu}}^B = \sum_{|k|>1} \frac{2\alpha}{\pi^2} \frac{1}{|k|^{1+s}} \sim O(1)$, as expected as the system belongs to the strict-extensivity class. A pseudo-code that sums up the reduction scheme for cluster generation can be found in Alg. 2.

The clock-method is also useful for the ECMC moves to deal with the $N$ cosine events. An event generated by the rate $\lambda_c^{i,i+k\hat{\tau},e}(\phi)$ can be interpreted as an event triggered by a bound rate $\lambda_k^B \ge \lambda_c^{i,i+k\hat{\tau},e}(\phi)$ and then resampled with probability $\frac{\lambda_c^{i,i+k\hat{\tau},e}(\phi)}{\lambda_k^B}$ as a true event. The tightest bound is $\lambda_{k \neq 0}^B = \frac{2\alpha}{\pi^2} \frac{1}{|k|^{1+s}}$ and $\lambda_0^B = 2g/\pi^2$. Using again the fact that the sum of Poisson processes of rates $\lambda_k^B$ is a Poisson process of rate $\lambda_{\text{EC}}^B = \sum_k \lambda_k^B$, simulating this process comes down to a thinning procedure [36, 38, 62, 71]:

* Sample the first bound event time $t_c$ following the Poisson process of rate $\lambda_{\text{EC}}^B$, i.e. $t_c = -\ln \nu / \lambda_{\text{EC}}^B$ with $\nu \sim \text{ran}(0, 1)$.

* Assign this bound event to an interaction at distance $k$ with probability $\frac{\lambda_c^B}{\lambda_{\text{EC}}^B}$ using a Walker table.

* Resample this bound event as a true event with probability $\lambda_c^{i,i+k\hat{\tau},e}(\phi)/\lambda_k^B$. If this is not a true event, go to step 1.

The complexity here identifies with the number of bound events one has to check before finding a true

event. This is roughly the inverse of the average resampling probability $\sum_k \frac{\lambda_c^B}{\lambda_{\text{EC}}^B} \frac{\lambda_c^{i,i+k\hat{\tau},e}(\phi)}{\lambda_k^B} \sim O(1)$ since the bound is tight. A pseudo-code that sums up the reduction scheme for the event times can be found in Alg. 3.

### D. Algorithm performance

To test the performance of our algorithm, we study the staggered magnetization $m = \cos(2\phi_i - 2\phi_j)$ introduced in Sec. IIC. Similarly to the situation in spin systems [41, 42], we expect this observable to be a good indicator of algorithmic performance as it quantifies all the fluctuations of $\phi$ around its average $\bar{\phi}_i$ and should capture the slowest relaxation dynamics, as the average itself $\bar{\phi}_i$ thermalizes quickly thanks to the ECMC and cluster moves. We test the four following algorithms: the Metropolis-Hastings algorithm (Met), the Metropolis-Hastings algorithm with cluster moves (Met-Clu) as described in Sec. III B, the event-chain Monte Carlo algorithm without cluster moves (ECMC), and the event-chain Monte Carlo algorithm with cluster moves (Clu-EC). We do not implement any complexity reduction scheme on the Met because the gain would be marginal since its moves are not persistent like in the ECMC. The al-
<figure><img src="image_6.png" /><figcaption>Figure 5. The averaged staggered magnetization $\langle m \rangle = \langle \overline{\cos(2\phi_i - 2\phi_i)} \rangle$ as a function of the system size $N = 8, 16, 32, 64, 128, 256, 512, 1024$ for a subohmic (left) and superohmic (right) bath. The curving lines correspond to the AFM while the straight lines are in the LL. The red regions are bounded by small-$N$ fits $\propto N^{-K_R}$ where $K_R$ gives the renormalized Luttinger parameter. The error bars are smaller than the marker size and were obtained by averaging over 40 runs with independent gaussian initial conditions $\phi_i \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0, 1)$.</figcaption></figure>

gorithms including cluster moves share their computational time equally between cluster moves and local updates (implemented through the Met or the ECMC). In order to draw a fair comparison between algorithms, we define the algorithmic time $t$ expressed in sweeps (i.e. $N^2$ operations) as increasing by 1 each time $N^2$ pairwise interactions are evaluated. The CPU time may not be perfectly reflected in the algorithmic time but we expect both to share the same scalings and the latter to be less sensitive to specific code implementation details.

We consider the autocorrelation function defined as

$$C_m(t) = \frac{\langle m(t)m(0) \rangle - \langle m \rangle^2}{\langle m^2 \rangle - \langle m \rangle^2}. \quad (44)$$

For all algorithms tested here, $C_m(t)$ decays exponentially and we extract its characteristic time scale through the integrated autocorrelation time $\tau_{\text{int}} = \frac{1}{2} \sum_{t=-\infty}^{+\infty} C_m(t)$. The results for all four algorithms are shown in Fig. 4. Deep in both phases, the ECMC and Clu-EC both exhibit $\tau_{\text{int}} \sim N^0$ while the Met and Clu-Met have $\tau_{\text{int}} \sim N^1$. We interpret this as coming from the complexity reduction that was implemented in the ECMC and the cluster algorithm. Deep in the AFM ($K = 0.1$), the algorithms with cluster moves are slightly slower than those without. This comes from the fact that the field $\phi$ gets

trapped around a single minimum, so any attempted cluster move will flip the entire field, resulting in a trivial move. At the phase transition, the Met has $\tau_{\text{int}} \sim N^3$ and the Clu-Met seems to have a similar large-size behavior although with a much better prefactor. This is traced back to the diffusive dynamics of the Met which accounts for a scaling $\sim N^2$ combined with the absence of complexity reduction which adds an extra $\sim N^1$. Next, the integrated correlation of the ECMC behaves as $\tau_{\text{int}} \sim N^2$ which signals diffusive dynamics with complexity reduction. Finally, the autocorrelation time of the Clu-EC shows evidence for a $\tau_{\text{int}} \sim N^{1.6}$ scaling, resulting in a clear reduction of the critical slowing down phenomenon, in addition to complexity reduction. It is very clear that the superior algorithm is the Clu-EC which, at the transition, is roughly a 1000 times faster than the Met for $N = 128$, allowing for simulations up to $N = 1024$ to be done in a matter of days on a standard computer cluster.

## IV. SIMULATION OF THE DISSIPATIVE XXZ SPIN CHAIN

We now use the ECMC algorithm with cluster moves (Clu-EC) to study the phase diagram of the dissipative XXZ spin chain introduced in Sec. II.
<figure><img src="image_7.png" /><figcaption>Figure 6. Left: phase diagram for a subohmic bath ($s = 0.5$). Right: phase diagram for a superohmic bath ($s = 1.5$). We distinguish two phases: an antiferromagnet (AFM) and a Luttinger liquid (LL). The numerical data are the dark dots. The error bars are shown but are often smaller than the marker size and were obtained using the procedure described in Fig. 7.</figcaption></figure>

A first approach consists of computing the averaged staggered magnetization $\langle m \rangle$ which, using its finite-size scalings (11,12), allows to distinguish both phases. However, this approach is not the most efficient as it requires exploring large system sizes for a limited precision in return. We therefore follow a second approach where the phase transition is located from the known value of the coupling $K_R$ at the transition, as is usually done for BKT phase transitions [67, 72, 73].

### A. Large system size determination of the phase diagram

The averaged staggered magnetization $\langle m \rangle$ can be computed up to large sizes to determine the location of the phase transition. In Fig. 5, we computed $\langle m \rangle$ for a subohmic ($s = 0.5$) and a superohmic ($s = 1.5$) bath up to sizes $N = 1024$. The data splits into two types of curves: those which in the log-log plot are lines and thus belong to the LL phase, and those which deviate from a perfect straight line and belong to the AFM. In the LL, fitting the curves gives access to the renormalized Luttinger parameter $K_R$ from the finite size scaling $\langle m \rangle_{LL} \propto N^{-K_R}$. Fitting the data above and below the transitions gives estimates for the critical value $K_R^c$ of $K_R$ at the transition. For $s = 0.5$, we find $K_R^c \in [0.74, 0.79]$, while for $s = 1.5$ we find $K_R^c \in [0.48, 0.54]$. While the lower bounds of these estimates are fully trustworthy as they correspond to curves bending and therefore in the AFM,

the upper bounds have to be considered with greater caution as they correspond to curves which remain straight up to $N = 1024$ with no guarantee that this shall remain at larger $N$. Nevertheless, these results seem to confirm the RG analysis in Ref. [22] stating that the transition lies at $K_R^c = \max(1 - s/2, 1/2)$, i.e. $K_R^c = 0.75$ for $s = 0.5$ and $K_R^c = 0.5$ for $s = 1.5$.

### B. Precise phase diagram

Having confirmed the validity of the RG prediction $K_R^c = \max(1 - s/2, 1/2)$ in the previous section, we can now use this knowledge to precisely determine the phase diagram. Since, in the LL, the finite-size scaling $\langle m \rangle_{LL} = (N_0/N)^{K_R}$ holds even for the smallest sizes [74], we can determine $K_R$ only from the system sizes up to $N = 128$, saving computational time and enabling the computations of more points on the phase diagram. For a given set of parameters, we extract $K_R$ by fitting the scaling of $\langle m \rangle$ for $N = 8, 16, 32, 64, 128$ and $N = 16, 32, 64, 128$ (see Fig. 7). At large $K$, since the system has a well-defined LL finite-size scaling even at small sizes, the fitting of $\langle m \rangle$ works well and the determination of $K_R$ comes with little error. At small $K$ and entering the AFM, the LL finite-size scaling does not work anymore, so extracting an apparent $K_R$ yields large error bars. We stress that these error bars reflect the physical breakdown of the LL phase and not a poor convergence of the algorithm. For the two curves shown in Fig. 7, one can directly read out the
<figure><img src="image_8.png" /><figcaption>Figure 7. Extraction of $K_R$ from the fitting of the staggered magnetization $m$ for sizes $N = 8, 16, 32, 64, 128$ ($N_{min} = 8$) and $N = 16, 32, 64, 128$ ($N_{min} = 16$). The fit works well at large $K$ (LL phase) while it does not at smaller $K$ (AFM phase). The inset plot shows the value of $K^c$ obtained for both curves. Since most of the error on the determination of $K^c$ comes from the small sizes, we keep the value of $K^c$ for $N_{min} = 16$ and use that for $N_{min} = 8$ to estimate the error, leading to $K^c = 0.866(6)$.</figcaption></figure>

critical microscopic coupling $K^c$ corresponding to $K_R^c = \max(1 - s/2, 1/2)$. Since the finite-size scaling gets more accurate at large sizes, we keep the value obtained for $N = 16, 32, 64, 128$ and compare it with that for $N = 8, 16, 32, 64, 128$ to get an estimate of the error (see the inset in Fig. 7). The same procedure can be used to extract the microscopic couplings corresponding to values of $K_R \ge K_R^c$. Varying the parameters and repeating the same analysis for each set of parameters leads to the phase diagrams shown in Fig. 6. As the spin-bath coupling strength $\alpha$ is increased, the AFM gets larger while the LL shrinks, which is in agreement with Ref. [22]. As a side note, the value of the slope of the AFM-LL phase boundary at $\alpha = 0$ does not agree with the RG predictions [22], but this was to be expected as this slope is non-universal and depends on the exact short-distance cutoff used, i.e. a lattice for our simulations and a real-space hard-core repulsion between operators for the RG computation.

V. DISCUSSION AND CONCLUSION

In this work, we have proposed an enhanced-ECMC method to study bosonized 1D quantum dissipative systems. The key idea is to work not on

the original microscopic Hamiltonian, be it fermions, bosons or spins, but on its bosonized formulation which provides a unified description of 1D quantum systems [2, 3, 6]. In the path integral formalism, a bosonized field theory is that of a simple scalar field and thus lends itself to a study by Monte-Carlo algorithms usually used to study classical systems. As a test model, we have considered a dissipative XXZ spin chain analytically studied in [22, 34] which is particularly costly to simulate due to the presence of long-range interactions, minima degeneracy and critical slowing down. This led us to design an efficient algorithm combining local, persistent and continuous moves from an event-chain Monte Carlo algorithm with global discrete moves from a cluster algorithm. Both move types are relevant to address the different dynamical bottlenecks and can be implemented in $O(1)$-complexity. We performed a detailed performance analysis and showed such moves result in a combined algorithm whose integrated autocorrelation time, expressed in sweeps, scales as $O(N^{1.6})$, to be compared to $O(N^3)$ for a standard Metropolis algorithm.

Analytical studies of the XXZ dissipative spin chain had pointed out the existence of two phases: a Luttinger liquid and an antiferromagnet separated by a BKT phase transition. Thanks to the enhanced-ECMC algorithm, we were able to first confirm the analytically predicted location of the phase transition, and then obtain a precise picture of the global phase diagram.

Compared to current state-of-the-art Quantum Monte Carlo (QMC) algorithms used to study dissipative 1D quantum systems, the enhanced-ECMC algorithm stands out as being relatively straightforward to implement while being able to simulate system sizes up to $N = 1024$, which matches or outperforms recent QMC studies of 1D dissipative quantum systems [24–26]. It is worth mentioning that, although we considered here a dissipative quantum system, this method can of course be applied to any bosonized system. However, even if fermionic systems can be described through bosonization, we note that such an algorithm does not solve the sign problem in 1D as the resulting bosonic action may be imaginary, as is the case for a spin-1 chain [9].

Acknowledgements

This work was made possible by Institut Pascal at Université Paris-Saclay with the support of the
program *Investissements d'avenir* ANR-11-IDEX-0003-01. We thank Federico Ghimenti for early discussions about this project. We also thank Laura Foini and Saptarshi Majumdar for fruitful discussions about the physics of quantum dissipative systems and Alexandre Martin for great advices on code implementation choices. O.B.-D. acknowledges the

support of the French ANR under the grant ANR-22-CMAS-0001 (*QuantTEdu-France* project). M.M. acknowledges the support of the French ANR under the grant ANR-20-CE46-0007 (*SuSa* project). A.R. acknowledges the support of the French ANR under the grant ANR-23-CE30-0031-04 (*DISCREEP* project).

## Appendix A: Restricting the configuration space

The action $S(\phi)$ defined in Eq. (14), the measure of fields $d\phi$, and all physical observables (see for instance Eqs. (4,10)) are invariant under the global discrete shift symmetry $\phi \to \phi + \pi$ because the field $\phi$ is a boson compactified on a circle of radius $\pi$ [53]. The symmetry group $\{\phi \to \phi + n\pi, n \in \mathbb{Z}\}$ having an infinite volume, the partition function is clearly divergent. What we actually want to sample are the equivalence classes $[\phi] = \{\phi' \in \mathbb{R}^{N^2} | \exists m \in \mathbb{Z}, \phi' = \phi + m\pi\}$ which can be represented by any of their members $\phi$ in a restricted configuration space [75]. For instance, one can work in the space

$$ \Omega_i = \{\phi \in \mathbb{R}^{N^2} | \phi_i \in [-\pi/2, \pi/2]\}. \qquad (A1) $$

The rest of this appendix first shows in Appendix A1 how this ensures that the probability distribution $P(\phi) \propto e^{-S(\phi)}$ is normalizable, and then in Appendix A2 why one can forget about the restricted configuration space $\Omega_i$ while performing a sampling routine.

### 1. Probability distribution normalization

This subsection aims at showing that the partition function $Z = \int_{\Omega_i} d\phi e^{-S(\phi)}$ is finite, making the associated probability distribution normalizable. Without loss of generality, the proof is done for the configuration space $\Omega_0$. We start by noticing that the action $S(\phi)$ is bounded from below by

$$ \begin{aligned} S(\phi) &\ge \frac{1}{2\pi K} \sum_{j=1}^{N} \sum_{i=1}^{N-1} (\phi_{(i,j)} - \phi_{(i+1,j)})^2 + \frac{1}{2\pi K} \sum_{j=1}^{N-1} (\phi_{(0,j)} - \phi_{(0,j+1)})^2 - \frac{g}{2\pi^2} N^2 - \frac{\alpha}{2\pi^2} N^2 \sum_{j=1}^{N-1} \mathcal{D}(j) \\ &= \frac{1}{2\pi K} \sum_{j=1}^{N} \sum_{i=1}^{N-1} n_{(i,j)}^2 + \frac{1}{2\pi K} \sum_{j=1}^{N-1} m_j^2 - \frac{gN^2}{2\pi^2} - \frac{\alpha N^2}{2\pi^2} \sum_{j=1}^{N-1} \mathcal{D}(j) \end{aligned} \qquad (A2) $$

where we have introduced $N^2 - N$ space-like bond variables $n_{(i,j)} = \phi_{(i,j)} - \phi_{(i+1,j)}$ and $N-1$ time-like bond variables $m_j = \phi_{(0,j)} - \phi_{(0,j+1)}$. One can make a change of variables from the field components $\{\phi_{(i,j)}\}$ to the bond variables $\{n_{(i,j)}\}$, $\{m_j\}$ and the constrained field component $\phi_0$. Since this change of variable has a unit Jacobian,

$$ Z \le \int_{-\pi/2}^{\pi/2} d\phi_0 \left[ \prod_{j=1}^{N} \prod_{i=1}^{N-1} \int_{-\infty}^{+\infty} dn_{(i,j)} e^{-\frac{1}{2\pi K} n_{(i,j)}^2} \right] \left[ \prod_{j=1}^{N} \int_{-\infty}^{+\infty} dm_j e^{-\frac{1}{2\pi K} m_j^2} \right] \exp \left( \frac{gN^2}{2\pi^2} + \frac{\alpha N^2}{2\pi^2} \sum_{j=1}^{N-1} \mathcal{D}(j) \right) $$

which is clearly convergent.

### 2. Switching configuration spaces while sampling

During the algorithm execution, the field $\phi$ might escape the configuration space $\Omega_i = \{\phi \in \mathbb{R}^{N^2} | \phi_i \in [-\pi/2, \pi/2]\}$. This can happen during a deterministic shift of the ECMC if the lifting variable is $(\hat{i}, \pm 1)$ or
during a cluster flip if the site $\hat{i}$ gets flipped. Following Ref. [54], a possible fix for the ECMC part would be to add a boundary Markov kernel to trigger a global shift $\phi \to \phi \pm \pi$ whenever the boundary of $\Omega_{\hat{i}}$ is reached. For the cluster part, a similar solution can be implemented by globally shifting the field as $\phi \to \phi + n\pi$ after each cluster flip with $n \in \mathbb{Z}$ chosen to maintain $\phi$ in $\Omega_{\hat{i}}$. However, because all observables are invariant under such a shift, we decide not to implement such eventually unnecessary boundary jumps.

## Appendix B: Global balance of the ECMC algorithm

A PDMP described by a generator $\mathcal{A}$ samples the probability distribution $\pi(\phi)$ if it is its stationary (or invariant) distribution, i.e.

$$ \sum_{v \in \mathcal{V}} \frac{1}{2N^2} \int_{\Omega_i} \mathcal{A} f(\phi, v) \pi(\phi) d\phi = 0. \qquad (B1) $$

In this appendix, we show, following the lines of [54], that the probability distribution $\pi(\phi) \propto e^{-S(\phi)}$ is indeed a stationary distribution of the generator

$$ \begin{aligned} \mathcal{A} f(\phi, (e, i)) = & \partial_{\phi_i} f(\phi, (e, i))e \\ & + \sum_{j \in \partial_{\phi_i}} [e \partial_{\phi_i} S_{\phi}^{i,j}(\phi)]_+ (f(\phi, (e, j)) - f(\phi, (e, i))) \\ & + \sum_{j \in \partial_{\phi_i}} [e \partial_{\phi_i} S_{\phi}^{i,j}(\phi)]_+ \sum_{k \in \{i, j\}, e'} \frac{[(-e') \partial_{\phi_k} S_{\phi}^{i,j}(\phi)]_+ (f(\phi, (e', j)) - f(\phi, (e, i)))}{\sum_{l \in \{i, j\}, e''} [(-e'') \partial_{\phi_l} S_{\phi}^{i,j}(\phi)]_+} \end{aligned} $$

where all sums with $e, e', e''$ are over $\{-1, 1\}$, $\partial_{\phi/c} i$ denotes the set of sites connected to $i$ through a quadratic/cosine interaction, and the remaining notations have been introduced in Sec. III A. We also introduce for later use the total quadratic and cosine interactions

$$ S_{\phi/c}(\phi) = \sum_{\langle i, j \rangle_{\phi/c}} S_{\phi/c}^{i, j}(\phi), \qquad (B2) $$

where $\langle i, j \rangle_{\phi/c}$ denotes pairs of sites connected through a quadratic/cosine interaction. Plugging the generator (B2) into Eq. (B1), one must show that the following expression vanishes

$$ \begin{aligned} & (\star) \sum_{v \in \mathcal{V}} \frac{1}{2N^2} \int_{\Omega_i} d\phi \pi(\phi) \partial_{\phi_i} f(\phi, (e, i))e \\ & (\star\star) + \sum_{v \in \mathcal{V}} \frac{1}{2N^2} \int_{\Omega_i} d\phi \pi(\phi) \sum_{j \in \partial_{\phi_i}} [e \partial_{\phi_i} S_{\phi}^{i,j}(\phi)]_+ (f(\phi, (e, j)) - f(\phi, (e, i))) \\ & (\star\star\star) + \sum_{v \in \mathcal{V}} \frac{1}{2N^2} \int_{\Omega_i} d\phi \pi(\phi) \sum_{j \in \partial_{\phi_i}} [e \partial_{\phi_i} S_{\phi}^{i,j}(\phi)]_+ \sum_{k \in \{i, j\}, e'} \frac{[(-e') \partial_{\phi_k} S_{\phi}^{i,j}(\phi)]_+ (f(\phi, (e', k)) - f(\phi, (e, i)))}{\sum_{l \in \{i, j\}, e''} [(-e'') \partial_{\phi_l} S_{\phi}^{i,j}(\phi)]_+} \end{aligned} \qquad (B3) $$
Let us begin with the quadratic term (★). Using the pairwise symmetry $[e\partial_{\phi_i} S_q^{i,j}]_+ = [(-e)\partial_{\phi_j} S_q^{i,j}]_+$ leads to

$$
\begin{align*}
(\star\star) &= \sum_{e,i} \frac{1}{2N^2} \int_{\Omega_i} d\phi \, \pi(\phi) \sum_{j \in \partial_q i} [e\partial_{\phi_i} S_q^{i,j}(\phi)]_+ (f(\phi, (e, j)) - f(\phi, (e, i))) \\
&= \sum_e \frac{1}{2N^2} \int_{\Omega_i} d\phi \, \pi(\phi) \left[ \sum_{j, i \in \partial_q j} [(-e)\partial_{\phi_j} S_q^{i,j}(\phi)]_+ f(\phi, (e, j)) - \sum_{i, j \in \partial_q} [e\partial_{\phi_i} S_q^{i,j}(\phi)]_+ f(\phi, (e, i)) \right] \\
&= \sum_{e,i} \frac{1}{2N^2} \int_{\Omega_i} d\phi \, \pi(\phi) f(\phi, (e, i)) \sum_{j \in \partial_q} [(-e)\partial_{\phi_i} S_q^{i,j}(\phi)]_+ - [e\partial_{\phi_i} S_q^{i,j}(\phi)]_+ \\
&= - \sum_{e,i} \frac{1}{2N^2} \int_{\Omega_i} d\phi \, \pi(\phi) f(\phi, (e, i)) \sum_{j \in \partial_q} e\partial_{\phi_i} S_q^{i,j}(\phi) \\
&= - \sum_{e,i} \frac{1}{2N^2} \int_{\Omega_i} d\phi \, \pi(\phi) f(\phi, (e, i)) e\partial_{\phi_i} S_q(\phi), \tag{B4}
\end{align*}
$$

where we recognized the total quadratic action $S_q(\phi)$ in the last line. The cosine term (★★★) is then computed by first ordering the sums as

$$
(\star\star\star) = \int_{\Omega_i} d\phi \, \pi(\phi) \frac{1}{2N^2} \sum_{e,i} \sum_{j \in \partial_c i} \sum_{k \in \{i,j\}, e'} [e\partial_{\phi_i} S_c^{i,j}(\phi)]_+ \frac{[(-e')\partial_{\phi_k} S_c^{i,j}(\phi)]_+ (f(\phi, (e', k)) - f(\phi, (e, i)))}{\sum_{l \in \{i,j\}, e''} [(-e'')\partial_{\phi_l} S_c^{i,j}(\phi)]_+} \tag{B5}
$$

Next, using the identity $\sum_i \sum_{j \in \partial_c i} g(i) = \sum_{\langle i,j \rangle_c} \sum_{m \in \{i,j\}} g(m)$, (★★★) becomes

$$
\begin{align*}
(\star\star\star) &= \int_{\Omega_i} d\phi \, \pi(\phi) \left[ \frac{1}{2N^2} \sum_{e,e'} \sum_{\langle i,j \rangle_c} \sum_{m \in \{i,j\}} \sum_{k \in \{i,j\}} [e\partial_{\phi_m} S_c^{i,j}(\phi)]_+ + \frac{[(-e')\partial_{\phi_k} S_c^{i,j}(\phi)]_+ f(\phi, (e', k))}{\sum_{l \in \{i,j\}, e''} [(-e'')\partial_{\phi_l} S_c^{i,j}(\phi)]_+} \right. \\
&\qquad \left. - \frac{1}{2N^2} \sum_{e,i} f(\phi, (e,i)) \sum_{j \in \partial_c i} [e\partial_{\phi_i} S_c^{i,j}(\phi)]_+ \sum_{k \in \{i,j\}, e'} \frac{[(-e')\partial_{\phi_k} S_c^{i,j}(\phi)]_+}{\sum_{l \in \{i,j\}, e''} [(-e'')\partial_{\phi_l} S_c^{i,j}(\phi)]_+} \right] \\
&= \int_{\Omega_i} d\phi \, \pi(\phi) \frac{1}{2N^2} \left[ \sum_{e'} \sum_{\langle i,j \rangle_c} \sum_{k \in \{i,j\}} f(\phi, (e', k)) [(-e')\partial_{\phi_k} S_c^{i,j}(\phi)]_+ \right. \\
&\qquad \left. - \sum_{e,i} f(\phi, (e,i)) \sum_{j \in \partial_c i} [e\partial_{\phi_i} S_c^{i,j}(\phi)]_+ \right] \\
&= \int_{\Omega_i} d\phi \, \pi(\phi) \frac{1}{2N^2} \left[ \sum_{e',k} f(\phi, (e', k)) \sum_{j \in \partial_c k} [(-e')\partial_{\phi_k} S_c^{k,j}(\phi)]_+ - \sum_{e,i} f(\phi, (e,i)) \sum_{j \in \partial_c i} [e\partial_{\phi_i} S_c^{i,j}(\phi)]_+ \right] \\
&= \int_{\Omega_i} d\phi \, \pi(\phi) \frac{1}{2N^2} \sum_{e,i} f(\phi, (e,i)) \sum_{j \in \partial_c i} [(-e)\partial_{\phi_i} S_c^{i,j}(\phi)]_+ - [e\partial_{\phi_i} S_c^{i,j}(\phi)]_+ \\
&= - \int_{\Omega_i} d\phi \, \pi(\phi) \frac{1}{2N^2} \sum_{e,i} f(\phi, (e,i)) \sum_{j \in \partial_c i} e\partial_{\phi_i} S_c^{i,j}(\phi) \\
&= - \sum_{e,i} \frac{1}{2N^2} \int_{\Omega_i} d\phi \, \pi(\phi) f(\phi, (e,i)) e\partial_{\phi_i} S_c(\phi). \tag{B6}
\end{align*}
$$
<figure><img src="image_9.png" /><figcaption>Figure 8. Integrated autocorrelation time $\tau_{int}$ for the staggered magnetization $m$ with respect to the parameter $p$ defined in (C1) and with the Clu-EC algorithm. The error bars were obtained by averaging over 40 runs with independent gaussian initial conditions $\phi_i \stackrel{i.i.d.}{\sim} \mathcal{N}(0, 1)$.</figcaption></figure>

Putting everything together yields

$$ \sum_{v \in V} \frac{1}{2N^2} \int_{\Omega_i} \mathcal{A} f(\phi, v) \pi(\phi) d\phi = \sum_{e, i} \frac{1}{2N^2} \int_{\Omega_i} d\phi \pi(\phi) e \left[ \partial_{\phi_i} f(\phi, (e, i)) - f(\phi, (e, i)) \partial_{\phi_i} (S_q(\phi) + S_c(\phi)) \right]. \quad (B7) $$

Finally, upon using $\pi(\phi) \propto e^{-(S_q+S_c)}$ one arrives at

$$ \begin{aligned} \sum_{v \in V} \frac{1}{2N^2} \int_{\Omega_i} \mathcal{A} f(\phi, v) \pi(\phi) \mu(v) d\phi dv &= \sum_{e, i} \frac{1}{2N^2} \int_{\Omega_i} d\phi e \left[ \pi(\phi) \partial_{\phi_i} f(\phi, (e, i)) + f(\phi, (e, i)) \partial_{\phi_i} \pi(\phi) \right] \\ &= \sum_{e, i} \frac{e}{2N^2} \int_{\Omega_i} d\phi \partial_{\phi_i} \left[ \pi(\phi) f(\phi, (e, i)) \right] = 0, \end{aligned} \quad (B8) $$

which finishes the proof.

## Appendix C: Choice of the cluster reflections

When performing cluster moves as defined in Sect. IIIB, one has the freedom to choose the reflection operator $R_n$. We argued that efficient moves are those which do not flip a field component $\phi_i$ too far from its original position. Starting a cluster from an initial site $i_0$, we thus propose to uniformly pick a reflection $R_n$ among those with an index

$$ n \in [2m-p, 2m-p+1, \dots, 2m+p-1, 2m+p], \quad (C1) $$

where $m = \lfloor \frac{2}{\pi}(\phi_{i_0} + \frac{\pi}{4}) \rfloor$ is such that $m\pi/2$ is the potential minima closest to $\phi_{i_0}$. To optimize the parameter $p$, we performed a series of runs at the phase transition for different values of $p$ and plotted the integrated autocorrelation time $\tau_{int}$ of the staggered magnetization $m$ in Fig. 8. The value $p=2$ seems to be the better one. Noticeably, this includes odd $(R_{2n+1})$ and even $(R_{2n})$ reflections, which was identified as a crucial performance factor for cluster algorithms applied to the sine-Gordon model [47] or the discrete Gaussian model [76].
Appendix D: Detailed balance of the cluster algorithm

The cluster algorithm described in Sec. IIIB samples the correct distribution $\pi(\phi) \propto e^{-S(\phi)}$ if it satisfies the detailed balance condition between two field configurations $\phi$ and $\phi'$

$$P_{\phi \to \phi'} A_{\phi \to \phi'} \pi(\phi) = P_{\phi' \to \phi} A_{\phi' \to \phi} \pi(\phi'), \quad (D1)$$

where $A_{\phi \to \phi'}$ is the probability of proposing the move $\phi \to \phi'$ and $P_{\phi \to \phi'}$ is that of accepting it. Since the cluster algorithm is rejection-free, Eq. (D1) should be satisfied with $P_{\phi \to \phi'} = P_{\phi' \to \phi} = 1$. To prove this, let us first consider the probability of adding site $j$ to the cluster while looking at site $i$ with the clock implementation. For a short-range neighbor, it is by definition given by $p_{i,j}(\phi)$ (see Eq. (39)). For a long-range neighbor, we explicit the Fukui-Todo implementation [67]. A long-range neighbor is added by creating $n$ bound events from $i$ with probability $f(n, \lambda_{\text{Clu}}^B) = e^{-\lambda_{\text{Clu}}^B} \lambda_{\text{Clu}}^{n/n!}$, attributing them to $i + k\hat{\tau}$ with probability $\lambda_k^B/\lambda_{\text{Clu}}^B$, and then accepting it with probability $p_{i,i+k\hat{\tau}}(\phi)/p_k^B$ if there is at least one bound event on $i + k\hat{\tau}$. This defines the probability $p'_{i,i+k\hat{\tau}}(\phi)$ of adding a long-range neighbor as

$$\begin{aligned}
p'_{i,i+k\hat{\tau}}(\phi) &= \sum_{n=1}^{+\infty} \frac{e^{-\lambda_{\text{Clu}}^B} \lambda_{\text{Clu}}^{n/n!}}{n!} \sum_{n_k=1}^{n} \binom{n}{n_k} \left(\frac{\lambda_k^B}{\lambda_{\text{Clu}}^B}\right)^{n_k} \left(1 - \frac{\lambda_k^B}{\lambda_{\text{Clu}}^B}\right)^{n-n_k} \frac{p_{i,i+k\hat{\tau}}(\phi)}{p_k^B} \\
&= \left(1 - e^{-\lambda_{\text{Clu}}^B}\right) \frac{p_{i,i+k\hat{\tau}}(\phi)}{p_k^B} \\
&= p_{i,i+k\hat{\tau}}(\phi),
\end{aligned} \quad (D2)$$

which shows that short-range and long-range neighbors are both added with probability $p_{i,j}(\phi)$.

Keeping the previous result in mind, a cluster $C$ that relates $\phi \to \phi'$ is built by, first, selecting the correct reflection $R_n$, and then generating a tree graph $G(C) = (C, E(C))$ that originates from an initial site $i_0 \in C$ and spans $C$. If one thinks of adding site $j$ to the cluster from site $i$ as activating the edge $\langle i, j \rangle$, the tree $G(C)$ is such that each edge $\langle i, j \rangle \in E(C)$ in this graph has to be accepted whereas all the edges on its border $\partial C = \{\langle i, j \rangle \mid j \in \partial i, i \in C, j \notin C\}$ have to be rejected. Formally, this means that

$$A(\phi \to \phi') = \sum_{i_0 \in C} P(i_0 \text{ is the initial node}) \sum_n P(R_n \text{ is the correct reflection } | i_0) \\ \times \sum_{E(C)} \prod_{\langle i,j \rangle \in E(C)} p_{i,j}(\phi) \prod_{\langle i,j \rangle \in \partial C} (1 - p_{i,j}(\phi)). \quad (D3)$$

The initial node being chosen uniformly over the lattice of size $N^2$, $P(i_0)$ is the initial node) = $1/N^2$. Furthermore, since the reflection index $n$ is conditioned on the initial node $i_0$, it can be written as $n = 2m + p$ with $m = \lfloor \frac{2}{\pi}\phi_{i_0} + \frac{1}{2} \rfloor$ and $p \in \{-2, -1, 0, 1, 2\}$. Therefore, $R_n$ is a reflection connecting $\phi$ and $\phi'$ if

$$n \frac{\pi}{2} - \phi_{i_0} = \phi'_{i_0} \Rightarrow p = \frac{2}{\pi} \phi_{i_0} + \frac{2}{\pi} \phi'_{i_0} - 2m. \quad (D4)$$

Writing $\phi_{i_0} = \frac{\pi}{2}(m+x)$ with $m = \lfloor \frac{2}{\pi}\phi_{i_0} + \frac{1}{2} \rfloor \in \mathbb{Z}$ and $|x| < \frac{1}{2}$ (resp. $\phi'_{i_0} = \frac{\pi}{2}(m'+x')$ with $m' = \lfloor \frac{2}{\pi}\phi_{i_0} + \frac{1}{2} \rfloor \in \mathbb{Z}$ and $|x'| < \frac{1}{2}$) leads to a reflection being valid if

$$p = x + x' + m' - m \Rightarrow \begin{cases} p = m' - m, \\ x + x' = 0, \end{cases} \quad (D5)$$

which, because $p$ is uniformly chosen in $\{-2, -1, 0, 1, 2\}$, implies

$$\sum_n P(R_n \text{ is the correct reflection } | i_0) = \frac{1}{5} P(x + x' = 0) \mathbf{1}_{\{-2, \dots, 2\}}(m' - m). \quad (D6)$$
The proposal probability $A(\phi \to \phi')$ is thus rewritten as

$$A(\phi \to \phi') = \sum_{i_0 \in C} \frac{1}{5N^2} P(x + x' = 0) \mathbf{1}_{\{-2, \dots, 2\}}(m' - m) \sum_{E(C)} \prod_{\langle i,j \rangle \in E(C)} p_{i,j}(\phi) \prod_{\langle i,j \rangle \in \partial C} (1 - p_{i,j}(\phi)). \quad (D7)$$

Putting everything together and using $p_{i,j}(\phi) = p_{i,j}(R_n\phi)$ leads to

$$\begin{aligned} \frac{A(\phi \to \phi')}{A(\phi' \to \phi)} &= \frac{\sum_{i_0 \in C} P(x + x' = 0) \mathbf{1}_{\{-2, \dots, 2\}}(m' - m) \sum_{E(C)} \prod_{\langle i,j \rangle \in E(C)} p_{i,j}(\phi) \prod_{\langle i,j \rangle \in \partial C} (1 - p_{i,j}(\phi))}{\sum_{i_0 \in C} P(x' + x = 0) \mathbf{1}_{\{-2, \dots, 2\}}(m - m') \sum_{E(C)} \prod_{\langle i,j \rangle \in E(C)} p_{i,j}(\phi') \prod_{\langle i,j \rangle \in \partial C} (1 - p_{i,j}(\phi'))} \\ &= \prod_{\langle i,j \rangle \in \partial C} \frac{1 - p_{i,j}(\phi)}{1 - p_{i,j}(\phi')}. \end{aligned} \quad (D8)$$

From the definition $p_{i,j}(\phi') = 1 - \exp(-[\Delta S_{\text{pair}}^{R_n}(\phi_i, \phi_j)]_+)$ and the fact that $\forall j \notin C, \phi_j = \phi'_j$, one arrives at

$$\begin{aligned} \frac{A(\phi \to \phi')}{A(\phi' \to \phi)} &= \prod_{\langle i,j \rangle \in \partial C} \frac{\exp(-[\Delta S_{\text{pair}}^{R_n}(\phi_i, \phi_j)]_+)}{\exp(-[\Delta S_{\text{pair}}^{R_n}(\phi'_i, \phi'_j)]_+)} \\ &= \exp\left(-\sum_{\langle i,j \rangle \in \partial C} \Delta S_{\text{pair}}^{R_n}(\phi_i, \phi_j)\right) \\ &= \frac{\pi(\phi')}{\pi(\phi)}, \end{aligned} \quad (D9)$$

which is the detailed balance condition (D1).

[1] F. D. M. Haldane, 'Luttinger liquid theory' of one-dimensional quantum fluids. I. Properties of the Luttinger model and their extension to the general 1D interacting spinless Fermi gas, Journal of Physics C: Solid State Physics 14, 2585 (1981).

[2] T. Giamarchi, Quantum physics in one dimension, International series of monographs on physics (Clarendon Press, Oxford, 2004).

[3] J. von Delft and H. Schoeller, Bosonization for beginners — refermionization for experts, Annalen der Physik 510, 225 (1998).

[4] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol, One dimensional bosons: From condensed matter systems to ultracold gases, Rev. Mod. Phys. 83, 1405 (2011).

[5] D. Sénéchal, A.-M. Tremblay, and C. Bourbonnais, Theoretical Methods for Strongly Correlated Electrons (Springer, New York, 2003).

[6] N. Dupuis, Field Theory of Condensed Matter and Ultracold Gases, Vol. 2 (unpublished) (WORLD SCIENTIFIC (EUROPE), 2025) see https://www.lptmc.jussieu.fr/users/dupuis.

[7] A. Altland and B. D. Simons, Condensed matter field theory (Cambridge university press, 2010).

[8] T. Giamarchi and H. J. Schulz, Anderson localization and interactions in one-dimensional metals, Phys. Rev. B 37, 325 (1988).

[9] H. J. Schulz, Phase diagrams and correlation exponents for quantum spin chains of arbitrary spin quantum number, Phys. Rev. B 34, 6372 (1986).

[10] T. Giamarchi, Mott transition in one dimension, Physica B: Condensed Matter 230-232, 975 (1997), proceedings of the International Conference on Strongly Correlated Electron Systems.

[11] R. Daviet and N. Dupuis, Flowing bosonization in the nonperturbative functional renormalization-group approach, SciPost Phys. 12, 110 (2022).

[12] A. O. Caldeira and A. J. Leggett, Influence of Dissipation on Quantum Tunneling in Macroscopic Systems, Phys. Rev. Lett. 46, 211 (1981).

[13] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Rev. Mod. Phys. 59, 1 (1987).

[14] A. Winter, H. Rieger, M. Vojta, and R. Bulla, Quantum Phase Transition in the Sub-Ohmic Spin-Boson
Model: Quantum Monte Carlo Study with a Continuous Imaginary Time Cluster Algorithm, Phys. Rev. Lett. 102, 030601 (2009).

[15] G. De Filippis, A. de Candia, L. M. Cangemi, M. Sassetti, R. Fazio, and V. Cataudella, Quantum phase transitions in the spin-boson model: Monte Carlo method versus variational approach à la Feynman, Phys. Rev. B 101, 180408 (2020).

[16] A. Schmid, Diffusion and localization in a dissipative quantum system, Phys. Rev. Lett. 51, 1506 (1983).

[17] P. Werner and M. Troyer, Efficient simulation of resistively shunted Josephson junctions, Phys. Rev. Lett. 95, 060201 (2005).

[18] N. Paris, L. Giacomelli, R. Daviet, C. Ciuti, N. Dupuis, and C. Mora, Resilience of the quantum critical line in the schmid transition, Phys. Rev. B 111, 064509 (2025).

[19] R. Citro, E. Orignac, and T. Giamarchi, Adiabatic-antiadiabatic crossover in a spin-Parierls chain, Phys. Rev. B 72, 024434 (2005).

[20] A. M. Lobos, A. Iucci, M. Müller, and T. Giamarchi, Dissipation-driven phase transitions in superconducting wires, Phys. Rev. B 80, 214515 (2009).

[21] E. Malatsetxebarria, Z. Cai, U. Schollwöck, and M. A. Cazalilla, Dissipative effects on the superfluid-to-insulator transition in mixed-dimensional optical lattices, Physical Review A 88, 10.1103/physrev.a.88.063630 (2013).

[22] O. Bouverot-Dupuis, S. Majumdar, A. Rosso, and L. Foini, Antiferromagnetic order enhanced by local dissipation, Phys. Rev. B 109, 205148 (2024).

[23] A. M. Lobos, A. Iucci, M. Müller, and T. Giamarchi, Dissipation-driven phase transitions in superconducting wires, Phys. Rev. B 80, 214515 (2009).

[24] M. Weber, F. F. Assaad, and M. Hohenadler, Directed-loop quantum monte carlo method for retarded interactions, Phys. Rev. Lett. 119, 097401 (2017).

[25] M. Weber, D. J. Luitz, and F. F. Assaad, Dissipation-induced order: The $s = 1/2$ quantum spin chain coupled to an ohmic bath, Phys. Rev. Lett. 129, 056402 (2022).

[26] A. L. S. Ribeiro, P. McClarty, P. Ribeiro, and M. Weber, Dissipation-induced long-range order in the one-dimensional Bose-Hubbard model, Phys. Rev. B 110, 115145 (2024).

[27] N. Laflorencie, G. Lemarié, and N. Macé, Chain breaking and kosterlitz-thouless scaling at the many-body localization transition in the random-field heisenberg spin chain, Physical Review Research 2, 10.1103/physrevresearch.2.042033 (2020).

[28] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).

[29] U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).

[30] N. Kawashima and K. Harada, Recent Developments of World-Line Monte Carlo Methods, Journal of the Physical Society of Japan 73, 1379 (2004), https://doi.org/10.1143/JPSJ.73.1379.

[31] K. Van Houcke, E. Kozik, N. Prokof'ev, and B. Svis tunov, Diagrammatic monte carlo, Physics Procedia 6, 95 (2010), computer Simulations Studies in Condensed Matter Physics XXI.

[32] R. G. Edwards and A. D. Sokal, Generalization of the fortuin-kasteleyn-swendsen-wang representation and monte carlo algorithm, Phys. Rev. D 38, 2009 (1988).

[33] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equation of state calculations by fast computing machines, The Journal of Chemical Physics 21, 1087 (1953), https://pubs.aip.org/aip/jcp/article-pdf/21/6/1087/18802390/1087_1_online.pdf.

[34] O. Bouverot-Dupuis, Mapping a dissipative quantum spin chain onto a generalized Coulomb gas, SciPost Phys. 17, 130 (2024).

[35] M. Michel, S. C. Kapfer, and W. Krauth, Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps, The Journal of Chemical Physics 140, 054116 (2014).

[36] M. Michel, X. Tan, and Y. Deng, Clock monte carlo methods, Phys. Rev. E 99, 010105 (2019).

[37] E. P. Bernard and W. Krauth, Two-Step Melting in Two Dimensions: First-Order Liquid-Hexatic Transition, Phys. Rev. Lett. 107, 155704 (2011).

[38] S. C. Kapfer and W. Krauth, Two-Dimensional Melting: From Liquid-Hexatic Coexistence to Continuous Transitions, Phys. Rev. Lett. 114, 035702 (2015).

[39] M. Isobe and W. Krauth, Hard-sphere melting and crystallization with event-chain Monte Carlo, The Journal of Chemical Physics 143, 084509 (2015).

[40] M. Klement, S. Lee, J. A. Anderson, and M. Engel, Newtonian Event-Chain Monte Carlo and Colli sion Prediction with Polyhedral Particles, Journal of Chemical Theory and Computation 17, 4686 (2021), https://doi.org/10.1021/acs.jctc.1c00311.

[41] M. Michel, J. Mayer, and W. Krauth, Event-chain monte carlo for classical continuous spin models, EPL (Europhysics Letters) 112, 20003 (2015).

[42] Y. Nishikawa, M. Michel, W. Krauth, and K. Hukushima, Event-chain algorithm for the Heisenberg model: Evidence for $z \simeq 1$ dynamic scaling, Phys. Rev. E 92, 063306 (2015).

[43] T. A. Kampmann, H.-H. Boltz, and J. Kierfeld, Monte carlo simulation of dense polymer melts using event chain algorithms, The Journal of Chemical Physics 143, 044105 (2015).

[44] J. Harland, M. Michel, T. A. Kampmann, and J. Kierfeld, Event-chain Monte Carlo algorithms for three-and many-particle interactions, EPL (Europhysics Letters) 117, 30001 (2017).

[45] F. Ghimenti, L. Berthier, and F. van Wijland, Irreversible monte carlo algorithms for hard disk glasses: From event-chain to collective swaps, Phys. Rev. of the Physical Society of Japan 73, 1379 (2004), https://doi.org/10.1143/JPSJ.73.1379.
Lett. 133, 028202 (2024).

[46] R. C. Brower and P. Tamayo, Embedded dynamics for $\phi^4$ theory, Phys. Rev. Lett. 62, 1087 (1989).

[47] M. Hasenbusch, M. Marcu, and K. Pinn, The sine gordon model: perturbation theory and cluster monte carlo, Physica A: Statistical Mechanics and its Applications 211, 255 (1994).

[48] https://plmlab.math.cnrs.fr/stoch-algo-phys/ecmc/enhanced-ecmc-bosonized.

[49] S. Majumdar, L. Foini, T. Giamarchi, and A. Rosso, Bath-induced phase transition in a Luttinger liquid, Phys. Rev. B 107, 165113 (2023).

[50] S. Majumdar, L. Foini, T. Giamarchi, and A. Rosso, Localization induced by spatially uncorrelated sub-ohmic baths in one dimension, Phys. Rev. B 108, 205138 (2023).

[51] This comes as no surprise since the isolated spin chain (i.e. with no bath) already exhibits a transition to an antiferromagnetic phase.

[52] U. Weiss, *Quantum Dissipative Systems*, 4th ed. (WORLD SCIENTIFIC, 2012) https://www.worldscientific.com/doi/pdf/10.1142/8334.

[53] D. Tong, Lecture Notes on Gauge Theory (2018).

[54] A. Monemvassitis, A. Guillin, and M. Michel, PDMP characterisation of event-chain monte carlo algorithms for particle systems, Journal of Statistical Physics 190, 66 (2023).

[55] P. Diaconis, S. Holmes, and R. Neal, Analysis of a nonreversible markov chain sampler, Ann. Appl. Probab. 10 (2000).

[56] M. Michel, A. Durmus, and S. Sénécal, Forward Event-Chain Monte Carlo: Fast Sampling by Randomness Control in Irreversible Markov Chains, Journal of Computational and Graphical Statistics 29, 689 (2020), https://doi.org/10.1080/10618600.2020.1750417.

[57] T. Guyon, A. Guillin, and M. Michel, Necessary and sufficient symmetries in Event-Chain Monte Carlo with generalized flows and application to hard dimers, The Journal of Chemical Physics 160, 024117 (2024).

[58] M. H. A. Davis, *Markov models and optimization*, Monographs on Statistics and Applied Probability, Vol. 49 (Chapman & Hall, London, 1993) pp. xiv+295.

[59] Z. Lei, W. Krauth, and A. C. Maggs, Event-chain monte carlo with factor fields, Phys. Rev. E 99, 043301 (2019).

[60] J. Lu and L. Wang, On explicit $L^2$-convergence rate estimate for piecewise deterministic Markov processes in MCMC algorithms, The Annals of Applied Probability 32, 1333 (2022).

[61] A. Eberle, A. Guillin, L. Hahn, F. Lörler, and M. Michel, Convergence of non-reversible Markov

processes via lifting and flow Poincaré inequality (2025), arXiv:2503.04238 [math.AP].

[62] P. A. W. Lewis and G. S. Shedler, Simulation of non-homogeneous poisson processes by thinning, Naval Research Logistics Quarterly 26, 403 (1979).

[63] R. H. Swendsen and J.-S. Wang, Nonuniversal critical dynamics in monte carlo simulations, Phys. Rev. Lett. 58, 86 (1987).

[64] N. Prokof'ev and B. Svistunov, Worm algorithms for classical statistical models, Phys. Rev. Lett. 87, 160601 (2001).

[65] U. Wolff, Collective monte carlo updating for spin systems, Phys. Rev. Lett. 62, 361 (1989).

[66] E. Luijten and H. Blöte, Monte Carlo method for spin models with long-range interactions, Int. J. Mod. Phys C6, 359 (1995).

[67] K. Fukui and S. Todo, Order-n cluster monte carlo method for spin systems with long-range interactions, Journal of Computational Physics 228, 2629 (2009).

[68] E. Flores-Sola, M. Weigel, R. Kenna, and B. Berche, Cluster Monte Carlo and dynamical scaling for long-range interactions, The European Physical Journal Special Topics 226, 581-594 (2017).

[69] A. J. Walker, An efficient method for generating discrete random variables with general distributions, ACM Trans. Math. Softw. 3, 253-256 (1977).

[70] G. Marsaglia, W. W. Tsang, and J. Wang, Fast generation of discrete random variables, Journal of Statistical Software 11, 1-11 (2004).

[71] A. Bouchard-Côté, S. J. Vollmer, and A. Doucet, The bouncy particle sampler: A nonreversible rejection-free markov chain monte carlo method, Journal of the American Statistical Association 113, 855 (2018).

[72] H. Weber and P. Minnhagen, Monte carlo determination of the critical temperature for the two-dimensional xy model, Phys. Rev. B 37, 5986 (1988).

[73] K. Harada and N. Kawashima. Universal jump in the helicity modulus of the two-dimensional quantum xy model, Phys. Rev. B 55, R11949 (1997).

[74] In fact, in one of the first Monte-Carlo studies of a BKT transition [72], the finite-size scaling relation was shown to be valid down to sizes $N = 3$, enabling a precise determination of the critical point from the system sizes $N = 3, 4, 5, 6, 7, 8$.

[75] This is exactly like fixing a gauge in electromagnetism to get rid of the gauge symmetry redundancy.

[76] H. G. Evertz, M. Hasenbusch, M. Marcu, K. Pinn, and S. Solomon, Stochastic cluster algorithms for discrete gaussian (sos) models, Physics Letters B 254, 185 (1991).