Dissipation-induced phase transitions in a frustrated four-spin plaquette spin-boson model

Abstract

We investigate a frustrated four-spin plaquette spin-boson model with competing nearest-neighbor and diagonal Ising couplings, where each spin is coupled to an independent bosonic bath. Combining a path-integral strong-coupling analysis with variational matrix product state simulations, we obtain the ground-state phase diagram. In the strong-dissipation limit we map the model onto a classical plaquette and derive analytic phase boundaries between ferromagnetic, Néel, and stripe phases. At intermediate dissipation we find a delocalized phase and two localized ordered phases with Néel and stripe character. We show that the localized phases are separated by a first-order line, while each is connected to the delocalized regime via a continuous (second-order) localization transition, and that these three boundaries meet at a quantum triple point. Analysis of spin correlations and reduced density matrices further reveals that entanglement concentrates on nearest-neighbor (diagonal) bonds in the Néel (stripe) phase, whereas in the delocalized regime intra-plaquette two-spin entanglement is strongly suppressed in favor of enhanced spin-bath correlations.

Graphical abstract

Keywords

spin-boson model / open quantum systems / dissipation-induced phase transition / quantum phase transitions / frustration / variational matrix product states / strong-coupling

1 Introduction

Quantum phase transitions (QPTs), driven by quantum rather than thermal fluctuations, provide a powerful framework for understanding how correlations and entanglement reorganize at zero temperature in strongly interacting systems [1–3]. In particular, light-matter platforms and open quantum impurity models have emerged as versatile arenas where strong coupling, geometry, and dissipation can be engineered and probed with high controllability [4, 20].

A paradigmatic example is the spin-boson model (SBM), which describes a two-level system coupled to a structured bosonic bath and captures the competition between coherent tunneling and environment-induced localization [5]. In its sub-Ohmic variants, the SBM exhibits a continuous QPT associated with the spontaneous breaking of a global {\mathbb{Z}}_2 symmetry, and has played a central role in the theory of dissipation in quantum systems. Multi-spin generalizations of the SBM naturally introduce competing spin-spin interactions on top of local spin-bath couplings, opening the door to collective ordering phenomena, frustration, and symmetry breaking in open many-body settings [6–9].

In parallel, few-body light−matter models based on the quantum Rabi and Dicke paradigms have been shown to emulate rich magnetic behavior when arranged in ring or plaquette geometries and subjected to photon hopping. A notable example is the quantum Rabi ring threaded by an artificial magnetic field, which can be mapped onto an effective spin model with XY exchange and Dzyaloshinskii−Moriya (DM) interactions [10]. By tuning the gauge phase of the photon hopping, this system realizes ferro-, antiferro-, and chiral magnetic phases in the light sector, and reveals superradiant transitions whose critical properties can be understood in terms of competing magnetic couplings and geometric frustration on the ring.

Building on this analogy, a quantum Rabi square with both nearest- and next-nearest-neighbor photon hopping has been proposed as a minimal optical realization of a frustrated J_1−J_2 spin model on a plaquette geometry [11]. In this model, the relative strength of the diagonal hopping plays a role analogous to that of a gauge phase, driving first-order transitions between antiferromagnetic and frustrated superradiant phases, in addition to the usual second-order transition between normal and superradiant regimes. This demonstrates that long-range cavity coupling can serve as an alternative to synthetic gauge fields for engineering competing orders and frustration in few-body light−matter systems.

These advances naturally raise the question of how analogous mechanisms manifest in dissipative spin-boson settings, where each spin is strongly coupled to its own continuum bath rather than to coherent cavity modes. In particular, it is not a priori clear whether frustration and competing orders can be stabilized or even enhanced when bath-induced localization suppresses tunneling, nor how the presence of independent sub-Ohmic environments modifies the balance between different ordering tendencies familiar from frustrated spin models. A square plaquette of four spins provides the simplest geometry in which nearest-neighbor and diagonal Ising couplings compete on equal footing, giving rise to ferromagnetic, Néel (checkerboard), and stripe-like configurations in the classical J_1−J_2 limit. Coupling each spin to an independent bosonic bath introduces an additional competition between spatial ordering and environment-induced localization, and offers a minimal setting in which to investigate how dissipation-assisted frustration emerges in an open quantum system.

In this work, we study a four-spin spin-boson model on a square plaquette, where the spins interact via nearest-neighbor and diagonal Ising-type couplings and are each coupled to an independent sub-Ohmic bosonic bath. In the strong-dissipation regime, integrating out the baths maps the model onto an effective classical J_1−J_2 plaquette, whose ground states realize ferromagnetic, Néel, and stripe orders with analytically accessible phase boundaries. Using the Variational Matrix Product State (VMPS) method, we go beyond this limit to obtain the full quantum phase diagram as a function of the spin-spin couplings and bath strength, and demonstrate how the numerical boundaries continuously approach their classical counterparts as the dissipation increases. We further characterize each phase in terms of spin correlations, reduced density matrices, and multipartite entanglement, thereby clarifying how frustration and strong environmental coupling jointly organize the ground state of this minimal open plaquette.

The remainder of the paper is organized as follows. Section 2 introduces the four-spin spin-boson Hamiltonian, the bath spectral density, and the VMPS approach used to obtain the ground state. Section 3 presents the strong-coupling analysis, the ground-state phase diagram, and the correlation and entanglement properties of the competing ordered phases. Finally, Section 4 summarizes our main results and discusses possible extensions.

2 Model and methodologies

2.1 Model

The generalized multi-spin-boson model studied in this work consists of four interacting spins located on a square plaquette, where each spin is strongly coupled to its own independent bosonic bath. The Hamiltonian of the full system can be written as

(1)H=H_{\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}+H_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}+H_{\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}\text{−}\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}.

The spin Hamiltonian is given by

(2)\begin{array}{rl}{H}_{\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}=& \frac{\mathrm{\Delta }}{2}\sum _{i=\mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D}}{\sigma }_x^{(i)}+J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}\sum _{⟨ij⟩}{\sigma }_z^{(i)}{\sigma }_z^{(j)}\\ & +J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}\sum _{⟨⟨ij⟩⟩}{\sigma }_z^{(i)}{\sigma }_z^{(j)},\end{array}

where \mathrm{\Delta } denotes the tunneling strength of each spin. The interaction J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} acts along the four nearest-neighbor edges of the plaquette, while J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} couples the two diagonal spin pairs. Depending on the sign of J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} and J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}, the spin interactions may favor parallel or antiparallel alignment, and the competition between edge and diagonal couplings naturally leads to frustration on the square geometry.

The bath Hamiltonian reads

(3)H_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}={\displaystyle \sum _{i=\mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D}}{\displaystyle \sum _k{\omega }_k{b}_{k,i}^{†}{b}_{k,i},}}

where b_{k,i}^{†} and b_{k,i} are the bosonic creation and annihilation operators for mode k in the bath coupled to spin i, and {\omega }_k is the corresponding mode frequency. The spin−bath interaction takes the form

(4)H_{\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}\text{−}\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}={\displaystyle \sum _{i=\mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D}}{\displaystyle \sum _k{\displaystyle \frac{1}{2}}{g}_k{\sigma }_z^{(i)}(b_{k,i}^{†}+b_{k,i}),}}

with g_k characterizing the coupling of each spin to the environmental modes. Its frequency dependence is encoded in the spectral density

(5)J(\omega )=\pi {\displaystyle \sum _k{g}_k^2\delta (\omega -{\omega }_k)=2\pi \alpha {\omega }_c^{1-s}{\omega }^s,}

where \alpha is the dimensionless dissipation strength, {\omega }_c is the cutoff frequency, and s is the bath exponent. The sub-Ohmic regime () enhances the weight of low-frequency modes and strengthens localization tendencies. This four-spin system exhibits two key symmetries. The first is a global {\mathbb{Z}}_2 symmetry, identical to that of the standard spin-boson model, under which the transformation {\sigma }_z^{(i)}\to -{\sigma }_z^{(i)} must be accompanied by b_{k,i}\to -b_{k,i} and b_{k,i}^{†}\to -b_{k,i}^{†} for all baths.

