# Supplemental Material for “Bath-engineering magnetic order in quantum spin
chains: An analytic mapping approach”

Brett Min<EMAIL_ADDRESS>Department of Physics and Centre for
Quantum Information and Quantum Control, University of Toronto, 60 Saint
George St., Toronto, Ontario, M5S 1A7, Canada  Nicholas Anto-Sztrikacs
Department of Physics and Centre for Quantum Information and Quantum Control,
University of Toronto, 60 Saint George St., Toronto, Ontario, M5S 1A7, Canada
Marlon Brenes  Department of Physics and Centre for Quantum Information and
Quantum Control, University of Toronto, 60 Saint George St., Toronto, Ontario,
M5S 1A7, Canada  Dvira Segal<EMAIL_ADDRESS>Department of Chemistry,
University of Toronto, 80 Saint George St., Toronto, Ontario, M5S 3H6, Canada
Department of Physics and Centre for Quantum Information and Quantum Control,
University of Toronto, 60 Saint George St., Toronto, Ontario, M5S 1A7, Canada

(August 28, 2024)

This Supplemental Material includes the following: In Sec. S1, we provide the
details on the effective Hamiltonian mapping (EFFH) technique and its utility
in studying the equilibrium state of a many-body system strongly coupled to
its environment. In Sec. S2, we go through a detailed derivation of extracting
the equilibrium state of a dissipative spin chain via this mapping technique.
Importantly, we compare results of the EFFH methods to the more accurate
numerical reaction coordinate technique, and validate the EFFH predictions. In
Sec. S3, we provide an alternative approach of obtaining the equilibrium state
of the dissipative spin chain via the general polaron transform. We compare
the two mapping approaches, the EFFH and the polaron, in Sec. S4. Finally, in
Sec. S5, we apply both the EFFH and the general polaron technique to the Ising
chain to clearly visualize the preferred magnetic order in various coupling
schemes depicted in the main text.

##  S1 Effective Hamiltonian Mapping

In this Section, we show how to transform a Hamiltonian into its “effective”
Hamiltonian model, such that in the new picture, the interaction of the system
with the surrounding bath(s) can be made weaker than it was in the original
model. Remarkably, this transformation, which is not exact, immediately
exposes the impact of system-bath couplings on the system in the form of
renormalizing and mixing parameters, and in generating bath-mediated
couplings.

The approach described here builds on Ref. Anto-Sztrikacs  _et al._ (2023a,
b), where it was exercised on impurity models only such as the spin-boson
model, a three level absorption refrigerator and a double quantum dot
thermoelectric power generator. It is enacted here for the first time on spin
chain models, and while varying the ‘range’ of the bath(s), whether being
global, local or of an intermediate locality range.

The EFFH method builds on a reaction coordinate mapping, extracting a
collective coordinate from the bath, followed by a polaron rotation of that
mode and its truncation. In Sec. S3 we discuss an alternative mapping approach
that enacts polaron transform directly on all modes in the bath.

Consider a many-body system described by the Hamiltonian $\hat{H}_{S}$ coupled
to a structured bosonic environment, modelled by an infinite set of harmonic
modes. The total Hamiltonian of the system, environment and their interaction
reads

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{I}$ (S1) $\displaystyle=$
$\displaystyle\hat{H}_{S}+\sum_{k}\nu_{k}\hat{c}_{k}^{\dagger}\hat{c_{k}}+\hat{S}\sum_{k}t_{k}\left(\hat{c}_{k}^{\dagger}+\hat{c}_{k}\right),$

where the set $\\{\hat{c}_{k}\\}$ are canonical bosonic operators, and
$\nu_{k}$ and $t_{k}$ are the frequencies and the environment-system couplings
between the $k$-th mode and the many-body system. $\hat{S}$ is an operator
with support over system degrees of freedom that dictates the nature of the
effect of the environment onto the system. Furthermore, the effect of the
environment is typically described via spectral density function, which may be
defined as $K(\omega)=\sum_{k}t^{2}_{k}\delta(\omega-\nu_{k})$.

A reaction-coordinate (RC) mapping Iles-Smith  _et al._ (2014); Anto-Sztrikacs
and Segal (2021) of Eq. (S1) leads to

$\displaystyle\hat{H}_{\rm{RC}}=\hat{H}_{S}+\Omega\hat{a}^{\dagger}\hat{a}+\lambda\hat{S}\left(\hat{a}^{\dagger}+\hat{a}\right)+\sum_{k}f_{k}\left(\hat{a}^{\dagger}+\hat{a}\right)\left(\hat{b}_{k}^{\dagger}+\hat{b}_{k}\right)+\sum_{k}\omega_{k}\hat{b}_{k}^{\dagger}\hat{b}_{k},$
(S2)

where $\\{\hat{a}\\}$ is a canonical bosonic operator for the reaction
coordinate with frequency $\Omega$. The parameter $\lambda$ now describes the
coupling between the original system $\hat{H}_{S}$ and the reaction coordinate
via the system operator $\hat{S}$. The enlarged system comprising
$\hat{H}_{S}$ and the reaction coordinate are now coupled to a residual bath
with the effect described by the new spectral function $K^{\rm
RC}(\omega)=\sum_{k}f^{2}_{k}\delta(\omega-\omega_{k})$. The harmonic modes of
the residual reservoir are described via the canonical operators
$\\{\hat{b}_{k}\\}$ with frequencies $\omega_{k}$, which may be expressed as
linear combinations of the original harmonic modes $\\{\hat{c}_{k}\\}$. We
note that both $\lambda$ and $\Omega$ follow from the original spectral
density function, via
$\lambda^{2}=\frac{1}{\Omega}\int_{0}^{\infty}{\rm{d}}\omega\;\omega
K(\omega)$ and
$\Omega^{2}=\frac{\int_{0}^{\infty}{\rm{d}}\omega\;\omega^{3}K(\omega)}{\int_{0}^{\infty}{\rm{d}}\omega\;\omega
K(\omega)}$ Nazir and Schaller (2018). Furthermore, the spectral density of
the redefine-residual bath is
$K^{\text{RC}}(\omega)=\frac{2\pi\lambda^{2}K(\omega)}{[\mathcal{P}\int\frac{K(\omega^{\prime})}{\omega^{\prime}-\omega}d\omega^{\prime}]^{2}+\pi^{2}K^{2}(\omega)}$
with $\mathcal{P}$ indicating a principal-value integral Nazir and Schaller
(2018); Iles-Smith  _et al._ (2014); Anto-Sztrikacs and Segal (2021).

Next, we employ the polaron transformation on the RC mode,
$\hat{\tilde{H}}=\hat{U}\hat{H}_{\rm{RC}}\hat{U}^{\dagger}$ with
$\hat{U}={\rm{exp}}\left[\frac{\lambda}{\Omega}(\hat{a}^{\dagger}-\hat{a})\hat{S}\right]$,
such that

$\displaystyle\hat{\tilde{H}}=\hat{\tilde{H}}_{S}+\Omega\hat{a}^{\dagger}\hat{a}-\frac{\lambda^{2}}{\Omega}\hat{S}^{2}+\sum_{k}f_{k}\left(\hat{a}^{\dagger}+\hat{a}-\frac{2\lambda}{\Omega}\hat{S}\right)\left(\hat{b}_{k}^{\dagger}+\hat{b}_{k}\right)+\sum_{k}\omega_{k}\hat{b}_{k}^{\dagger}\hat{b}_{k}.$
(S3)

Note that the polaron transform was enacted on a single mode, the reaction
coordinate (which represents a collective mode of the original bath) and the
system Hamiltonian.

We now perform an approximation by considering only the ground-state level of
the reaction coordinate, by projecting $\hat{\tilde{H}}$ onto this manifold
via the projector $\hat{Q}_{0}=|0\rangle\langle 0|$, to obtain Anto-Sztrikacs
_et al._ (2023a)

$\displaystyle\hat{H}^{\rm eff}$
$\displaystyle=\hat{Q}_{0}\hat{\tilde{H}}_{S}\hat{Q}_{0}-\frac{\lambda^{2}}{\Omega}\hat{S}^{2}-\hat{S}\sum_{k}\frac{2\lambda
f_{k}}{\Omega}\left(\hat{b}_{k}^{\dagger}+\hat{b}_{k}\right)+\sum_{k}\omega_{k}\hat{b}_{k}^{\dagger}\hat{b}_{k}.$
(S4)

Crucially, Eq. (S4) has the same form as the original total Hamiltonian
$\hat{H}$ from Eq. (S1) with the new effective system Hamiltonian defined as

$\displaystyle\hat{H}^{\rm
eff}_{S}=\hat{Q}_{0}\hat{\tilde{H}}_{S}\hat{Q}_{0}-\frac{\lambda^{2}}{\Omega}\hat{S}^{2}.$
(S5)

We also provide a useful identity that allows a quick extraction of the
effective system Hamiltonian, generalizing Ref. Anto-Sztrikacs  _et al._
(2023a):

$\hat{H}^{\text{eff}}_{S}=e^{-(\lambda^{2}/2\Omega^{2})\hat{S}^{2}}\left(\sum^{\infty}_{n=0}\frac{\lambda^{2n}}{\Omega^{2n}n!}\hat{S}^{n}\hat{H}_{S}\hat{S}^{n}\right)e^{-(\lambda^{2}/2\Omega^{2})\hat{S}^{2}}-\frac{\lambda^{2}}{\Omega}\hat{S}^{2}.$
(S6)

