#  Photoinduced non-reciprocal magnetism

Ryo Hanai<EMAIL_ADDRESS>Center for Gravitational Physics and
Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto 606-8502, Japan  Daiki Ootsuki
<EMAIL_ADDRESS>Graduate School of Human and Environmental
Studies, Kyoto University, Sakyo-ku, Kyoto 606-8501, Japan  Rina Tazai
<EMAIL_ADDRESS>Center for Gravitational Physics and Quantum
Information, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto
606-8502, Japan

###### Abstract

Out of equilibrium, the action-reaction symmetry of the interactions is often
broken, leading to the emergence of various collective phenomena with no
equilibrium counterparts. Although ubiquitous in classical active systems,
implementing such non-reciprocal interactions in solid-state systems has
remained challenging, as the known quantum schemes require precise control
over the system on a single-site level. Here, we propose a novel dissipation-
engineering protocol to induce non-reciprocal interactions in solid-state
platforms with light, which we expect to be achievable with state-of-the-art
experimental techniques. Focusing on magnetic metals for concreteness, we show
microscopically that a light injection that introduces the decay channel to a
virtually excited state gives rise to non-reciprocal interactions between
localized spins. One can even realize a situation where spin A tries to align
with spin B but the B tries the opposite, resulting in a chase-and-runaway
dynamics. Applying our scheme to layered ferromagnets, we show that a non-
reciprocal phase transition from a static to a many-body time-dependent chiral
phase emerges. Our work paves the way to bring solid-state systems to the
realm of non-reciprocal science, providing yet another possibility to control
quantum matter with light.

Figure 1:  Light control of magnetic interactions and magnetism via
dissipation. (a) In the absence of light, the interaction between the local
spins (thick blue arrow) is reciprocal. The spins are driven towards the
equilibrium configuration [alignment in the ferromagnetic case illustrated
here] through a magnetic friction called the Gilbert damping (green arrows).
(b) When the light is tuned to a frequency $h\nu$ that selectively activates
the red spin [in the way illustrated in Fig. 2], the light-induced torque
(pink arrows) acts on the activated spin. As a result, the two spins
effectively interact non-reciprocally, where the active (quiescent) spin tries
to anti-align (align) with the opponent’s spin. (c) Phase diagram of a layered
ferromagnet under light injection, determined by our meanfield description
(Eq. (8)). Here, the two ferromagnetic layers (A and B) are separated by a
non-magnetic metal and the laser is injected to introduce active dissipation
to the B layer at rate $\gamma_{\rm B}$, making the interlayer magnetic
interactions non-reciprocal. When the light is off ($\gamma_{\rm B}=0$), the
magnetization of the two layers aligns for ferromagnetic interlayer
interaction $j_{\rm AB}(=j_{\rm BA})>0$ (blue region). As the light-induced
dissipation is turned on $\gamma_{\rm B}>0$, the system exhibits a phase
transition to an antialigned configuration (red) at $\gamma_{\rm
B}\simeq\alpha_{\rm B}|g_{\rm B}|$ and a non-reciprocal phase transition to a
time-dependent chiral phase [where the two magnetizations exhibit many-body
chase-and-runaway dynamics](cyan). [See text and Fig. 4 for further details.]
We set the intralayer interaction and Gilbert damping of the layer A (B) as
$j_{\rm AA}=10{\rm meV}(j_{\rm BB}=5{\rm meV}$) and $\alpha_{\rm
A}=0.1(\alpha_{\rm B}=2\times 10^{-3})$, respectively. The sd coupling
strength for B is $g_{\rm B}=-10{\rm meV}$, the filling $n=1$, and the
temperature is $k_{\rm B}T=9{\rm meV}$. The dashed lines are the phase
boundary at a lower temperature $k_{\rm B}T=5{\rm meV}$.

In equilibrium, the free energy minimization principle states that all
interactions between constituents must have action-reaction symmetry. However,
this constraint is no longer present once the system is out of equilibrium [1,
2]. In fact, non-reciprocal interactions are ubiquitous in Nature: the brain
is composed of inhibitory and exhibitory neurons that non-reciprocally
interact [3, 4, 5, 6]; the predator chases the prey and the prey runs away
because their interaction is asymmetric; colloids immersed in a
chemically/optically active media exhibit non-reciprocal interactions [7, 8,
9, 10]. Recent studies revealed that such non-reciprocal interactions
fundamentally affect the collective properties of many-body systems [11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. A prominent example
is the emergence of non-reciprocal phase transitions [15, 16, 17, 18, 19, 20,
21, 22, 23, 25, 24, 27, 28], where a time-dependent phase that exhibits a
collective and persistent chase-and-runaway motion between macroscopic
quantities arises. Its critical point is characterized by the coalescence of a
collective mode to the Nambu-Goldstone mode [17, 15, 20, 21, 22, 23, 25]
(instead of merely being degenerate as in the conventional cases), which is
strictly forbidden in equilibrium.

Non-reciprocal interactions are not restricted to the classical systems
mentioned above. Quantum systems may also exhibit non-reciprocity [29, 30, 28]
by carefully tailoring dissipation microscopically [31, 32]. Indeed, such non-
reciprocity has been implemented in synthetic quantum systems such as cold
atoms [33], optomechanics [34], and circuit QED [35]. However, these schemes
require fine control of dissipation and gauge flux at a single-site/plaquette
level, imposing a challenge to realize non-reciprocal interactions in solid-
state systems.

Figure 2:  Dissipation engineering RKKY interactions in magnetic metals. (a)
Schematic description of magnetic metals. It is composed of localized spins
(thick blue arrows) and conduction electrons (yellow spheres with arrows
attached). The conduction electrons form a Fermi sea up to the Fermi energy
$\varepsilon_{\rm F}$. The localized electrons are modeled as a two-level
system with energy $\varepsilon_{d}$ with an on-site Coulomb interaction $U$.
The two types of electrons couple through a spin exchange coupling called the
sd coupling. Here, the sd spin exchange coupling arises through the second-
order process where the conduction electron tunnels into the localized orbital
to virtually excite to a double-occupied state and back, which involves spin-
flip (the bottom panel). Using the conduction electrons as a medium, a
magnetic interaction between the localized spins arises (the so-called RKKY
interaction). As these perturbative processes lower the energy, the ground
state configuration is the state that activates these processes the most. (b)
Our dissipation-engineering scheme is to inject light at a frequency $h\nu$
resonant with the double-occupied state and a higher-level state (say, an
unoccupied higher-level f-orbital state at the energy $\varepsilon_{f,\rm B}$
that quickly dissipates with rate $\Gamma_{f,\rm B}$ in the figure). Since the
decay occurs only when the localized-conduction electron exchange process
activates (which lowers the energy), this process selectively destroys the
energetically favored states, giving rise to a light-induced torque (pink
arrow in Fig. 1(b)) that applies opposite to Gilbert damping (green arrows in
Fig. 1(a),(b)). When this decay is turned on only to spin B but not A, the
light-induced torque is applied only to B spins and hence a non-reciprocal
interaction emerges. The lost electron is immediately resupplied from the
surrounding environment.

In this work, we propose a novel dissipation-engineering scheme to realize
non-reciprocal interactions in solid-state systems with light (Fig. 1(a),(b)),
which we expect to be achievable with state-of-the-art experimental
techniques. We show that interactions between localized spins in magnetic
metals mediated by conduction electrons, known as the Ruderman-Kittel-Kasuya-
Yosida (RKKY) interaction [36, 37, 38], can be engineered to switch from
ferromagnetic to antiferromagnetic (or vice versa) with light. This is
achieved by injecting light at an appropriately tuned frequency that
introduces the decay to a virtually excited state, which selectively destroys
the energetically favored state (Fig. 2). By tuning the light frequency to
activate only some fraction of the spins [where the other spins are off-
resonant from the radiated light], one can engineer non-reciprocal
interactions where spin A tries to align with spin B but spin B tries the
opposite, resulting in chase and runaway dynamics, see Fig. 1(b). We estimate
the injection power needed for their emergence is within reach of the current
techniques and the heating effect should be minimal. These results are derived
microscopically by combining the projection method generalized to open quantum
systems [39, 40, 41, 42, 43, 44] and Keldysh theory [45, 46, 47]. By applying
this scheme to a layered ferromagnet, we show that the light-induced
interlayer non-reciprocal interaction triggers a non-reciprocal phase
transition [15] to a time-dependent chiral phase where the magnetization of
the two layers exhibit a collective chase and runaway motion (Fig. 1(c)).

We expect our scheme to be applicable to a wide class of quantum materials
such as Mott insulating phases in strongly correlated electrons [48, 49],
multi-band superconductivity [50, 51], and optical phonon-mediated
superconductivity that arises e.g. in ${\rm SrTiO_{3}}$ [52, 53]. This broad
applicability is anticipated because our proposal does not rely on any
properties specific to magnetic metals as long as the interaction is mediated
via a virtual high-energy state.

The effect of non-reciprocal interactions on the collective properties of
many-body systems is currently under heavy investigation in many different
disciplines of science, ranging from active soft matter [54, 55, 8, 56, 57, 9,
10, 58, 15, 59], photonics [7], robotics [60], living matter [61, 62], open
quantum systems [29, 30, 31, 63, 27, 28], ecology [64, 65, 66, 67, 68, 69],
and neuroscience [3, 4, 5, 6], to sociology [70]. Our work serves as a
starting point to study quantum non-reciprocal science in solid-state
platforms, offering an alternative route to control quantum matter with light
[71, 72, 73, 20, 74, 75, 76, 77, 78, 79, 80].

##  I The dissipation engineering scheme

To illustrate our idea to create non-reciprocal interactions with light, let
us briefly review the spin-spin interactions in magnetic metals [81], which is
the focus of this work. [See Fig. 2(a).] These materials comprise localized
spins (responsible for magnetic properties) and conduction electrons that are
free to move (responsible for metallic properties), which couple through the
spin-exchange coupling called the sd coupling [82, 41]. The sd coupling arises
through the exchange of electrons, where the conduction electrons tunnel to
the localized orbital to virtually excite the localized electron to a double-
occupied state (see the “speech bubble” in Fig. 2(a)). The localized spins
interact indirectly through this sd coupling using the conduction electrons as
the mediating field, giving rise to magnetic interactions known as the RKKY
interactions [36, 37, 38].

We propose to dissipation engineer the RKKY interaction with light. Our idea
is to inject light that has its energy $h\nu$ tuned to selectively turn on the
tunneling from the double-occupied state to a higher-energy state that quickly
relaxes (Fig. 2(b)). This introduces a finite lifetime to the virtual state,
directly affecting the properties of the sd coupling and hence the RKKY
interaction. The lost electron is immediately resupplied such that localized
electrons are always singly occupied.

We show below from explicit microscopic calculations that this light-induced
dissipation gives rise to a torque (pink arrows in Fig. 1(b)), which
interestingly applies in the direction opposite to the energetically favored
configuration. This light-induced torque competes with the magnetic friction
called the Gilbert damping [83, 84] (green arrows) that relaxes the system to
the ground state. When we choose the frequency of the light $h\nu$ in such a
way that it only resonates with a portion of the spins (red spins in Fig.
1(b)), this light-induced torque only applies to those activated spins, giving
rise to non-reciprocal interactions.

Why does our dissipation scheme give rise to torque that applies in the
opposite direction from the Gilbert damping? This can be understood from the
general picture we provide below. Our light-injection scheme only dissipates
the virtual state of the second-order process illustrated in Fig. 2(b). Hence,
the decay channel is turned on only when this process is activated; the
configuration that activates the process more is the one that is likely to
decay faster. Note crucially that these second-order processes always lower
the energy in equilibrium because, according to the second-order perturbation
theory, the energy change due to this process is given by $\Delta
E=\sum_{m}|v_{i,m}|^{2}/(E_{i}-E_{m})<0$, where $E_{i(m)}$ is the energy of
the initial (intermediate) state of the unperturbed system and $v_{i,m}$ is
the matrix element between these states. (Note that $E_{i}<E_{m}$.) For
example, the sd exchange coupling is antiferromagnetic because the above
process can only be activated when conduction and localized spins are
antialigned due to the Pauli-blocking effect. Similarly, the sign of RKKY
interactions is determined by which configuration the perturbative processes
activate the most frequently. Therefore, the energetically favorable state is
the state that experiences the strongest decay. This results in dissipative
interaction in the direction opposite from the friction towards the ground
state arising due to the coupling with the surrounding environment such as the
substrate, bulk phonons, etc. We remark that similar physics was discussed in
Refs. [85, 86] in the context of cold atomic systems, where they also
dissipation-engineered sign-reversal of (reciprocal) interactions; but see
also Methods for their crucial differences arising from the absence of
surrounding environment.

As the above scenario does not rely on features specific to magnetic systems,
we expect our scheme to be equally relevant to a wide class of quantum
materials [48, 49, 50, 51, 52, 53].

##  II Quantum master equation

Our goal from here on is to develop a formalism that allows us to describe the
effective spin dynamics of the open quantum system with dissipation induced by
continuous light injection. In this work, we consider magnetic metals modeled
by the Anderson impurity model [87, 81] (Fig 2(a)), described by the
Hamiltonian $\hat{H}=\hat{H}_{0}+\hat{H}_{cd}$,

$\displaystyle\hat{H}_{0}$ $\displaystyle=$
$\displaystyle\sum_{\bm{k},\sigma=\uparrow,\downarrow}\varepsilon_{\bm{k}}\hat{c}^{\dagger}_{\bm{k}\sigma}\hat{c}_{\bm{k}\sigma}$
$\displaystyle+$
$\displaystyle\sum_{a}\Big{[}\sum_{\sigma=\uparrow,\downarrow}\varepsilon_{d,a}\hat{d}^{\dagger}_{\sigma,a}\hat{d}_{\sigma,a}+U_{a}\hat{d}^{\dagger}_{\uparrow,a}\hat{d}_{\uparrow,a}\hat{d}^{\dagger}_{\downarrow,a}\hat{d}_{\downarrow,a}\Big{]},$
$\displaystyle\hat{H}_{cd}$ $\displaystyle=$
$\displaystyle\sum_{a,\sigma}\big{[}v_{a}\hat{d}_{\sigma,a}^{\dagger}\hat{c}_{\bm{R}_{a}\sigma}+{\rm
h.c.}\big{]}.$ (2)

Here, conduction electrons are modeled as free electrons, where
$\hat{c}_{\bm{k},\sigma}$ is the fermionic annihilation operator of the
conduction electrons with spin $\sigma=\uparrow,\downarrow$ and momentum
$\bm{k}$ and
$\hat{c}_{\bm{r},\sigma}=\sum_{\bm{k}}e^{i\bm{k}\cdot\bm{r}}\hat{c}_{\bm{k},\sigma}$
is its Fourier transform. $\varepsilon_{\bm{k}}$ is the kinetic energy. The
conduction electrons are assumed to be large enough to be always in thermal
equilibrium, where its distribution is given by the Fermi distribution
function $f(\varepsilon_{\bm{k}})=[e^{(\varepsilon_{\bm{k}}-\varepsilon_{\rm
F})/(k_{\rm B}T)}+1]^{-1}$ ($\varepsilon_{\rm F}$ is the Fermi energy and
$k_{\rm B}$ is the Boltzmann constant) at low temperature $k_{\rm
B}T\ll\varepsilon_{\rm F}$. The localized electrons at site $a$ are modeled as
a two-level system that has energy $\varepsilon_{d,a}$ and an on-site Coulomb
repulsion $U_{a}(>0)$, where $\hat{d}_{\sigma,a}$ is a fermionic annihilation
operator at site $a$. $v_{a}$ is the conduction-localized electron mixing (c-d
mixing), and $\bm{R}_{a}$ is the position of the localized electron. The
energy level of localized electrons $\varepsilon_{d,a}$ sit below the Fermi
energy $\varepsilon_{\rm F}$ but the double-occupied state
$\varepsilon_{d,a}+U_{a}$ is above it. The resultant localized electrons are
always singly occupied. Below, we consider systems with strong enough Coulomb
repulsion $U_{a}\gg v_{a}$ that justifies the perturbative treatment of the
c-d mixing $v_{a}$ (See SI Sec. III.C for a more precise condition).

As mentioned previously, we introduce a decay channel to localized electrons
in the double-occupied state by injecting a laser that couples this state to a
higher-energy state that quickly relaxes with the rate $\Gamma_{f,a}$ (see
Fig. 2(b)). This dissipative process can be safely regarded as a Markov
process as long as the higher-energy state decays fast enough compared to its
re-population rate, which is true in the range of interest (see Methods for
details). Such Markovian open quantum systems are generally described by the
Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master equation [39, 40] (also
called the Lindblad master equation; see e.g. Ref. [88] and SI Sec. I for a
brief review), given by (where $\hat{\rho}$ is the reduced system density
operator),

$\displaystyle\partial_{t}\hat{\rho}$ $\displaystyle=$
$\displaystyle-i[\hat{H},\hat{\rho}]+\sum_{a,\sigma}\kappa_{a}{\mathcal{D}}[\hat{d}_{\sigma,a}\hat{P}_{\uparrow\downarrow}^{a}]\hat{\rho},$
(3)

for our system. The dissipator
${\mathcal{D}}[\hat{L}]\hat{\rho}=\hat{L}\hat{\rho}\hat{L}^{\dagger}-\frac{1}{2}\\{\hat{L}^{\dagger}\hat{L},\hat{\rho}\\}$
makes the time evolution non-unitary. Here, $\hat{P}_{\uparrow\downarrow}^{a}$
is a projection operator onto the double-occupied state at site $a$, turning
on a decay channel at the rate $\kappa_{a}$ only when site $a$ is double-
occupied.

We wish to derive the localized spin dynamics in the presence of light-induced
dissipation $\kappa_{a}>0$. In the equilibrium limit $\kappa_{a}\rightarrow
0$, a standard procedure to analyze the Anderson impurity model (Eq. (II) and
(2)) is to map the localized electrons in the fermionic picture to localized
spins, which is performed by projecting out the virtual excited states that
have fast oscillations [41, 81]. This incorporates the second-order process in
terms of $v_{a}$ illustrated in Fig. 2(a). Here, we perform the same procedure
in spirit but employ it to the GKSL master equation (3) [89, 90, 42, 43, 44],
see Supplemental Information (SI) Sec. I-III for details. It yields,

$\displaystyle\partial_{t}\hat{\rho}=-i[\hat{H}_{\rm sd},\hat{\rho}]$ (4)
$\displaystyle+\sum_{a}\Big{[}\gamma_{a}{\mathcal{D}}\big{[}\sum_{\sigma}\hat{d}_{\sigma,a}^{\dagger}\hat{c}_{\bm{R}_{a},\sigma}\hat{P}_{\rm
s}^{a}\big{]}+\sum_{\sigma}\kappa_{a}{\mathcal{D}}[\hat{d}_{\sigma,a}\hat{P}_{\uparrow\downarrow}^{a}]\Big{]}\hat{\rho}.$

The first term on the right-hand side describes the coherent dynamics that
have an identical form to those found in equilibrium and the second and the
third are the light-induced dissipative terms. The sd Hamiltonian [82]
$\hat{H}_{\rm sd}=-(1/2)\sum_{a}g_{a}\hat{P}_{\rm
s}^{a}[\hat{\bm{\tau}}(\bm{R}_{a})\cdot\hat{\bm{S}}_{a}]\hat{P}_{\rm s}^{a}$
describes the spin exchange coupling between the conduction and localized
spins, where $\hat{P}_{\rm s}^{a}$ is the projection operator to singly-
occupied localized electron states at site $a$. Here,
$(\hat{\bm{S}}_{a})_{i}=\sum_{\sigma,\sigma^{\prime}}\hat{d}_{\sigma,a}^{\dagger}(\sigma_{i})_{\sigma\sigma^{\prime}}\hat{d}_{\sigma^{\prime},a}$
is the localized spins, $\bm{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})$ are
the Pauli matrices,
$\hat{\bm{\tau}}(\bm{R}_{a})=\sum_{\sigma,\sigma^{\prime}}(\hat{c}_{\bm{R}_{a},\sigma^{\prime}}^{\dagger}\bm{\sigma}_{\sigma^{\prime},\sigma}\hat{c}_{\bm{R}_{a},\sigma})$
is the conduction spin at position $\bm{R}_{a}$, and
$g_{a}\simeq-|v_{a}|^{2}[(\varepsilon_{d,a}+U_{a}-\varepsilon_{\rm
F})^{-1}+(\varepsilon_{\rm F}-\varepsilon_{d,a})^{-1}]<0$ (where
$\kappa_{a}\ll U_{a}$ is assumed) is the sd coupling strength that is
antiferromagnetic. (Note that $\varepsilon_{d,a}<\varepsilon_{\rm
F}<\varepsilon_{d,a}+U_{a}$.) As usual, we have assumed that only the
excitations near the Fermi surface are responsible. We have also ignored the
impurity potential of conduction electrons, as they play a minor role.