The second important symmetry is associated with the square geometry. In addition to the trivial relabeling of the four spins, the model possesses a nontrivial diagonal-flip operation: flipping one of the two diagonal spin pairs (for definiteness, spins \mathrm{B} and \mathrm{D}) together with the corresponding bath transformation (b_{k,\mathrm{B}},b_{k,\mathrm{D}})\to-(b_{k,\mathrm{B}},b_{k,\mathrm{D}}) leaves both the bath Hamiltonian and the spin-bath coupling invariant, while changing J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} to -J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} without affecting the diagonal interaction J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}. As a result, the Hamiltonians at \pm J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} are exactly unitarily equivalent. Introducing the unitary operator

(6)U\equiv {\sigma }_x^{\mathrm{B}}{\sigma }_x^{\mathrm{D}}{P}_{\mathrm{B}}{P}_{\mathrm{D}},P_i\equiv \exp (\mathrm{i}\pi {\displaystyle \sum _k{b}_{k,i}^{†}{b}_{k,i}),}

where i=\mathrm{B},\mathrm{D}, one has

(7)H(-J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}},J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}})=UH(J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}},J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}})U^{†},

where {\sigma }_x^{\mathrm{B}}{\sigma }_x^{\mathrm{D}} flips {\sigma }_z^{\mathrm{B}(\mathrm{D})}\to -{\sigma }_z^{\mathrm{B}(\mathrm{D})} and the bath-parity operators implement b_{k,\mathrm{B}(\mathrm{D})}\to -b_{k,\mathrm{B}(\mathrm{D})}. This exact equivalence implies that the spectrum and all phase boundaries are strictly symmetric under J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}\to -J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} for identical baths. This equivalence relies on the fact that all four baths are identical and are independently coupled to their respective spins with the same spectral density and coupling strength. Consequently, all phase boundaries for can be obtained from the J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}>0 data by this mapping, and in our VMPS calculations we therefore restrict ourselves to J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}\ge 0; the case will be discussed by exploiting this symmetry in Section 3.

These two symmetries play a central role in determining the quantum phases of the model: the global {\mathbb{Z}}_2 symmetry is spontaneously broken at the delocalized-localized transition, while the competition between J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} and J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} leads to frustration and symmetry breaking among different spin configurations on the plaquette.

2.2 Methodologies

To obtain the ground state of the Hamiltonian, we apply the Variational Matrix Product States (VMPS) method, which provides an efficient representation for one-dimensional systems with short-range couplings owing to the entanglement area law [12, 13]. Following the standard treatment of the spin-boson model [14, 15], each continuous bath spectral density is logarithmically discretized into N bosonic modes (chain sites) labeled by n=0,1,\cdots ,N-1, with a higher resolution in the low-frequency region that dominates the sub-Ohmic critical behavior. The discretization scheme and numerical settings (including the logarithmic parameter \mathrm{\Lambda } and the construction of the discretized modes) are given in Appendix B.

Each discretized bath is then mapped to a one-dimensional chain using a standard tridiagonalization (Lanczos) procedure applied to the discretized star Hamiltonian, which is mathematically equivalent to an orthogonal-polynomial-based unitary transformation [16, 17], as illustrated schematically in Fig. 1. The mapping

(8)a_{n,i}={\displaystyle \sum _k{U}_{nk}{b}_{k,i},}

transforms the original star geometry [Fig. 1(a)] into a semi-infinite tight-binding chain [Fig. 1(b)], in which each spin couples only to the first site of its own chain. After the transformation, the Hamiltonian becomes

(9)H^{\prime }=H_{\mathrm{S}}+{\displaystyle \sum _{i=A}^{\mathrm{D}}(H_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h},i}^{\prime }+H_{\mathrm{i}\mathrm{n}\mathrm{t},i}^{\prime }),}

with

(10)H_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h},i}^{\prime }={\displaystyle \sum _{n=0}^{N-1}{\omega }_n{a}_{n,i}^{†}{a}_{n,i}+{\displaystyle \sum _{n=0}^{N-2}{t}_n(a_{n,i}^{†}{a}_{n+1,i}+a_{n+1,i}^{†}{a}_{n,i})}},

(11)H_{\mathrm{i}\mathrm{n}\mathrm{t},i}^{\prime }={\displaystyle \frac{c_0}{2}}{\sigma }_z^{(i)}(a_{0,i}+a_{0,i}^{†}).

To represent the four mapped bath chains within a single MPS, we adopt an interlaced site ordering, as sketched in Fig. 1(b): the chains associated with baths \mathrm{A} and \mathrm{B} are interlaced on the left part of the MPS, while those for baths \mathrm{C} and \mathrm{D} are interlaced on the right, and the four spins are grouped into two-site blocks (AB and CD) in the center. Importantly, this is purely a numerical site layout used in the VMPS/OBB implementation and does not modify the physical Hamiltonian. Its only effect is that the nearest-neighbor hopping terms along each individual Wilson chain become non-nearest-neighbor (typically next-nearest-neighbor) with respect to the chosen MPS ordering; these couplings are therefore included explicitly during the variational optimization.

To efficiently represent the large bosonic displacements along each chain, we employ the optimized boson basis (OBB) method [14], which constructs dynamically displaced local bases adapted to the instantaneous expectation values of {\hat{x}}_{n,i}=\frac{1}{\sqrt{2}}(a_{n,i}+a_{n,i}^{†}). In practice, each bosonic site is first truncated in the bare Fock basis to N_{\mathrm{p}\mathrm{h}} states, and then compressed into an optimized local basis of dimension d via the OBB transformation. This substantially reduces the required local Hilbert-space dimension (typically d\sim 8) without loss of accuracy near criticality. Convergence with respect to the truncation parameters (N_{\mathrm{p}\mathrm{h}},d) and the MPS bond dimension D is discussed in Appendix A.

After implementing VMPS with the OBB representation, we obtain the converged ground-state wave function |{\psi }_{\mathrm{G}\mathrm{S}}⟩ of the full four-spin-four-bath system. From this state, we compute the ground-state energy E_{\mathrm{G}\mathrm{S}}, the local magnetizations ⟨{\sigma }_z^{(i)}⟩, and the bosonic excitation numbers ⟨a_{n,i}^{†}{a}_{n,i}⟩, which together characterize the delocalized regime and the localized ordered regimes with Néel and stripe character.

In addition to these local observables, we will also make use of reduced spin density matrices in order to resolve the entanglement structure on the plaquette. In particular, we extract the two-spin density matrices {\rho }_{\mathrm{A}\mathrm{B}} and {\rho }_{\mathrm{A}\mathrm{C}} and the four-spin reduced state {\rho }_S={\mathrm{T}\mathrm{r}}_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}(|{\psi }_{\mathrm{G}\mathrm{S}}⟩⟨{\psi }_{\mathrm{G}\mathrm{S}}|). From these we compute Wootters concurrences on the nearest-neighbor bond \mathrm{A}\mathrm{B} and the diagonal bond \mathrm{A}\mathrm{C}, the von Neumann entropy of the four-spin block, and classical weights in the Néel and stripe sectors defined via simple projectors onto the corresponding Ising configurations. The explicit definitions and physical interpretation of these diagnostics are presented in Section 3.3, where they are used to analyze how frustration and dissipation reorganize quantum correlations between intra-plaquette bonds and the environmental degrees of freedom.

3 Results and discussion

3.1 Strong-coupling limit: Effective classical plaquette

As summarized in Appendix C, integrating out the independent bosonic baths generates long-range interactions along the imaginary-time direction. In the strong-dissipation regime (\alpha \gg 1), the bath-induced kernel strongly suppresses spin flips and pins each Ising worldline to an almost time-independent configuration,

(12)s_i(\tau )\simeq s_i=\pm 1,i=\mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D}.

Within this \alpha \to \mathrm{\infty } freezing approximation (i.e., neglecting rare spin flips), the bath contribution reduces, to leading order, to a configuration-independent energy shift. Consequently, the competition between different plaquette patterns is controlled primarily by the Ising couplings J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} and J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}, which yields the simple analytic boundaries quoted below. The effective classical energy therefore becomes

(13)E_{\mathrm{c}\mathrm{l}}(\{s_i\})=J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}{\displaystyle \sum _{⟨ij⟩}{s}_i{s}_j+J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}{\displaystyle \sum _{⟨⟨ij⟩⟩}{s}_i{s}_j+E_0(\alpha ,s),}}