Let us now discuss some key aspects of the transformed Hamiltonian described
in Eq. (S4). First, note the coupling of the system to the modes of the
residual bath is re-scaled $f_{k}\to 2\lambda f_{k}/\Omega$. If this residual
coupling can be made a perturbative parameter, strong coupling effects can be
studied as a function of $\lambda$. This constitutes a specific form of a
Markovian embedding, where strong-coupling effects become embedded onto the
effective system Hamiltonian. Furthermore, the spectral density of the
effective bath will be further dressed by the RC parameters:
$K^{\text{eff}}(\omega)=\frac{4\lambda^{2}}{\Omega^{2}}K^{\text{RC}}(\omega)$.
Worthy of note, for the case of an open system initially coupled to a bosonic
bath with Brownian form,

$\displaystyle
K(\omega)=\frac{4\gamma\Omega^{2}\lambda^{2}\omega}{(\omega^{2}-\Omega^{2})^{2}+(2\pi\gamma\Omega\omega)^{2}},$
(S7)

the effective spectral density function maps to
$K^{\rm{eff}}(\omega)=\frac{4\lambda^{2}}{\Omega^{2}}\gamma\omega$, where
$\gamma$ initially played the role of a width parameter, now represents the
residual coupling strength between the residual bath and open system. However,
the mapping is not limited to the Brownian function. In fact, one can treat
Eq. (S6) as a general effective Hamiltonian emerging from Eq. (S1) with
$\Omega$ as the frequency of the bath and $\lambda$ as a measure for the
system-bath coupling energy. From Eq. (S6) we note that strong coupling
effects (large $\lambda$) are now embedded in the effective system
Hamiltonian, while the coupling to the residual bath can be made weak, by
assuming $\gamma\ll 1$ in the above mentioned Brownian example. We also note
that the coupling operator $\hat{S}$ has intricate consequences in relation to
environmental effects. Indeed, while thermodynamic equilibration between a
system and environment is a phenomenon that occurs irrespective of the
microscopic details of the nature of the coupling between each subsystem; at
strong coupling, on the other hand, it has been established that microscopic
details are of the essence Cresser and Anders (2021); Chiu  _et al._ (2022).

The effective Hamiltonian treatment allows one to understand equilibrium
thermodynamics in many-body systems coupled to thermal environments. Through
the Markovian embedding of the reaction-coordinate mapping, we conjecture that
the equilibrium state of the many-body system is a Gibbs state of the
effective Hamiltonian

$\displaystyle\hat{\rho}^{\rm eff}_{S}=e^{-\beta\hat{H}^{\rm eff}_{S}}/Z^{\rm
eff}.$ (S8)

Here, $Z^{\rm eff}=\text{Tr}\left[e^{-\beta\hat{H}^{\rm eff}_{S}}\right]$ is
the partition function of the effective system. It should be noted that at
sufficiently strong $\lambda$, this state does not coincide with the thermal
state of the isolated system, $\hat{H}_{S}$, and strong-coupling effects may
be observed in the many-body equilibrium state. Indeed, $\hat{H}^{\rm
eff}_{S}$ already contains terms that depend on $\lambda$ that generate
interactions among all the spins composing the many-body system. The
equilibrium state in Eq. (S8) was analyzed in Refs. Anto-Sztrikacs  _et al._
(2023b, a) for the spin-boson model, where we demonstrated that it provided an
excellent approximation for equilibrium expectation values from the weak to
the strong coupling limit. Moreover, it can be shown that the equilibrium
state of the EFFH method is an exact result for the spin boson model in the
ultrastrong coupling regime Anto-Sztrikacs  _et al._ (2023a).

##  S2 Application of the EFFH mapping to the Dissipative Heisenberg chain

In this Section, we apply the effective Hamiltonian mapping as described in
Sec. S1 to the general dissipative spin-chain $\hat{H}_{S}$ given by the
following Hamiltonian:

$\hat{H}_{S}=\sum_{i=1}^{N}\Delta_{i}\hat{\sigma}^{z}_{i}+\sum_{\alpha}\sum^{N-1}_{i=1}J_{\alpha}\hat{\sigma}^{\alpha}_{i}\hat{\sigma}^{\alpha}_{i+1},$
(S9)

where $\Delta_{i}>0$ represents the spin splitting of the $i$th qubit,
$J_{\alpha}>0$ is the interaction strength between neighbouring spins in the
$\alpha=\\{x,y,z\\}$ direction. In the following subsections, we will examine
the structure of the equilibrium state under four different dissipative
coupling schemes which corresponds to Fully-global S2.1, Fully-local S2.2,
Half-and-half S2.3, and Pairwise S2.4 coupling schemes.

As for notation, we denote the original spin sites by $i$ and $j$, and index
baths with $n$. However, once the mapping is operated and depending on the
type of the bath, we use notation as convenient and clarify it in the text.

In the global model only one bath couples to the system. In the fully local
model, the number of independent baths equals the number of sites. In the
half-and-half model, the chain couples to two independent baths, while in the
pairwise model there are $N/2$ independent baths, with $N$ spins in the chain.
These models are depicted in Fig. 1 in the Main text.

###  S2.1 Fully-global coupling model and Benchmarking

We imagine the entire spin chain to be coupled to a single bosonic bath
described by the following Hamiltonian,

$\hat{H}_{\text{glob}}=\hat{H}_{S}+\hat{S}_{\text{glob}}\sum_{k}t_{k}\left(\hat{c}^{\dagger}_{k}+\hat{c}_{k}\right)+\sum_{k}\nu_{k}\hat{c}^{\dagger}_{k}\hat{c}_{k}.$
(S10)

Here, $\hat{H}_{S}$ is given by Eq. (S9) and
$\hat{S}_{\text{glob}}=\sum^{N}_{i=1}\hat{\sigma}^{x}_{i}$; $i$ is an index of
spin sites in the chain. We apply the EFFH framework to the global spin-chain
following the framework outlines in Sec. S1, extracting a single reaction
coordinate from the global bath, polaron transforming it, then truncting it. A
crucial point in the derivation is that
$\hat{S}_{\text{glob}}^{2}=N\hat{I}+2\sum_{i<j}\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{j}$,
which implies that we generate new interactions between spins in the system
due to the influence of the global reservoir. Following through the steps
outlined in the previous section S1, we arrive at the following effective
system Hamiltonian,

$\displaystyle\hat{H}^{\rm
eff}_{\text{glob},S}=\sum_{i=1}^{N}\tilde{\Delta}_{i}\hat{\sigma}^{i}_{z}+\sum_{i=1}^{N-1}(J_{x}\hat{\sigma}^{x}_{i}\hat{\sigma}_{i+1}^{x}+\tilde{J}_{y}\hat{\sigma}_{i}^{y}\hat{\sigma}_{i+1}^{y}+\tilde{J}_{z}\hat{\sigma}^{z}_{i}\hat{\sigma}^{z}_{i+1})-\frac{\lambda^{2}}{\Omega}N\hat{I}-\frac{2\lambda^{2}}{\Omega}\sum_{i<j}\hat{\sigma}_{i}^{x}\hat{\sigma}_{j}^{x},$
(S11)

where, $\tilde{\Delta}_{i}=\Delta_{i}e^{-\frac{2\lambda^{2}}{\Omega^{2}}}$ and

$\displaystyle\tilde{J}_{y}$
$\displaystyle=\frac{J_{y}}{2}\left(1+e^{-8\lambda^{2}/\Omega^{2}}\right)+\frac{J_{z}}{2}\left(1-e^{-8\lambda^{2}/\Omega^{2}}\right),$
$\displaystyle\tilde{J}_{z}$
$\displaystyle=\frac{J_{z}}{2}\left(1+e^{-8\lambda^{2}/\Omega^{2}}\right)+\frac{J_{y}}{2}\left(1-e^{-8\lambda^{2}/\Omega^{2}}\right).$
(S12)

After the mapping, the effective system weakly couples to the residual bath,
thus we can calculate equilibrium expectation values using the equilibrium
state in Eq. (S8). We do not write here explicitly the residual bath’s
Hamiltonian and the EFFH system-bath interaction term, which are included in
Eq. (S4), since they do not play a role in the equilibrium state of the
effective model.

Before showing any calculations, already at the level of the Hamiltonian
described in Eq. (S11), we can see that the spins should manifest a crossover
in the $x$ direction with increasing $\lambda$, from the anti-ferromagnetic
(AFM) alignment of spins due to $J_{x}>0$, to a FM order. This transition
takes place due to the last negative term in the effective system’s
Hamiltonian, Eq. (S11). Importantly, since this long-range term scales as
$N^{2}$, with $N$ the total number of spins, the value of $\lambda$ at which
the AFM to FM transition occurs will keep shifting to lower values with
increasing the size $N$. This is in contrast to the fully-connected chain
(with $J_{x}$ connecting all spins) where a fixed critical point takes place,
at a specific $\lambda$. The fully-connected model is described in the main
text.