The second term is the emergent correlated dissipation with the rate
$\gamma_{a}\simeq\kappa_{a}|v_{a}|^{2}/(\varepsilon_{d,a}+U_{a}-\varepsilon_{\rm
F})^{2}$ that arises from the interplay between the strong correlation effect
and the light-induced decay (see Fig. 5 in Methods). This term induces
dissipative tunneling of electrons from the conduction band to the localized
orbital when the electron $a$ is singly occupied. The third term describes the
decay of electrons from the double-occupied state at a much faster rate than
the correlated tunneling ($\kappa_{a}\gg\gamma_{a}$), driving the system
immediately back to the singly occupied state.

The emergent correlated dissipation (the second term in Eq. (4)) already
captures the underlying mechanism of sign inversion of interactions described
in the previous section. This can be seen from the localized spin dynamics
that are derived from the master equation (4),

$\displaystyle\dot{{\bm{S}}}_{a}=g_{a}\left<\hat{\bm{S}}_{a}\times\hat{\bm{\tau}}(\bm{R}_{a})\right>-\gamma_{a}n{\bm{S}}_{a}+\gamma_{a}{\bm{\tau}}(\bm{R}_{a}),$
(5)

where $\langle\hat{O}\rangle={\rm tr}[\hat{\rho}\hat{O}]$ is the expectation
value,
$\bm{S}_{a}=\langle\hat{\bm{S}}_{a}\rangle,\bm{\tau}(\bm{R}_{a})=\left<\hat{\bm{\tau}}(\bm{R}_{a})\right>$,
and $n$ is the filling of the conduction electron. The first and second terms
describe the coherent dynamics arising from sd interaction and the light-
induced decay of the dipole moment of the localized spin, respectively. The
third term is the emergent dissipative torque, which drives the localized spin
toward alignment with the conduction spins. This is the opposite of what is
expected from energetics, where the sd coupling is antiferromagnetic
$g_{a}<0$. We note that it was crucial for this light-induced dissipator in
Eq. (4) to have the form ${\mathcal{D}}[\sum_{\sigma}\hat{o}_{\sigma}]$
instead of $\sum_{\sigma}{\mathcal{D}}[\hat{o}_{\sigma}]$ for the dissipative
torque to appear, see SI Sec. III.C for details.

##  III Landau-Lifshitz-Gilbert equation with light-induced interactions

The remaining task to obtain the effective interaction between the localized
spins [36, 37, 38] (i.e. the dissipation-engineered RKKY interaction) is to
integrate out the conduction electron degrees of freedom that we regard as a
(non-Markovian) bath. We perform this by mapping the master equation (4) to a
Keldysh action [46] that allows us to utilize field-theoretic approaches,
taking into account non-adiabatic effects from the Fermi statistics of the
conduction electrons. Once the Keldysh action is obtained, we integrate out
the conduction electrons’ degrees of freedom within the second-order
perturbation in terms of sd coupling $g_{a}$ and light-induced correlated
decay rate $\gamma_{a}$ under the gradient approximation (i.e., Markov
approximation plus a first-order non-Markovian correction). We extract the
localized spin dynamics in the saddle-point approximation from the obtained
reduced Keldysh action consisting only of localized electrons’ degrees of
freedom. We detail the procedure in SI Sec. IV.

The obtained semiclassical localized spin dynamics are

$\displaystyle\dot{\bm{S}}_{a}=\sum_{b(\neq
a)}J_{ab}\bm{S}_{a}\times\bm{S}_{b}-\alpha_{a}\bm{S}_{a}\times\dot{\bm{S}}_{a}$
$\displaystyle\ \ \ \ \ -\gamma_{a}n\bm{S}_{a}-\sum_{b(\neq
a)}\Omega_{ab}\bm{S}_{b},$ (6)

which is one of the main results of this work. We emphasize that all terms,
including the Gilbert damping term, are obtained microscopically.

The first two terms on the right-hand side are the Landau-Lifshitz-Gilbert
equation [83, 84], while the last two terms are terms generated through our
controlled dissipation. The first term gives rise to the coherent precession
motion around the effective magnetic field $\bm{B}_{\rm eff}=\sum_{b(\neq
a)}J_{ab}(R_{ab})\bm{S}_{b}$, where $J_{ab}(=J_{ba})$ is the RKKY interaction
strength given by $J_{ab}(R_{ab})=-9\pi[(g_{a}g_{b})/\varepsilon_{\rm
F}]n^{2}F(2k_{\rm F}R_{ab})$ in the case of parabolic dispersion
$\varepsilon_{\bm{k}}=\hbar^{2}\bm{k}^{2}/(2m)$ in three spatial dimensions,
where $m$ is the conduction electron mass, $k_{\rm F}=\sqrt{2m\varepsilon_{\rm
F}}/\hbar$ is the Fermi momentum, $R_{ab}=|\bm{R}_{a}-\bm{R}_{b}|$ is the
inter-spin distance, and $F(x)=[-x\cos x+\sin x]/x^{4}$, and $\hbar$ is the
Dirac constant. The second term, obtained as the first-order correction to the
Markov approximation [91], is the Gilbert damping term describing the spins’
magnetic friction [83, 84]. This drives the system toward the ground state
configuration when combined with the first term (green arrows in Fig. 1).
Here, $\alpha_{a}$ is the Gilbert damping rate, which in the parabolic
dispersion case reads $\alpha_{a}=(9\pi^{2}/4)n^{2}(g_{a}/\varepsilon_{\rm
F})^{2}$. The relaxation rate for such a process is $\gamma_{\rm
Gil}^{ab}\sim\alpha_{a}|J_{ab}|$.

Figure 3:  Typical energy scales in different physical pictures. The energy
scales in (left panel) localized electrons immersed in conduction electrons
picture described by Eq. (3), (middle panel) localized spins immersed in
conduction electrons picture described by Eq. (4), and (right panel) in the
interacting localized spin picture (where the conduction electrons are
integrated out) described by Eq. (III). These pictures map from one to the
other via the projection method (SI Sec. III) and Keldysh theory (SI Sec. IV).
Each panel lists the energy scales in the coherent (Hamiltonian) and
dissipative dynamics. The equilibration occurs at the timescale set by the
energy scale in the Hamiltonian multiplied by the Gilbert damping rate
$\alpha_{a}$. This competes with the light-induced dissipation and the sign
reversal of interaction occurs when the latter exceeds the former.

The light-induced dissipative interactions
$\Omega_{ab}(R_{ab})=(\gamma_{a}/|g_{a}|)J_{ab}(R_{ab})(\simeq(\kappa_{a}/U_{a})J_{ab}(R_{ab}))$
compete with this equilibration dynamics (pink arrows in Fig. 1(b)). In
addition to the unavoidable decay with rate $\gamma_{a}$ of the dipole moment
described by the third term, our light induces effective interactions that
cannot be written as the derivative of the energy function. (We briefly note
that we have ignored the contribution from self-dissipative interaction
$\Omega_{aa}$, which merely renormalizes the decay rate $\gamma_{a}$.) This
dissipative interaction drives the system towards a configuration opposite
from the ground state configuration. When $\kappa_{a}\ \raise
1.29167pt\hbox{$>$}\kern-8.00003pt\lower 3.01385pt\hbox{$\sim$}\
\alpha_{a}U_{a}$, this light-induced contribution
($\Omega_{ab}=\kappa_{a}/U_{a}\cdot J_{ab}$) exceeds the Gilbert damping
($\gamma_{\rm Gil}^{ab}\sim\alpha_{a}J_{ab}$), causing the effective
interaction to change its sign. When this sign flip occurs to one of the spins
but not the other, non-reciprocal interactions with effective opposite signs
emerge, resulting in chase-and-runaway dynamics, the situation illustrated in
Fig. 1(b).

In Methods, we estimate the pumping power $P$ required to achieve this regime
as

$\displaystyle P\ \raise 1.29167pt\hbox{$>$}\kern-8.00003pt\lower
3.01385pt\hbox{$\sim$}\ \alpha_{a}\frac{2\pi U_{a}\nu
m_{0}\epsilon_{0}c}{e^{2}}\Gamma_{f,a}$ (7)

using a Lorentz oscillator model, where $e$ is the electron charge, $m_{0}$ is
the electron mass, $\epsilon_{0}$ is the vaccuum dielectric constant, and $c$
is the speed of light. Setting the typical values $U_{a}\sim 1{\rm
eV},\alpha_{a}\sim 10^{-2},\Gamma_{f,a}\sim 10{\rm meV},h\nu\sim 1{\rm eV}$,
the required pump power is $P\ \raise 1.29167pt\hbox{$>$}\kern-8.00003pt\lower
3.01385pt\hbox{$\sim$}\ 10^{8}{\rm W/cm^{2}}$. Not only is this pump power
achievable, e.g. with Raman lasers with pulse duration of $O(10{\rm ns})$ [92]
or even with a steady-state resource [93], but the heating effect should be
minimal for magnetic metals. For instance, in Ref. [78], it was reported that
the sample did not demagnetize until the pump power exceeded $P\sim
10^{10}{\rm W/cm^{2}}$. For convenience for the readers, we have summarized
the typical energy scales in Fig. 3.

Figure 4:  Non-reciprocal phase transitions in light-activated ferromagnetic
layers. (a) Layered ferromagnet composed of A and B layers separated by a non-
magnetic metal exposed of light injection that activates the B layer. The
interlayer RKKY interaction is mediated by the itinerant electrons in the non-
magnetic metal layer, which is modified by light. In particular, the
magnetization in the A layer aligns with the rate $\alpha_{\rm A}j_{\rm AB}$
for ferromagnetic interlayer interactions $j_{\rm AB}>0$, while the B layer
may align or antialign with layer A depending on the sign of the modified
effective interaction $\alpha_{\rm B}j_{\rm AB}-\Omega_{\rm AB}$.
$\varphi_{\rm A(B)}$ is the orientation direction of the magnetization
$\bm{m}_{\rm A(B)}$ of A(B)-layer ferromagnet. (b)-(d) Different phases
arising in this system. (b) Aligned phase ($\Delta\varphi=\varphi_{\rm
A}-\varphi_{\rm B}=0$) realized in the equilibrium limit $\gamma_{\rm B}=0$.
(c) Antialigned phase ($\Delta\varphi=\pi$) at $\gamma_{\rm B}=0.022{\rm
meV}$. (d) Chiral phase ($\Delta\varphi\neq 0,\pi$) at $\gamma_{\rm
B}=0.03{\rm meV}$. (e) The orientation angle difference $\Delta\varphi$ and
(f) the magnitude of the magnetization $|\bm{m}_{\rm A,B}|$ as a function of
$\gamma_{\rm B}$. We set $j_{\rm AB}=5{\rm meV}$ and the other parameters are
the same as those used in Fig. 1(c). The simulations were run from random
initial conditions.

##  IV Non-reciprocal phase transitions

So far, we have shown from microscopic calculations that one can have non-
reciprocal interactions with light. An intriguing question is how such non-
reciprocal interaction affects the collective properties of many-body systems.
In this work, at the outset, we consider a simple setup illustrated in Fig.
4(a), where two layers of ferromagnets (that we label A and B) sandwich a non-
magnetic metal, and light that activates only layer B spins is exposed to the
system. In the absence of light, this is a typical MRAM device. The itinerant
electrons in the non-magnetic metal layer mediate the interlayer RKKY
interaction [94, 95, 96, 79]. The light injection induces additional torque to
the B layer applying oppositely from conventional RKKY interaction, giving
rise to non-reciprocal interlayer interaction (Fig. 4(a)).

We wish to predict the magnetization dynamics of this system. Unfortunately,
deriving the governing equation of the collective magnetization dynamics from
a microscopic model is a highly non-trivial task that requires a beyond-
saddle-point approximation we employed above [97]. Although deriving such a
coarse-grained description from microscopics is an important challenge, here,
we take a phenomenological approach below.

We make the following observations: (a) In the absence of light, the magnetism
$\bm{m}_{a={\rm A,B}}$ should converge to the known equilibrium value. (b)
Since the relaxation towards this state occurs through the Gilbert damping,
their relaxation time for the $a$-layer is expected to be $\tau_{\rm
Gil}^{\alpha}=O([\alpha_{a}\sum_{b}j_{ab}]^{-1})$, where $j_{ab}(=j_{ba})$ is
the interaction strength $J_{ab}$ multiplied by the number of spins a given
spin couples to. (c) When light is injected into the B layer, the light-
induced torque $\Omega_{\rm BA}$ and a decay with the rate $\gamma_{\rm B}$ is
introduced to B layer magnetization $\bm{m}_{\rm B}$. This brings us to
propose the following phenomenological meanfield description:

$\displaystyle\dot{\bm{m}}_{\rm A}=\alpha_{\rm A}\bigg{[}-k_{\rm
B}T\bm{m}_{\rm A}+\Big{[}1-\frac{(\bm{h}_{\rm eff}^{\rm A})^{2}}{3(k_{\rm
B}T)^{2}}\Big{]}\bm{h}_{\rm eff}^{\rm A}\bigg{]},$ (8a)
$\displaystyle\dot{\bm{m}}_{\rm B}=\alpha_{\rm B}\bigg{[}-k_{\rm
B}T\bm{m}_{\rm B}+\Big{[}1-\frac{(\bm{h}_{\rm eff}^{\rm B})^{2}}{3(k_{\rm
B}T)^{2}}\Big{]}\bm{h}_{\rm eff}^{\rm B}\bigg{]}$ $\displaystyle\ \ \ \ \ \ \
-\gamma_{\rm B}n\bm{m}_{\rm B}-\Omega_{\rm BA}\bm{m}_{\rm A},$ (8b)

where $\bm{h}_{\rm eff}^{a}=\sum_{b={\rm A,B}}j_{ab}\bm{m}_{b}$ is the
effective field applied to $\bm{m}_{a}$ and $\Omega_{\rm BA}=(\gamma_{\rm
B}/|g_{\rm B}|)j_{\rm AB}$ is the light-induced torque. For simplicity, we
have assumed that the system is close enough to the disordered-ordered
transition point that the Ginzburg-Landau expansion is justified and the
anisotropy is strong enough that $z$-component of the magnetization vanishes
$m_{a}^{z}=0$. In the absence of light $\gamma_{\rm B}=0$, the steady state
$\dot{\bm{m}}_{a}=0$ reproduces the known result from the Weiss theory. This
is of the general form introduced in Ref. [15].

Figure 1(c) shows the phase diagram obtained by simulating Eq. (8). In the
absence of the light injection $\gamma_{\rm B}=0$, unsurprisingly, the
magnetization orientation of the two layers aligns $\Delta\varphi=\varphi_{\rm
A}-\varphi_{\rm B}=0$ (Fig. 4(b)), where the orientation of the magnetism is
defined by $\bm{m}_{a}=|\bm{m}_{a}|(\cos\varphi_{a},\sin\varphi_{a})$. As one
increases the laser power that increases $\gamma_{\rm B}$, the light-induced
torque $\Omega_{\rm BA}$ weakens the ferromagnetic interaction, until it swaps
the sign at $\Omega_{\rm BA}\ \raise 1.29167pt\hbox{$>$}\kern-8.00003pt\lower
3.01385pt\hbox{$\sim$}\ \alpha_{\rm B}j_{\rm AB}$ or $\gamma_{\rm B}\ \raise
1.29167pt\hbox{$>$}\kern-8.00003pt\lower 3.01385pt\hbox{$\sim$}\ \alpha_{\rm
B}|g_{\rm B}|$ (see Fig. 4(e)). This causes a transition from aligned
$\Delta\varphi=0$ to antialigned configuration $\Delta\varphi=\pi$ (blue thin
line in Fig. 1(c) and Fig. 4(c)). Remarkably, the B layer completely
demagnetizes at the transition point $|\bm{m}_{\rm B}|=0$, while the A layer
is still ferromagnetic $\bm{m}_{\rm A}\neq 0$ (Fig. 4(f)) even though the
interlayer coupling is still present.

As $\gamma_{\rm B}$ is further increased, the system exhibits a non-reciprocal
phase transition [15, 26] to a time-dependent chiral phase (see Fig. 1(c))
exhibiting a many-body chase-and-runaway motion (Fig. 4(d)). The parity
spontaneously breaks in this chiral phase, where the relative orientation
angle converges to a state $\Delta\varphi=(\neq 0,\pi)$ that is not invariant
under the parity operation $\Delta\varphi\rightarrow-\Delta\varphi$. This
cannot be understood from the Landau theory [15], where the critical point is
characterized by the coalescence of the collective modes to the Nambu-
Goldstone mode [17, 20]. Interestingly, there are regions where this symmetry-
broken phase expands when increasing the temperature, which is opposite from
what is conventionally expected (see the dashed lines in Fig. 1(c) that show
the phase boundary at a lower temperature). This is a signature of order-by-
disorder phenomena discussed in Ref. [26], where a direct analogy between the
geometrically frustrated systems and non-reciprocal matter was drawn.

##  V Discussion

In summary, we have proposed a scheme to dissipation-engineer non-reciprocal
interactions with light. We showed microscopically that the light injection to
magnetic metals that introduces decay of a virtually excited state induces
non-reciprocal interaction between localized spins. Applying this method to
layered ferromagnets, we showed that a non-reciprocal phase transition to a
time-dependent chiral phase emerges [15, 26]. The pump intensity required to
achieve this is estimated to be within reach of the current experimental
techniques.

Recent works showed that non-reciprocal phase transitions exhibit
unconventional critical phenomena associated with anomalously giant
fluctuations [17], fluctuation-induced first-order transition [20], and
diverging entropy production [21, 22, 23], due to the unique feature that the
criticality is driven by the coalescence of the modes to the Nambu-Goldstone
mode. Furthermore, non-reciprocal interactions have been shown to give rise to
emergent features such as odd elasticity [11, 12] (anti-symmetric part of the
static elastic modulus tensor), long-ranged order in two spatial dimensions
[13, 14] (in an apparent violation of the Hohenberg-Mermin-Wagner theorem),
and phenomena analogous to those occurring in geometrically frustrated systems
[26]. Our dissipation-engineering scheme may allow exploring these non-
reciprocal physics in solid-state platforms.

Methods

## Estimation of the required power

Below, we estimate the required laser power $P$ to realize the sign-inversion
of the interactions, which occurs when the decay rate of the double-occupied
state $\kappa_{a}$ exceeds $\alpha_{a}U_{a}$ (see main text and Fig. 3). Our
scheme considers the situation where the double-occupied (at energy
$\varepsilon_{d,a}+U_{a}$) and higher-level states (at energy
$\varepsilon_{f,a}$) are coupled through the injected laser. We assume that
the higher-level state is localized and dissipates with the rate
$\Gamma_{f,a}$, so one can model it with a Lorentz oscillator model [98].