where E_0(\alpha ,s) is a global constant that is irrelevant for phase competition. Equation (13) shows that the strong-coupling localized limit is exactly equivalent to a classical frustrated four-spin plaquette controlled solely by (J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}},J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}).

With the spin labeling as in Fig. 1(a), the plaquette admits three inequivalent ordered patterns (up to global inversion and point-group operations):

(i) Ferromagnetic (FM). All spins align, (s_{\mathrm{A}},s_{\mathrm{B}},s_{\mathrm{C}},s_{\mathrm{D}})=(+,+,+,+). All four nearest-neighbor bonds and both diagonals contribute with the same sign, yielding

(14)E_{\mathrm{F}\mathrm{M}}=4J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}+2J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}+E_0.

(ii) Néel (checkerboard). Nearest neighbors alternate around the square, e.g., (s_A,s_B,s_C,s_D)=(+,-,+,-). The four edge bonds are antiparallel while the two diagonal pairs are parallel, so that

(15)E_{\mathrm{N}\text{é}\mathrm{e}\mathrm{l}}=-4J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}+2J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}+E_0.

(iii) Stripe. Spins align in two adjacent pairs, e.g., (s_{\mathrm{A}},s_{\mathrm{B}},s_{\mathrm{C}},s_{\mathrm{D}})=(+,+,-,-) (horizontal stripes) or (+,-,-,+) (vertical stripes). In both cases, the four nearest-neighbor contributions cancel pairwise while both diagonals are antiparallel, giving

(16)E_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p}\mathrm{e}}=-2J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}+E_0.

The horizontal and vertical stripe states are degenerate.

The analytical phase boundaries follow from simple energy comparisons. Comparing Eqs. (14) and (15) gives

(17)E_{\mathrm{F}\mathrm{M}}=E_{\mathrm{N}\text{é}\mathrm{e}\mathrm{l}}\iff J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}=0,

which separates FM order () from Néel order (J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}>0) at fixed J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}. Comparing the stripe energy with the Néel (FM) energy yields

(18)E_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p}\mathrm{e}}=E_{\mathrm{N}\text{é}\mathrm{e}\mathrm{l}}\iff J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}=J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}},J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}>0,

(19)

Equivalently, the stripe phase is stabilized for sufficiently strong positive diagonal coupling,

(20)J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}>|J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}|,

while for the sign of J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} selects between Néel and ferromagnetic order, in agreement with Fig. 2. All three boundaries meet at (J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}},J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}})=(0,0), which is the classical triple point of the plaquette.

In the strong-dissipation regime, the phase boundaries are essentially independent of the bath exponent s and the coupling strength \alpha. More precisely, in the \alpha \to \mathrm{\infty } freezing limit (neglecting rare spin-flip events), the baths mainly suppress the tunneling dynamics and shift the absolute ground-state energy, while the competition between different static plaquette configurations is governed by J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} and J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} and thus fixes the phase boundaries. Subleading quantum-fluctuation corrections at finite \alpha may lead to small quantitative shifts but do not affect the qualitative phase structure discussed here.

3.2 Phase diagram

We now turn to the ground-state phase diagram of the four-spin spin-boson model. Figure 3(a) shows the phase diagram in the (\alpha ,J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}) plane for the sub-Ohmic bath exponent s=0.3, fixed tunneling \mathrm{\Delta }=0.1, and nearest-neighbor coupling J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}=0.02. Three distinct regions are identified: (i) a delocalized phase in which all local magnetizations vanish, ⟨{\sigma }_i^z⟩=0 for i=\mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D}; (ii) a localized Néel-ordered phase at smaller J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}, in which nearest neighbors alternate while the two diagonal pairs remain parallel; and (iii) a localized stripe-ordered phase at sufficiently large positiveJ_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}, where two neighboring spins align while the opposite pair anti-aligns, realizing a (+,+,-,-) pattern up to point-group symmetry. These ordered phases continuously connect to the Néel and stripe configurations of the strong-coupling classical plaquette discussed in Section 3.1 and Fig. 2: as \alpha increases, bath-induced localization suppresses tunneling and the numerical phase boundaries approach the strong-coupling lines J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}=0 and J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}=\pm J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}. The complementary phase diagram in the (\alpha ,J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}) plane for fixed J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}=0.02 is shown in Fig. 3(b) and exhibits the same three phases.

For a fixed J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}, increasing the dissipation strength \alpha drives a localization transition out of the delocalized phase. For J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} below the stripe threshold the localized phase is Néel ordered, whereas for sufficiently large positive J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} the system instead develops stripe order. In both cases the local magnetizations ⟨{\sigma }_i^z⟩ emerge continuously from zero at the boundary of the delocalized phase. Consistently, the ground-state energy remains smooth and no discontinuity is resolved in \mathrm{\partial }{E}_{\mathrm{G}\mathrm{S}}/\mathrm{\partial }\alpha at the transition within our numerical resolution, supporting continuous (second-order) quantum phase transitions driven by bath-induced localization, in close analogy with the standard {\mathbb{Z}}_2 transition of the single-spin spin-boson model.

To characterize the competition between Néel and stripe order, we consider a vertical cut at fixed \alpha \simeq 0.026 [dashed line in Fig. 3(a)]. Along this cut, as shown in Fig. 3(c), the four magnetizations ⟨{\sigma }_i^z⟩ form two symmetry-related pairs in the ordered regimes. For J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}\lesssim 0.016 the system realizes a Néel pattern, with nearest neighbors anti-aligned and the two diagonal spins parallel. In an intermediate window of J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} the magnetizations collapse to numerically vanishing values, signaling a reentrant delocalized regime. Here “reentrant” means that upon increasing J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} at fixed \alpha, the system leaves a localized ordered phase, enters a delocalized window, and then re-enters a localized ordered phase (with stripe character at larger positive J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}). For J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}\gtrsim 0.027 the magnetizations bifurcate again into a stripe configuration, in which a neighboring spin pair points up while the opposite pair points down (up to symmetry).

The corresponding ground-state energy E_{\mathrm{G}\mathrm{S}}(J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}) and its numerical derivative \mathrm{\partial }{E}_{\mathrm{G}\mathrm{S}}/\mathrm{\partial }{J}_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} are shown in Fig. 3(d). The two delocalized-to-ordered boundaries along this cut are continuous, consistent with the smooth onset of ⟨{\sigma }_i^z⟩ and the absence of a kink in E_{\mathrm{G}\mathrm{S}}. By contrast, the Néel-stripe boundary is accompanied by a sharp kink in E_{\mathrm{G}\mathrm{S}} and a pronounced jump-like feature in \mathrm{\partial }{E}_{\mathrm{G}\mathrm{S}}/\mathrm{\partial }{J}_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}, which signals a first-order transition between the two distinct localized orders. The three boundaries separating the delocalized, Néel, and stripe regions meet at a three-phase junction in the (\alpha ,J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}) plane (a quantum triple point in this few-body setting), where the system is simultaneously susceptible to localization and to the competition between the two types of spatial order. Figures 3(c) and (d) correspond to the same cut of Fig. 3(a) and therefore share the same fixed value of \alpha.

Finally, we briefly comment on the role of the sign of J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}. As discussed in Section 2, a diagonal-flip transformation on spins \mathrm{B} and \mathrm{D}, combined with the corresponding bath-parity operation, maps the Hamiltonian at \pm J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} onto one another via the exact unitary equivalence H(-J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}},J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}})=UH(J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}},J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}})U^{†} [Eq. (7)]. This implies that the spectrum, and hence the phase boundaries in the (\alpha ,J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}) and (\alpha ,J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}) planes, are strictly symmetric under J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}\to -J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} (for identical baths). At the level of ordered states, the mapping interchanges the Néel and ferromagnetic configurations, while leaving the stripe and delocalized regimes invariant up to a relabeling of sites. Consequently, the phase diagram for J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}=-0.02 is obtained from Fig. 3(a) by replacing the Néel region with a ferromagnetic region, and Fig. 3(b) only displays the J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}>0 half of the (\alpha ,J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}) plane; the part follows by reflection about J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}=0. Along the cut shown in Figs. 3(c) and (d), the same symmetry implies that the magnetization curves for are obtained from Fig. 3(c) by reflecting all ⟨{\sigma }_i^z⟩ about the horizontal axis, while E_{\mathrm{G}\mathrm{S}}(J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}) and its derivative remain unchanged.