To validate these predictions, based on the EFFH equilibrium state
approximation, in Figure S1, we show a comparison between the structure factor
$S_{\alpha}=\frac{1}{N^{2}}\sum_{ij}\langle\hat{\sigma}^{\alpha}_{i}\hat{\sigma}^{\alpha}_{j}\rangle$
for $\alpha=x,y,z$ of a general Heisenberg $N=10$-site chain computed from the
EFFH and the numerically accurate reaction coordinate (RC) method Iles-Smith
_et al._ (2014); Anto-Sztrikacs and Segal (2021). We choose generic
parameters, $\Delta_{i}=0$, $J_{x}=0.77$, $J_{y}=1.23$, $J_{z}=0.89$, and
$\Omega=8$. For RC simulations, we used a Brownian bath for the original
Hamiltonian, which is converted to an Ohmic bath in the RC picture.
Importantly, we observe in Figure S1 an excellent agreement between the
equilibrium state EFFH and the RC numerical method. We also find that the
accuracy of the EFFH method becomes better as $N$ increases. Note that
prediction of the RC numerical method were benchmarked against numerically
exact Hierarchical Equation of Motion simulations in Refs. Iles-Smith  _et
al._ (2016); Anto-Sztrikacs  _et al._ (2023b). In the first row of Fig. S3 we
further present a heat map of this crossover using a super-Ohmic spectral
function as an example, with details explained in Sec. S5.

Figure S1: Comparison of the expectation values of the spin structure factor
$S_{\alpha}=\frac{1}{N^{2}}\sum_{ij}\langle\hat{\sigma}^{\alpha}_{i}\hat{\sigma}^{\alpha}_{j}\rangle$
for $\alpha=\\{x,y,z\\}$, computed from both the numerical $\hat{\rho}^{\rm
RC}_{S}$ and the analytical $\hat{\rho}^{\rm eff}_{S}$ for different system
sizes $N$ of the anisotropic Heisenberg chain. The parameters of the
calculation are $J_{x}=0.77$, $J_{y}=1.23$, $J_{z}=0.89$ and $\Omega=8$. The
main panels display the results for $T=1$, while the insets show results for
$T=0.1$. With the increasing coupling strength $\lambda$, we find a crossover
from $S_{\alpha}=0$ to $S_{\alpha}\neq 0$, that becomes more steeper at lower
temperature. The limiting value for $S_{x}$ becomes $1$ as the entire chain
prefers a ferromagnetic alignment in the $x$ direction while the limiting
value for both $S_{y}$ and $S_{z}$ tends to $1/N$ as the entire spin-spin
correlations in both $y$ and $z$ directions are completely suppressed except
the autocorrelators. We also note that the EFFH results show better agreement
with the numerically accurate RC method at large $N$.

###  S2.2 Fully-local coupling model

Next, we move on to a fully-local coupling scheme where each spin couples to
their own independent bosonic bath. The total Hamiltonian of the model is

$\hat{H}_{\text{local}}=\hat{H}_{S}+\sum^{N}_{n=1}\hat{S}_{\text{local,n}}\sum_{k}t_{n,k}\left(\hat{c}^{\dagger}_{n,k}+\hat{c}_{n,k}\right)+\sum_{n,k}\nu_{n,k}\hat{c}^{\dagger}_{n,k}\hat{c}_{n,k}.$
(S13)

Here, $\hat{H}_{S}$ is again given by Eq. (S9) while
$\hat{S}_{\text{local},n}=\hat{\sigma}^{x}_{n}$. That is, unlike the fully-
global coupling scheme studied in the previous subsection S2.1,
$\hat{S}_{\text{local},n}$ acts on the $n$th individual site. Here, the index
$n$ that counts the baths is identical to the site index $i$. Upon proceeding
with the EFFH mapping, we extract a reaction coordinate from each independent
bath, polaron transform each of these collective modes and truncate them. The
results is that there are no bath-induced spin-spin correlations via the
$\hat{S}^{2}_{\text{local},n}$ term. The effective system Hamiltonian instead
transforms to

$\displaystyle\hat{H}^{\rm{eff}}_{\text{local},S}=$
$\displaystyle\sum_{n=1}^{N}\tilde{\Delta}_{n}\hat{\sigma}^{z}_{n}+\sum_{i=1}^{N-1}\left[J_{x}\hat{\sigma}^{x}_{n}\hat{\sigma}^{x}_{n+1}+\tilde{J}_{y}^{(n,n+1)}\hat{\sigma}^{y}_{n}\hat{\sigma}^{y}_{n+1}+\tilde{J}_{z}^{(n,n+1)}\hat{\sigma}^{z}_{n}\hat{\sigma}^{z}_{n+1}\right],$
(S14)

where now the transformed parameters map as
$\tilde{\Delta}_{n}=\Delta_{n}e^{-\frac{2\lambda_{n}^{2}}{\Omega_{n}^{2}}}$,
$\tilde{J}_{y}^{(n,n+1)}=J_{y}e^{-\frac{2\lambda_{n}^{2}}{\Omega_{n}^{2}}}e^{-\frac{2\lambda_{n+1}^{2}}{\Omega_{n+1}^{2}}}$,
and
$\tilde{J}_{z}^{(n,n+1)}=J_{z}e^{-\frac{2\lambda_{n}^{2}}{\Omega_{n}^{2}}}e^{-\frac{2\lambda_{n+1}^{2}}{\Omega_{n+1}^{2}}}$,
while $J_{x}$ stays intact. We comment that the interactions with each
reservoir can be made distinct, which offers flexibility in engineering the
effective system, as the values of the parameters on each site can be tuned by
the coupling.

Furthermore, in contrast to the global scheme, there are only nearest-
neighbour effects within this local bath picture. As the coupling to the baths
increases, we expect the model to go from a Heisenberg chain to a quantum
Ising model. This is due to the suppression of both of the $J_{y}$ and $J_{z}$
interaction terms with increasing coupling to the bath. In the second row of
Fig. S3 we present a heat map of this crossover, with details explained in
Sec. S5.

###  S2.3 Half-and-half coupling model

We now demonstrate that by tuning the locality of the baths we can engineer
different magnetic orderings in a spin chain. One simple example is a chain of
$N$ sites, in which the left-half of the chain couples to one bath and the
right-half of the chain couples to a separate bath. This model is described by
the following Hamiltonian

$\hat{H}_{\text{half}}=\hat{H}_{S}+\sum_{n}\hat{S}_{\text{half},n}\sum_{k}t_{n,k}\left(\hat{c}^{\dagger}_{n,k}+\hat{c}_{n,k}\right)+\sum_{n,k}\nu_{n,k}\hat{c}^{\dagger}_{n,k}\hat{c}_{n,k},$
(S15)

where $\hat{H}_{S}$ is again given by Eq. (S9), while $n=\\{L,R\\}$, with
$\hat{S}_{\text{half},L}=\sum^{N/2}_{i=1}\hat{\sigma}^{x}_{i}$, and
$\hat{S}_{\text{half},R}=\sum^{N}_{N/2+1}\hat{\sigma}^{x}_{i}$. Crucially, now
there are two different coupling depending on which half of the chain the spin
is coupled to. As a result, the EFFH approach yields an intermediate scheme
between the local and global approaches. For instance, in analogy with the
global case, we see the generation of bath-induced long-ranged spin-spin
interaction term via
$\hat{S}_{n}^{2}=\frac{N}{2}\hat{I}+2\sum_{i<j}\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{j}$
(where indices $i$ and $j$ should run appropriately depending on $n$, whether
it is the left or right half). On top of it, we expect a suppression of both
$J_{y}$ and $J_{z}$ spin-spin interaction for the boundary spins at the right
(left)-end of the left (right) half of the chain.

Technically, we extract a reaction coordinate mode from each of the two baths.
The polaron transformation involved in the mapping is made of two
transformations that commute with each other,

$\displaystyle\hat{U}={\rm{exp}}\left[\frac{\lambda_{L}}{\Omega_{L}}(\hat{a}_{L}^{\dagger}-\hat{a}_{L})\sum_{i=1}^{N/2}\hat{\sigma}^{x}_{i}\right]{\rm{exp}}\left[\frac{\lambda_{R}}{\Omega_{R}}(\hat{a}_{R}^{\dagger}-\hat{a}_{R})\sum_{i=N/2}^{N}\hat{\sigma}^{x}_{i}\right]=\hat{U}_{L}\hat{U}_{R},$
(S16)

where $\hat{a}^{\dagger}_{n}\leavevmode\nobreak\ (\hat{a}_{n})$ corresponds to
the creation (annihilation) operator of the RC extracted from left ($n=L$) or
right $(n=R)$ bath. We emphasize the necessity of
$[\hat{U}_{L},\hat{U}_{R}]=0$, which allows the EFFH method to proceed with
minor modifications, similarly to the previous two examples. We note that the
transformation for sites $i=1$ to $i=N/2-1$ and $i=N/2+2$ to $i=N$ will follow
analogously to the global bath, while terms at boundary, sites $i=N/2$ and
$i=N/2+1$ yield a transformation akin to the local bath. Using the results
from the previous two sections, we find that the effective Hamiltonian of the
system can be decomposed into a portion for the left bath (spin $1$ to $N/2$),
the right bath (spin $N/2+1$ to $N$) and the boundary (spin $N/2$ and $N/2+1$)
as

$\displaystyle\hat{H}^{\rm eff}_{\text{half},S}=\hat{H}^{\rm
eff}_{L,S}+\hat{H}^{\rm eff}_{\text{bound},S}+\hat{H}^{\rm eff}_{R,S}.$ (S17)

It follows that from the mapping, for each partition of the Hamiltonian we
obtain