Figure 5:  Light-injection induced dissipation and their energy scales. (a)
The double-occupied (at the energy $\varepsilon_{d,a}+U_{a}$) and the higher
energy states (at the energy $\varepsilon_{f,a}$) are coupled by the injection
of a resonant laser $h\nu=\varepsilon_{f,a}-(\varepsilon_{d,a}+U_{a})$. The
higher-energy state dissipates with the rate $\Gamma_{f,a}$. The localized
electrons are typically in a single-occupied state but may virtually excite
once in a while to a double-occupied state via the c-d mixing $v_{a}$. Note
that no electrons decay in the absence of c-d mixing $v_{a}=0$ because they
are always in the single-occupied state. (b) Localized spin picture obtained
after projecting out the double-occupied states (see Eq. (5)). In this
picture, one finds that the effective transfer rate from the localized
electron to the higher-energy state is given by
$\gamma_{a}\approx\kappa_{a}|v_{a}|^{2}/U_{a}^{2}$. We require
$\gamma_{a}\ll\Gamma_{f,a}$ to justify the Markov approximation.

When a laser with the pump power $P$ is injected into the material, the
dissipation causes the energy loss of the laser intensity if the system is in
a double-occupied state. The lost energy density per unit time and volume $W$
is given by

$\displaystyle
W=\frac{1}{2}\epsilon_{0}\omega\chi^{\prime\prime}(\omega)|E|^{2}=\frac{\omega\chi^{\prime\prime}(\omega)}{c}P,$
(9)

where $\omega=2\pi\nu$ is the laser frequency. Here, we have expressed the
pump power $P=\frac{1}{2}c\epsilon_{0}|E|^{2}$ in terms of speed of light $c$,
vacuum dielectric constant $\epsilon_{0}$, and electric field $E$. The
absorption susceptibility $\chi^{\prime\prime}(\omega)$ is computed according
to the Lorentz oscillator model as

$\displaystyle\chi^{\prime\prime}(\omega)=\frac{ne^{2}}{\epsilon_{0}m_{0}}\frac{\omega\Gamma_{f,a}}{(\omega^{2}_{0}-\omega^{2})^{2}+\omega^{2}\Gamma^{2}_{f,a}}\simeq\frac{ne^{2}}{\epsilon_{0}m_{0}}\frac{1}{\omega_{0}\Gamma_{f,a}}.$
(10)

Here, $\omega_{0}$ is the resonant frequency (which, in our case, corresponds
to $\hbar\omega_{0}=\varepsilon_{f,a}-(\varepsilon_{d,a}+U_{a})$), $n$ is the
number of electrons per unit volume, and $m_{0}$ is the electron mass. In the
second equality, we have set the laser frequency to be on resonance
$h\nu=\hbar\omega=\hbar\omega_{0}$.

The decay rate of the double-occupied state $\kappa_{a}$ per electron is
estimated as

$\displaystyle\kappa_{a}=\frac{W}{n\cdot\omega_{0}}.$ (11)

This needs to be above $\alpha_{a}U_{a}$ to achieve the regime for showing
laser-induced magnetization switching. This condition is given by,

$\displaystyle\kappa_{a}=\frac{\chi^{\prime\prime}(\omega_{0})P}{n\cdot c}\
\raise 1.29167pt\hbox{$>$}\kern-8.00003pt\lower 3.01385pt\hbox{$\sim$}\
\alpha_{a}U_{a}.$ (12)

This yields the condition,

$\displaystyle P\ \raise 1.29167pt\hbox{$>$}\kern-8.00003pt\lower
3.01385pt\hbox{$\sim$}\ \frac{\alpha_{a}U_{a}n\cdot
c}{\chi^{\prime\prime}(\omega_{0})}=\alpha_{a}\frac{U_{a}\omega_{0}m_{0}c\epsilon_{0}}{e^{2}}\Gamma_{f,a},$
(13)

shown in Eq. (7).

## Justification of Markov approximation

For the Markov approximation to be valid, the relaxation rate of the bath must
be much faster than the timescale of the system dynamics. We argue here that
this is likely to be justified in the range of interest at realistic
parameters for magnetic metals.

In our setup, the double-occupied state couples to the higher-level state with
the decay rate $\Gamma_{f,a}$ via a laser injection tuned to be resonant with
the two states. The higher-level state can be regarded as our “external bath”
in the context of open quantum systems. As explained above, this process gives
rise to the decay rate $\kappa_{a}$ once the site $a$ is double occupied (see
Fig. 5(a)). Note crucially that, as illustrated in Fig. 5(b), this is
different from the rate at which the relevant system pumps an electron to this
high-energy state because the relevant system transfers into a double-occupied
state only once in a while when the conduction electron tunnels to the
localized orbit. For example, when the c-d mixing is absent $v_{a}=0$, there
would be no electron transfer from the relevant system to the higher-energy
state. The relevant transfer rate of a localized electron from the relevant
system to the higher-energy state is estimated to be
$\gamma_{a}\simeq\kappa_{a}|v_{a}|^{2}/U_{a}^{2}$, which we have derived in
Eq. (4).

For the Markov approximation to be valid, the dissipation rate of this higher
energy state (i.e. the “external bath”) $\Gamma_{f,a}$ must be much faster
than the supply rate to this state $\gamma_{a}\ll\Gamma_{f,a}$. This is
because when the higher-level state is occupied, the Pauli blocking effect
would suppress the decay. A slow relaxation of the occupancy of the state
would lead to a non-Markovian effect.

It should be relatively easy to satisfy this Markov condition at the regime of
interest. We are interested in the regime where we see the sign-reversal of
the RKKY interactions, which happens when the dissipation rate exceeds
$\kappa_{a}\ \raise 1.29167pt\hbox{$>$}\kern-8.00003pt\lower
3.01385pt\hbox{$\sim$}\ \alpha_{a}U_{a}$ or $\gamma_{a}\ \raise
1.29167pt\hbox{$>$}\kern-8.00003pt\lower 3.01385pt\hbox{$\sim$}\
\alpha_{a}|g_{a}|$ (see the discussion above Eq. (7)). In the case
$\alpha_{a}=10^{-2}$ and $|g_{a}|=10{\rm meV}$, this sets the condition,
$\gamma_{a}\ \raise 1.29167pt\hbox{$>$}\kern-8.00003pt\lower
3.01385pt\hbox{$\sim$}\ 0.1{\rm meV}$. This required dissipation rate is less
than the typical value of linewidth $\Gamma_{f,a}$, satisfying the
justification condition for the Markov approximation,
$\gamma_{a}\ll\Gamma_{f,a}$.

It is worth summarizing the requirements for our scheme to work. Putting the
conditions in one line,

$\displaystyle\alpha_{a}U_{a}\frac{|v_{a}|^{2}}{U_{a}^{2}}\ \raise
1.29167pt\hbox{$<$}\kern-8.00003pt\lower 3.01385pt\hbox{$\sim$}\
\frac{Pe^{2}}{\epsilon_{0}m_{0}c}\frac{1}{\omega_{0}\Gamma_{f,a}}\frac{|v_{a}|^{2}}{U_{a}^{2}}=\gamma_{a}\ll\Gamma_{f,a}.$
(14)

This shows that the required pump power $P$ to achieve sign inversion becomes
less by choosing the higher-energy state that has a longer lifetime (i.e.
smaller $\Gamma_{f,a}$) as shown in Eq. (7) but small $\Gamma_{f,a}$ makes it
more difficult to satisfy the Markovian condition $\gamma_{a}\ll\Gamma_{f,a}$.
We remark that a smaller c-d mixing $|v_{a}|/U_{a}$ helps satisfy the
Markovian condition $\gamma_{a}\ll\Gamma_{f,a}$ without modifying the
condition for the sign inversion, which physically makes sense because small
$v_{a}$ makes the double-occupied state rarer and hence rarer for the
electrons to escape from the relevant system.

## Comparison to cold atom experiments

It is interesting to compare our proposal relevant to solid-state systems to
the recent cold atom experiment of dissipative Fermi Hubbard model [85, 86],
where they demonstrated a dynamic sign reversal of interactions. In their
experiment, similar to our situation, they introduced a controlled decay
channel that is activated only when the sites are double occupied, causing
both atoms to decay whenever they are on the same site. They demonstrated that
this engineered dissipation decreases the anti-ferromagnetic correlation
(present in the ground state) and increases the ferromagnetic correlation,
again analogous to our sign reversal of effective interactions. We briefly
note, however, that the possibility of implementing non-reciprocal
interactions was not explored in their work.

In contrast to our solid-state case, where the coupling to the environment is
unavoidable, their cold atomic systems are almost perfectly isolated from the
environment other the decay to the vacuum they purposely introduced. This
fundamental difference critically impacts the resulting dynamics. First, in
our proposal, the surrounding environment (i.e. the conduction band)
immediately compensates for the lost electrons such that localized orbitals
are always singly occupied, while their atomic system only has a loss channel.
As a result, their dissipation stops activating when the atoms stop colliding,
resulting in a strong initial state dependence on the final configuration.
This is in stark contrast to ours, where no initial-state dependence is
present and even exhibits persistent time-dependent to a collective chase and
runaway phase due to non-reciprocal interaction. Second, the Gilbert damping
(present in our system) is absent in the cold atomic systems, as there is no
environment where the atoms can dissipate their spin angular momentum. The
competition between the friction (that drives the system toward equilibrium)
and the light-induced dissipative interaction is hence a unique feature of our
proposal for solid-state systems.

Acknowledgement

RH was supported by a Grant-in-Aid for Research Activity Start-up from JSPS in
Japan (No. 23K19034) and the National Research Foundation (NRF) funded by the
Ministry of Science of Korea (Grant No. RS-2023-00249900). DO was supported by
a Grant-in-Aid for Early-Career Scientists from JSPS in Japan (No. 21K13882).
RT was supported by a Grant-in-Aid for Early-Career Scientists (No. 22K14003)
and for Research Activity Start-up (No. 20K22328) from JSPS in Japan.

SUPPLEMENTARY INFORMATION

This Supplemental Information (SI) provides full details of our formalism.

In Sec. I, we give a brief review of the formalism of open quantum systems,
focusing on those used in this work. We provide theoretical methods to treat
Markovian (Sec. I.1) and non-Markovian (Sec. I.2) baths.

In Sec. II, we introduce our model.

In Sec. III, we derive Eq. (4) in the main text, which describes the effective
localized spin dynamics immersed in conduction electrons with light-induced
decay. In the derivation, we use the results from the projection method
introduced in Sec. I.1.2.

In Sec. IV, we derive the Landau-Lifshitz-Gilbert equation with light-induced
interactions (Eq. (6) in the main text), which is our main result of this
work. We use the Keldysh formalism introduced in Sec. I.2 in the derivation.

##  I  Mini-review: Theoretical description of open quantum systems

This section briefly reviews the theoretical description of open quantum
systems. Open quantum systems are systems where the relevant system is coupled
to baths that model the effect of the surrounding environment. When the
relevant system is (weakly) attached to a single thermal bath, the system
reaches a thermal equilibrium with the bath. In contrast, when the relevant
system is attached to two (or more) baths, an energy current may constantly
flow from one bath to the system that dissipates to the other baths, driving
the state into a nonequilibrium steady state. The latter is the situation we
will be interested in in this work.

The full microscopic information of the total system composed of the relevant
system plus the baths, described by the total Hamiltonian $\hat{H}_{\rm
tot}=\hat{H}_{s}+\hat{H}_{b}+\hat{H}_{sb}$, is encoded in the density operator
of the full system $\hat{\rho}_{\rm tot}$. Here, $\hat{H}_{s(b)}$ is the
relevant system (bath) Hamiltonian that is composed of operators acting on the
relevant system (bath), and $\hat{H}_{sb}$ is the system-bath coupling. The
bath Hamiltonian $\hat{H}_{b}$ may contain multiple baths, in which case
drives the system into a nonequilibrium steady state. The dynamics are
governed by Liouville’s equation,

$\displaystyle\partial_{t}\hat{\rho}_{\rm tot}=-i[\hat{H}_{\rm
tot},\hat{\rho}_{\rm tot}].$ (S1)

We will be interested in the properties of the reduced-density matrix
$\hat{\rho}={\rm Tr}_{\rm B}\hat{\rho}_{\rm tot}$, where the baths’ degrees of
freedom are traced out. When the relevant system is weakly attached to a
single thermal bath, the system generically reaches an equilibrium state,
which is fully characterized by eigenenergy $E_{n}$ and eigenstates
$\lvert\phi_{n}\rangle$ of the system Hamiltonian $\hat{H}_{s}$
($\hat{H}_{s}\lvert\phi_{n}\rangle=E_{n}\lvert\phi_{n}\rangle$), thanks to the
equilibrium statistical mechanics that tells us that the reduced density
operator $\hat{\rho}$ reaches the canonical ensemble (when attached to a
thermal bath at the temperature $T$),

$\displaystyle\hat{\rho}_{\rm eq}=\frac{1}{Z}e^{-\hat{H}_{s}/(k_{\rm
B}T)}=\frac{1}{Z}\sum_{n}e^{-E_{n}/(k_{\rm
B}T)}\lvert\phi_{n}\rangle\langle\phi_{n}\rvert.$ (S2)

Here, $Z={\rm tr}[e^{-\hat{H}_{s}/(k_{\rm B}T)}]$ is the partition function,
where ${\rm tr}[\cdots]$ is the trace over the relevant system degrees of
freedom, and $k_{\rm B}$ is the Boltzmann constant. Therefore, once we compute
the eigenenergy and the eigenstates, the density operator $\hat{\rho}$ is
fully determined and it is not necessary to compute Liouville’s equation (S1)
directly.

However, for systems attached to multiple baths, the system converges to non-
equilibrium systems, where such a general relation is unknown. Therefore, one
is required to directly analyze the dynamics of the density operator
$\hat{\rho}={\rm Tr}_{\rm B}\hat{\rho}_{\rm tot}$. This is indeed what we will
do in the following.

Roughly speaking, there are two types of baths, i.e. Markovian and non-
Markovian baths. Markovian baths are baths where they relax fast enough such
that the evolution of the relevant system is determined solely by the state of
at the instantaneous time; memory effects are absent. Non-Markovian baths are
the baths where such memory effects cannot be ignored. As we explain in Secs.
III and IV, our dissipation-engineering scheme requires treating both types of
baths and therefore, we give a short review of their theoretical description
of both cases.

This section is organized as follows. In Sec. I.1, we review the theoretical
treatment of the quantum master equation relevant to systems attached to
Markovian baths. After introducing the Gorini-Kossakowski-Sudarshan-Lindblad
(GKSL) master equation [39, 40], which is the most general description of
Markov dynamics of the reduced density operator $\hat{\rho}$, we describe the
general properties of the GKSL master equation in Sec. I.1. We then introduce
a projection method we use to derivate Eq. (4) in the main text. In Sec. I.2,
we introduce the so-called Keldysh formalism, which allows us to treat non-
Markovian open systems on a field-theoretical basis. This will be used to
derive Eq. (6) in the main text.

###  I.1 Markov environment: the GKSL master equation

In this section, we will be concerned with situations where the system is
weakly attached to a bath that relaxes fast enough such that the relevant
system dynamics can be regarded as Markovian. One strategy to treat this
system is to start from Liouville’s equation (S1) and trace out the bath
degrees of freedom. One can microscopically derive the governing equation for
reduced density matrix $\hat{\rho}(t)$ dynamics, e.g., by employing a second-
order Born-Markov approximation in terms of system-bath coupling strength (see
e.g., [99]). Here, we take a more phenomenological approach, where we take
advantage of the property that in Markovian dynamics, by definition, the
density matrix $\hat{\rho}(t+\Delta t)$ is determined solely by transforming
the density matrix $\hat{\rho}(t)$ [i.e., $\hat{\rho}(t+\Delta t)$ does not
depend on the density matrix before the time $t$] and that the density
operator must fulfill certain mathematical properties [99, 88].

Let us denote such a transformation as $\hat{\rho}(t+\Delta
t)={\mathcal{V}}\hat{\rho}(t)$, where ${\mathcal{V}}$ operates on the density
matrix $\hat{\rho}$. [Operators that act on operators are called
superoperators.] There are several properties that this operation must
satisfy, arising from demanding that the reduced operator $\hat{\rho}(t+\Delta
t)$ retains all the properties density operators must satisfy. The first
property is the trace-preservation, i.e., ${\rm tr}[\hat{\rho}(t+\Delta
t)]={\rm tr}[{\mathcal{V}}\hat{\rho}(t)]={\rm tr}[\hat{\rho}(t)](=1)$, which
is required because the eigenvalues $p_{n}$ of the density operator
$\hat{\rho}=\sum_{n}p_{n}\lvert n\rangle\langle n\rvert$ describes the
probability of being in the state $\lvert n\rangle$, which must sum to unity
$\sum_{n}p_{n}=1$.

The second property is a property called complete positivity. Positivity is
the property that all the eigenvalues $p_{n}$ are positive, which is again
required because $p_{n}$ is the probability and must be positive. The
operation ${\mathcal{V}}$ is required to satisfy a stronger condition called
complete positivity, where not only the reduced density operator $\hat{\rho}$
but the total system $\hat{\rho}_{\rm tot}$ stays positive after the operation
$[{\mathcal{V}}\otimes{\mathcal{I}}]\hat{\rho}_{\rm tot}$, where
${\mathcal{V}}$ operates on the relevant system and ${\mathcal{I}}$, the
identity superoperator, acts on the bath degrees of freedom. See e.g., [99]
for details. The map that satisfies both these properties is called the CPTP
map.

Under this condition, it has been shown [39, 40] that this operation must have
the form (see e.g. Ref. [99] for the proof),

$\displaystyle\hat{\rho}(t+\Delta
t)={\mathcal{V}}\hat{\rho}(t)=\hat{\rho}(t)-i[\hat{H},\hat{\rho}(t)]\Delta
t+\sum_{\ell}\Big{[}\hat{L}_{\ell}\hat{\rho}(t)\hat{L}_{\ell}^{\dagger}-\frac{1}{2}\\{\hat{L}_{\ell}^{\dagger}\hat{L}_{\ell},\hat{\rho}(t)\\}\Big{]}\Delta
t+O((\Delta t)^{2})$ (S3)

or in the $\Delta t\rightarrow 0$ limit,

$\displaystyle\partial_{t}\hat{\rho}$ $\displaystyle=$
$\displaystyle{\mathcal{L}}\hat{\rho}=-i[\hat{H},\hat{\rho}]+\sum_{\ell}\kappa_{\ell}{\mathcal{D}}[\hat{L}_{\ell}]\hat{\rho}.$
(S4)

This is called the GKSL or Lindblad master equation [39, 40]. Here, the
superoperator ${\mathcal{L}}$ is called the Lindbladian. $\hat{H}$ is the
Hamiltonian of the system which determines the coherent dynamics. The
dissipator,
${\mathcal{D}}[\hat{L}]\hat{\rho}=\hat{L}\hat{\rho}\hat{L}^{\dagger}-\frac{1}{2}\\{\hat{L}^{\dagger}\hat{L},\hat{\rho}\\}$,
is the term that makes the time evolution non-unitary. The so-called jump
operator $\hat{L}_{\ell}$ consists of operators acting on the system, which
describes how the environment interacts with the system with the dissipation
rate $\kappa_{\ell}$ in channel $\ell$. For instance, suppose the system
consists of bosonic/fermionic particles and one takes $\hat{L}_{\ell}=\hat{a}$
(where $\hat{a}$ is an annihilation operator of a particle in the system) and
$\kappa_{\ell}=\kappa$. This corresponds to the case where the environment
induces one-body loss from the system with the loss rate $\kappa$. When one
takes $\hat{L}=\hat{a}^{\dagger}$, on the other hand, it describes the
situation where the particles are pumped into the system. Taking
$\hat{L}=\hat{a}^{2}$ means two-body loss, etc.

####  I.1.1  General properties of the GKSL master equation

Here we review some general properties of the GKSL master equation (S4)
relevant to our study. To formulate them, it is useful to express the GKSL
master equation in terms of the vectorized form to make explicit connections
to the more familiar Schrödinger equation. The vectorization (Choi-
Jamiolkowski isomorphism) is chosen to be performed by the replacement [88],