3.3 Entanglement structure

To elucidate how frustration and dissipation reorganize quantum correlations on the plaquette, we analyze several entanglement diagnostics along the vertical cut at fixed dissipation strength \alpha \simeq 0.026 in Fig. 3(a), with J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}=0.02, \mathrm{\Delta }=0.1, and s=0.3. From the converged VMPS ground state |{\psi }_{\mathrm{G}\mathrm{S}}⟩ we form the reduced four-spin density matrix

(21){\rho }_S={\mathrm{T}\mathrm{r}}_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}\left(|{\psi }_{\mathrm{G}\mathrm{S}}⟩⟨{\psi }_{\mathrm{G}\mathrm{S}}|\right),

and obtain two-spin reduced states {\rho }_{\mathrm{A}\mathrm{B}} and {\rho }_{\mathrm{A}\mathrm{C}} by tracing out the remaining spins. Here and below, product states such as |\uparrow \downarrow \uparrow \downarrow ⟩ are understood in the {\sigma }^z basis and ordered as (\mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D}).

Diagnostics. The entanglement between two spins is quantified by the Wootters concurrence C(\rho ) for two qubits [18]. We denote by C_{\mathrm{A}\mathrm{B}} the concurrence of {\rho }_{\mathrm{A}\mathrm{B}} (a nearest-neighbor bond) and by C_{\mathrm{A}\mathrm{C}} the concurrence of {\rho }_{\mathrm{A}\mathrm{C}} (a diagonal bond). To characterize the coupling between the plaquette and its environment we use the von Neumann entropy of the four-spin block,

(22)S_{4\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}=-\mathrm{T}\mathrm{r}\left({\rho }_S\ln {\rho }_S\right),

with the natural logarithm.

To monitor the competition between distinct classical ordering tendencies we also define configuration weights

(23)p_{\mu }=\mathrm{T}\mathrm{r}(P_{\mu }{\rho }_S),\mu =\mathrm{N}\text{é}\mathrm{e}\mathrm{l},\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p}\mathrm{e},

where P_{\mathrm{N}\text{é}\mathrm{e}\mathrm{l}} and P_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p}\mathrm{e}} project onto the subspaces spanned by the corresponding sets of Ising configurations. Explicitly, we use the six classical plaquette states

(24)\begin{array}{rl}& |N_1⟩=|\uparrow \downarrow \uparrow \downarrow ⟩,|N_2⟩=|\downarrow \uparrow \downarrow \uparrow ⟩,\\ & |S_1⟩=|\uparrow \uparrow \downarrow \downarrow ⟩,|S_2⟩=|\downarrow \downarrow \uparrow \uparrow ⟩,\\ & |S_3⟩=|\uparrow \downarrow \downarrow \uparrow ⟩,|S_4⟩=|\downarrow \uparrow \uparrow \downarrow ⟩,\end{array}

and define

(25)P_{\mathrm{N}\text{é}\mathrm{e}\mathrm{l}}={\displaystyle \sum _{m=1}^2|N_m⟩⟨N_m|,P_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p}\mathrm{e}}={\displaystyle \sum _{m=1}^4|S_m⟩⟨S_m|.}}

Note that p_{\mathrm{N}\text{é}\mathrm{e}\mathrm{l}}+p_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p}\mathrm{e}}\le 1 in general, since {\rho }_S may also have weight in other spin configurations (including local defects) as well as coherences between different sectors.

Bond entanglement and environment entanglement.Figure 4(a) shows that deep in the Néel phase (J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}\lesssim 0.016) the concurrence is concentrated on the nearest-neighbor bond: C_{\mathrm{A}\mathrm{B}}>0 while C_{\mathrm{A}\mathrm{C}}\approx 0. In the stripe phase (J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}\gtrsim 0.027) the situation is reversed: C_{\mathrm{A}\mathrm{C}}>0 whereas C_{\mathrm{A}\mathrm{B}}\approx 0, reflecting that the dominant correlations shift from nearest-neighbor to diagonal bonds when the underlying spin pattern switches from checkerboard (Néel) to stripe order. In the intermediate region 0.016\lesssim J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}\lesssim 0.027, where the phase diagram exhibits a reentrant delocalized regime, both concurrences are strongly suppressed and become numerically indistinguishable from zero.

The four-spin entropy S_{4\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}} in Fig. 4(b) complements this picture. In the ordered Néel and stripe phases, S_{4\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}} remains relatively small, indicating that the plaquette is only weakly entangled with the baths when a single ordering tendency dominates. Upon entering the reentrant delocalized window, S_{4\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}} increases rapidly and develops a broad maximum. This signals enhanced plaquette–bath entanglement concomitant with the quenching of two-spin entanglement within the plaquette, consistent with the notion that strong dissipation can transfer quantum correlations from intra-plaquette bonds to environmental degrees of freedom.

Competition between Néel and stripe sectors. The configuration weights in Fig. 4(c) provide a direct view of how the dominant classical sector evolves. For small J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} the Néel manifold dominates, p_{\mathrm{N}\text{é}\mathrm{e}\mathrm{l}}\gg p_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p}\mathrm{e}}, consistent with the finite C_{\mathrm{A}\mathrm{B}}. Across the reentrant delocalized window, weight is gradually transferred from the Néel to the stripe manifold, and for J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}\gtrsim 0.027 the stripe sector becomes dominant, p_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p}\mathrm{e}}>p_{\mathrm{N}\text{é}\mathrm{e}\mathrm{l}}, in one-to-one correspondence with the emergence of diagonal concurrence C_{\mathrm{A}\mathrm{C}}. The reentrant regime thus appears as an intermediate region in which neither ordering manifold fully dominates, and bond entanglement is quenched while plaquette–bath entanglement is enhanced.

Microscopic structure at the center of the reentrance. To clarify the structure of the reentrant delocalized regime, we inspect a representative point near its center, (\alpha ,J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}})\simeq (0.026,0.02). We restrict {\rho }_S to the eight-dimensional subspace spanned by the Néel, stripe, and ferromagnetic configurations \{|N_{1,2}⟩,|S_{1-4}⟩,|F_{1,2}⟩\} with |F_1⟩=|\uparrow \uparrow \uparrow \uparrow ⟩ and |F_2⟩=|\downarrow \downarrow \downarrow \downarrow ⟩, and introduce the associated cat/anti-cat basis

(26)\begin{array}{rl}& |{\varphi }_N^{\pm }⟩=(|N_1⟩\pm |N_2⟩)/\sqrt{2},|{\varphi }_F^{\pm }⟩=(|F_1⟩\pm |F_2⟩)/\sqrt{2},\\ & |{\varphi }_{S,h}^{\pm }⟩=(|S_1⟩\pm |S_2⟩)/\sqrt{2},|{\varphi }_{S,v}^{\pm }⟩=(|S_3⟩\pm |S_4⟩)/\sqrt{2}.\end{array}

The corresponding diagonal weights w_{\mu }=⟨{\varphi }_{\mu }|{\rho }_S|{\varphi }_{\mu }⟩ are displayed in Fig. 5 (top). We stress that these eight states span only a subspace of the full sixteen-dimensional four-spin Hilbert space; the remaining weight resides in configurations with local defects and in additional spin-bath dressing, and also contributes to the entropy S_{4\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}. Within the cat/anti-cat subspace, the three symmetric cat states |{\varphi }_N^{+}⟩, |{\varphi }_{S,h}^{+}⟩, and |{\varphi }_{S,v}^{+}⟩ carry comparable and largest weights (w_{N^{+}}\simeq w_{S_h^{+}}\simeq w_{S_v^{+}}\approx 0.13), while the antisymmetric cats and the ferromagnetic cats are subleading. Together with the residual weight outside the subspace, this indicates that in the reentrant regime {\rho }_S is best viewed as a strongly mixed state built from multiple competing Néel- and stripe-like cat configurations dressed by substantial plaquette–bath entanglement, rather than being dominated by a single ordered cat state.