$\displaystyle\hat{H}^{\rm
eff}_{\text{half},L}=\sum_{i=1}^{N/2}\tilde{\Delta}_{i}\hat{\sigma}^{z}_{i}+\sum_{i=1}^{N/2-1}J_{x}\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{i+1}+\tilde{J}_{y}\hat{\sigma}^{y}_{i}\hat{\sigma}^{y}_{i+1}+\tilde{J}_{z}\hat{\sigma}^{z}_{i}\hat{\sigma}^{z}_{i+1}-\frac{N}{2}\frac{\lambda_{L}^{2}}{\Omega_{L}}\hat{I}-\frac{2\lambda_{L}^{2}}{\Omega_{L}}\sum_{i<j=1}^{N/2}\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{j}.$
(S18)

We note that the spin gets suppress as
$\tilde{\Delta}_{i}=\Delta_{i}e^{-\frac{2\lambda_{L}^{2}}{\Omega_{L}^{2}}}$,
and the interactions via
$\tilde{J}_{y(z)}=\frac{J_{y(z)}}{2}\left(1+e^{-8\lambda_{L}^{2}/\Omega_{L}^{2}}\right)+\frac{J_{z(y)}}{2}\left(1-e^{-8\lambda_{L}^{2}/\Omega_{L}^{2}}\right)$.
Similarly, the right portion of the chain is modifed as

$\displaystyle\hat{H}^{\rm
eff}_{\text{half},R}=\sum_{i=N/2+1}^{N}\tilde{\Delta}_{i}\hat{\sigma}^{z}_{i}+\sum_{i=N/2+1}^{N}J_{x}\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{i+1}+\tilde{J}_{y}\hat{\sigma}^{y}_{i}\hat{\sigma}^{y}_{i+1}+\tilde{J}_{z}\hat{\sigma}^{z}_{i}\hat{\sigma}^{z}_{i+1}-\frac{N}{2}\frac{\lambda_{R}^{2}}{\Omega_{R}}\hat{I}-\frac{2\lambda_{R}^{2}}{\Omega_{R}}\sum_{i<j=N/2+1}^{N}\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{j}.$
(S19)

Similarly to the left bath, we have that the spin gets suppress as
$\tilde{\Delta}_{i}=\Delta_{i}e^{-\frac{2\lambda_{R}^{2}}{\Omega_{R}^{2}}}$,
and the interactions via
$\tilde{J}_{y(z)}=\frac{J_{y(z)}}{2}\left(1+e^{-8\lambda_{R}^{2}/\Omega_{R}^{2}}\right)+\frac{J_{z(y)}}{2}\left(1-e^{-8\lambda_{R}^{2}/\Omega_{R}^{2}}\right)$.
Lastly, the two terms at the boundary are represented, with the
renormalization already taken into account as

$\displaystyle\hat{H}^{\rm
eff}_{\text{half},\text{bound}}=J_{x}\hat{\sigma}^{x}_{N/2}\hat{\sigma}^{x}_{N/2+1}+J_{y}e^{-\frac{2\lambda_{L}^{2}}{\Omega_{L}^{2}}}e^{-\frac{2\lambda_{R}^{2}}{\Omega_{R}^{2}}}\hat{\sigma}^{y}_{N/2}\hat{\sigma}^{y}_{N/2+1}+J_{z}e^{-\frac{2\lambda_{L}^{2}}{\Omega_{L}^{2}}}e^{-\frac{2\lambda_{R}^{2}}{\Omega_{R}^{2}}}\hat{\sigma}^{z}_{N/2}\hat{\sigma}^{z}_{N/2+1}.$
(S20)

Due to the bath-induced ferromagnetic interaction, we expect to observe a
transition of spin alignments in the $x$ direction from an AFM to FM order,
for spins coupled to a common bath. In contrast, the boundary spin retains
their AFM interaction. This will lead to a domain wall structure at
sufficiently strong coupling. We present these transition, from an AFM order
to a FM with a domain wall, in the third row of Fig. S3 as discussed in Sec.
S5.

###  S2.4 Pairwise coupling model

In this last example, we couple each pairs of spins to a common bath. In the
model, we take every odd site of the chain and its nearest neighbour to the
right and couple them to the same bath. The total model is described by the
following Hamiltonian

$\hat{H}_{\text{pair}}=\hat{H}_{S}+\sum^{N/2}_{n=1}\hat{S}_{\text{pair},n}\sum_{k}t_{n,k}\left(\hat{c}^{\dagger}_{n,k}+\hat{c}_{n,k}\right)+\sum_{n,k}\nu_{n,k}\hat{c}^{\dagger}_{n,k}\hat{c}_{n,k}.$
(S21)

Here, $\hat{S}_{\text{pair},n}=\hat{\sigma}^{x}_{2n-1}+\hat{\sigma}^{x}_{2n}$;
$n\in\\{1,\dots,N/2\\}$ is the bath index.

$\displaystyle\hat{H}^{\text{eff}}_{\text{pair},S}=$
$\displaystyle\sum^{N/2}_{n=1}\left(\bar{\Delta}_{2n-1}\hat{\sigma}^{z}_{2n-1}+\bar{\Delta}_{2n}\hat{\sigma}^{z}_{2n}\right)-\sum^{N/2}_{n=1}\frac{\lambda^{2}_{n}}{\Omega_{n}}\hat{S}^{2}_{\text{pair},n}+\sum_{\alpha}\sum^{N/2}_{n=1}\left(\bar{J}_{\alpha,2n-1}\hat{\sigma}^{\alpha}_{2n-1}\hat{\sigma}^{\alpha}_{2n}+\bar{J}_{\alpha,2n}\hat{\sigma}^{\alpha}_{2n}\hat{\sigma}^{\alpha}_{2n+1}\right).$
(S22)

The new parameters are now defined as follows:
$\bar{\Delta}_{2n-1}=\Delta_{2n-1}\exp(-2\lambda^{2}_{n}/\Omega^{2}_{n})$,
$\bar{\Delta}_{2n}=\Delta_{2n}\exp(-2\lambda^{2}_{n}/\Omega^{2}_{n})$,
$\bar{J}_{x,2n-1}=\bar{J}_{x,2n}=J_{x}$,

$\displaystyle\bar{J}_{y(z),2n-1}=\frac{J_{y(z)}}{2}\left(1+e^{-\frac{8\lambda^{2}_{n}}{\Omega^{2}_{n}}}\right)+\frac{J_{z(y)}}{2}\left(1-e^{-\frac{8\lambda^{2}_{n}}{\Omega^{2}_{n}}}\right),$
(S23)
$\displaystyle\bar{J}_{y(z),2n}=J_{y(z)}\exp\left(-2\left(\frac{\lambda^{2}_{2n-1}}{\Omega^{2}_{2n-1}}+\frac{\lambda^{2}_{2n}}{\Omega^{2}_{2n}}\right)\right).$

That is, at sufficiently strong coupling, the spin-spin interactions at the
boundary of each two baths, corresponding to even site indices $(2n)$, in
$y(z)$ directions will be suppressed, the result of a local-coupling. On the
other hand we again generate FM interaction terms in the $x$ direction for
spins coupled to a common bath. Overall, we expect to observe an extended Neel
order where the anti-ferromagnetic wavelength doubles. The development of this
extended Neel state is presented in the fourth row of Fig. S3 and is further
exemplified in the main text.

##  S3 General Polaron Mapping and its application to the Dissipative
Heisenberg chain

In this section, we provide an alternative approach of obtaining the effective
Hamiltonian via the polaron transform enacted on all modes in the bath(s). The
system’s Hamiltonian is again a spin chain given by

$\hat{H}_{S}=\sum_{i=1}^{N}\Delta_{i}\hat{\sigma}^{z}_{i}+\sum_{\alpha}\sum^{N-1}_{i=1}J_{\alpha}\hat{\sigma}^{\alpha}_{i}\hat{\sigma}^{\alpha}_{i+1}.$
(S24)

where $\Delta_{i}>0$ represents the spin splitting of the $i$th qubit,
$J_{\alpha}>0$ is the interaction strength between neighbouring spins in the
$\alpha=\\{x,y,z\\}$ direction. We now derive the effective Hamiltonian using
an alternative approach to that presented in Sec. S1 for the four models
considered above, fully global S3.1, fully local S3.2, half-and-half S3.3 and
pairwise coupling S3.4. The effective Hamiltonian created by the polaron
mapping here is completely parallel to the result of the reaction-coordinate
based EFFH method of Sec. S2. We compare and discuss the pros and cons of
these two mapping methods in Sec. S4.

###  S3.1 Fully-global coupling model

For the fully-global coupling scheme, recall we work with the following total
Hamiltonian,

$\displaystyle\hat{H}_{\text{glob}}$ $\displaystyle=$
$\displaystyle\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{I}$ (S25) $\displaystyle=$
$\displaystyle\hat{H}_{S}+\sum_{k}\nu_{k}\hat{c}^{\dagger}_{k}\hat{c}_{k}+\sum^{N}_{n=1}\hat{\sigma}^{x}_{n}\sum_{k}t_{n,k}\left(\hat{c}^{\dagger}_{k}+\hat{c}_{k}\right),$

where here we use $n$ as the index for the spin sites. Note that we can
consider here a rather general case, with spins coupled at different strength
to the common bath. Following the procedure described in Ref. Anto-Sztrikacs
_et al._ (2023b), we perform a series of $n=1,...,N$ polaron transform via the
unitary