$\displaystyle\lvert i\rangle\langle j\rvert\rightarrow\lvert
j\rangle\otimes\lvert i\rangle.$ (S5)

For example, an operator $\hat{A}=\sum_{i,j}a_{i,j}\lvert i\rangle\langle
j\rvert$ is vectorized as

$\displaystyle{\rm vec}(\hat{\rho})=\sum_{i,j}a_{i,j}\lvert
j\rangle\otimes\lvert i\rangle.$ (S6)

Their dual $\hat{A}^{\dagger}=\sum_{i,j}a_{i,j}^{*}\lvert j\rangle\langle
i\rvert$ is vectorized as

$\displaystyle[{\rm vec}(\hat{A})]^{\dagger}=\sum_{i,j}a_{i,j}^{*}\langle
j\rvert\otimes\langle i\rvert.$ (S7)

This choice of vectorization is useful, thanks to the following properties.
First, the Hilbert-Schmidt inner product
$(\hat{A},\hat{B})={\rm{tr}}[{\hat{A}^{\dagger}\hat{B}}]$ can be expressed as

$\displaystyle(\hat{A},\hat{B})$ $\displaystyle=$
$\displaystyle{\rm{tr}}[{\hat{A}^{\dagger}\hat{B}}]=\sum_{i,j,i^{\prime},j^{\prime}}a_{i,j}^{*}b_{i^{\prime},j^{\prime}}{\rm
tr}[(\lvert i\rangle\langle j\rvert)^{\dagger}(\lvert i^{\prime}\rangle\langle
j^{\prime}\rvert)]=\sum_{i,j}a_{i,j}^{*}b_{i,j}$ (S8) $\displaystyle=$
$\displaystyle\sum_{i,j,i^{\prime},j^{\prime}}(a_{j,i}^{*}\langle
i\rvert\otimes\langle j\rvert)(b_{i^{\prime},j^{\prime}}\lvert
j^{\prime}\rangle\otimes\lvert i^{\prime}\rangle)=[{\rm
vec}(\hat{A})]^{\dagger}{\rm vec}(\hat{B})$

which is just the inner product of the two vectors. Another useful property is
the relation

$\displaystyle{\rm
vec}(\hat{A}\hat{B}\hat{C})=(\hat{C}^{\mathrm{T}}\otimes\hat{A}){\rm
vec}({\hat{B}}),$ (S9)

which can be shown similarly by decomposing the matrix in terms of the chosen
eigenbasis. This relation implies (where $\hat{I}=\sum_{i}\lvert
i\rangle\otimes\langle i\rvert$ is the identity operator)

$\displaystyle{\rm vec}(\hat{A}\hat{\rho}\hat{B})$ $\displaystyle=$
$\displaystyle(\hat{B}^{\mathrm{T}}\otimes\hat{A}){\rm vec}(\hat{\rho})$ (S10)
$\displaystyle{\rm vec}(\hat{A}\hat{\rho})$ $\displaystyle=$
$\displaystyle{\rm vec}(\hat{A}\hat{\rho}\hat{I})=(\hat{I}\otimes\hat{A}){\rm
vec}(\hat{\rho})$ (S11) $\displaystyle{\rm vec}(\hat{\rho}\hat{A})$
$\displaystyle=$ $\displaystyle{\rm
vec}(\hat{I}\hat{\rho}\hat{A})=(\hat{A}^{\mathrm{T}}\otimes\hat{I}){\rm
vec}(\hat{\rho}),$ (S12)

allowing us to express the GKSL master equation (S4) in the vectorized form,

$\displaystyle\partial_{t}\rm{vec}(\hat{\rho})=\breve{\mathcal{L}}\rm{vec}(\hat{\rho}),$
(S13)

where

$\displaystyle\breve{\mathcal{L}}=-i(\hat{I}\otimes\hat{H}-\hat{H}^{\mathrm{T}}\otimes\hat{I})+\sum_{\ell}\kappa_{\ell}\Big{[}\hat{L}_{\ell}^{*}\otimes\hat{L}_{\ell}-\frac{1}{2}\hat{I}\otimes\hat{L}_{\ell}^{\dagger}\hat{L}_{\ell}-\frac{1}{2}(\hat{L}_{\ell}^{\dagger}\hat{L}_{\ell})^{\mathrm{T}}\otimes\hat{I}\Big{]}.$
(S14)

Equation (S13) formally has a similar form to the Schrödinger equation of
quantum mechanics. However, unlike in the Schrödinger equation, the operator
$\check{\mathcal{L}}$ that governs the dynamics is non-Hermitian. In those
cases, it is useful to characterize the $\check{\mathcal{L}}$ with both the
right and left eigenvectors ($\hat{r}_{n}$ and $\hat{l}_{n}$, respectively),

$\displaystyle\breve{\mathcal{L}}{\rm{vec}}(\hat{r}_{n})=\lambda_{n}{\rm{vec}}(\hat{r}_{n}),\qquad\breve{\mathcal{L}}^{\dagger}{\rm{vec}}(\hat{l}_{n})=\lambda_{n}^{*}{\rm{vec}}(\hat{l}_{n}),$
(S15)

where the eigenvalues $\lambda_{n}$ are sorted as
$\lambda_{N-1}\leq\lambda_{N-2}\leq...\leq\lambda_{2}\leq\lambda_{1}\leq\lambda_{0}$.
Note that in Hermitian cases, $\hat{r}_{n}=\hat{l}_{n}$ but this does not
generally hold otherwise. In the operator form, these read

$\displaystyle{\mathcal{L}}\hat{r}_{n}=\lambda_{n}\hat{r}_{n},\qquad\hat{l}_{n}^{\dagger}{\mathcal{L}}=\hat{l}_{n}^{\dagger}\lambda_{n},$
(S16)

which form a biorthogonal and complete basis,

$\displaystyle{\rm vec}(\hat{l}_{n})^{\dagger}{\rm
vec}(\hat{r}_{m})=\delta_{nm},\qquad\sum_{n}{\rm vec}(\hat{r}_{n}){\rm
vec}(\hat{l}_{n})^{\dagger}=\hat{I}.$ (S17)

Here, we used the property that an inner product is properly defined (Eq.
(S8)). One can then decompose the density operator into eigenmodes of
$\breve{\mathcal{L}}$:

$\displaystyle{\rm vec}(\hat{\rho}(t))=e^{\breve{\mathcal{L}}t}{\rm
vec}(\hat{\rho}(0))=\sum_{n}e^{\lambda_{n}t}{\rm vec}(\hat{r}_{n}){\rm
vec}(\hat{l}_{n})^{\dagger}\hat{\rho}(0),$ (S18)

which is equivalent to

$\displaystyle\hat{\rho}(t)=e^{{\mathcal{L}}t}\hat{\rho}(0)=\sum_{n}e^{\lambda_{n}t}\hat{r}_{n}{\rm
tr}[\hat{l}_{n}^{\dagger}\hat{\rho}(0)].$ (S19)

As one sees from here, $-{\rm Re}\lambda_{n}$ characterizes the decay rate of
the mode $n$ and ${\rm Im}\lambda_{n}$ characterizes its oscillation
frequency.

We remark that the GKSL master equation (S4) is guaranteed to have at least
one eigenstate with a vanishing eigenvalue $\lambda_{0}=0$, corresponding to
the steady state $\partial_{t}\hat{\rho}_{\rm ss}=0$. This can be shown from
the trace-preserving property of the GKSL master equation that the dynamics
assure ${\rm tr}\hat{\rho}(t)=1$ at arbitrary times,

$\displaystyle\partial_{t}{\rm tr}\hat{\rho}={\rm
tr}[{\mathcal{L}}\hat{\rho}]=0.$ (S20)

In the vectorized GKSL master equation, this reads

$\displaystyle 0={\rm vec}(\hat{I})^{\dagger}\breve{\mathcal{L}}{\rm
vec}(\hat{\rho}).$ (S21)

Here, we have used the relation ${\rm tr}(\hat{A})={\rm
vec}(\hat{I})^{\dagger}{\rm vec}(\hat{A})$ where ${\rm
vec}(\hat{I})=\sum_{i}\lvert i\rangle\otimes\lvert i\rangle$. For this to be
true for arbitrary $\hat{\rho}$,

$\displaystyle 0={\rm vec}(I)^{\dagger}\breve{\mathcal{L}}$ (S22)

must be satisfied. This assures the presence of a left eigenstate given by

$\displaystyle\hat{l}_{0}=\hat{I},$ (S23)

with zero eigenvalue $\lambda_{0}=0$. The right eigenstate is the state
satisfying

$\displaystyle{\mathcal{L}}(\hat{r}_{0})=0,$ (S24)

which is nothing but a steady state $\hat{r}_{0}=\hat{\rho}_{ss}$. Using this
property, we find

$\displaystyle\hat{\rho}(t)=\hat{\rho}_{ss}+\sum_{n\neq
0}e^{\lambda_{n}t}\hat{r}_{n}{\rm tr}[\hat{l}_{n}^{\dagger}\hat{\rho}(0)].$
(S25)

Here, the coefficient in front of $\hat{\rho}_{ss}$ is unity ${\rm
tr}[\hat{l}_{0}^{\dagger}\hat{\rho}(0)]=1$, which can be shown from
orthogonality relation (S17) that implies ${\rm
tr}[\hat{l}_{m}\hat{r}_{n}]=\delta_{n,m}$ and $\hat{l}_{0}=\hat{I}$. One also
finds that the $n\neq 0$ right eigenstates are traceless, from ${\rm
tr}[\hat{r}_{n}]={\rm tr}[\hat{l}_{0}^{\dagger}\hat{r}_{n}]=0~{}(n\neq 0)$.

####  I.1.2 Projection method

There are many situations in condensed matter physics and quantum optics where
there is a separation of timescales in the problem. In such cases, it is often
useful to project out the higher-energy degrees of freedom to obtain an
effective low-energy description [41, 89, 90]. These approaches have been
proven useful in describing, e.g., electrons in a potential perturbed by
impurities can be described by an effective $k\cdot p$ theory [100]; the
Fermi-Hubbard model at the half-filling can be mapped to an anti-
ferromagnetically interacting spin model [101]; the Anderson impurity model
can be mapped to an sd exchange coupling model [41] that has successfully
described the Kondo effect [82].

Here, we extend this concept to an open quantum system, which will be used in
Sec. III. We develop here a projection method that is consistent with the so-
called Lindblad perturbation theory [42, 43, 44], which is similar to the
generalized Schrieffer-Wolff transformation technique [102, 103, 104] in
spirit. We consider a Lindbladian that can be split into the non-perturbative
and perturbative part
${\mathcal{L}}={\mathcal{L}}_{0}+\epsilon{\mathcal{L}}_{1}$ (where
$\epsilon=1$ is a book-keeping constant). The non-perturbative Lindbladian
master equation,

$\displaystyle\partial_{t}{\rm{vec}}(\hat{\rho}^{(0)})=\breve{\mathcal{L}}_{0}\rm{vec}(\hat{\rho}^{(0)})$
(S26)

is assumed to have a separation of timescales. We label a set of eigenvectors
$\hat{r}_{n}[n\in\mathfrak{s(f)}]$ with small [in magnitude] eigenvalues
[$|\lambda_{n\in\mathfrak{s}}|<|\lambda_{n\in\mathfrak{f}}|$] as slow (fast)
modes. We split the density operator into slow and fast part
$\hat{\rho}^{(0)}=\hat{\rho}_{s}^{(0)}+\hat{\rho}_{f}^{(0)}$. Here,
$\hat{\rho}_{s}^{(0)}={\mathcal{P}}\hat{\rho}^{(0)}=\sum_{n\in\mathfrak{s}}c_{n}\hat{r}_{n}^{(0)}$
is the density operator spanned projected to the slow variable space (where
${\mathcal{P}}\hat{\rho}=\sum_{n\in\mathfrak{s}}\hat{r}_{n}^{(0)}{\rm
tr}[\hat{l}_{n}^{(0)\dagger}\hat{\rho}]$ is a projection operator onto the
slow space) and
$\hat{\rho}_{f}^{(0)}={\mathcal{Q}}\hat{\rho}^{(0)}=({\mathcal{I}}-{\mathcal{P}})\hat{\rho}^{(0)}=\sum_{n\in\mathfrak{f}}c_{n}\hat{r}_{n}^{(0)}$
is the fast variable space (where
${\mathcal{Q}}\hat{\rho}=\sum_{n\in\mathfrak{f}}\hat{r}_{n}^{(0)}{\rm
tr}[\hat{l}_{n}^{(0)\dagger}\hat{\rho}]$ is a projection operator onto the
fast variable space and ${\mathcal{I}}$ is an identity superoperator). The
non-perturbative GKSL master equation
$\partial_{t}\hat{\rho}_{s}^{(0)}+\partial_{t}\hat{\rho}_{f}^{(0)}={\mathcal{L}}_{0}\hat{\rho}_{s}^{(0)}+{\mathcal{L}}_{0}\hat{\rho}_{f}^{(0)}$
reads

$\displaystyle\partial_{t}\begin{pmatrix}\breve{\mathcal{P}}{\rm
vec}(\hat{\rho}^{(0)})\\\ \breve{\mathcal{Q}}{\rm
vec}(\hat{\rho}^{(0)})\end{pmatrix}=\begin{pmatrix}\breve{\mathcal{P}}\breve{\mathcal{L}}_{0}&0\\\
0&\breve{\mathcal{Q}}\breve{\mathcal{L}}_{0}\end{pmatrix}\begin{pmatrix}\breve{\mathcal{P}}{\rm
vec}(\hat{\rho}^{(0)})\\\ \breve{\mathcal{Q}}{\rm
vec}(\hat{\rho}^{(0)})\end{pmatrix}$ (S27)

in the vectorized representation, where we used the property that these two
spaces do not talk to each other via ${\mathcal{L}}_{0}$. Here, we have
introduced a projection operator in the vectorized representation, ${\rm
vec}({\mathcal{P}}\hat{\rho})\equiv\breve{\mathcal{P}}{\rm vec}(\hat{\rho})$
and similarly for $\breve{\mathcal{Q}}$.

Let us now add back the perturbation ${\mathcal{L}}_{1}$. In such a case, we
find (Note that
$\breve{\mathcal{L}}_{sf},\breve{\mathcal{L}}_{fs}=O(\epsilon)$)

$\displaystyle\partial_{t}\begin{pmatrix}\breve{\mathcal{P}}{\rm
vec}(\hat{\rho})\\\ \breve{\mathcal{Q}}{\rm vec}(\hat{\rho})\end{pmatrix}$
$\displaystyle\approx$
$\displaystyle\begin{pmatrix}\partial_{t}\breve{\mathcal{P}}{\rm
vec}(\hat{\rho})\\\
0\end{pmatrix}=\begin{pmatrix}\breve{\mathcal{L}}_{ss}&\breve{\mathcal{L}}_{sf}\\\
\breve{\mathcal{L}}_{fs}&\breve{\mathcal{L}}_{ff}\end{pmatrix}\begin{pmatrix}\breve{\mathcal{P}}{\rm
vec}(\hat{\rho})\\\ \breve{\mathcal{Q}}{\rm vec}(\hat{\rho})\end{pmatrix}$
(S28) $\displaystyle=$
$\displaystyle\begin{pmatrix}\breve{\mathcal{P}}\breve{\mathcal{L}}&\epsilon\breve{\mathcal{P}}\breve{\mathcal{L}}_{1}\\\
\epsilon\breve{\mathcal{Q}}\breve{\mathcal{L}}_{1}&\breve{\mathcal{Q}}\breve{\mathcal{L}}\end{pmatrix}\begin{pmatrix}\breve{\mathcal{P}}{\rm
vec}(\hat{\rho})\\\ \breve{\mathcal{Q}}{\rm vec}(\hat{\rho})\end{pmatrix}$

where in the second equality, we have done an adiabatic approximation for the
fast modes (i.e., the fast modes relax fast enough such that the fast
variables adiabatically follow the slow variables). This leads to

$\displaystyle\breve{\mathcal{Q}}{\rm
vec}(\hat{\rho})\approx-(\breve{\mathcal{L}}_{ff})^{-1}\breve{\mathcal{L}}_{fs}\breve{\mathcal{P}}{\rm
vec}(\hat{\rho}),$ (S29)

therefore,

$\displaystyle\partial_{t}\breve{\mathcal{P}}{\rm
vec}(\hat{\rho})=\breve{\mathcal{L}}_{\rm eff}\breve{\mathcal{P}}{\rm
vec}(\hat{\rho})=[\breve{\mathcal{L}}_{ss}-\breve{\mathcal{L}}_{sf}(\breve{\mathcal{L}}_{ff})^{-1}\breve{\mathcal{L}}_{fs}]\breve{\mathcal{P}}{\rm
vec}(\hat{\rho})$ (S30)

or

$\displaystyle\breve{\mathcal{L}}_{\rm eff}$ $\displaystyle=$
$\displaystyle\breve{\mathcal{P}}\breve{\mathcal{L}}\breve{\mathcal{P}}-\epsilon^{2}\breve{\mathcal{P}}\breve{\mathcal{L}}_{1}\breve{\mathcal{Q}}(\breve{\mathcal{L}})^{-1}\breve{\mathcal{Q}}\breve{\mathcal{L}}_{1}\breve{\mathcal{P}}$
(S31) $\displaystyle=$
$\displaystyle\breve{\mathcal{P}}\bigg{[}\breve{\mathcal{L}}_{0}+\epsilon\breve{\mathcal{L}}_{1}-\epsilon^{2}\sum_{m\in\mathfrak{f}}\frac{\breve{\mathcal{L}}_{1}{\rm
vec}(\hat{r}_{m}^{(0)}){\rm
vec}(\hat{l}_{m}^{(0)})^{\dagger}\breve{\mathcal{L}}_{1}}{\lambda_{m}^{(0)}}\bigg{]}\breve{\mathcal{P}}+O(\epsilon^{3})$

The matrix element of the effective Lindbladian $({\mathcal{L}}_{\rm
eff})_{n_{l}\in\mathfrak{s},n_{r}\in\mathfrak{s}}\equiv(\hat{l}_{n_{l}\in\mathfrak{s}}^{(0)},{\mathcal{L}_{\rm
eff}\hat{r}_{n_{r}\in\mathfrak{s}}^{(0)}})$, which contains all the
information about the system in the slow dynamics of interest, reads

$\displaystyle({\mathcal{L}}_{\rm eff})_{n_{l},n_{r}}$ $\displaystyle\equiv$
$\displaystyle(\hat{l}_{n_{l}}^{(0)},{\mathcal{L}_{\rm
eff}\hat{r}_{n_{r}}^{(0)}})={\rm
vec}(\hat{l}_{n_{l}}^{(0)})^{\dagger}\breve{\mathcal{L}}_{\rm eff}{\rm
vec}(\hat{r}_{n_{r}}^{(0)})$ (S32) $\displaystyle=$
$\displaystyle({\mathcal{L}}_{0})_{n_{l},n_{r}}+\epsilon({\mathcal{L}}_{1})_{n_{l},n_{r}}-\epsilon^{2}\sum_{m\in\mathfrak{f}}\frac{{\rm
tr}[\hat{l}_{n_{l}}^{(0)\dagger}{\mathcal{L}}_{1}\hat{r}_{m}^{(0)}]{\rm
tr}[\hat{l}_{m}^{(0)\dagger}{\mathcal{L}}_{1}\hat{r}_{n_{r}}^{(0)}]}{\lambda_{m}^{(0)}}+O(\epsilon^{3}).$

Equation (S32) is the central relation that will be used below to derive the
effective spin-exchange dynamics between the localized and conduction
electrons (Eq. (4) in the main text).

###  I.2 Non-Markov environment: Keldysh theory