Bell-basis structure and vanishing concurrence. To connect this mixed cat structure to the suppression of bond entanglement, we transform {\rho }_{\mathrm{A}\mathrm{B}} and {\rho }_{\mathrm{A}\mathrm{C}} into the Bell basis \{|{\mathrm{\Phi }}^{\pm }⟩,|{\mathrm{\Psi }}^{\pm }⟩\} defined by

|{\mathrm{\Phi }}^{\pm }⟩={\displaystyle \frac{|\uparrow \uparrow ⟩\pm |\downarrow \downarrow ⟩}{\sqrt{2}}},|{\mathrm{\Psi }}^{\pm }⟩={\displaystyle \frac{|\uparrow \downarrow ⟩\pm |\downarrow \uparrow ⟩}{\sqrt{2}}}.

As shown in Fig. 5 (bottom), both reduced states are nearly Bell diagonal. For Bell-diagonal two-qubit states, the concurrence is determined solely by the largest Bell population and becomes nonzero only if that population exceeds 1/2. Here all Bell populations remain below 1/2, implying C_{\mathrm{A}\mathrm{B}}\approx C_{\mathrm{A}\mathrm{C}}\approx 0, in agreement with Fig. 4(a). We therefore conclude that the reentrant delocalized regime is characterized by strong plaquette–bath entanglement and a highly mixed competition between Néel- and stripe-like cat configurations, which together wash out residual two-spin entanglement on both nearest-neighbor and diagonal bonds.

Remark on . By the diagonal-flip unitary equivalence discussed earlier (a local {\sigma }_x flip on spins \mathrm{B} and \mathrm{D} accompanied by the corresponding bath parity transformations), the Hamiltonians with \pm J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} are related by a unitary map. Since concurrence and S_{4\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}} are invariant under local unitaries on the spin degrees of freedom, the curves C_{\mathrm{A}\mathrm{B}}(J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}), C_{\mathrm{A}\mathrm{C}}(J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}), and S_{4\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}(J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}) along the same vertical cut are unchanged when J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} is replaced by -J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}. What changes is the identification of the dominant classical manifold: the mapping interchanges the Néel and ferromagnetic sectors, while leaving the stripe sector invariant up to relabeling. Accordingly, the sector weights for are obtained from Fig. 4(c) by replacing the Néel configurations \{|N_{1,2}⟩\} (and cats |{\varphi }_N^{\pm }⟩) with their ferromagnetic counterparts \{|F_{1,2}⟩\} (and |{\varphi }_F^{\pm }⟩).

4 Conclusions and outlook

We have studied a frustrated four-spin plaquette spin-boson model in which each spin is coupled to an independent sub-Ohmic bath and the spins interact via nearest-neighbor and diagonal Ising-type couplings. This setting embeds the competition between Néel and stripe ordering on a square plaquette into a dissipative multi-impurity environment with a global {\mathbb{Z}}_2 symmetry. Using a strong-coupling path-integral analysis together with VMPS simulations performed in the full spin-bath Hilbert space, we have clarified how frustration and dissipation jointly organize the ground state and its correlation structure.

In the strong-dissipation limit, integrating out the baths generates long-ranged temporal correlations that suppress tunneling and freeze the spin worldlines, reducing the problem to an effective classical plaquette with nearest-neighbor and diagonal Ising couplings. The resulting ferromagnetic, Néel, and stripe phases are separated by simple analytic boundaries J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}=0 and J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}=\pm J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}, which are independent of the bath parameters. Consistent with this benchmark, our VMPS phase diagrams approach these straight lines as \alpha increases and, at intermediate dissipation, reveal a delocalized phase separating two localized ordered phases with Néel and stripe character. The transitions out of the delocalized phase are continuous and correspond to spontaneous breaking of the global {\mathbb{Z}}_2 symmetry, while the Néel-stripe transition is first order and meets the delocalized boundary at a triple point.

We emphasize that the present system is a minimal four-spin (zero-dimensional) plaquette coupled to independent bosonic baths. Accordingly, the “phase transitions” discussed here should be understood as well-defined dissipation-driven ground-state transitions (and corresponding changes in order parameters and reduced density matrices) in a few-body open quantum system, rather than thermodynamic long-range ordering in the strict N\to \mathrm{\infty } sense. Establishing true long-range order requires extended lattices or arrays in the thermodynamic limit. Nevertheless, the plaquette provides a minimal building block that captures the microscopic mechanism by which dissipation enhances localization and selects competing Néel/stripe patterns, and it naturally motivates future studies of coupled plaquettes and quasi-one-dimensional arrays.

Beyond identifying phase boundaries, we have analyzed how entanglement is redistributed across the plaquette and its environment. In the Néel (stripe) phase, two-spin entanglement is predominantly carried by a nearest-neighbor (diagonal) bond, while the four-spin entropy S_{4\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}} remains relatively small, indicating comparatively weak plaquette–bath entanglement. In the reentrant delocalized regime, both bond concurrences are strongly suppressed whereas S_{4\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}} develops a broad maximum, showing that independent baths absorb a substantial fraction of the total quantum correlations. A decomposition of the reduced spin density matrix into cat and anti-cat superpositions of classical Néel, stripe, and ferromagnetic configurations further indicates that this regime is best viewed as a strongly mixed ensemble of competing cat-like patterns dressed by pronounced spin−bath entanglement and configuration defects, rather than a single ordered cat state.

More broadly, our work illustrates how structured environments can both stabilize frustration-induced ordering and reshuffle quantum correlations between intra-plaquette bonds and environmental degrees of freedom. The combination of VMPS (with optimized boson bases), strong-coupling benchmarks, and reduced-density-matrix diagnostics developed here can be extended in several directions, including larger frustrated clusters and quasi-one-dimensional arrays of plaquettes, as well as settings with phase-dependent couplings or partially shared baths. These extensions provide a natural route to explore the interplay of dissipation-driven localization with geometric frustration, and, in dynamical or nonequilibrium variants, with chiral response and nonreciprocal transport.

5 Appendix A: Convergence analysis of the VMPS method

In this Appendix, we present a detailed analysis of the numerical convergence of the variational matrix product state (VMPS) method employed in this work. VMPS is particularly well suited for strongly correlated quantum systems with an effectively one-dimensional structure. In the present four-spin plaquette spin-boson model, the star-to-chain mapping transforms each independent bosonic bath into a Wilson chain. These four chains are subsequently arranged into a single MPS geometry with an interlaced ordering, and the spin degrees of freedom are grouped into two-site blocks (AB and CD), as sketched in Fig. 1(b). This nontrivial one-dimensional embedding makes careful convergence checks especially important, in particular in parameter regimes with enhanced entanglement.

The accuracy of the VMPS calculations is controlled by several numerical parameters. The MPS bond dimension D sets the maximum number of Schmidt states retained across each bond. For the bosonic degrees of freedom, each Wilson-chain site is first truncated in the bare Fock basis to N_{\mathrm{p}\mathrm{h}} states, and subsequently compressed using the optimized boson basis (OBB) with local dimension d. Finally, we denote by N the Wilson-chain length kept for each bath after logarithmic discretization and star-to-chain mapping. Increasing N reduces discretization and finite-chain effects, while increasing (D,d,N_{\mathrm{p}\mathrm{h}}) systematically controls the MPS truncation and local bosonic truncation errors.

To illustrate the convergence behavior, we fix the microscopic parameters to representative values used throughout the main text, namely bath exponent s=0.3, dissipation strength \alpha =0.026, diagonal coupling J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}=0.02, and tunneling amplitude \mathrm{\Delta }=0.1. We then study the stripe order parameter across a cut close to the phase boundary. Using the site-resolved local magnetizations m_i\equiv ⟨{\sigma }_i^z⟩ (i=\mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D}), we define

(A1)\begin{array}{c}{m}_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p}\mathrm{e}}=\left|{\displaystyle \frac{m_{\mathrm{A}}+m_{\mathrm{B}}-m_{\mathrm{C}}-m_{\mathrm{D}}}{4}}\right|.\end{array}

Figure A1(a) shows the order parameter obtained for several sets of (D,d,N_{\mathrm{p}\mathrm{h}},N), as listed in the legend. Away from the phase boundary, the results converge rapidly and are essentially indistinguishable already for moderate bond dimensions. The largest deviations appear in the vicinity of the transition around J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}\simeq 0.015, where quantum fluctuations and bath-induced correlations are most pronounced. Importantly, even in this critical regime the curves systematically approach a stable result as the numerical parameters are increased, indicating that the phase boundary and the order-parameter behavior are well converged on the scale relevant for this work.