$\hat{W}_{n}=\exp(-i\hat{\sigma}^{x}_{n}\hat{B}_{n}/2)\quad\text{where}\quad\hat{B}_{n}=2i\sum_{k}\frac{g_{n,k}}{\nu_{k}}\left(\hat{c}^{\dagger}_{k}-\hat{c}_{k}\right).$
(S26)

This transformation is referred to as “full-polaron” if the variational
parameters $\\{g_{n,k}\\}$ are simply set to $\\{t_{n,k}\\}$, the original
system-reservoir couplings. If, instead, the optimal values for
$\\{g_{n,k}\\}$ are obtained by minimizing the Gibbs-Bogoliubov-Feynman upper
bound on the free energy, the transformation is called “variational” Xu and
Cao (2016). After applying consecutive $\hat{W}_{n}$, we subtract and add
appropriate bath averaged quantities to obtain the total transformed
Hamiltonian in the following form:

$\hat{H}^{\text{pol}}=\hat{H}^{\text{pol}}_{S}+\hat{H}^{\text{pol}}_{B}+\hat{H}^{\text{pol}}_{I}.$
(S27)

Here, the new system Hamiltonian $\hat{H}^{\text{pol}}_{S}$ is now given by

$\displaystyle\hat{H}^{\text{pol}}_{S}=$
$\displaystyle\sum^{N}_{n=1}E^{(n)}_{0}\hat{I}+\sum^{N}_{n=1}\Delta_{n}\langle\hat{\mathcal{C}}_{n}\rangle\hat{\sigma}^{z}_{n}+\sum^{N-1}_{n=1}\left(J_{x}-2E^{(n)}_{I}\right)\hat{\sigma}^{x}_{n}\hat{\sigma}^{x}_{n+1}-2\sum_{n+1<m}E^{(nm)}_{I}\hat{\sigma}^{x}_{n}\hat{\sigma}^{x}_{m}$
(S28) $\displaystyle+$
$\displaystyle\sum^{N-1}_{n=1}\Big{(}\left[J_{y}\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle+J_{z}\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle\right]\hat{\sigma}^{y}_{n}\hat{\sigma}^{y}_{n+1}+\left[J_{z}\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle+J_{y}\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle\right]\hat{\sigma}^{z}_{n}\hat{\sigma}^{z}_{n+1}\Big{)}$

where we define $\hat{\mathcal{C}}_{n}=\cos(\hat{B}_{n})$ and
$\hat{\mathcal{S}}_{n}=\sin(\hat{B}_{n})$ (definition for $\hat{B}_{n}$ is
given in Eq. (S26)). The averages $\langle\dots\rangle$ are done over the
unchanged bath Hamiltonian
$\hat{H}^{\text{pol}}_{B}=\sum_{k}\nu_{k}\hat{c}^{\dagger}_{k}\hat{c}_{k}$
with inverse temperature $\beta$. The renormalized system parameters now
depend on the bath spectral function and temperature. Explicit forms are given
as follows

$\displaystyle E^{(n)}_{0}=$
$\displaystyle\sum_{n,k}\frac{g_{n,k}(g_{n,k}-2t_{n,k})}{\nu_{k}}$ (S29)
$\displaystyle E^{(n)}_{I}=$
$\displaystyle\sum_{k}\left[\left(\frac{g_{n,k}(t_{n+1,k}-g_{n+1,k}/2)}{\nu_{k}}\right)+\left(\frac{g_{n+1,k}(t_{n,k}-g_{n,k}/2)}{\nu_{k}}\right)\right]$
$\displaystyle E^{(nm)}_{I}=$
$\displaystyle\sum_{k}\left[\left(\frac{g_{n,k}(t_{m,k}-g_{m,k}/2)}{\nu_{k}}\right)+\left(\frac{g_{m,k}(t_{n,k}-g_{n,k}/2)}{\nu_{k}}\right)\right]$
$\displaystyle\langle\hat{\mathcal{C}}_{n}\rangle=$
$\displaystyle\Bigg{\langle}\cos\left(2i\sum_{k}\frac{g_{n,k}}{\nu_{k}}\left(\hat{c}^{\dagger}_{k}-\hat{c}_{k}\right)\right)\Bigg{\rangle}=\exp\left(-2\sum_{k}\frac{g^{2}_{n,k}}{\nu^{2}_{k}}\coth\left(\frac{\beta\nu_{k}}{2}\right)\right)$
$\displaystyle\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle=$
$\displaystyle\Bigg{\langle}\cos\left(2i\sum_{k}\frac{g_{n,k}}{\nu_{k}}\left(\hat{c}^{\dagger}_{k}-\hat{c}_{k}\right)\right)\cos\left(2i\sum_{k}\frac{g_{n+1,k}}{\nu_{k}}\left(\hat{c}^{\dagger}_{k}-\hat{c}_{k}\right)\right)\Bigg{\rangle}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left[\exp\left(-2\sum_{k}\frac{(g_{n,k}+g_{n+1,k})^{2}}{\nu^{2}_{k}}\coth\left(\frac{\beta\nu_{k}}{2}\right)\right)+\exp\left(-2\sum_{k}\frac{(g_{n,k}-g_{n+1,k})^{2}}{\nu^{2}_{k}}\coth\left(\frac{\beta\nu_{k}}{2}\right)\right)\right]$
$\displaystyle\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle=$
$\displaystyle\Bigg{\langle}\sin(2i\sum_{k}\frac{g_{n,k}}{\nu_{k}}\left(\hat{c}^{\dagger}_{k}-\hat{c}_{k}\right))\sin\left(2i\sum_{k}\frac{g_{n+1,k}}{\nu_{k}}\left(\hat{c}^{\dagger}_{k}-\hat{c}_{k}\right)\right)\Bigg{\rangle}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\left[\exp\left(-2\sum_{k}\frac{(g_{n,k}+g_{n+1,k})^{2}}{\nu^{2}_{k}}\coth\left(\frac{\beta\nu_{k}}{2}\right)\right)-\exp\left(-2\sum_{k}\frac{(g_{n,k}-g_{n+1,k})^{2}}{\nu^{2}_{k}}\coth\left(\frac{\beta\nu_{k}}{2}\right)\right)\right]$
$\displaystyle\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{C}}_{n+1}\rangle=$
$\displaystyle\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{S}}_{n+1}\rangle=0,\,\,\,\langle\hat{\mathcal{S}}_{n}\rangle=\Bigg{\langle}\sin\left(2i\sum_{k}\frac{g_{n,k}}{\nu_{k}}\left(\hat{c}^{\dagger}_{k}-\hat{c}_{k}\right)\right)\Bigg{\rangle}=0.$

Note that $E^{(n)}_{I}$ and $E^{(nm)}_{I}$ are nearest-neighbour and long-
ranged bath induced spin-spin interactions, respectively, that we separate in
the equation for the convenience of presentation. In the weak coupling limit,
we find $E^{(n)}_{0}\rightarrow
0,\langle\hat{\mathcal{C}}_{n}\rangle\rightarrow 1,E^{(n)}_{I}\rightarrow
0,\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle\rightarrow
1,E^{(nm)}_{I}\rightarrow
0,\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle\rightarrow 0$.
Hence, we recover the correct weak-coupling limit. The new system-bath
interaction Hamiltonian is given by

$\displaystyle\hat{H}^{\text{pol}}_{I}=$
$\displaystyle\sum_{n}\Delta_{n}\Big{(}\left(\hat{\mathcal{C}}_{n}-\langle\hat{\mathcal{C}}_{n}\rangle\right)\hat{\sigma}^{z}_{n}-\hat{\mathcal{S}}_{n}\hat{\sigma}^{y}_{n}\Big{)}+\sum_{n,k}\hat{\sigma}^{x}_{n}\left(t_{n,k}-g_{n,k}\right)\left(\hat{b}^{\dagger}_{k}+\hat{b}_{k}\right)$
(S30) $\displaystyle+$
$\displaystyle\sum_{n}\Big{[}J_{y}\left(\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}-\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle\right)+J_{z}\left(\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}-\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle\right)\Big{]}\hat{\sigma}^{y}_{n}\hat{\sigma}^{y}_{n+1}$
$\displaystyle+$
$\displaystyle\sum_{j}\Big{[}J_{z}\left(\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}-\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle\right)+J_{y}\left(\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}-\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle\right)\Big{]}\hat{\sigma}^{z}_{n}\hat{\sigma}^{z}_{n+1}$
$\displaystyle+$
$\displaystyle\sum_{n}\Big{[}\left(J_{y}\hat{\mathcal{C}}_{n}\hat{\mathcal{S}}_{n+1}-J_{z}\hat{\mathcal{S}}_{n}\hat{\mathcal{C}}_{n+1}\right)\hat{\sigma}^{y}_{n}\hat{\sigma}^{z}_{n+1}+\left(J_{y}\hat{\mathcal{S}}_{n}\hat{\mathcal{C}}_{n+1}-J_{z}\hat{\mathcal{C}}_{n}\hat{\mathcal{S}}_{n+1}\right)\hat{\sigma}^{z}_{n}\hat{\sigma}^{y}_{n+1}\Big{]}.$

Here, the bath averaged quantities are identical to the expressions summarized
in Eq. (S29). If we now consider a symmetric model, with the full-polaron
mapping ($J_{y}=J_{z}$ and $g_{n,k}=t_{n,k}=t_{k}$), the interaction
Hamiltonian reduces to