The previous section discussed open quantum systems whose dynamics are assumed
to be Markovian. However, there are situations where such assumptions are
inappropriate, where the baths do not relax fast enough compared to the
relevant system and the dynamics of the baths play a role as a result. A
prominent example is a many-body system attached to a thermal bath. In this
case, the relevant system generically approaches an equilibrium state, where
general properties such as the fluctuation-dissipation theorem [105] should be
satisfied. However, these relations cannot be satisfied when we neglect the
non-Markovianity arising from Fermi/Bose distribution functions [106]. As we
will be interested in situations where we continuously tune from equilibrium
to non-equilibrium, a formalism that includes such effects is desired.

This subsection briefly reviews the Keldysh approach that can incorporate such
effects. (For further details, we refer to e.g., Refs. [45, 46, 47].) In this
approach, we describe the dynamics of the density matrix $\hat{\rho}$ in terms
of a path integral, which enables us to use field-theoretical techniques such
as perturbative diagrammatic calculations, renormalization group methods, and
makes the connection to the semi-classical limit explicit. This will prove
useful for obtaining the equation of motion of observables and will be used in
the analysis performed in Sec. IV.

####  I.2.1 Keldysh formalism for closed systems

Let us start by formulating the Keldysh theory for closed systems described by
the Hamiltonian $\hat{H}$. The Liouville equation reads,

$i\partial_{t}\hat{\rho}=[\hat{H},\hat{\rho}],$ (S33)

which has a formal solution,

$\displaystyle\hat{\rho}(t)=\hat{U}(t,t_{0})\hat{\rho}(t_{0})\hat{U}^{-1}(t,t_{0}),$
(S34)

with the time evolution operator given by,

$\displaystyle\hat{U}(t,t_{0})=e^{-i\hat{H}(t-t_{0})}.$ (S35)

We set the initial state to be a thermal state at temperature $T$,

$\displaystyle\hat{\rho}(t=t_{0})=\frac{1}{Z_{0}}\sum_{n}e^{-E_{n}/(k_{\rm
B}T)}\lvert n\rangle\langle n\rvert.$ (S36)

Here, $E_{n}$ is an eigenenergy, $\lvert n\rangle$ is the eigenstate, and
$Z_{0}=\sum_{n}e^{-E_{n}/(k_{\rm B}T)}$. We define the Keldysh partition
function as,

$\displaystyle Z\equiv{\rm tr}\hat{\rho}(t_{f})={\rm
tr}[\hat{U}(t_{f},t_{0})\hat{\rho}(t_{0})\hat{U}^{-1}(t_{f},t_{0})].$ (S37)

Although this quantity itself is trivially computed as unity $Z=1$, the
functional integral form of this quantity serves as a generating function that
gives us useful macroscopic information we will be interested in.

To illustrate how this formalism works, we will focus on a simple, single-mode
fermion system $\hat{H}=\varepsilon_{0}\hat{\psi}^{\dagger}\hat{\psi}$, where
$\hat{\psi}$ is a fermionic annihilation operator. Similarly to the path
integral formalism in quantum mechanics, we expand the time evolution operator
into a product of infinitesimally short time interval
$\delta_{t}=(t_{f}-t_{0})/N$ (with $N\rightarrow\infty$) and replace the
annihilation operators $\hat{\psi}_{j}$ ($j$ labels discretized time
$t_{j}=j\delta_{t}$) to a Grassmann variable in the fermionic case $\psi_{i}$
with coherent states

$\displaystyle\hat{\psi}_{j}\lvert\psi_{j}\rangle=\psi_{j}\lvert\psi_{j}\rangle,\qquad\langle\psi_{j}\rvert\hat{\psi}_{j}^{\dagger}=\langle\psi_{j}\rvert\psi_{j}^{*}.$
(S38)

by inserting the completeness relation

$\displaystyle\hat{1}=\int
d\psi_{j}d\psi_{j}^{*}e^{-\psi_{j}^{*}\psi_{j}}\lvert\psi_{j}\rangle\langle\psi_{j}\rvert.$
(S39)

We obtain,

$\displaystyle Z$ $\displaystyle=$ $\displaystyle
Z_{0}^{-1}\int\prod_{n}d\bar{\psi}_{+,n}d\psi_{+,n}d\bar{\psi}_{-,n}d\psi_{-,n}$
(S40) $\displaystyle\times$ $\displaystyle
e^{\sum_{j}[\bar{\psi}_{j,+}(\psi_{j-1,+}-\psi_{j,+}-i\delta_{t}\varepsilon_{0}\psi_{j-1,+})+\bar{\psi}_{j,-}(\psi_{j-1,-}-\psi_{j,-}+i\delta_{t}\varepsilon_{0}\psi_{j-1,-})]}e^{-\bar{\psi}_{1,+}\psi_{1,+}-\bar{\psi}_{1,-}\psi_{1,-}+\zeta\bar{\psi}_{1,+}\psi_{N,-}+\bar{\psi}_{1,-}\psi_{N,+}}$
$\displaystyle\equiv$
$\displaystyle\int{\mathcal{D}}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]e^{iS[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]},$

where we have used

$\displaystyle\langle\psi_{+,n+1}|\psi_{+,n}\rangle=e^{\bar{\psi}_{+,n+1}\psi_{+,n}},\qquad\langle-\psi_{-,n}|-\psi_{-,n+1}\rangle=e^{\bar{\psi}_{-,n}\psi_{-,n+1}}$
(S41)

and
$Z_{0}^{-1}=\langle\psi_{+}(t_{0})\rvert\hat{\rho}(t_{0})\lvert-\psi_{-}(t_{0})\rangle$.
Here, the action has the form $S=\bar{\psi}G^{-1}\psi$, where we have
introduced a set of coherent states for both the time evolution of ket ($+$)
and bra ($-$), collectively denoted as

$\displaystyle\psi=\begin{pmatrix}\psi_{1,+}&\cdots&\psi_{N,+}&\psi_{1,-}&\cdots&\psi_{N,-}\end{pmatrix}^{\rm
T},\qquad\bar{\psi}=\begin{pmatrix}\psi_{1,+}&\cdots&\psi_{N,+}&\psi_{1,-}&\cdots&\psi_{N,-}\end{pmatrix}$
(S42)

with $\zeta=-e^{-\varepsilon_{0}/(k_{\rm B}T)}$ and the kernel is given by,

$\displaystyle G^{-1}=\begin{pmatrix}(G^{-1})_{++}&(G^{-1})_{+-}\\\
(G^{-1})_{-+}&(G^{-1})_{--}\end{pmatrix}$ (S43)

with

$\displaystyle(G^{-1})_{++}$ $\displaystyle=$ $\displaystyle
i\begin{pmatrix}1&0&0&\cdots&0\\\ -a_{+}&1&0&\cdots&0\\\
0&-a_{+}&1&\cdots&0\\\ 0&0&\ddots&\ddots&0\\\
0&0&0&-a_{+}&1\end{pmatrix},\qquad(G^{-1})_{--}=i\begin{pmatrix}1&0&0&\cdots&0\\\
-a_{-}&1&0&\cdots&0\\\ 0&-a_{-}&1&\cdots&0\\\ 0&0&\ddots&\ddots&0\\\
0&0&0&-a_{-}&1\end{pmatrix},$ (S44) $\displaystyle(G^{-1})_{+-}$
$\displaystyle=$ $\displaystyle i\begin{pmatrix}0&0&0&\cdots&-\zeta\\\
0&0&0&\cdots&0\\\ 0&0&0&\cdots&0\\\ 0&0&\ddots&\ddots&0\\\
0&0&0&0&0\end{pmatrix},\qquad(G^{-1})_{-+}=i\begin{pmatrix}0&0&0&\cdots&1\\\
0&0&0&\cdots&0\\\ 0&0&0&\cdots&0\\\ 0&0&\ddots&\ddots&0\\\
0&0&0&0&0\end{pmatrix},$ (S45)

where $a_{\pm}=1\mp i\varepsilon_{0}\delta_{t}$. The action in the continuous-
time limit $N\rightarrow\infty$ reads,

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int
dtdt^{\prime}\begin{pmatrix}\bar{\psi}_{+}(t)&\bar{\psi}_{-}(t)\end{pmatrix}\begin{pmatrix}G_{++}(t,t^{\prime})&G_{+-}(t,t^{\prime})\\\
G_{-+}(t,t^{\prime})&G_{--}(t,t^{\prime})\end{pmatrix}^{-1}\begin{pmatrix}\psi_{+}(t^{\prime})\\\
\psi_{-}(t^{\prime})\end{pmatrix}$ (S46) $\displaystyle=$ $\displaystyle\int
dt\begin{pmatrix}\bar{\psi}_{+}(t)&\bar{\psi}_{-}(t)\end{pmatrix}\begin{pmatrix}i\partial_{t}-\varepsilon_{0}&0\\\
0&-(i\partial_{t}-\varepsilon_{0})\end{pmatrix}\begin{pmatrix}\psi_{+}(t)\\\
\psi_{-}(t)\end{pmatrix}.$

Taking the inverse of $G^{-1}$, one obtains in the continuous-time limit
$N\rightarrow\infty$ (see e.g. Ref. [45] for derivation),

$\displaystyle G_{+-}(t,t^{\prime})$
$\displaystyle=if(\varepsilon_{0})e^{-i\varepsilon_{0}(t-t^{\prime})},$ (S47a)
$\displaystyle G_{-+}(t,t^{\prime})$
$\displaystyle=i(1+f(\varepsilon_{0}))e^{-i\varepsilon_{0}(t-t^{\prime})},$
(S47b) $\displaystyle G_{++}(t,t^{\prime})$
$\displaystyle=\theta(-t)G_{+-}(t,t^{\prime})+\theta(t)G_{-+}(t,t^{\prime}),$
(S47c) $\displaystyle G_{--}(t,t^{\prime})$
$\displaystyle=\theta(t)G_{+-}(t,t^{\prime})+\theta(-t)G_{-+}(t,t^{\prime}).$
(S47d)

where $f(\omega)=1/(e^{\omega/(k_{\rm B}T)}+1)$ is the Fermi distribution
function and $\theta(t)$ is a step function. Their Fourier transformation is
given by,

$\displaystyle G_{+-}(\omega)$ $\displaystyle=2\pi
i\delta(\omega-\varepsilon_{0})f(\varepsilon_{0}),$ (S48a) $\displaystyle
G_{-+}(\omega)$ $\displaystyle=-2\pi
i\delta(\omega-\varepsilon_{0})(1-f(\varepsilon_{0})),$ (S48b) $\displaystyle
G_{++}(\omega)$
$\displaystyle=\frac{1-f(\varepsilon_{0})}{\omega-\varepsilon_{0}+i\delta}+\frac{f(\varepsilon_{0})}{\omega-\varepsilon_{0}-i\delta},$
(S48c) $\displaystyle G_{--}(\omega)$
$\displaystyle=-\frac{f(\varepsilon_{0})}{\omega-\varepsilon_{0}+i\delta}-\frac{1-f(\varepsilon_{0})}{\omega-\varepsilon_{0}-i\delta}.$
(S48d)

These are related to the various two-point functions as [45, 46, 47],

$\displaystyle G_{+-}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-i\left<\psi_{+}(t)\bar{\psi}_{-}(t^{\prime})\right>=-i\int{\mathcal{D}}(\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-})\psi_{+}(t)\bar{\psi}_{-}(t^{\prime})e^{iS}=i\left<\hat{\psi}^{\dagger}(t^{\prime})\hat{\psi}(t)\right>$
(S49) $\displaystyle G_{-+}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-i\left<\psi_{-}(t)\bar{\psi}_{+}(t^{\prime})\right>=-i\int{\mathcal{D}}(\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-})\psi_{-}(t)\bar{\psi}_{+}(t^{\prime})e^{iS}=-i\left<\hat{\psi}(t)\hat{\psi}^{\dagger}(t^{\prime})\right>$
(S50) $\displaystyle G_{++}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-i\left<\psi_{+}(t)\bar{\psi}_{+}(t^{\prime})\right>=-i\int{\mathcal{D}}(\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-})\psi_{+}(t)\bar{\psi}_{+}(t^{\prime})e^{iS}=-i\left<{\mathcal{T}}\hat{\psi}(t)\hat{\psi}^{\dagger}(t^{\prime})\right>$
(S51) $\displaystyle G_{--}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-i\left<\psi_{-}(t)\bar{\psi}_{-}(t^{\prime})\right>=-i\int{\mathcal{D}}(\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-})\psi_{-}(t)\bar{\psi}_{-}(t^{\prime})e^{iS}=-i\left<\bar{\mathcal{T}}\hat{\psi}(t)\hat{\psi}^{\dagger}(t^{\prime})\right>$
(S52)

where ${\mathcal{T}}$ and $\bar{\mathcal{T}}$ are the time-ordered and anti-
time-ordered operators, respectively.

It is often useful to change the variables via a “Wick rotation”,

$\displaystyle\psi_{c}=\frac{1}{\sqrt{2}}(\psi_{+}+\psi_{-}),\qquad\psi_{q}=\frac{1}{\sqrt{2}}(\psi_{+}-\psi_{-}),$
(S53)

to transform the action into,

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int
dtdt^{\prime}\begin{pmatrix}\bar{\psi}_{q}(t)&\bar{\psi}_{c}(t)\end{pmatrix}\begin{pmatrix}G_{qq}(t,t^{\prime})&G_{qc}(t,t^{\prime})\\\
G_{cq}(t,t^{\prime})&G_{cc}(t,t^{\prime})\end{pmatrix}^{-1}\begin{pmatrix}\psi_{q}(t^{\prime})\\\
\psi_{c}(t^{\prime})\end{pmatrix}$ (S54)

where the kernel in this representation is related to the correlation
functions and response functions as,

$\displaystyle G^{\rm K}(t,t^{\prime})\equiv G_{cc}(t,t^{\prime})$
$\displaystyle=-i\left<\psi_{c}(t)\bar{\psi}_{c}(t^{\prime})\right>=-i\int{\mathcal{D}}(\psi_{q},\bar{\psi}_{q},\psi_{c},\bar{\psi}_{c})\psi_{c}(t)\bar{\psi}_{c}(t^{\prime})e^{iS}=-i\left<[\hat{\psi}(t),\hat{\psi}^{\dagger}(t^{\prime})]\right>$
(S55a) $\displaystyle G^{\rm R}(t,t^{\prime})\equiv G_{cq}(t,t^{\prime})$
$\displaystyle=-i\left<\psi_{c}(t)\bar{\psi}_{q}(t^{\prime})\right>=-i\int{\mathcal{D}}(\psi_{q},\bar{\psi}_{q},\psi_{c},\bar{\psi}_{c})\psi_{c}(t)\bar{\psi}_{q}(t^{\prime})e^{iS}=-i\theta(t-t^{\prime})\left<\\{\hat{\psi}(t),\hat{\psi}^{\dagger}(t^{\prime})\\}\right>$
(S55b) $\displaystyle G^{\rm A}(t,t^{\prime})\equiv G_{qc}(t,t^{\prime})$
$\displaystyle=-i\left<\psi_{q}(t)\bar{\psi}_{c}(t^{\prime})\right>=-i\int{\mathcal{D}}(\psi_{q},\bar{\psi}_{q},\psi_{c},\bar{\psi}_{c})\psi_{q}(t)\bar{\psi}_{c}(t^{\prime})e^{iS}=i\theta(t^{\prime}-t)\left<\\{\hat{\psi}(t),\hat{\psi}^{\dagger}(t^{\prime})\\}\right>$
(S55c) $\displaystyle 0=G_{qq}(t,t^{\prime})$
$\displaystyle=-i\left<\psi_{q}(t)\bar{\psi}_{q}(t^{\prime})\right>=-i\int{\mathcal{D}}(\psi_{q},\bar{\psi}_{q},\psi_{c},\bar{\psi}_{c})\psi_{q}(t)\bar{\psi}_{q}(t^{\prime})e^{iS}.$
(S55d)

Here, $G^{\rm R},G^{\rm A}$, and $G^{\rm K}$ are called the retarded,
advanced, and Keldysh Green’s functions, respectively. Their Fourier transform
is given by,

$\displaystyle G^{\rm R}(\omega)$
$\displaystyle=\frac{1}{\omega-\varepsilon_{0}+i\delta},$ (S56a)
$\displaystyle G^{\rm A}(\omega)$
$\displaystyle=\frac{1}{\omega-\varepsilon_{0}-i\delta},$ (S56b)
$\displaystyle G^{\rm K}(\omega)$ $\displaystyle=-\pi
i(1-2f(\omega))\delta(\omega-\varepsilon_{0}).$ (S56c)

Physically, $G^{\rm R}(\omega)[=G^{\rm A*}(\omega)]$ describes the response to
an external field, while $G^{\rm K}(\omega)$ describes the correlation
function. One can check that the fluctuation-dissipation relation for fermions
[105]

$\displaystyle G^{\rm K}(\omega)=(1-2f(\omega))[G^{\rm R}(\omega)-G^{\rm
A}(\omega)]$ (S57)

is satisfied, as it should in general closed equilibrium systems.

It is worth pointing out that the action Eq. (S46) (or Eq. (S54)) vanishes
when one puts $\psi_{+}=\psi_{-}$ (or equivalently, $\psi_{q}=0$). This so-
called causality structure originates from the trace-preserving property of
the governing dynamics [45, 46, 47]. This can be seen by noting that,
replacing $\psi_{+}$ to $\psi_{-}$ implies moving the time evolution operator
$\hat{U}(t_{f},t_{0})$ that time-evolves the ket in the Keldysh partition
function (S37) to the bra-dynamics, using the cyclic property of the trace.

It is straightforward to generalize this result to many-body systems. An
example relevant to our study is a free electron gas described by the
Hamiltonian
$\hat{H}=\sum_{\bm{k}}\varepsilon_{\bm{k}}\hat{c}^{\dagger}_{\bm{k}}\hat{c}_{\bm{k}}$.
The action is given by,

$\displaystyle S=\int
dtdt^{\prime}\sum_{\bm{k}}\begin{pmatrix}\bar{\psi}_{\bm{k}}^{+}(t)&\bar{\psi}_{\bm{k}}^{-}(t)\end{pmatrix}\begin{pmatrix}G_{\bm{k}}^{++}(t,t^{\prime})&G_{\bm{k}}^{+-}(t,t^{\prime})\\\
G_{\bm{k}}^{-+}(t,t^{\prime})&G_{\bm{k}}^{--}(t,t^{\prime})\end{pmatrix}^{-1}\begin{pmatrix}\psi_{\bm{k}}^{+}(t^{\prime})\\\
\psi_{\bm{k}}^{-}(t^{\prime})\end{pmatrix}$ (S58)

where

$\displaystyle\begin{pmatrix}G_{0}^{++}(\bm{k},\omega)&G_{0}^{+-}(\bm{k},\omega)\\\
G_{0}^{-+}(\bm{k},\omega)&G_{0}^{--}(\bm{k},\omega)\end{pmatrix}=\begin{pmatrix}\big{[}\frac{1-f_{\bm{k}}}{\omega-\varepsilon_{\bm{k}}+i\delta}+\frac{f_{\bm{k}}}{\omega-\varepsilon_{\bm{k}}-i\delta}\big{]}&2\pi
i\delta(\omega-\varepsilon_{\bm{k}})f_{\bm{k}}\\\ -2\pi
i\delta(\omega-\varepsilon_{\bm{k}})(1-f_{\bm{k}})&\big{[}-\frac{f_{\bm{k}}}{\omega-\varepsilon_{\bm{k}}+i\delta}-\frac{1-f_{\bm{k}}}{\omega-\varepsilon_{\bm{k}}-i\delta}\big{]}\end{pmatrix}.$
(S59)

####  I.2.2 Open quantum system attached to a thermal bath

Keldysh formalism is useful for formulating open quantum systems. To
illustrate this, we will show an example of a single-mode fermion attached to
a thermal bath modeled by a collection of free fermions. The Hamiltonian
$\hat{H}=\hat{H}_{s}+\hat{H}_{b}+\hat{H}_{sb}$ is given by,