To further assess convergence at the level of entanglement structure, we analyze the bond entanglement entropy

(A2)\begin{array}{c}S(m)=-{\displaystyle \sum _i{\lambda }_i^2(m)\log \left({\lambda }_i^2(m)\right),}\end{array}

where {\lambda }_i(m) are the Schmidt coefficients associated with bond m. In practice, the sum is evaluated using the retained Schmidt spectrum of dimension D at bond m. Figure A1(b) displays S(m) as a function of the bond index along the MPS chain for J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}=0.015, chosen close to the transition point. Here, m=0 corresponds to the bond separating the two central spin blocks, while positive and negative indices label bonds extending into the interlaced bosonic chains.

In the interlaced ordering, the central bond m=0 separates the two spin blocks (AB and CD) and thus mainly probes correlations internal to the plaquette, whereas the neighboring bonds (m=\pm 1) already cut between a spin block and a part of the surrounding bath degrees of freedom. Consistent with the enhanced system–bath entanglement along this cut [cf. the broad maximum of S_{4\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}} in Fig. 4(b)] and the suppressed two-spin concurrences in Fig. 4(a), the maximal bond entropy is therefore found at m=\pm 1 rather than at m=0.

A pronounced peak in the entanglement entropy is observed near the central region, reflecting the strong correlations induced by frustration and dissipation. As the numerical parameters are increased, the entanglement profiles rapidly converge and essentially collapse onto a single curve, demonstrating that the VMPS ansatz captures the relevant entanglement structure across the entire chain.

Finally, Fig. A1(c) shows the decay of the squared Schmidt coefficients {\lambda }_k^2 at the bond where S(m) is maximal (here m=\pm 1), plotted as {\log }_{10}({\lambda }_k^2). The rapid decay of {\lambda }_k^2 indicates that only a small number of Schmidt states contribute significantly. Consequently, the discarded weight (truncation error) at bond m,

(A3)\begin{array}{c}{ϵ}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{c}}(m)=1-{\displaystyle \sum _{i=1}^{\mathrm{D}}{\lambda }_i^2(m),}\end{array}

remains very small (typically below {10}^{-6}) for the largest parameter sets shown. In Fig. A1(c) we evaluate {ϵ}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{c}} at the same bond where S(m) is maximal. This provides strong evidence that the numerical truncation errors are well controlled.

We emphasize that convergence near phase transitions is intrinsically more demanding due to the enhanced entanglement and the growth of the relevant correlation scale near criticality. For this reason, while the present data reliably locate phase boundaries and characterize the distinct phases, we refrain from extracting critical exponents from finite-D calculations.

6 Appendix B: Logarithmic discretization and star-to-chain mapping

In this Appendix, we summarize the logarithmic discretization of the bath spectral density and the subsequent star-to-chain mapping used to obtain the Wilson-chain parameters in the VMPS simulations. Our implementation follows the improved logarithmic discretization scheme introduced in Ref. [14] (including its supplementary material), which reduces discretization artifacts and yields a weaker dependence of critical properties on the discretization parameter.

6.1 Logarithmic discretization and star representation

For each spin i=\mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D}, the independent bosonic bath is characterized by the spectral density (with cutoff {\omega }_c)

(B1)

We discretize the interval (0,{\omega }_c) into logarithmic bins using a discretization parameter \mathrm{\Lambda }>1 and (for simplicity) choose z=1. We do not perform z-averaging in the present work; the residual \mathrm{\Lambda }-dependence is mitigated by the improved discretization procedure. The bin boundaries are chosen as

(B2)\begin{array}{c}{\omega }_n^{-}={\omega }_c{\mathrm{\Lambda }}^{-n},{\omega }_n^{+}={\omega }_c{\mathrm{\Lambda }}^{-(n+1)},n=0,1,\cdots ,\end{array}

so that the n-th interval corresponds to \omega \in ({\omega }_n^{+},{\omega }_n^{-}). Within each bin, the continuum bath is approximated by a single effective mode, leading to the star representation

(B3)\begin{array}{c}{H}_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}^{(i)}={\displaystyle \sum _{n=0}^{N-1}{\xi }_n{b}_{n,i}^{†}{b}_{n,i},H_{\mathrm{i}\mathrm{n}\mathrm{t}}^{(i)}={\displaystyle \frac{{\sigma }_z^{(i)}}{2}}{\displaystyle \sum _{n=0}^{N-1}{\gamma }_n\left(b_{n,i}+b_{n,i}^{†}\right)}},\end{array}

where {\xi }_n and {\gamma }_n are the effective mode frequencies and couplings, respectively.

A convenient and widely used choice is to match the first two moments of J(\omega ) in each bin. Defining

(B4)\begin{array}{c}{\gamma }_n^2={\displaystyle \frac{1}{\pi }}{\displaystyle {\int }_{{\omega }_n^{+}}^{{\omega }_n^{-}}J(\omega )\mathrm{d}\omega ,{\xi }_n={\displaystyle \frac{{\displaystyle {\int }_{{\omega }_n^{+}}^{{\omega }_n^{-}}\omega J(\omega )\mathrm{d}\omega }}{{\displaystyle {\int }_{{\omega }_n^{+}}^{{\omega }_n^{-}}J(\omega )\mathrm{d}\omega }}},}\end{array}

we obtain for the power-law form (B1)

(B5)\begin{array}{c}{\gamma }_n^2={\displaystyle \frac{2\alpha {\omega }_c^{1-s}}{1+s}}\left[({\omega }_n^{-}{)}^{1+s}-({\omega }_n^{+}{)}^{1+s}\right],\end{array}

(B6)\begin{array}{c}{\xi }_n=\frac{1+s}{2+s}{\displaystyle \frac{({\omega }_n^{-}{)}^{2+s}-({\omega }_n^{+}{)}^{2+s}}{({\omega }_n^{-}{)}^{1+s}-({\omega }_n^{+}{)}^{1+s}}}.\end{array}

In the improved discretization scheme of Ref. [14], the above moment-matching is refined by incorporating the logarithmic weight within each interval such that low-frequency contributions are represented more accurately; in practice, this modifies the effective parameters ({\xi }_n,{\gamma }_n) and significantly reduces discretization artifacts already at moderate \mathrm{\Lambda }. Our numerical coefficients are computed following the explicit prescriptions in Ref. [14].

6.2 Star-to-chain mapping and Wilson-chain parameters

The discretized star Hamiltonian (B3) is mapped to a one-dimensional chain via a standard tridiagonalization (Lanczos) procedure applied to the logarithmically discretized star representation. This mapping is mathematically equivalent to an orthogonal-polynomial-based unitary transformation, and we mention the orthogonal-polynomial connection only in this equivalence sense [16, 17]. In the present (real) bosonic bath problem, the transformation matrix can be chosen real, so that “unitary” reduces to an orthogonal matrix.

We first define the collective bath mode that couples directly to the spin,

(B7)\begin{array}{c}{c}_0^2={\displaystyle \sum _{n=0}^{N-1}{\gamma }_n^2,a_{0,i}={\displaystyle \frac{1}{c_0}}{\displaystyle \sum _{n=0}^{N-1}{\gamma }_n{b}_{n,i},}}\end{array}

so that the interaction becomes

(B8)\begin{array}{c}{H}_{\mathrm{i}\mathrm{n}\mathrm{t}}^{(i)}={\displaystyle \frac{c_0}{2}}{\sigma }_z^{(i)}\left(a_{0,i}+a_{0,i}^{†}\right).\end{array}

The remaining orthonormal chain modes a_{m,i} (m=1,2,\cdots ,N-1) are generated by the same transformation a_{m,i}={\sum }_n{U}_{mn}{b}_{n,i} that tridiagonalizes the bath Hamiltonian. In the resulting chain representation, each bath reads

(B9)\begin{array}{c}{H}_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}^{(i)}={\displaystyle \sum _{m=0}^{N-1}{\omega }_m{a}_{m,i}^{†}{a}_{m,i}+{\displaystyle \sum _{m=0}^{N-2}{t}_m\left(a_{m,i}^{†}{a}_{m+1,i}+a_{m+1,i}^{†}{a}_{m,i}\right)}},\end{array}