$\hat{H}^{\text{pol}}_{I}=\sum_{n}\Delta_{n}\Big{(}\left(\hat{\mathcal{C}}_{n}-\langle\hat{\mathcal{C}}_{n}\rangle\right)\hat{\sigma}^{z}_{n}-\hat{\mathcal{S}}_{n}\hat{\sigma}^{y}_{n}\Big{)},$
(S31)

which will be indeed weak when $\Delta_{n}$ is small. As a result, we can
safely approximate the equilibrium steady-state of the dissipative chain to be

$\hat{\rho}_{S}^{\text{pol}}=e^{-\beta\hat{H}^{\text{pol}}_{S}}/Z^{\text{pol}}_{S},$
(S32)

from which various steady-state observable can be computed.
$Z^{\text{pol}}_{S}$ is the partition function of the model.

Next, we express the coefficients in the Hamiltonian using the continuous
spectral density function,
$K_{n}(\omega)=\sum_{k}t^{2}_{n,k}\delta(\omega-\nu_{k})$. The full-polaron
mapping is known to be most accurate for the super-Ohmic spectral function
given by

$\displaystyle
K_{n}(\omega)=\alpha_{n}\frac{\omega^{3}}{\omega^{2}_{c}}e^{-\omega/\omega_{c}},$
(S33)

where $\alpha_{n}$ is a dimensionless coupling strength for n${}^{\text{th}}$
spin and $\omega_{c}$ is the cutoff (or characteristic) frequency of the bath.
Using this model, we are able to represent all the parameters in Eq. (S29) in
terms these parameters, and the bath’s temperature. Specifically,

$\displaystyle E^{(n)}_{0}$ $\displaystyle=-2\omega_{c}\sum_{n}\alpha_{n}$
(S34) $\displaystyle E^{(n)}_{I}$
$\displaystyle=2\omega_{c}\sqrt{\alpha_{n}\alpha_{n+1}}$ $\displaystyle
E^{(nm)}_{I}$ $\displaystyle=2\omega_{c}\sqrt{\alpha_{n}\alpha_{m}}$
$\displaystyle\langle\hat{\mathcal{C}}_{n}\rangle$
$\displaystyle=\exp\left(-2\int^{\infty}_{0}d\omega\frac{K_{n}(\omega)}{\omega^{2}}\coth\left(\frac{\beta\omega}{2}\right)\right)$
$\displaystyle\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle$
$\displaystyle=\frac{1}{2}\left[\exp\left(-2\int
d\omega\frac{K_{n}(\omega)+2\sqrt{K_{n}(\omega)K_{n+1}(\omega)}+K_{n+1}(\omega)}{\omega^{2}}\coth\left(\frac{\beta\omega}{2}\right)\right)\right]$
$\displaystyle+\frac{1}{2}\left[\exp\left(-2\int
d\omega\frac{K_{n}(\omega)-2\sqrt{K_{n}(\omega)K_{n+1}(\omega)}+K_{n+1}(\omega)}{\omega^{2}}\coth\left(\frac{\beta\omega}{2}\right)\right)\right]$
$\displaystyle\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle=$
$\displaystyle-\frac{1}{2}\left[\exp\left(-2\int
d\omega\frac{K_{n}(\omega)+2\sqrt{K_{n}(\omega)K_{n+1}(\omega)}+K_{n+1}(\omega)}{\omega^{2}}\coth\left(\frac{\beta\omega}{2}\right)\right)\right]$
$\displaystyle+\frac{1}{2}\left[\exp\left(-2\int
d\omega\frac{K_{n}(\omega)-2\sqrt{K_{n}(\omega)K_{n+1}(\omega)}+K_{n+1}(\omega)}{\omega^{2}}\coth\left(\frac{\beta\omega}{2}\right)\right)\right].$

In the special case when all the spins couple to the bath in an identical
manner,

$\displaystyle E_{I}=2\omega_{c}\alpha$ (S35)

is identified as the bath-induced spin-spin interaction energy (with a factor
of 2 further appearing in Eq. (S28)). As expected, this interaction energy
depends on properties of the bath: Its characteristic frequency and its
coupling parameter to the system. Interestingly, it does not depend on the
temperature of the bath; temperature dependence of $E_{I}$ could show up in
the variational polaron treatment Anto-Sztrikacs  _et al._ (2023b) and when
considering other types of baths, e.g., spin baths.

###  S3.2 Fully local coupling model

We now study the Heisenberg model with a fully-local coupling scheme, that is,
with individual separate baths coupled to each site. The total Hamiltonian is

$\hat{H}_{\text{loca}}=\hat{H}_{S}+\sum^{N}_{n=1}\hat{S}_{\text{loca,n}}\sum_{k}t_{n,k}\left(\hat{c}^{\dagger}_{n,k}+\hat{c}_{n,k}\right)+\sum_{n,k}\nu_{n,k}\hat{c}^{\dagger}_{n,k}\hat{c}_{n,k}$
(S36)

where $\hat{S}_{\text{loca},n}=\hat{\sigma}^{x}_{n}$. After consecutive
polaron transformation via the unitary

$\hat{W}_{n}=\exp(-i\hat{\sigma}^{x}_{n}\hat{B}_{n}/2)\quad\text{where}\quad\hat{B}_{n}=2i\sum_{k}\frac{g_{n,k}}{\nu_{k}}\left(\hat{c}^{\dagger}_{n,k}-\hat{c}_{n,k}\right).$
(S37)

we obtain the polaron-transformed Hamiltonian

$\hat{H}^{\text{pol}}=\hat{H}^{\text{pol}}_{S}+\hat{H}^{\text{pol}}_{B}+\hat{H}^{\text{pol}}_{I},$
(S38)

where the system’s Hamiltonian is

$\hat{H}^{\text{pol}}_{S}=\sum^{N}_{n=1}\Delta_{n}\langle\hat{\mathcal{C}}_{n}\rangle+\sum^{N-1}_{n=1}\left(J_{x}\hat{\sigma}^{x}_{n}\hat{\sigma}^{x}_{n+1}+J_{y}\langle\hat{\mathcal{C}}_{n}\rangle\langle\hat{\mathcal{C}}_{n+1}\rangle\hat{\sigma}^{y}_{n}\hat{\sigma}^{y}_{n+1}+J_{z}\langle\hat{\mathcal{C}}_{n}\rangle\langle\hat{\mathcal{C}}_{n+1}\rangle\hat{\sigma}^{z}_{n}\hat{\sigma}^{z}_{n+1}\right).$
(S39)

Recall that averages are done with respect to the state of the baths, which
are given, as before, by

$\hat{H}^{\text{pol}}_{B}=\sum_{n,k}\nu_{n,k}\hat{c}^{\dagger}_{n,k}\hat{c}_{n,k},$
(S40)

and the system-bath interaction Hamiltonian is

$\displaystyle\hat{H}^{\text{pol}}_{I}=$
$\displaystyle\sum^{N}_{n=1}\Delta_{n}\left(\left[\hat{\mathcal{C}}_{n}-\langle\hat{\mathcal{C}}_{n}\rangle\right]\hat{\sigma}^{z}_{n}-\hat{\mathcal{S}}_{n}\hat{\sigma}^{y}_{n}\right)$
(S41)
$\displaystyle\sum^{N-1}_{n=1}\left[J_{y}\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}-J_{y}\langle\hat{\mathcal{C}}_{n}\rangle\langle\hat{\mathcal{C}}_{n+1}\rangle+J_{z}\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\right]\hat{\sigma}^{y}_{n}\hat{\sigma}^{y}_{n+1}$
$\displaystyle+$
$\displaystyle\sum^{N-1}_{n=1}\left[J_{z}\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}-J_{z}\langle\hat{\mathcal{C}}_{n}\rangle\langle\hat{\mathcal{C}}_{n+1}\rangle+J_{y}\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\right]\hat{\sigma}^{z}_{n}\hat{\sigma}^{z}_{n+1}$
$\displaystyle+$
$\displaystyle\sum^{N-1}_{n=1}\Big{(}\left[J_{y}\hat{\mathcal{C}}_{n}\hat{\mathcal{S}}_{n+1}-J_{z}\hat{\mathcal{S}}_{n}\hat{\mathcal{C}}_{n+1}\right]\hat{\sigma}^{y}_{n}\hat{\sigma}^{z}_{n+1}+\left[J_{y}\hat{\mathcal{S}}_{n}\hat{\mathcal{C}}_{n+1}-J_{z}\hat{\mathcal{S}}_{n}\hat{\mathcal{C}}_{n+1}\right]\hat{\sigma}^{z}_{n}\hat{\sigma}^{y}_{n+1}\Big{)}.$

The expressions in this Hamiltonian are given by Eq. (S29), before enacting
the thermal average. The important difference is that each site $n$ here
couples to a distinct bath, and the baths are uncorrelated, leading to
$\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle=0$ and
$\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle=\langle\hat{\mathcal{C}}_{n}\rangle\langle\hat{\mathcal{C}}_{n+1}\rangle$.

###  S3.3 Half-half coupling model

Next, we move on to the half-and-half coupling scheme. Here, an $N$-site chain
is coupled to two reservoirs, $L$ and $R$, according to the following
Hamiltonian

$\hat{H}_{\text{half}}=\hat{H}_{S}+\sum_{n}\hat{S}_{\text{half,n}}\sum_{k}t_{n,k}\left(\hat{c}^{\dagger}_{n,k}+\hat{c}_{n,k}\right)+\sum_{n,k}\nu_{n,k}\hat{c}^{\dagger}_{n,k}\hat{c}_{n,k}.$
(S42)