$\displaystyle\hat{H}_{s}=\varepsilon_{0}\hat{\psi}^{\dagger}\psi,\qquad\hat{H}_{b}=\sum_{\mu}\omega_{\mu}\hat{\phi}_{\mu}^{\dagger}\phi_{\mu},\qquad\hat{H}_{sb}=\sum_{\mu}[v_{\mu}\hat{\phi}_{\mu}^{\dagger}\hat{\psi}+{\rm
h.c.}],$ (S60)

where $\hat{\phi}_{\mu}$ and $\omega_{\mu}$ are the fermionic annihilation
operator and the energy of bath fermions at mode $\mu$, respectively. The
Keldysh partition function can be obtained for this system as,

$\displaystyle
Z=\int{\mathcal{D}}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-},\phi_{+},\bar{\phi}_{+},\phi_{-},\bar{\phi}_{-}]e^{iS_{\rm
tot}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-},\phi_{+},\bar{\phi}_{+},\phi_{-},\bar{\phi}_{-}]},$
(S61)

where the total action is composed of three parts $S_{\rm
tot}=S_{0}+S_{b}+S_{sb}$, which are given by,

$\displaystyle S_{0}$ $\displaystyle=$ $\displaystyle\int
dtdt^{\prime}\begin{pmatrix}\bar{\psi}_{+}(t)&\bar{\psi}_{-}(t)\end{pmatrix}\begin{pmatrix}G_{++}(t,t^{\prime})&G_{+-}(t,t^{\prime})\\\
G_{-+}(t,t^{\prime})&G_{--}(t,t^{\prime})\end{pmatrix}^{-1}\begin{pmatrix}\psi_{+}(t^{\prime})\\\
\psi_{-}(t^{\prime})\end{pmatrix},$ (S62) $\displaystyle S_{b}$
$\displaystyle=$ $\displaystyle\sum_{\mu}\int
dtdt^{\prime}\begin{pmatrix}\bar{\phi}_{+}^{\mu}(t)&\bar{\phi}_{-}^{\mu}(t)\end{pmatrix}\begin{pmatrix}B_{++}^{\mu}(t,t^{\prime})&B_{+-}^{\mu}(t,t^{\prime})\\\
B_{-+}^{\mu}(t,t^{\prime})&B_{--}^{\mu}(t,t^{\prime})\end{pmatrix}^{-1}\begin{pmatrix}\phi_{+}^{\mu}(t^{\prime})\\\
\phi_{-}^{\mu}(t^{\prime})\end{pmatrix},$ (S63) $\displaystyle S_{sb}$
$\displaystyle=$ $\displaystyle\sum_{\mu}\int
dt[(v_{\mu}\bar{\phi}_{+}(t)\psi_{+}(t)+{\rm
c.c.})-(v_{\mu}\bar{\phi}_{-}(t)\psi_{-}(t)+{\rm c.c.})].$ (S64)

Here, $\phi_{\pm},\bar{\phi}_{\pm}$ are the Grassmann variables of the baths,
$G_{s,s^{\prime}}(s,s^{\prime}=+,-)$ is given by Eq. (S47), and

$\displaystyle\begin{pmatrix}B_{++}^{\mu}(t,t^{\prime})&B_{+-}^{\mu}(t,t^{\prime})\\\
B_{-+}^{\mu}(t,t^{\prime})&B_{--}^{\mu}(t,t^{\prime})\end{pmatrix}$
$\displaystyle=\begin{pmatrix}i\theta(-t)f(\omega_{\mu})e^{-i\varepsilon_{0}(t-t^{\prime})}+i\theta(t)(1+f(\omega_{\mu}))e^{-i\omega_{\mu}(t-t^{\prime})}&if(\omega_{\mu})e^{-i\omega_{\mu}(t-t^{\prime})}\\\
i(1+f(\omega_{\mu}))e^{-i\omega_{\mu}(t-t^{\prime})}&i\theta(t)f(\omega_{\mu})e^{-i\omega_{\mu}(t-t^{\prime})}+i\theta(-t)(1+f(\omega_{\mu}))e^{-i\omega_{\mu}(t-t^{\prime})}\end{pmatrix}$

are the various two-point correlation functions of the baths.

Since the total action $S_{\rm tot}$ is quadratic in terms of the bath
Grassmann variables $\phi_{\pm}$, one can analytically integrate out the bath
degrees of freedom as,

$\displaystyle
Z=\int{\mathcal{D}}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]e^{iS_{\rm
eff}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]},$ (S66)

where

$\displaystyle e^{iS_{\rm
eff}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]}$ $\displaystyle=$
$\displaystyle
e^{iS_{0}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]}\int{\mathcal{D}}[\phi_{+},\bar{\phi}_{+},\phi_{-},\bar{\phi}_{-}]e^{iS_{b}[\phi_{+},\bar{\phi}_{+},\phi_{-},\bar{\phi}_{-}]+S_{sb}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-},\phi_{+},\bar{\phi}_{+},\phi_{-},\bar{\phi}_{-}]}$
(S67) $\displaystyle=$ $\displaystyle
e^{iS_{0}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]}e^{i\Delta
S[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]}.$

The contribution arising from the attachment of the bath $\Delta S$ is
computed as,

$\displaystyle\Delta S[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]$
$\displaystyle=$ $\displaystyle\sum_{\mu}|v_{\mu}|^{2}\int
dtdt^{\prime}\begin{pmatrix}\bar{\psi}_{+}(t)&\bar{\psi}_{-}(t)\end{pmatrix}\begin{pmatrix}B_{++}^{\mu}(t,t^{\prime})&B_{+-}^{\mu}(t,t^{\prime})\\\
B_{-+}^{\mu}(t,t^{\prime})&B_{--}^{\mu}(t,t^{\prime})\end{pmatrix}\begin{pmatrix}\psi_{+}(t^{\prime})\\\
\psi_{-}(t^{\prime})\end{pmatrix},$ (S68) $\displaystyle=$
$\displaystyle\sum_{\mu}|v_{\mu}|^{2}\int\frac{d\omega}{2\pi}\begin{pmatrix}\bar{\psi}_{q}(\omega)&\bar{\psi}_{c}(\omega)\end{pmatrix}\begin{pmatrix}0&\frac{1}{\omega-\omega_{\mu}-i\delta}\\\
\frac{1}{\omega-\omega_{\mu}+i\delta}&-\pi
i(1-2f(\omega))\delta(\omega-\omega_{\mu})\end{pmatrix}\begin{pmatrix}\psi_{q}(\omega)\\\
\psi_{c}(\omega)\end{pmatrix},$ $\displaystyle\simeq$
$\displaystyle\int\frac{d\omega}{2\pi}\begin{pmatrix}\bar{\psi}_{q}(\omega)&\bar{\psi}_{c}(\omega)\end{pmatrix}\begin{pmatrix}0&i\frac{\gamma(\omega)}{2}\\\
-i\frac{\gamma(\omega)}{2}&i\frac{\gamma(\omega)}{2}(1-2f(\omega))\end{pmatrix}\begin{pmatrix}\psi_{q}(\omega)\\\
\psi_{c}(\omega)\end{pmatrix}$

Here, in the second equality, we have moved to the quantum-classical
representation Eq. (S53). In the third, we have ignored the contribution from
the real part of the propagators (the “Lamb shift”) and have introduced a
dissipation rate
$\gamma(\omega)=\pi\sum_{\mu}|v_{\mu}|^{2}\delta(\omega-\omega_{\mu})$, where
we have taken the spectrum of the bath fermions to be a continuum.

The effective action $S_{\rm eff}=S_{0}+\Delta S$ is therefore summarized as,

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int
dtdt^{\prime}\begin{pmatrix}\bar{\psi}_{q}(\omega)&\bar{\psi}_{c}(\omega)\end{pmatrix}\begin{pmatrix}0&G^{\rm
A}(\omega)\\\ G^{\rm R}(\omega)&G^{\rm
K}(\omega)\end{pmatrix}^{-1}\begin{pmatrix}\psi_{q}(\omega)\\\
\psi_{c}(\omega)\end{pmatrix},$ (S69)

where the correlation functions and response functions are given by,

$\displaystyle G^{\rm R}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{\omega-\varepsilon_{0}+i\gamma(\omega)},$ (S70)
$\displaystyle G^{\rm A}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{\omega-\varepsilon_{0}-i\gamma(\omega)},$ (S71)
$\displaystyle G^{\rm K}(\omega)$ $\displaystyle=$
$\displaystyle-i\pi(1-2f(\omega))\frac{1}{\pi}\frac{\gamma(\omega)}{(\omega-\varepsilon_{0})^{2}+\gamma^{2}(\omega)}.$
(S72)

Compared with the closed system counterpart (Eq. (S56)), the delta-functions
of the spectrum in Eq. (S56) is replaced by a Lorentzian. This shows that the
coupling to the bath gives rise to dissipation in the relevant system.

Remarkably, the fluctuation-dissipation theorem (Eq. (S57)) is still satisfied
with this attachment to the thermal bath. This illustrates how the Keldysh
theory rightly describes how an attachment to a thermal bath thermalizes the
relevant system. (In fact, this can be understood as a result of a certain
symmetry of the Keldysh action that is respected in the thermal equilibrium,
see Ref. [106].) We remark that this could not have been captured in the GKSL
formalism discussed in Sec. I.1 [106], where the Markovian approximation
ignores the dynamical effects arising from the Fermi distribution function.

####  I.2.3 Keldysh technique applied to the GKSL master equation

So far, we have considered the case where the relevant system is attached to a
single bath, reproducing the results known from the equilibrium statistical
mechanics. The real advantage of using the Keldysh approach is that it can
treat nonequilibrium systems arising from attaching to multiple baths.

In our study, we will be concerned with a system (localized electrons)
attached to two baths. One of the baths is Markovian (the higher-energy state)
and the other is non-Markovian (conduction electrons). As we have reviewed in
Sec. I.1, once we integrated out the Markovian bath, the reduced density
operator $\hat{\rho}$ is governed by the GKSL equation (S4), reproduced below
for convenience:

$\displaystyle\partial_{t}\hat{\rho}$ $\displaystyle=$
$\displaystyle{\mathcal{L}}\hat{\rho}=-i[\hat{H},\hat{\rho}]+\sum_{\ell}\kappa_{\ell}{\mathcal{D}}[\hat{L}_{\ell}]\hat{\rho},$
(S73)

with
${\mathcal{D}}[\hat{L}]\hat{\rho}=\hat{L}\hat{\rho}\hat{L}^{\dagger}-\frac{1}{2}\\{\hat{L}^{\dagger}\hat{L},\hat{\rho}\\}$.
Therefore, we will have in mind the situation where the Hamiltonian $\hat{H}$
is composed of the relevant system and the non-Markovian bath (but the
formalism below is general).

Following Refs. [46, 47], we formulate a Keldysh path integral formalism for
the GKSL equation to treat this type of open system. As done in the previous
subsections, we define the Keldysh action

$\displaystyle Z$ $\displaystyle=$ $\displaystyle{\rm
tr}[\hat{\rho}(t_{f})]={\rm
tr}\big{[}e^{{\mathcal{L}}(t_{f}-t_{0})}\hat{\rho}(t_{0})\big{]}={\rm
tr}\Big{[}\lim_{N\rightarrow\infty}\prod_{n=0}^{N}\hat{\rho}_{n}\Big{]}$ (S74)

where we have discretized the GKSL equation (S73) as

$\displaystyle\hat{\rho}_{n+1}=(\hat{1}+\delta_{t}{\mathcal{L}})[\hat{\rho}_{n}]+O(\delta_{t}^{2}).$
(S75)

with $\delta_{t}=(t_{f}-t_{0})/N$ and defined
$\hat{\rho}_{0}\equiv\hat{\rho}(t_{0})$.

Let us restrict ourselves, for simplicity, to a single-mode fermion system. As
before, we expand the Keldysh partition function into coherent states. This
can be computed once
$\langle\psi_{+,n+1}\rvert\hat{\rho}_{n+1}\lvert-\psi_{-,n+1}\rangle=\langle\psi_{+,n+1}\rvert(\hat{1}+\delta_{t}{\mathcal{L}})[\hat{\rho}_{n}]\lvert-\psi_{-,n+1}\rangle$
is calculated. Using

$\displaystyle\langle\psi_{+,n+1}\rvert{\mathcal{L}}\hat{\rho}_{n}\lvert-\psi_{-,n+1}\rangle$
(S76) $\displaystyle=\int
d\bar{\psi}_{+,n}d\psi_{+,n}d\bar{\psi}_{-,n}d\psi_{-,n}e^{-\bar{\psi}_{+,n}\psi_{+,n}-\bar{\psi}_{-,n}\psi_{-,n}}\langle\psi_{+,n}\rvert\hat{\rho}_{n}\lvert-\psi_{-,n+1}\rangle\langle\psi_{+,n+1}\rvert{\mathcal{L}}[\lvert\psi_{+,n}\rangle\langle-\psi_{-,n}\rvert]\lvert-\psi_{-,n+1}\rangle$
$\displaystyle=\int
d\bar{\psi}_{+,n}d\psi_{+,n}d\bar{\psi}_{-,n}d\psi_{-,n}e^{-\bar{\psi}_{+,n}\psi_{+,n}-\bar{\psi}_{-,n}\psi_{-,n}}\langle\psi_{+,n}\rvert\hat{\rho}_{n}\lvert-\psi_{-,n}\rangle$
$\displaystyle\times$
$\displaystyle\Bigg{[}(-i)\Big{[}\langle\psi_{+,n+1}\rvert\hat{H}\lvert\psi_{+,n}\rangle\langle-\psi_{-,n}|-\psi_{-,n+1}\rangle-\langle\psi_{+,n+1}|\psi_{+,n}\rangle\langle-\psi_{-,n}\rvert\hat{H}\lvert-\psi_{-,n+1}\rangle\Big{]}$
$\displaystyle+\sum_{\ell}\kappa_{\ell}\Big{[}\langle\psi_{+,n+1}\rvert\hat{L}_{\ell}\lvert\psi_{+,n}\rangle\langle-\psi_{-,n}\rvert\hat{L}_{\ell}^{\dagger}\lvert-\psi_{-,n+1}\rangle$
$\displaystyle\ \ \ \ \ \ \ \ \ \
-\frac{1}{2}\Big{[}\langle\psi_{+,n+1}\rvert\hat{L}_{\ell}^{\dagger}\hat{L}_{\ell}\lvert\psi_{+,n}\rangle\langle-\psi_{-,n}|-\psi_{-,n+1}\rangle+\langle\psi_{+,n+1}|\psi_{+,n}\rangle\langle-\psi_{-,n}\rvert\hat{L}_{\ell}^{\dagger}\hat{L}_{\ell}\lvert-\psi_{-,n+1}\rangle\Big{]}\Bigg{]}$
$\displaystyle=\int
d\bar{\psi}_{+,n}d\psi_{+,n}d\bar{\psi}_{-,n}d\psi_{-,n}e^{-(\bar{\psi}_{+,n}-\bar{\psi}_{+,n+1})\psi_{+,n}-\bar{\psi}_{-,n}(\psi_{-,n}-\psi_{-,n+1})}\langle\psi_{+,n}\rvert\hat{\rho}_{n}\lvert-\psi_{-,n}\rangle$
$\displaystyle\times$
$\displaystyle\Bigg{[}(-i)\Big{[}H_{+}(\psi_{+,n},\bar{\psi}_{+,n})-H_{-}(\psi_{-,n},\bar{\psi}_{-,n})\Big{]}+\sum_{\ell}\kappa_{\ell}\Big{[}L_{+,n,\ell}(-\bar{L}_{-,n,\ell})-\frac{1}{2}\big{(}\bar{L}_{+,n,\ell}L_{+,n,\ell}+\bar{L}_{-,n,\ell}L_{-,n,\ell}\big{)}\Big{]}\Bigg{]}$

where

$\displaystyle
H_{+}(\psi_{+,n},\bar{\psi}_{+,n})=\langle\psi_{+,n}\rvert\hat{H}\lvert\psi_{+,n}\rangle,\quad
H_{-}(\psi_{-,n},\bar{\psi}_{-,n})=\langle-\psi_{-,n}\rvert\hat{H}\lvert-\psi_{-,n}\rangle,$
(S77) $\displaystyle
L_{+,n,\ell}=\langle\psi_{+,n}\rvert\hat{L}_{\ell}\lvert\psi_{+,n}\rangle,\qquad
L_{-,n,\ell}=\langle-\psi_{-,n}\rvert\hat{L}_{\ell}\lvert-\psi_{-,n}\rangle,$
(S78)

This yields,

$\displaystyle\langle\psi_{+,n+1}\rvert\hat{\rho}_{n+1}\lvert-\psi_{-,n+1}\rangle$
$\displaystyle=$ $\displaystyle\int
d\bar{\psi}_{+,n}d\psi_{+,n}d\bar{\psi}_{-,n}d\psi_{-,n}e^{-(\bar{\psi}_{+,n}-\bar{\psi}_{+,n+1})\psi_{+,n}-\bar{\psi}_{-,n}(\psi_{-,n}-\psi_{-,n+1})}e^{\delta_{t}{\mathcal{L}}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]}$

with

$\displaystyle{\mathcal{L}}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]$
$\displaystyle=$
$\displaystyle(-i)\Big{[}H_{+}(\psi_{+,n},\bar{\psi}_{+,n})-H_{-}(\psi_{-,n},\bar{\psi}_{-,n})\Big{]}$
(S80) $\displaystyle+$
$\displaystyle\sum_{\ell}\kappa_{\ell}\Big{[}L_{+,n,\ell}(-\bar{L}_{-,n,\ell})-\frac{1}{2}\big{(}\bar{L}_{+,n,\ell}L_{+,n,\ell}+\bar{L}_{-,n,\ell}L_{-,n,\ell}\big{)}\Big{]}$

Here, we have assumed that the Hamiltonian $\hat{H}$, the dissipators
$\hat{L}_{\ell}$, and $\hat{L}_{\ell}^{\dagger}\hat{L}_{\ell}$ are normal
ordered. We note that there are subtleties in this normal ordering procedure,
where the causality structure may be broken [107]. This can be mitigated by
introducing an infinitesimal time difference that enforces the ordering of the
operators, see e.g. Appendix A in Ref. [47] for details. Taking the continuous
limit, the Keldysh partition function can be computed as

$\displaystyle
Z=\int{\mathcal{D}}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]e^{iS[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]},$
(S81)

with the Keldysh action $S$ given by,

$\displaystyle S[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]$
$\displaystyle=$ $\displaystyle\int
dt\Big{[}\bar{\psi}_{+}(t)(i\partial_{t}-H_{+}(t))\psi_{+}(t)-\bar{\psi}_{-}(t)(i\partial_{t}-H_{-}(t))\psi_{+}(t)$
(S82) $\displaystyle+$
$\displaystyle\sum_{\ell}\kappa_{\ell}\big{[}L_{+,\ell}(t)\bar{L}_{-,\ell}(t)-\frac{1}{2}(\bar{L}_{+,\ell}(t_{+\delta})L_{+,\ell}(t_{-\delta})+\bar{L}_{-,\ell}(t_{-\delta})L_{-,\ell}(t_{+\delta}))\big{]}\Big{]}$

Here, we have introduced $t_{\pm\delta}=t\pm 0^{+}$ to ensure the causality
structure of the action, i.e.,
$S[\psi_{+}=\psi_{-},\bar{\psi}_{+}=\bar{\psi}_{-}]=0$. Equation (S82) that
gives the relation between the GKSL equation (S73) and the Keldysh action $S$
is the central relation used in Sec. IV.

##  II The model

Equipped with the theoretical tools introduced in Sec. I, we formulate below
our dissipation-engineering scheme to control local spin-spin interactions in
magnetic metals. This section summarizes our model.