which is the form used in the main text.

Operationally, ({\omega }_m,t_m) can be obtained from the Lanczos recursion applied to the star Hamiltonian: starting from the normalized vector |v_0⟩\propto {\sum }_n{\gamma }_n|n⟩, one constructs an orthonormal basis \{|v_m⟩\} such that the diagonal matrix \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\xi }_n) is represented in this basis by a symmetric tridiagonal matrix with diagonal entries {\omega }_m and off-diagonal entries t_m. This procedure is numerically stable and guarantees that the first site m=0 carries the entire spin-bath coupling c_0 defined in Eq. (B7).

Throughout this work, the chain length N is chosen sufficiently large such that observables of interest (ground-state energy, local magnetizations, and entanglement diagnostics) are converged with respect to further increases of N (see Appendix A). The discretization parameter \mathrm{\Lambda } and the cutoff {\omega }_c are fixed consistently across all calculations; throughout this work we use \mathrm{\Lambda }=2. All four baths are taken to be identical and are discretized using the same (\mathrm{\Lambda },{\omega }_c,N) throughout.

7 Appendix C: Bath-induced long-range interactions and effective two-dimensional description

In this Appendix we outline how integrating out the independent bosonic baths generates effective long-range interactions in imaginary time and leads to a natural two-dimensional (2D) description of the multi-spin-boson model. This framework explains why the four-spin plaquette with independent baths exhibits spin configurations that closely resemble those of frustrated 2D Ising lattices, despite the small number of spins and purely local spin−bath couplings.

7.1 Path-integral representation and bath integration

We start from the generic multi-spin-boson Hamiltonian

(C1)\begin{array}{c}H=H_S+{\displaystyle \sum _{i,k}{\omega }_k{b}_{k,i}^{†}{b}_{k,i}+{\displaystyle \sum _{i,k}{\displaystyle \frac{1}{2}}{g}_k{\sigma }_z^{(i)}(b_{k,i}+b_{k,i}^{†}),}}\end{array}

where H_S denotes the spin Hamiltonian, including the tunneling term \mathrm{\Delta } and the Ising interactions J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} and J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} on the plaquette, while the index i labels the four spins \mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D}. The baths are independent and characterized by a spectral density

(C2)\begin{array}{c}J(\omega )=\pi {\displaystyle \sum _k{g}_k^2\delta (\omega -{\omega }_k)=2\pi \alpha {\omega }_c^{1-s}{\omega }^s,}\end{array}

with dissipation strength \alpha, cutoff frequency {\omega }_c, and bath exponent s.

The equilibrium partition function at inverse temperature \beta reads

(C3)\begin{array}{c}Z(\beta )=\mathrm{T}\mathrm{r}{\mathrm{e}}^{-\beta H}.\end{array}

Introducing a path-integral representation in imaginary time \tau \in [0,\beta ], and working in the {\sigma }_z basis, we denote by s_i(\tau )=\pm 1 the path of spin i, corresponding to the eigenvalues of {\sigma }_z^{(i)}(\tau ). In terms of these paths and the bosonic fields b_{k,i}(\tau ), the partition function can be written schematically as [5, 19]

(C4)\begin{array}{c}Z={\displaystyle \int \mathcal{D}[s_i(\tau )]\mathcal{D}[b,b^{†}]\exp [-S_S[s]-S_B[b]-S_{SB}[s,b]]},\end{array}

where S_S[s] arises from H_S, S_B[b] from the free bosons, and S_{SB}[s,b] from the spin−bath coupling.

Since the baths are Gaussian (free oscillators with linear coupling), the bosonic fields can be integrated out exactly. Performing the Gaussian functional integrals over b_{k,i} and b_{k,i}^{†} yields an influence functional that depends only on the spin paths [5],

(C5)\begin{array}{c}{\displaystyle \int \mathcal{D}[b,b^{†}]{\mathrm{e}}^{-S_B[b]-S_{SB}[s,b]}={\mathrm{e}}^{-S_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}[s]}.}\end{array}

Thus the partition function becomes

(C6)\begin{array}{c}Z={\displaystyle \int \mathcal{D}[s_i(\tau )]{\mathrm{e}}^{-S_{\mathrm{e}\mathrm{f}\mathrm{f}}[s]},S_{\mathrm{e}\mathrm{f}\mathrm{f}}[s]=S_S[s]+S_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}[s].}\end{array}

For each spin i, the bath contribution can be written in the standard form (up to an additive constant independent of the spin path) [5, 19]

(C7)\begin{array}{c}{S}_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}^{(i)}[s_i]={\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^{\beta }\mathrm{d}\tau {\displaystyle {\int }_0^{\beta }\mathrm{d}{\tau }^{\prime }K(\tau -{\tau }^{\prime })[s_i(\tau )-s_i({\tau }^{\prime }){]}^2,}}\end{array}

where the memory kernel K(\tau ) is determined by J(\omega ) via

(C8)\begin{array}{c}K(\tau )={\displaystyle {\int }_0^{\mathrm{\infty }}{\displaystyle \frac{\mathrm{d}\omega }{\pi }}{\displaystyle \frac{J(\omega )}{{\omega }^2}}{\displaystyle \frac{\cosh [\omega ({\displaystyle \frac{\beta }{2}}-|\tau |)]}{\sinh ({\displaystyle \frac{\beta \omega }{2}})}}.}\end{array}

Summing over all spins,

(C9)\begin{array}{c}{S}_{\text{bath}}[s]={\displaystyle \sum _i{S}_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}^{(i)}[s_i].}\end{array}

In the sub-Ohmic regime and at low temperatures (\beta {\omega }_c\gg 1), the kernel shows a characteristic long-time power-law tail,

(C10)

more precisely for \beta {\omega }_c\gg 1 and time separations 1/{\omega }_c\ll |\tau |\ll \beta. The long-time behavior is controlled solely by the low-frequency scaling J(\omega )\propto {\omega }^s and is insensitive to the detailed ultraviolet cutoff shape.

7.2 Imaginary-time discretization and effective 2D classical model

To make contact with a classical spin model, we discretize the imaginary-time interval into N_{\tau } slices:

(C11)\begin{array}{c}{\tau }_{\ell }=\ell \mathrm{\Delta }\tau ,\ell =0,1,\cdots ,N_{\tau }-1;\mathrm{\Delta }\tau =\beta /N_{\tau },\end{array}

and define Ising variables

(C12)\begin{array}{c}{s}_{i,\ell }\equiv s_i({\tau }_{\ell })=\pm 1,\end{array}

with periodic boundary condition s_{i,N_{\tau }}\equiv s_{i,0}. (We use N_{\tau } to avoid confusion with the Wilson-chain length N employed in the VMPS calculations.)

The contribution from the spin Hamiltonian H_S yields, after a Trotter decomposition, a short-range action S_S[s] that is local or nearest-neighbor in the discrete time index \ell. For the Ising couplings at fixed time slice, one obtains terms of the form

(C13)\begin{array}{c}{S}_J[s]\approx -\mathrm{\Delta }\tau {\displaystyle \sum _{\ell }[J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}{\displaystyle \sum _{⟨ij⟩}{s}_{i,\ell }{s}_{j,\ell }+J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}{\displaystyle \sum _{⟨⟨ij⟩⟩}{s}_{i,\ell }{s}_{j,\ell }],}}}\end{array}

where ⟨ij⟩ and ⟨⟨ij⟩⟩ denote nearest-neighbor and diagonal pairs on the plaquette, respectively. The tunneling term \mathrm{\Delta }{\sum }_i{\sigma }_x^{(i)} produces nearest-neighbor couplings along the imaginary-time direction,

(C14)\begin{array}{c}{S}_{\mathrm{\Delta }}[s]=-{\displaystyle \sum _{i,\ell }{K}_{\mathrm{\Delta }}{s}_{i,\ell }{s}_{i,\ell +1},}\end{array}

where K_{\mathrm{\Delta }} is a known function of \mathrm{\Delta }\mathrm{\Delta }\tau arising from the Trotter decomposition of the transverse-field term and is strictly increasing with \mathrm{\Delta } for fixed \mathrm{\Delta }\tau.