The system $\hat{H}_{S}$ is again given by Eq. (S9), $n=\\{L,R\\}$. As for the
coupling operators, we use
$\hat{S}_{\text{half,L}}=\sum^{N/2}_{n=1}\hat{\sigma}^{x}_{n}$, and
$\hat{S}_{\text{half},R}=\sum^{N}_{N/2+1}\hat{\sigma}^{x}_{n}$. That is, the
left half of the chain is coupled to bath $L$ while the right half is coupled
to bath $R$. We can already infer the effect of the half-and-half coupling
scheme. That is, the system develops a long-range spin-spin interaction in the
$x$ direction for all spins coupled to a common bath. On the other hand, the
$J_{y(z)}$ interaction at the boundary between the two segments is expected to
be suppressed. To observe the effect of the transformation more clearly, we
break $\hat{H}_{S}$ into three sectors,

$\displaystyle\hat{H}_{S}=\hat{H}_{S,L}+\hat{H}_{S,\text{boundary}}+\hat{H}_{S,R}.$
(S43)

Here,

$\displaystyle\hat{H}_{S,L}=$
$\displaystyle\sum^{N/2}_{n=1}\Delta_{n}\hat{\sigma}^{z}_{n}+\sum_{\alpha}\sum^{N/2-1}_{n=1}J_{\alpha}\hat{\sigma}^{\alpha}_{n}\hat{\sigma}^{\alpha}_{n+1}$
(S44) $\displaystyle\hat{H}_{S,\text{boundary}}=$
$\displaystyle\sum_{\alpha}J_{\alpha}\hat{\sigma}^{\alpha}_{N/2}\hat{\sigma}^{\alpha}_{N/2+1}$
$\displaystyle\hat{H}_{S,R}=$
$\displaystyle\sum^{N}_{n=N/2+1}\Delta_{n}\hat{\sigma}^{z}_{n}+\sum_{\alpha}\sum^{N-1}_{n=N/2+1}J_{\alpha}\hat{\sigma}^{\alpha}_{n}\hat{\sigma}^{\alpha}_{n+1}$

Note that the $\hat{H}_{S,L}$ ($\hat{H}_{S,R}$) is the left (right) half of
the chain that are coupled to a common bath $L$ ($R$). The boundary term only
involves the spin-spin interaction at the right most end of chain $L$ and the
left-most end of chain $R$. The polaron transformed Hamiltonian is

$\hat{H}^{\text{pol}}_{S}=\hat{H}^{\text{pol}}_{S,L}+\hat{H}^{\text{pol}}_{S,\text{boundary}}+\hat{H}^{\text{pol}}_{S,R}$
(S45)

where

$\displaystyle\hat{H}^{\text{pol}}_{S,L}=$
$\displaystyle\sum^{N/2}_{n=1}E^{(n)}_{0}\hat{I}+\sum^{N/2}_{n=1}\Delta_{n}\langle\hat{\mathcal{C}}_{n}\rangle\hat{\sigma}^{z}_{n}+\sum^{N/2-1}_{n=1}\left(J_{x}-2E^{(n)}_{I}\right)\hat{\sigma}^{x}_{n}\hat{\sigma}^{x}_{n+1}-\sum_{|n-m|>1}E^{(nm)}_{I}\hat{\sigma}^{x}_{n}\hat{\sigma}^{x}_{m}$
(S46) $\displaystyle+$
$\displaystyle\sum^{N/2-1}_{n=1}\Big{(}\left[J_{y}\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle+J_{z}\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle\right]\hat{\sigma}^{y}_{n}\hat{\sigma}^{y}_{n+1}+\left[J_{z}\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle+J_{y}\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle\right]\hat{\sigma}^{z}_{n}\hat{\sigma}^{z}_{n+1}\Big{)},$
$\displaystyle\hat{H}^{\text{pol}}_{S,R}=$
$\displaystyle\sum^{N}_{n=N/2+1}E^{(n)}_{0}\hat{I}+\sum^{N}_{n=N/2+1}\Delta_{n}\langle\hat{\mathcal{C}}_{n}\rangle\hat{\sigma}^{z}_{n}+\sum^{N-1}_{n=N/2+1}\left(J_{x}-2E^{(n)}_{I}\right)\hat{\sigma}^{x}_{n}\hat{\sigma}^{x}_{n+1}-\sum_{|n-m|>1}E^{(nm)}_{I}\hat{\sigma}^{x}_{n}\hat{\sigma}^{x}_{m}$
(S47) $\displaystyle+$
$\displaystyle\sum^{N-1}_{n=N/2+1}\Big{(}\left[J_{y}\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle+J_{z}\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle\right]\hat{\sigma}^{y}_{n}\hat{\sigma}^{y}_{n+1}+\left[J_{z}\langle\hat{\mathcal{C}}_{n}\hat{\mathcal{C}}_{n+1}\rangle+J_{y}\langle\hat{\mathcal{S}}_{n}\hat{\mathcal{S}}_{n+1}\rangle\right]\hat{\sigma}^{z}_{n}\hat{\sigma}^{z}_{n+1}\Big{)},$

and the boundary term is

$\hat{H}^{\text{pol}}_{S,\text{bound}}=J_{x}\hat{\sigma}^{x}_{N/2}\hat{\sigma}^{x}_{N/2+1}+J_{y}\langle\hat{\mathcal{C}}_{N/2}\rangle\langle\hat{\mathcal{C}}_{N/2+1}\rangle\hat{\sigma}^{y}_{N/2}\hat{\sigma}^{y}_{N/2+1}+J_{z}\langle\hat{\mathcal{C}}_{N/2}\rangle\langle\hat{\mathcal{C}}_{N/2+1}\rangle\hat{\sigma}^{z}_{N/2}\hat{\sigma}^{z}_{N/2+1}$
(S48)

Note that at strong coupling, the $L$ and $R$ segments of the chain interact
with each other via the $\sigma^{x}_{N/2}\sigma^{x}_{N/2+1}$ coupling at the
boundary thus preferring the anti-ferromagnetic exchange $J_{x}>0$. In
contrast, the rest of the spins will aligned themselves in their respective
baths in the $x$ direction due to the FM term over dominating at strong
coupling.

Expressions for the renormalized parameters are given by Eq. (S29). The
system-bath interaction Hamiltonian can be trivially obtained by combining Eq.
(S30) and Eq. (S41). The bath Hamiltonians remain unaffected by the polaron
transform.

###  S3.4 Pairwise coupling model

Finally, we examine the pairwise coupling scheme given by the following
Hamiltonian

$\hat{H}_{\text{pair}}=\hat{H}_{S}+\sum^{N/2}_{n=1}\hat{S}_{\text{pair,n}}\sum_{k}t_{n,k}\left(\hat{c}^{\dagger}_{n,k}+\hat{c}_{n,k}\right)+\sum_{n,k}\nu_{n,k}\hat{c}^{\dagger}_{n,k}\hat{c}_{n,k}$
(S49)

where $\hat{S}_{\text{pair},n}=\hat{\sigma}^{x}_{2n-1}+\hat{\sigma}^{x}_{2n}$;
$n\in\\{1,\dots,N/2\\}$ is the bath index. Here, every odd spin will be
coupled to a common bath together with the even numbered spin to its right.
The resulting system’s Hamiltonian after the polaron transformation is

$\displaystyle\hat{H}^{\text{pol}}_{\text{pair},S}=$
$\displaystyle\sum^{N/2}_{n=1}\left(\Delta_{2n-1}\langle\hat{\mathcal{C}}_{2n-1}\rangle\hat{\sigma}^{z}_{2n-1}+\Delta_{2n}\langle\hat{\mathcal{C}}_{2n}\rangle\hat{\sigma}^{z}_{2n}\right)$
(S50) $\displaystyle+$
$\displaystyle\sum^{N/2}_{n=1}J_{x}\hat{\sigma}^{x}_{2n}\hat{\sigma}^{x}_{2n+1}+\sum^{N/2}_{n=1}\left(J_{x}-2E^{(2n-1,2n)}_{I}\right)\hat{\sigma}^{x}_{2n-1}\hat{\sigma}^{x}_{2n}$
$\displaystyle+$
$\displaystyle\sum^{N/2}_{n=1}\left[\left[J_{y}\langle\hat{\mathcal{C}}_{2n-1}\hat{\mathcal{C}}_{2n}\rangle+J_{z}\langle\hat{\mathcal{S}}_{2n-1}\hat{\mathcal{S}}_{2n}\rangle\right]\hat{\sigma}^{y}_{2n-1}\hat{\sigma}^{y}_{2n}+\left[J_{z}\langle\hat{\mathcal{C}}_{2n-1}\hat{\mathcal{C}}_{2n}\rangle+J_{y}\langle\hat{\mathcal{S}}_{2n-1}\hat{\mathcal{S}}_{2n}\rangle\right]\hat{\sigma}^{z}_{2n-1}\hat{\sigma}^{z}_{2n}\right]$
$\displaystyle+$
$\displaystyle\sum^{N/2-1}_{n=1}\left[J_{y}\langle\hat{\mathcal{C}}_{2n}\rangle\langle\hat{\mathcal{C}}_{2n+1}\rangle\hat{\sigma}^{y}_{2n}\hat{\sigma}^{y}_{2n+1}+J_{z}\langle\hat{\mathcal{C}}_{2n}\rangle\langle\hat{\mathcal{C}}_{2n+1}\rangle\hat{\sigma}^{z}_{2n}\hat{\sigma}^{z}_{2n+1}\right].$