As illustrated schematically in Fig. 2(a) in the main text, this system is
composed of localized electrons immersed in a Fermi sea of conduction
electrons. The localized electrons are responsible for magnetic properties,
while the conduction electrons are responsible for metallic properties. To
manipulate the effective spin-spin interactions between the localized sites
(Fig. 2(b) in the main text), we inject a laser light that gives rise to
dissipation to the double-occupied state at site $a$ at the rate $\kappa_{a}$.
By tuning the frequency of the laser to be resonant with a higher-level state
(say, upper-level $f$-orbital state) that has a linewidth $\Gamma_{f,a}$,
(Fig. 2(c) in the main text), we turn on tunneling from the double-occupied
state to the higher-level state ($a$ labels the localized electron sites).
This process can be regarded as a Markov process, as long as the dissipation
rate of the higher-level state $\Gamma_{f,a}$ (which is the timescale of the
memory of the bath) is larger than its supply rate. A detailed discussion on
this point is provided in Methods. As explained in Sec. I.1, the time
evolution of reduced system density operator $\hat{\rho}$ of such Markov
system is generally governed by the GKSL master equation. In our system, it
reads

$\displaystyle\partial_{t}\hat{\rho}$ $\displaystyle=$
$\displaystyle{\mathcal{L}}\hat{\rho}$ (S83) $\displaystyle=$
$\displaystyle-i[\hat{H},\hat{\rho}]+\sum_{a,\sigma}\kappa_{a}{\mathcal{D}}[\hat{d}_{\sigma,a}\hat{P}_{\uparrow\downarrow}^{a}]\hat{\rho}+\sum_{a,\sigma}\delta_{a}{\mathcal{D}}[\hat{d}_{\sigma,a}^{\dagger}\hat{P}_{\varnothing}^{a}]\hat{\rho}.$

The first term describes the coherent dynamics of the magnetic metals and the
second term is the light-induced dissipation. We have added the third term
phenomenologically to ensure the steady state of localized electrons is singly
occupied at each site. (We have omitted the third term in the main text, as
they only play a role in ensuring the state is in the singly occupied state.)

The magnetic metal is modeled by the Anderson impurity model
$\hat{H}=\hat{H}_{0}+\hat{H}_{cd}$ consisting of two parts. The unperturbed
part $\hat{H}_{0}=\hat{H}_{c0}+\hat{H}_{d0}$ describes the conduction
electrons and the localized electrons, given respectively by,

$\displaystyle\hat{H}_{c0}$ $\displaystyle=$
$\displaystyle\sum_{\bm{k},\sigma=\uparrow,\downarrow}\varepsilon_{\bm{k}}\hat{c}^{\dagger}_{\bm{k}\sigma}\hat{c}_{\bm{k}\sigma},$
(S84) $\displaystyle\hat{H}_{d0}$ $\displaystyle=$
$\displaystyle\sum_{a}\Big{[}\sum_{\sigma=\uparrow,\downarrow}\varepsilon_{d,a}\hat{d}^{\dagger}_{\sigma,a}\hat{d}_{\sigma,a}+U_{a}\hat{d}^{\dagger}_{\uparrow,a}\hat{d}_{\uparrow,a}\hat{d}^{\dagger}_{\downarrow,a}\hat{d}_{\downarrow,a}\Big{]}.$
(S85)

Here, the conduction electrons are modeled as free electrons, where
$\hat{c}_{\bm{k},\sigma}$ is a fermionic annihilation operator of conduction
electrons with momentum $\bm{k}$ and spin $\sigma=\uparrow,\downarrow$ and
$\varepsilon_{\bm{k}}$ is the kinetic energy. We assume that the conduction
electrons’ system size is large enough to ensure they are always in thermal
equilibrium. Its distribution obeys the Fermi distribution function
$f(\varepsilon_{\bm{k}})=[e^{(\varepsilon_{\bm{k}}-\varepsilon_{\rm
F})/(k_{\rm B}T)}+1]^{-1}$ that is characterized by temperature $T$ and the
Fermi energy $\varepsilon_{\rm F}$. We will focus on low-temperature regime
$k_{\rm B}T\ll\varepsilon_{\rm F}$. The localized electrons are modeled as a
two-level system with energy $\varepsilon_{d,a}$ and an on-site Coulomb
repulsion $U_{a}>0$ at site $a$. Here, $\hat{d}_{a,\sigma}$ is a fermionic
annihilation operator of localized electrons at site $a$. These two types of
electrons couple through the second part of the Hamiltonian $\hat{H}$, namely
the c-d mixing $\hat{H}_{cd}$, where the electrons in the conduction band
tunnel at rate $v_{a}$ to the localized orbital and vice versa, described by

$\displaystyle\hat{H}_{cd}$ $\displaystyle=$
$\displaystyle\sum_{a,\sigma}\big{[}v_{a}\hat{d}_{\sigma,a}^{\dagger}\hat{c}_{\bm{R}_{a}\sigma}+{\rm
h.c.}\big{]}=\sum_{a}\sum_{\bm{k},\sigma}\big{[}v_{a}e^{i\bm{k}\cdot\bm{R}_{a}}\hat{d}_{\sigma,a}^{\dagger}\hat{c}_{\bm{k}\sigma}+{\rm
h.c.}\big{]},$ (S86)

Here,
$\hat{c}_{\bm{r},\sigma}=\sum_{\bm{k}}e^{i\bm{k}\cdot\bm{r}}\hat{c}_{\bm{k},\sigma}$
is the Fourier transform of $\hat{c}_{\bm{k},\sigma}$ and $\bm{R}_{a}$ is the
position of the local electron at site $a$. In what follows, we restrict our
interest to systems where the Coulomb repulsion is strong enough ($U_{a}\gg
v_{a}$) to treat $\hat{H}_{cd}$ as a perturbation.

The second term in the GKSL master equation (S83) describes the light-induced
loss of an electron in the double-occupied state. The localized electron at
site $a$ can take four states: it can be vacant
$\lvert\varnothing\rangle_{a}$, singly occupied with
$\uparrow(\downarrow)$-spin
$\lvert\uparrow(\downarrow)\rangle_{a}=\hat{d}_{a,\uparrow(\downarrow)}^{\dagger}\lvert\varnothing\rangle$,
or double-occupied occupied
$\lvert\uparrow\downarrow\rangle_{a}=\hat{d}_{a,\uparrow}^{\dagger}\hat{d}_{a,\downarrow}^{\dagger}\lvert\varnothing\rangle_{a}$.
$\hat{P}_{\uparrow\downarrow}^{a}$ at site $a$ projects the state to the
double-occupied state $\lvert\uparrow\downarrow\rangle_{a}$, enforcing the
loss to activate only when site $a$ is double-occupied. It will become clear
in the following that this selective decay of localized electrons effectively
modifies the interactions between the localized spins.

The third term ($\delta_{a}\rightarrow 0^{+}$) is phenomenologically added to
drive the system towards a singly occupied state ($\lvert\uparrow\rangle_{a}$
or $\lvert\downarrow\rangle_{a}$), where $\hat{P}_{\varnothing}^{a}$ is the
projection operator to the vacant state $\lvert\varnothing\rangle$ at site
$a$.

##  III Effective localized spin dynamics immersed in conduction electrons
with light-induced correlated dissipation

To proceed, we resort to the projection method introduced in Sec. I.1.2 to
perturbatively derive the effective low-energy dynamics given by Eq. (4) in
the main text. We split the Lindbladian
${\mathcal{L}}={\mathcal{L}}_{0}+{\mathcal{L}}_{1}$ in Eq. (S83) to non-
perturbative (${\mathcal{L}}_{0}$) and perturbative part (${\mathcal{L}}_{1}$)
as follows:

$\displaystyle{\mathcal{L}}_{0}\hat{\rho}$ $\displaystyle=$
$\displaystyle-i[\hat{H}_{c0}+\hat{H}_{d0},\hat{\rho}]+\sum_{a,\sigma}\kappa_{a}{\mathcal{D}}[\hat{d}_{\sigma,a}\hat{P}_{\uparrow\downarrow}^{a}]\hat{\rho}+\sum_{a,\sigma}\delta_{a}{\mathcal{D}}[\hat{d}_{\sigma,a}^{\dagger}\hat{P}_{\varnothing}^{a}]\hat{\rho},$
(S87) $\displaystyle{\mathcal{L}}_{1}\rho$ $\displaystyle=$
$\displaystyle-i[\hat{H}_{cd},\hat{\rho}].$ (S88)

The perturbative part ${\mathcal{L}}_{1}$, i.e. the c-d mixing, gives a
contribution at most $O(v_{a}/U_{a})$, which is assumed to be small. In this
section, we will use the relation Eq. (S32) to derive the low-energy effective
Lindbladian of our system (Eq. (S83)) within the second-order perturbation in
${\mathcal{L}}_{1}$, which we reproduce below for convenience:

$\displaystyle({\mathcal{L}}_{\rm eff})_{n_{l},n_{r}}$ $\displaystyle\equiv$
$\displaystyle(\hat{l}_{n_{l}}^{(0)},{\mathcal{L}_{\rm
eff}\hat{r}_{n_{r}}^{(0)}})={\rm
tr}[\hat{l}_{n_{l}}^{(0)\dagger}{\mathcal{L}_{\rm eff}\hat{r}_{n_{r}}^{(0)}}]$
(S89) $\displaystyle=$ $\displaystyle{\rm
tr}[\hat{l}_{n_{l}}^{(0)\dagger}{\mathcal{L}_{0}\hat{r}_{n_{r}}^{(0)}}]+{\rm
tr}[\hat{l}_{n_{l}}^{(0)\dagger}{\mathcal{L}_{1}\hat{r}_{n_{r}}^{(0)}}]-\sum_{m\in\mathfrak{f}}\frac{{\rm
tr}[\hat{l}_{n_{l}}^{(0)\dagger}{\mathcal{L}}_{1}\hat{r}_{m}^{(0)}]{\rm
tr}[\hat{l}_{m}^{(0)\dagger}{\mathcal{L}}_{1}\hat{r}_{n_{r}}^{(0)}]}{\lambda_{m}^{(0)}}+O(({\mathcal{L}}_{1})^{3}).$

Here, $\hat{r}_{n_{r}}^{(0)}$ and $\hat{l}_{n_{l}}^{(0)}$ are the right and
left eigenstates, respectively, that form a basis of the slow variable space
we define below. $(\hat{A},\hat{B})={\rm tr}[\hat{A}^{\dagger}\hat{B}]$ is the
Hilbert-Schmidt inner product (see Eq. (S8)). In the third term, the sum is
taken over the fast degrees of freedom, also defined below.

###  III.1 Characterization of non-perturbative Lindbladian
${\mathcal{L}}_{0}$

As is clear from the expression of Eq. (S89), the first step to derive the
effective low-energy Lindbladian is to characterize the (left) right
eigenstates $\hat{r}_{n}^{(0)}(\hat{l}_{n}^{(0)})$ and eigenvalues of the
unperturbed Lindbladian ${\mathcal{L}}_{0}$ (Eq. (S87)). This can be done
exactly since the nonperturbative Lindbladian
${\mathcal{L}}_{0}={\mathcal{L}}_{c0}+\sum_{a}{\mathcal{L}}_{d0,a}$ is
composed of free conduction electrons and localized electrons that are
completely decoupled, where

$\displaystyle{\mathcal{L}}_{c0}\hat{\rho}$ $\displaystyle=$
$\displaystyle-i\big{[}\sum_{\bm{k},\sigma}\varepsilon_{\bm{k}}\hat{c}_{\bm{k},\sigma}^{\dagger}\hat{c}_{\bm{k},\sigma},\hat{\rho}\big{]},$
(S90) $\displaystyle{\mathcal{L}}_{d0,a}\hat{\rho}$ $\displaystyle=$
$\displaystyle-i\big{[}\big{(}\sum_{\sigma}\varepsilon_{d,a}\hat{d}^{\dagger}_{\sigma,a}\hat{d}_{\sigma,a}+U_{a}\hat{d}^{\dagger}_{\uparrow,a}\hat{d}_{\uparrow,a}\hat{d}^{\dagger}_{\downarrow,a}\hat{d}_{\downarrow,a}\big{)},\hat{\rho}]+\sum_{\sigma}\kappa_{a}{\mathcal{D}}[\hat{d}_{\sigma,a}\hat{P}_{\uparrow\downarrow}^{a}]\hat{\rho}+\sum_{\sigma}\delta_{a}{\mathcal{D}}[\hat{d}_{\sigma,a}^{\dagger}\hat{P}_{\varnothing}^{a}]\hat{\rho}.$
(S91)

The eigenvectors can be expressed as a direct product, i.e.
$\hat{r}_{n}^{(0)}=\hat{r}_{n}^{c(0)}\otimes\prod_{a}\hat{r}_{a,n}^{d(0)}$ and
$\hat{l}_{n}^{(0)}=\hat{l}_{n}^{c(0)}\otimes\prod_{a}\hat{l}_{a,n}^{d(0)}$.
Here, $\hat{r}_{n}^{c(0)}$ and $\hat{l}_{n}^{c(0)}$ are the right and left
eigenstates, respectively, of the (non-perturbative) conduction electron
Lindbladian ${\mathcal{L}}_{c0}$ that satisfies,

$\displaystyle{\mathcal{L}}_{c0}\hat{r}_{n}^{c(0)}=\lambda_{n}^{c(0)}\hat{r}_{n}^{c(0)},\qquad{\mathcal{L}}_{c0}^{\dagger}\hat{l}_{n}^{c(0)}=\lambda_{n}^{c(0)*}\hat{l}_{n}^{c(0)},$
(S92)

with an eigenvalue $\lambda_{n}^{c(0)}$. Similarly, $\hat{r}_{a,n}^{d(0)}$ and
$\hat{l}_{a,n}^{d(0)}$ are the right and left eigenstates, respectively, with
an eigenvalue $\lambda_{a,n}^{d(0)}$ for the localized electron part
${\mathcal{L}}_{d0,a}$

$\displaystyle{\mathcal{L}}_{d0,a}\hat{r}_{a,n}^{d(0)}=\lambda_{a,n}^{d(0)}\hat{r}_{a,n}^{d(0)},\qquad{\mathcal{L}}_{d0,a}^{\dagger}\hat{l}_{a,n}^{d(0)}=\lambda_{a,n}^{d(0)*}\hat{l}_{a,n}^{d(0)}.$
(S93)

The total eigenvalue $\lambda_{n}^{(0)}$ is given by their sum
$\lambda_{n}^{(0)}=\lambda_{n}^{c(0)}+\sum_{a}\lambda_{a,n}^{d(0)}$.

####  III.1.1 Characterization of ${\mathcal{L}_{d0,a}}$

Let us first characterize ${\mathcal{L}_{d0,a}}$. It is helpful to point out
that the superoperator
${\mathcal{N}_{a,\sigma}^{d}}\hat{\rho}=[\hat{d}_{a,\sigma}^{\dagger}\hat{d}_{a,\sigma},\hat{\rho}]$
commutes with ${\mathcal{L}_{d0,a}}$, implying that the two superoperators
share a set of eigenvectors. The eigenvector of the superoperator
${\mathcal{N}_{a,\sigma}^{d}}$ is the outer product of Fock states with
corresponding eigenvalue $m_{\uparrow}(m_{\downarrow})=0,\pm 1$. Here, the
eigenvalue $m_{\uparrow}(m_{\downarrow})=0,\pm 1$ is given by the number of
$\sigma=\uparrow(\downarrow)$-spin electrons in ket minus those in bra.
Therefore, the right and left eigenvectors of ${\mathcal{L}_{d0,a}}$ can be
cast into a block-diagonal form
${\mathcal{L}_{d0,a}}=\otimes_{m_{\uparrow},m_{\downarrow}}{\mathcal{L}_{d0,a,(m_{\uparrow},m_{\downarrow})}}$,
where $m_{\sigma}$ is the eigenvalue of ${\mathcal{N}}_{a,\sigma}^{d}$. The
simultaneous right and left eigenstates of ${\mathcal{L}_{d0,a}}$ and
${\mathcal{N}}_{a,\sigma}^{d}$ are given respectively by

$\displaystyle{\mathcal{L}}_{d0,a,(m_{\uparrow},m_{\downarrow})}\hat{r}_{a,\ell,(m_{\uparrow},m_{\downarrow})}^{d(0)}=\lambda_{a,\ell,(m_{\uparrow},m_{\downarrow})}^{d(0)}\hat{r}_{a,\ell,(m_{\uparrow},m_{\downarrow})}^{d(0)},\qquad{\mathcal{L}}_{d0,a,(m_{\uparrow},m_{\downarrow})}^{\dagger}\hat{l}_{a,\ell,(m_{\uparrow},m_{\downarrow})}^{d(0)}=\lambda_{a,\ell,(m_{\uparrow},m_{\downarrow})}^{d(0)*}\hat{l}_{a,\ell,(m_{\uparrow},m_{\downarrow})}^{d(0)}.$
(S94)

For example, $(m_{\uparrow},m_{\downarrow})=(0,0)$ sector has the right and
left eigenmode of the form

$\displaystyle\hat{r}_{a,\ell,m_{\uparrow}=0,m_{\downarrow}=0}^{d(0)}$
$\displaystyle=$ $\displaystyle
r_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\varnothing}^{(0)}\lvert\varnothing\rangle\langle\varnothing\rvert+r_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\uparrow}^{(0)}\lvert\uparrow\rangle\langle\uparrow\rvert$
(S95) $\displaystyle+$ $\displaystyle
r_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\downarrow}^{(0)}\lvert\downarrow\rangle\langle\downarrow\rvert+r_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\uparrow\downarrow}^{(0)}\lvert\uparrow\downarrow\rangle\langle\uparrow\downarrow\rvert,$
$\displaystyle\hat{l}_{a,\ell,m_{\uparrow}=0,m_{\downarrow}=0}^{d(0)}$
$\displaystyle=$ $\displaystyle
l_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\varnothing}^{(0)}\lvert\varnothing\rangle\langle\varnothing\rvert+l_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\uparrow}^{(0)}\lvert\uparrow\rangle\langle\uparrow\rvert$
(S96) $\displaystyle+$ $\displaystyle
l_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\downarrow}^{(0)}\lvert\downarrow\rangle\langle\downarrow\rvert+l_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\uparrow\downarrow}^{(0)}\lvert\uparrow\downarrow\rangle\langle\uparrow\downarrow\rvert.$

Similarly, the right and left eigenvectors have the form

$\displaystyle\hat{r}_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0)}^{d(0)}$
$\displaystyle=$ $\displaystyle
r_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0),\varnothing}^{(0)}\lvert\uparrow\rangle\langle\varnothing\rvert+r_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0),\downarrow}^{(0)}\lvert\uparrow\downarrow\rangle\langle\downarrow\rvert,$
(S97) $\displaystyle\hat{l}_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0)}^{d(0)}$
$\displaystyle=$ $\displaystyle
l_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0),\varnothing}^{(0)}\lvert\uparrow\rangle\langle\varnothing\rvert+l_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0),\downarrow}^{(0)}\lvert\uparrow\downarrow\rangle\langle\downarrow\rvert,$
(S98)

for $(m_{\uparrow},m_{\downarrow})=(1,0)$,

$\displaystyle\hat{r}_{a,\ell,(m_{\uparrow}=-1,m_{\downarrow}=0)}^{d(0)}$
$\displaystyle=$ $\displaystyle
r_{a,\ell,(m_{\uparrow}=-1,m_{\downarrow}=0),\uparrow}^{(0)}\lvert\varnothing\rangle\langle\uparrow\rvert+r_{a,\ell,(m_{\uparrow}=-1,m_{\downarrow}=0),\uparrow\downarrow}^{(0)}\lvert\downarrow\rangle\langle\uparrow\downarrow\rvert,$
(S99)
$\displaystyle\hat{l}_{a,\ell,(m_{\uparrow}=-1,m_{\downarrow}=0)}^{d(0)}$
$\displaystyle=$ $\displaystyle
l_{a,\ell,(m_{\uparrow}=-1,m_{\downarrow}=0),\uparrow}^{(0)}\lvert\varnothing\rangle\langle\uparrow\rvert+l_{a,\ell,(m_{\uparrow}=-1,m_{\downarrow}=0),\uparrow\downarrow}^{(0)}\lvert\downarrow\rangle\langle\uparrow\downarrow\rvert,$
(S100)