The bath contribution (C7) becomes, upon discretization,

(C15)\begin{array}{c}{S}_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}[s]\simeq {\displaystyle \frac{1}{2}}{\displaystyle \sum _i{\displaystyle \sum _{\ell \ne {\ell }^{\prime }}{\tilde{K}}_{\ell -{\ell }^{\prime }}[s_{i,\ell }-s_{i,{\ell }^{\prime }}{]}^2,}}\end{array}

with

(C16)\begin{array}{c}{\tilde{K}}_{\ell -{\ell }^{\prime }}\approx K({\tau }_{\ell }-{\tau }_{{\ell }^{\prime }})\mathrm{\Delta }{\tau }^2\propto {\displaystyle \frac{\alpha }{|\ell -{\ell }^{\prime }{|}^{1+s}}}\text{for large}|\ell -{\ell }^{\prime }|.\end{array}

The effective action (C6) thus takes the schematic form

(C17)\begin{array}{c}\begin{array}{rl}{S}_{\mathrm{e}\mathrm{f}\mathrm{f}}[s]=& S_J[s]+S_{\mathrm{\Delta }}[s]+S_{\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{h}}[s]\\ \approx & -\sum _{\ell }({\tilde{J}}_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}\sum _{⟨ij⟩}{s}_{i,\ell }{s}_{j,\ell }+{\tilde{J}}_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}\sum _{⟨⟨ij⟩⟩}{s}_{i,\ell }{s}_{j,\ell })\\ & -\sum _{i,\ell }{K}_{\mathrm{\Delta }}{s}_{i,\ell }{s}_{i,\ell +1}+\frac{1}{2}\sum _{i,\ell \ne {\ell }^{\prime }}{\tilde{K}}_{\ell -{\ell }^{\prime }}(s_{i,\ell }-s_{i,{\ell }^{\prime }}{)}^2,\end{array}\end{array}

where {\tilde{J}}_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}}=\mathrm{\Delta }\tau J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} and {\tilde{J}}_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}=\mathrm{\Delta }\tau J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}} are the discretized spatial couplings at each Trotter slice.

Equation (C7) shows that the original quantum four-spin-boson model is mapped onto a classical Ising model defined on a two-dimensional lattice spanned by the spin index i=\mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D} and the imaginary-time slice index \ell =0,\cdots ,N_{\tau }-1 [2, 3]. The spatial direction (i-index) carries nearest-neighbor and diagonal couplings inherited from J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}} and J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}, while the imaginary-time direction carries both short-range couplings (from \mathrm{\Delta }) and long-range power-law couplings (from the bath kernel {\tilde{K}}_{\ell -{\ell }^{\prime }}). The partition function can be written as

(C18)\begin{array}{c}Z(\beta )={\displaystyle \sum _{\{s_{i,\ell }=\pm 1\}}\exp [-S_{\mathrm{e}\mathrm{f}\mathrm{f}}(\{s_{i,\ell }\})],}\end{array}

which is the standard Boltzmann weight of a 2D classical spin system.

In the zero-temperature limit (\beta \to \mathrm{\infty }), and in particular in the strong-coupling (worldline-freezing) regime where rare spin flips are suppressed, the dominant contributions to Z are controlled by spin paths that minimize the effective action S_{\mathrm{e}\mathrm{f}\mathrm{f}}[s] to leading order. The bath-induced long-range kernel plays the role of an emergent interaction that links identical spins across different Trotter slices. This mechanism explains why independent baths can stabilize spin patterns that closely resemble frustrated configurations on genuine two-dimensional lattices: the temporal stacking effectively endows the four-spin unit cell with extended correlations. Consequently, phases such as stripe-like or checkerboard arrangements arise naturally as classical minima of the effective action in sub-Ohmic environments.

The same classical description also clarifies the behavior of reduced density matrices and entanglement diagnostics in the strong-coupling limit. For \alpha \to \mathrm{\infty } the bath kernel pins the Ising paths to almost time-independent values s_{i,\ell }\simeq s_i=\pm 1, so that the low-energy sector is spanned by classical plaquette configurations \{|s_{\mathrm{A}}{s}_{\mathrm{B}}{s}_{\mathrm{C}}{s}_{\mathrm{D}}⟩\} that minimize the classical plaquette energy E_{\mathrm{c}\mathrm{l}}(\{s_i\}). Strong dephasing from the independent baths drives the reduced four-spin density matrix {\rho }_S towards a form that is, to leading order, diagonal in this classical basis. Even if a small transverse field \mathrm{\Delta } admixes weak coherent superpositions within (nearly) degenerate classical manifolds, strong dephasing from the independent baths suppresses the corresponding off-diagonal elements in {\rho }_S. As a result, the two-spin reduced density matrices approach a form that is diagonal in the {\sigma }^z product basis in the localized (freezing) limit, implying asymptotically vanishing Wootters concurrences for fixed (J_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}},J_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}}) as \alpha \to \mathrm{\infty }.

References

[1] J. Rogel-Salazar, Quantum Phase Transitions, 2nd Ed., Cambridge University Press, 2012

[2] M. Vojta , Quantum phase transitions, Rep. Prog. Phys. 66(12), 2069 (2003)

[3] S. Sachdev , Quantum phase transitions, Phys. World 12(4), 33 (1999)

[4] Y. J. Yan , J. Jin , R. X. Xu , and X. Zheng , Dissipation equation of motion approach to open quantum systems, Front. Phys. (Beijing) 11(4), 110306 (2016)

[5] 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), 1 (1987)

[6] M. Vojta , N. H. Tong , and R. Bulla , Quantum phase transitions in the sub-ohmic spin-boson model: Failure of the quantum-classical mapping, Phys. Rev. Lett. 94(7), 070604 (2005)

[7] R. Bulla , N. H. Tong , and M. Vojta , Numerical renormalization group for bosonic systems and application to the sub-ohmic spin-boson model, Phys. Rev. Lett. 91(17), 170601 (2003)

[8] 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(3), 030601 (2009)

[9] Q. H. Chen , Y. Y. Zhang , T. Liu , and K. L. Wang , Numerically exact solution to the finite-size Dicke model, Phys. Rev. A 78(5), 051801 (2008)

[10] D. F. Padilla , H. Pu , G. J. Cheng , and Y. Y. Zhang , Understanding the quantum Rabi ring using analogies to quantum magnetism, Phys. Rev. Lett. 129(18), 183602 (2022)

[11] Y. Xu , F. X. Sun , Q. He , H. Pu , and W. Zhang , Quantum phase transition in a quantum Rabi square with next-nearest-neighbor hopping, Phys. Rev. A 110(2), 023702 (2024)

[12] F. Verstraete , V. Murg , and J. I. Cirac , Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Adv. Phys. 57(2), 143 (2008)

[13] A. Weichselbaum , F. Verstraete , U. Schollwöck , J. I. Cirac , and J. von Delft , Variational matrix-product-state approach to quantum impurity models, Phys. Rev. B 80(16), 165117 (2009)

[14] C. Guo , A. Weichselbaum , J. von Delft , and M. Vojta , Critical and strong-coupling phases in one- and two-bath spin-boson models, Phys. Rev. Lett. 108(16), 160401 (2012)

[15] A. J. Dunnett and A. W. Chin , Efficient bond-adaptive approach for finite-temperature open quantum dynamics using the one-site time-dependent variational principle for matrix product states, Phys. Rev. B 104(21), 214302 (2021)

[16] A. W. Chin , Á. Rivas , S. F. Huelga , and M. B. Plenio , Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials, J. Math. Phys. 51(9), 092109 (2010)

[17] F. A. Y. N. Schröder and A. W. Chin , Simulating open quantum dynamics with time-dependent variational matrix product states: Towards microscopic correlation of environment dynamics and reduced system evolution, Phys. Rev. B 93(7), 075105 (2016)

[18] W. K. Wootters , Entanglement of formation of an arbitrary state of two qubits, Phys. Rev. Lett. 80(10), 2245 (1998)

[19] U. Weiss, Quantum Dissipative Systems, World Scientific, 2012

[20] Y. A. Yan and J. Shao , Stochastic description of quantum Brownian dynamics, Front. Phys. (Beijing) 11(4), 110309 (2016)

RIGHTS & PERMISSIONS

Higher Education Press