The new system-bath interaction Hamiltonian can be obtained by combining Eq.
(S30) and Eq. (S41).

##  S4 Discussion of the two mapping approaches

We summarize the two mapping approaches. In the EFFH method, we first extract
a collective reaction coordinate mode, polaron transform this specific mode
and truncate it. This approach was exercised on spin chains in Sec. S2. A
polaron approach on the full bath was detailed in Sec. S3. Notably, the two
mapping approaches provide the same form for the system’s Hamiltonian. This
can be clearly observed for, e.g., the global coupling scheme by comparing the
system terms in Eq. (S11) to the system’s Hamiltonian in Eq. (S28). Similarly,
in the local-bath scheme we get Eq. (S14) corresponding to Eq. (S39). In the
half-and-half coupling model, Eq. (S17) matches Eq. (S43). In the pairwise
coupling model, Eq. (S21) corresponds to Eq. (S50). The significance of this
agreement is that we predict that equilibrium expectation values of these
system’s Hamiltonian will display the same physics. Moreover, the prediction
on the expected magnetic order hold for any spectral function—so long as one
can confirm that in the EFFH model, the residual coupling is weak, and that
the polaron approach provides a good approximation to exact results. Here, we
tested the EFFH model assuming a Brownian spectral function for the bath, and
the polaron mapping assuming a super-Ohmic function.

As for the system-bath interaction Hamiltonian, the EFFH and polaron mapping
methods yield different terms. Most notably, the interaction Hamiltonian under
the polaron mapping are more complex bringing more terms. Whether relaxation
dynamics to equilibrium or steady state will be noticeably different in this
two methods remains a topic of future studies.

The advantages of the EFFH method is the simplicity of derivation compared to
the polaron method. However, the polaron approach here was general in the
sense that different spins could couple with different strength to the same
bath. In contrast, the EFFH mapping needs to be further extended to handle
such an effect. Here the difficulty would lie in building a RC that balances
different coupling strength to the different spins.

##  S5 Application: The Transverse field Ising chain

In this section, we study the magnetic order in the Ising model. First, we
consider the case of an Ising chain coupled to a global bath, then look at
other coupling schemes.

We turn the general Heisenberg model into an Ising chain by setting
$J_{y}=J_{z}=0$. This allows us to clearly analyze and visualize the preferred
magnetic order under system-bath coupling. The starting Hamiltonian is

$\hat{H}^{\text{Ising}}_{S}=\sum^{N}_{i=1}\Delta_{i}\hat{\sigma}^{z}_{i}+\sum^{N-1}_{i=1}J_{x}\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{i+1}.$
(S51)

In the global-coupling scheme, the effective system’s Hamiltonian, obtained
from either the EFFH treatment or the polaron approach, is given by

$\displaystyle\hat{H}^{\text{eff}}_{\text{glob},S}=$
$\displaystyle\sum^{N}_{i=1}\tilde{\Delta}_{i}\hat{\sigma}^{z}_{i}+\sum^{N-1}_{i=1}\tilde{J}_{x}\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{i+1}-\frac{\lambda^{2}}{\Omega}\hat{S}^{2}_{\text{glob}}$
(S52) $\displaystyle\hat{H}^{\text{pol}}_{\text{glob},S}=$
$\displaystyle\sum^{N}_{n=1}E^{(n)}_{0}\hat{I}+\sum^{N}_{n=1}\Delta_{n}\langle\hat{\mathcal{C}}_{n}\rangle\hat{\sigma}^{z}_{n}+\sum^{N-1}_{n=1}\left(J_{x}-2E^{(n)}_{I}\right)\hat{\sigma}^{x}_{n}\hat{\sigma}^{x}_{n+1}-2\sum_{n+1<m}E^{(nm)}_{I}\hat{\sigma}^{x}_{n}\hat{\sigma}^{x}_{m}.$

It is obtained by setting $J_{y}=J_{z}=0$ in Eq. (S11) and Eq. (S28). The most
important aspect of this model is the development of bath-induced long range
spin-spin coupling with magnitude $\frac{2\lambda^{2}}{\Omega}$ in the
language of the EFFH mapping and $4\omega_{c}\alpha$ in the polaron super-
Ohmic bath mapping case, see Eq. (S35).

Figure S2: Structure factor
$S_{x}=\frac{1}{N^{2}}\sum_{ij}\langle\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{j}\rangle$
for a globally-coupled transverse field Ising chain. We consider models
corresponding to (a) super-Ohmic and (b) Brownian baths. The structure factor
indicates a clear crossover from an AFM to FM order due to the coupling to the
bath with results displayed for different chain lengths, $N=4,6,8,10$ sites,
and temperatures; dark to light corresponding to $T=0.2,0.1,0.05$. Insets
$(a_{1})$ and $(a_{2})$ correspond to spin correlations,
$\langle\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{j}\rangle$, at $\alpha=0$ and
$\alpha=0.3$, respectively for $N=10$ with $T=\Delta=0.1$ and $J_{x}=1$.
Insets $(b_{1})$ and $(b_{2})$ correspond to
$\langle\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{j}\rangle$ at $\lambda=0$ and
$\lambda=3$ respectively for $N=10$ with $T=\Delta=0.1$ and $J_{x}=1$.

In Figure S2, we plot the structure factor
$S_{x}=\frac{1}{N^{2}}\sum_{ij}\langle\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{j}\rangle$
as a function of coupling strength to the bath. The nature of the bath affects
the renormalization of parameters and the bath-generated spin coupling
$E_{I}$. As a demonstration, we assume a super-Ohmic bath spectral function in
panel (a) of Figure S2. Panel (b) of that figure could correspond to any
spectral function, as long as the residual coupling is weak. Here we associate
it as an example with a Brownian spectral function expressed in Eq. (S7).
Regardless of the bath spectral density, the structure function evinces that
the preferred magnetic order in the Ising chain with the original AFM
interaction goes into an FM order as the coupling to the bath increases.
Results are displayed for different chain sizes and at different temperatures.
Furthermore, we display as a heat map spin correlations in the $x$ direction
at both weak and strong coupling, exemplifying the AFM and the FM order,
respectively.

Figure S3: Spin-spin correlations in the $x$ direction for the transverse
field Ising chain coupled to super-Ohmic bath(s). Rows represents fully-
global, fully-local, half-and-half, and pairwise coupling schemes from top to
bottom corresponding to Fig. 1 in the Main text. Columns represent five
different coupled strengths to clearly demonstrate the crossover of the
preferred magnetic alignment for each coupling schemes (with the specific
coupling strengths indicated in the title). Parameters are chosen to be
$N=10$, $J_{x}=1$, $T=\Delta=0.1$.

Finally, in Figure S3, we present a comprehensive picture of bath-induced
magnetic order in the transverse field Ising chain with different schemes for
the bath coupling. We present the spin-spin correlations,
$\langle\hat{\sigma}_{i}\hat{\sigma}_{j}\rangle$, and use parameters
corresponding to a super-Ohmic bath spectral density function.

We study the four models of Fig. 1 (Main) and present the spin-spin
correlations in the $x$ direction in Fig. S3. The first row corresponds to the
global-bath model. Here, the system shows the transition from an AFM to a FM
order with increasing coupling $\alpha$. The second row corresponds to the
local bath case. Here, the interaction with the bath suppresses the energy
splitting thus turning the model into the zero-field Ising model. The third
row shows results for the half-and-half coupling model. Here, we observe the
development of a domain wall between the two segments, each showing internally
an FM order. The last row displays the development of an extended Neel phase
in the pairwise coupling model. Here each pair of spins align in the same
direction, but the interaction between every pair is still anti-ferromagnetic.
This case is described in details in Fig. 3 of the Main text.

## References

  * Anto-Sztrikacs  _et al._ (2023a) N. Anto-Sztrikacs, A. Nazir, and D. Segal, PRX Quantum 4, 020307 (2023a).
  * Anto-Sztrikacs  _et al._ (2023b) N. Anto-Sztrikacs, B. Min, M. Brenes, and D. Segal, Phys. Rev. B 108, 115437 (2023b).
  * Iles-Smith  _et al._ (2014) J. Iles-Smith, N. Lambert, and A. Nazir, Phys. Rev. A 90, 032114 (2014).
  * Anto-Sztrikacs and Segal (2021) N. Anto-Sztrikacs and D. Segal, New Journal of Physics 23, 063036 (2021).
  * Nazir and Schaller (2018) A. Nazir and G. Schaller, “The reaction coordinate mapping in quantum thermodynamics,” in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_ (Springer International Publishing, 2018) pp. 551–577.
  * Cresser and Anders (2021) J. D. Cresser and J. Anders, Phys. Rev. Lett. 127, 250601 (2021).
  * Chiu  _et al._ (2022) Y.-F. Chiu, A. Strathearn, and J. Keeling, Phys. Rev. A 106, 012204 (2022).
  * Iles-Smith  _et al._ (2016) J. Iles-Smith, A. G. Dijkstra, N. Lambert, and A. Nazir, The Journal of Chemical Physics 144, 044110 (2016).
  * Xu and Cao (2016) D. Xu and J. Cao, Front. Phys. 11, 110308 (2016).