for $(m_{\uparrow},m_{\downarrow})=(-1,0)$, and so on.

Let us start with the analysis of the eigenvectors of
$(m_{\uparrow},m_{\downarrow})=(0,0)$-sector, which is often called the
$T_{1}$-modes that physically correspond to the population dynamics of the
states. In this sector, one finds by plugging the form Eqs. (S95) and (S96)
into Eqs. (S94),

$\displaystyle\bm{M}_{a}\bm{r}_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0)}^{(0)}$
$\displaystyle=$
$\displaystyle\lambda_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0)}^{(0)}\bm{r}_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0)}^{(0)},\qquad\bm{r}_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0)}^{(0)}=\begin{pmatrix}r_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\varnothing}^{(0)}\\\
r_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\uparrow}^{(0)}\\\
r_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\downarrow}^{(0)}\\\
r_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\uparrow\downarrow}^{(0)}\end{pmatrix},$
(S101)
$\displaystyle(\bm{l}_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0)}^{(0)})^{\mathsf{T}}\bm{M}_{a}$
$\displaystyle=$
$\displaystyle\lambda_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0)}^{(0)}(\bm{l}_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0)}^{(0)})^{\mathsf{T}},\qquad\bm{l}_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0)}^{(0)}=\begin{pmatrix}l_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\varnothing}^{(0)}\\\
l_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\uparrow}^{(0)}\\\
l_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\downarrow}^{(0)}\\\
l_{a,\ell,(m_{\uparrow}=0,m_{\downarrow}=0),\uparrow\downarrow}^{(0)}\end{pmatrix},$
(S102)

with

$\displaystyle\bm{M}_{a}=\begin{pmatrix}-2\delta_{a}&0&0&0\\\
\delta_{a}&0&0&\kappa_{a}\\\ \delta_{a}&0&0&\kappa_{a}\\\
0&0&0&-2\kappa_{a}\end{pmatrix}.$ (S103)

Note that the Hamiltonian $\hat{H}$ gives no contribution to this Fock-space
diagonal subspace, where the dynamical matrix $\bm{M}_{a}$ contains no
contribution from $\hat{H}$. This is due to the property that $\hat{H}$ is
diagonal in the Fock space.

The first four rows of Table 1 report the eigenvalues and the right and left
eigenvectors in the $(m_{\uparrow},m_{\downarrow})=(0,0)$ sector, computed by
diagonalizing $\bm{M}_{a}$. Note that the first two eigenvectors are
degenerate with zero eigenvalues. This is because both spin configurations of
singly occupied state ($\lvert\uparrow\rangle$ and $\lvert\downarrow\rangle$)
do not experience any dissipation in ${\mathcal{L}}_{d0,a}$. We also note that
this zero-eigenvalue degenerate space of the left eigenvectors contains the
identity operator $\hat{I}$, as is expected from the general result (Eq.
(S23)).

We now move on to the analysis of other sectors
$(m_{\uparrow},m_{\downarrow})\neq(0,0)$ (often called the $T_{2}$-modes).
Noting that the jump term
$\hat{d}_{\sigma,a}\hat{P}_{\uparrow\downarrow}\hat{r}_{a,\ell,(m_{\uparrow},m_{\downarrow})}\hat{P}_{\uparrow\downarrow}\hat{d}_{\sigma,a}^{\dagger}=0~{}((m_{\uparrow},m_{\downarrow})\neq(0,0))$
does not contribute in this sector because of the presence of the projection
operator $\hat{P}_{\uparrow\downarrow}$, the eigenvalues and eigenstates are
given as a simple beamsplitter operator, see Table 1. This can be readily
checked. For example, for the $(m_{\uparrow},m_{\downarrow})=(1,0)$ sector,
${\mathcal{L}}_{d0}\hat{r}_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0)}$ can be
compute as,

$\displaystyle{\mathcal{L}}_{d0}(r_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0),\varnothing}^{(0)}\lvert\uparrow\rangle\langle\varnothing\rvert+r_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0),\downarrow}^{(0)}\lvert\uparrow\downarrow\rangle\langle\downarrow\rvert)$
$\displaystyle=\Big{(}-i\varepsilon_{d,a}+\frac{\delta_{a}}{2}\Big{)}r_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0),\varnothing}^{(0)}\lvert\uparrow\rangle\langle\varnothing\rvert+\Big{(}-i(\varepsilon_{d,a}+U_{a})-\frac{\kappa_{a}}{2}\Big{)}r_{a,\ell,(m_{\uparrow}=1,m_{\downarrow}=0),\downarrow}^{(0)}\lvert\uparrow\downarrow\rangle\langle\downarrow\rvert$
(S104)

which is diagonal, giving the eigenvalues and eigenstates listed in the
seventh and ninth row of Table 1.

As mentioned earlier, while the $T_{1}$-modes (the
$(m_{\uparrow},m_{\downarrow})=(0,0)$ sector) only involve dissipation, the
$T_{2}$-modes (the $(m_{\uparrow},m_{\downarrow})\neq(0,0)$ sector) involves
coherence between the two different states. The oscillation frequencies of the
$T_{2}$-modes are the difference between the energies of two states, which are
of order $O(\varepsilon_{d,a},U_{a})$. We regard the first six (the last six)
eigenvalues on Table 1 as the slow (fast) modes. We will perturbatively
project out the fast modes to obtain the effective dynamics including the
effect of the c-d mixing ${\mathcal{L}}_{1}$.

Table 1: Eigenvalues and eigenstates of ${\mathcal{L}}_{d0,a}$ | Eigenvalue $\lambda_{a,n}^{d(0)}$ | Right eigenstate $\hat{r}_{a,n}^{d(0)}$ | Left eigenstate $\hat{l}_{a,n}^{d(0)}$
---|---|---|---
slow | $0$ | $\frac{1}{2}[\lvert\uparrow\rangle\langle\uparrow\rvert+\lvert\downarrow\rangle\langle\downarrow\rvert]$ | $\hat{I}$
$0$ | $\frac{1}{2}[\lvert\uparrow\rangle\langle\uparrow\rvert-\lvert\downarrow\rangle\langle\downarrow\rvert]$ | $\frac{1}{2}[\lvert\uparrow\rangle\langle\uparrow\rvert-\lvert\downarrow\rangle\langle\downarrow\rvert]$
$-2\delta_{a}$ | $\lvert\uparrow\rangle\langle\uparrow\rvert+\lvert\downarrow\rangle\langle\downarrow\rvert-2\lvert\varnothing\rangle\langle\varnothing\rvert$ | $-\frac{1}{2}\lvert\varnothing\rangle\langle\varnothing\rvert$
$-2\kappa_{a}$ | $\lvert\uparrow\downarrow\rangle\langle\uparrow\downarrow\rvert-\frac{1}{2}\lvert\uparrow\rangle\langle\uparrow\rvert-\frac{1}{2}\lvert\downarrow\rangle\langle\downarrow\rvert$ | $\lvert\uparrow\downarrow\rangle\langle\uparrow\downarrow\rvert$
$0$ | $\lvert\uparrow\rangle\langle\downarrow\rvert$ | $\lvert\uparrow\rangle\langle\downarrow\rvert$
$0$ | $\lvert\downarrow\rangle\langle\uparrow\rvert$ | $\lvert\downarrow\rangle\langle\uparrow\rvert$
fast | $-i\varepsilon_{d,a}-\frac{\delta_{a}}{2}$ | $\lvert\sigma\rangle\langle\varnothing\rvert$ | $\lvert\sigma\rangle\langle\varnothing\rvert$
$i\varepsilon_{d,a}-\frac{\delta_{a}}{2}$ | $\lvert\varnothing\rangle\langle\sigma\rvert$ | $\lvert\varnothing\rangle\langle\sigma\rvert$
$-i(\varepsilon_{d,a}+U_{a})-\frac{\kappa_{a}}{2}$ | $\lvert\uparrow\downarrow\rangle\langle\sigma\rvert$ | $\lvert\uparrow\downarrow\rangle\langle\sigma\rvert$
$i(\varepsilon_{d,a}+U_{a})-\frac{\kappa_{a}}{2}$ | $\lvert\sigma\rangle\langle\uparrow\downarrow\rvert$ | $\lvert\sigma\rangle\langle\uparrow\downarrow\rvert$
$-i(2\varepsilon_{d,a}+U_{a})-\frac{\kappa_{a}}{2}$ | $\lvert\uparrow\downarrow\rangle\langle\varnothing\rvert$ | $\lvert\uparrow\downarrow\rangle\langle\varnothing\rvert$
$i(2\varepsilon_{d,a}+U_{a})-\frac{\kappa_{a}}{2}$ | $\lvert\varnothing\rangle\langle\uparrow\downarrow\rvert$ | $\lvert\varnothing\rangle\langle\uparrow\downarrow\rvert$

####  III.1.2 Characterization of ${\mathcal{L}_{c0}}$

Let us next analyze ${\mathcal{L}}_{c0}$. As noted earlier, we assume that the
system size of the conduction electrons is large enough such that they are
always in the equilibrium state. Assuming further that the system is at low
temperature $k_{\rm B}T\ll\varepsilon_{\rm F}$, the steady state
(${\mathcal{L}}_{c0}\hat{\rho}_{c,ss}^{(0)}=0$) is approximately given by a
Fermi sphere $\hat{\rho}_{c,ss}^{(0)}=\lvert F\rangle\langle F\rvert$, where

$\displaystyle\lvert F\rangle=\prod_{\varepsilon_{\bm{k}}<\varepsilon_{\rm
F}}\prod_{\sigma=\uparrow,\downarrow}\hat{c}^{\dagger}_{\bm{k},\sigma}\lvert
0\rangle,$ (S105)

and $\lvert 0\rangle$ is a vacuum state. Since ${\mathcal{L}}_{c0}$ is
composed of free electrons, one can obtain all the eigenvalues and eigenstates
by adding or subtracting electrons from this state to both ket and bra space.
For example, by creating or annihilating an electron above or below the Fermi
sphere, given respectively by,

$\displaystyle\lvert\bm{k}\sigma\rangle\equiv\hat{c}_{\bm{k},\sigma}^{\dagger}\lvert
F\rangle~{}(\varepsilon_{\bm{k}}>\varepsilon_{\rm
F}),\qquad\lvert\overline{\bm{k}\sigma}\rangle\equiv\hat{c}_{\bm{k},\sigma}\lvert
F\rangle~{}(\varepsilon_{\bm{k}}<\varepsilon_{\rm F})$ (S106)

then

$\displaystyle{\mathcal{L}}_{c0}(\lvert\bm{k},\sigma\rangle\langle F\rvert)$
$\displaystyle=-i\varepsilon_{\bm{k}}(\lvert\bm{k}\sigma\rangle\langle
F\rvert),$ (S107a) $\displaystyle{\mathcal{L}}_{c0}(\lvert
F\rangle\langle\bm{k},\sigma\rvert)$
$\displaystyle=+i\varepsilon_{\bm{k}}(\lvert
F\rangle\langle\bm{k},\sigma\rvert),$ (S107b)
$\displaystyle{\mathcal{L}}_{c0}(\lvert\overline{\bm{k}\sigma}\rangle\langle
F\rvert)$
$\displaystyle=+i\varepsilon_{\bm{k}}(\lvert\overline{\bm{k}\sigma}\rangle\langle
F\rvert),$ (S107c) $\displaystyle{\mathcal{L}}_{c0}(\lvert
F\rangle\langle\overline{\bm{k}\sigma}\rvert)$
$\displaystyle=-i\varepsilon_{\bm{k}}(\lvert
F\rangle\langle\overline{\bm{k}\sigma}\rvert),$ (S107d)

gives a set of eigenvalues and right eigenstates of ${\mathcal{L}}_{c0}$. The
left eigenstate is identical to the right eigenstate, which is immediate from
the property ${\mathcal{L}}_{c0}^{\dagger}={\mathcal{L}}_{c0}^{*}$. Similarly,
by further adding one more particle or a hole to these states, one obtains
further sets of eigenvalues and right eigenstates of ${\mathcal{L}}_{c0}$,

$\displaystyle{\mathcal{L}}_{c0}(\lvert\bm{k}\sigma,\overline{\bm{k}^{\prime}\sigma^{\prime}}\rangle\langle
F\rvert)$
$\displaystyle=-i(\varepsilon_{\bm{k}}-\varepsilon_{\bm{k}^{\prime}})(\lvert\bm{k}\sigma,\overline{\bm{k}^{\prime}\sigma^{\prime}}\rangle\langle
F\rvert),$ (S108a) $\displaystyle{\mathcal{L}}_{c0}(\lvert
F\rangle\langle\bm{k},\sigma,\overline{\bm{k}^{\prime}\sigma^{\prime}}\rvert)$
$\displaystyle=+i(\varepsilon_{\bm{k}}-\varepsilon_{\bm{k}^{\prime}})(\lvert
F\rangle\langle\bm{k},\sigma,\overline{\bm{k}^{\prime}\sigma^{\prime}}\rvert),$
(S108b)
$\displaystyle{\mathcal{L}}_{c0}(\lvert\bm{k}\sigma\rangle\langle\bm{k}^{\prime}\sigma^{\prime}\rvert)$
$\displaystyle=-i(\varepsilon_{\bm{k}}-\varepsilon_{\bm{k}^{\prime}})(\lvert\bm{k}\sigma\rangle\langle\bm{k}^{\prime}\sigma^{\prime}\rvert),$
(S108c)
$\displaystyle{\mathcal{L}}_{c0}(\lvert\overline{\bm{k}\sigma}\rangle\langle\overline{\bm{k}^{\prime}\sigma^{\prime}}\rvert)$
$\displaystyle=+i(\varepsilon_{\bm{k}}-\varepsilon_{\bm{k}^{\prime}})(\lvert\overline{\bm{k}\sigma}\rangle\langle\overline{\bm{k}^{\prime}\sigma^{\prime}}\rvert),$
(S108d)

etc., where

$\displaystyle\lvert\bm{k}\sigma,\overline{\bm{k}^{\prime}\sigma^{\prime}}\rangle\equiv\hat{c}_{\bm{k}\sigma}^{\dagger}\hat{c}_{\bm{k}^{\prime}\sigma^{\prime}}\lvert
F\rangle~{}(\varepsilon_{\bm{k}}>\varepsilon_{\rm
F},\varepsilon_{\bm{k}^{\prime}}<\varepsilon_{\rm F}).$ (S109)

One can continue the same procedure to obtain all sets of eigenvalues and
eigenstates.

Note crucially that, at the low-temperature regime we are interested in, the
energy of the created/removed electrons from the Fermi sphere are generically
near the Fermi surface, i.e.,
$\varepsilon_{\bm{k}},\varepsilon_{\bm{k}^{\prime}}\sim\varepsilon_{\rm F}$.
Therefore, the eigenvalues of the eigenstates with different total charges in
the ket and bra space (such as those in Eq. (S107)) are at least
$O(\varepsilon_{\rm F})$, while those with the same charge in ket and bra
space (such as those in Eq. (S108)) are generically small. We will regard the
latter (former) as the slow (fast) modes in the next subsection.

###  III.2 Derivation of Eq. (4) in the main text

Having derived the right and left eigenstates of the nonperturbative
Lindbladian ${\mathcal{L}}_{0}$, we are now in the position to derive the
effective low-energy master equation (Eq. (4) in the main text) that
incorporates the effect of c-d mixing ${\mathcal{L}}_{1}$ within the second-
order perturbation. In what follows, we will explicitly compute the matrix
element of the low-energy effective Lindbladian $({\mathcal{L}}_{\rm
eff,sd})_{n_{l},n_{r}}=(\hat{l}_{n_{l}},{\mathcal{L}}_{\rm
eff,sd}\hat{r}_{n_{r}}^{(0)})={\rm
tr}[\hat{l}_{n_{l}}^{(0)\dagger}{\mathcal{L}}_{\rm
eff,sd}\hat{r}_{n_{r}}^{(0)}]$ (see Eq. (S89)).

We will mainly be interested in the singly-occupied case, where we take the
right eigenstate $\hat{r}_{n_{r}}^{(0)}$ as a singly occupied state for the
localized electrons and a Fermi sphere for the conduction electron, given by
($\sigma_{+},\sigma_{-}=\uparrow,\downarrow$)

$\displaystyle\hat{r}_{n_{r}}^{(0)}=\hat{r}_{(\sigma_{+};F),(\sigma_{-};F)}^{(0)}=\lvert\sigma_{+}\rangle\lvert
F\rangle\langle\sigma_{-}\rvert\langle F\rvert.$ (S110)

This state has a zero eigenvalue in the unperturbed Lindbladian, i.e.,
${\mathcal{L}}_{0}\hat{r}^{(0)}_{(\sigma_{+};F),(\sigma_{-};F)}=0$. For
simplicity, we have temporarily dropped the localized electron label $a$ and
will focus on a single impurity problem (as the extension to the multiple
impurity case is straightforward). We will discuss the matrix elements that
involves $T_{1}$ modes with vacant
$\lvert\varnothing\rangle\langle\varnothing\rvert$ and double occupied state
$\lvert\uparrow\downarrow\rangle\langle\uparrow\downarrow\rvert$ later. The
left eigenstate $\hat{l}_{n_{l}}^{(0)}$ is taken from the slow variable space
defined in the previous section, see Table 1 and the final paragraph of Sec.
III.1.2. One can readily check that the zeroth order contribution in Eq. (S89)
vanishes, i.e., $({\mathcal{L}}_{\rm
eff,sd})_{n_{l},[(\sigma_{+};F),(\sigma_{-};F)]}^{(0)}={\rm
tr}[\hat{l}_{n_{l}}^{(0)\dagger}{\mathcal{L}}_{0}\hat{r}_{(\sigma_{+};F),(\sigma_{-};F)}^{(0)}]=0$.

For later convenience, we find it useful to introduce the c-d mixing
operators,

$\displaystyle\hat{V}_{d\Leftarrow
c}^{a\bm{k}\sigma}=v_{a}e^{i\bm{k}\cdot\bm{R}_{a}}\hat{d}_{\sigma,a}^{\dagger}\hat{c}_{\bm{k}\sigma},\qquad\hat{V}_{c\Leftarrow
d}^{a\bm{k}\sigma}=v_{a}^{*}e^{-i\bm{k}\cdot\bm{R}_{a}}\hat{c}_{\bm{k}\sigma}^{\dagger}\hat{d}_{\sigma,a}(=(\hat{V}_{d\Leftarrow
c}^{a\bm{k}\sigma})^{\dagger}),$ (S111)

that simplifies the notation of perturbative Lindbladian to

$\displaystyle{\mathcal{L}}^{(1)}\rho$ $\displaystyle=$
$\displaystyle-i[\hat{H}_{cd},\hat{\rho}]=-i\sum_{a}\sum_{\bm{k},\sigma}\big{[}\hat{V}_{d\Leftarrow
c}^{a\bm{k}\sigma}+\hat{V}_{c\Leftarrow d}^{a\bm{k}\sigma},\hat{\rho}\big{]}.$
(S112)

The state $\hat{r}^{(0)}_{\sigma_{+},\sigma_{-};F}$ is perturbed with
${\mathcal{L}}^{(1)}$ as,

$\displaystyle{\mathcal{L}}^{(1)}\hat{r}_{(\sigma_{+};F),(\sigma_{-};F)}^{(0)}$
$\displaystyle=$
$\displaystyle-i\sum_{\bm{k},\sigma}\big{[}(\hat{V}_{d\Leftarrow
c}^{a\bm{k}\sigma}+\hat{V}_{c\Leftarrow
d}^{a\bm{k}\sigma})(\lvert\sigma_{+}\rangle\lvert
F\rangle\langle\sigma_{-}\rvert\langle F\rvert)-(\lvert\sigma_{+}\rangle\lvert
F\rangle\langle\sigma_{-}\rvert\langle F\rvert)(\hat{V}_{d\Leftarrow