Quantum Zeno Effect in Noisy Integrable Quantum Circuits for Impurity Models

Abstract

We theoretically study the open quantum system dynamics (in the Trotterized limit) of integrable quantum circuits in the presence of onsite dephasing noise with a spin-1/2 “impurity” interacting at the edge. Using a combination of Bethe Ansatz (BA) and exact diagonalization (ED), we study the dynamics of both the bulk and the impurity for the XXX (Heisenberg) and the XX qubit chains in the presence and absence of bulk noise. In the absence of noise, we show that the impurity exhibits two distinct phases, the bound mode phase where the impurity keeps oscillating in time, and the Kondo phase where it decays with Kondo time $t_{K}$. Turning on the bulk dephasing noise, we find for the two models that in the long time limit in both regimes the quantum Zeno effect takes place where the dynamics of the impurity magnetization slows down as the noise strength $\gamma$ increases. The impurity magnetization in the bound mode regime shows the opposite effect, decaying faster as the noise strength increases for short times ($t\ll 1/\gamma$). We show that the bulk KPZ dynamics of the XXX model is converted to diffusive dynamics as in the XX case studied before by V. Alba, driving both systems to the Zeno effect slowing down the impurity dynamics in the long time limit.

Introduction: The rise of noisy intermediate scale quantum (NISQ) [1, 2, 3] devices in a wide range of physical platforms including superconducting qubits [4, 5, 6], trapped ions [7, 8, 9, 10], neutral atoms [11, 12, 13], and circuit quantum electrodynamics [14, 15, 16, 17, 18], has made understanding the effects of noise on these devices and on quantum many-body systems (QMBS) in general a question of paramount importance [19, 20, 21, 22]. Simulating large generic QMBS on a classical computer is extremely difficult  [23, 24, 25] as the size of the Hilbert space grows exponentially with the size of the system. To circumvent this challenge, quantum computers can be utilized to study QMBS [24, 25] (as they are built to hold such an exponentially large Hilbert space). This can be achieved by Trotterizing the continuous time evolution of a QMBS as a quantum circuit [26, 27]. While, generically solving a quantum circuit remains an equally challenging problem, integrable quantum circuits [28, 29, 30, 31, 32, 33], in particular, based on the $R$-matrices that solve the Yang-Baxter equation  [34, 35], can simulate the unitary evolution of initial states under integrable Hamiltonians [36, 37, 38, 39]. Integrability allows more control and insights into the behavior of these complex systems garnering significant attention recently  [40, 41, 42].

All quantum devices currently accessible are under a significant influence of noise [43, 44, 45, 46, 47, 48, 49, 50, 51, 52]. Noise affects the dynamics and transport properties of quantum circuits resulting in dramatic changes [19, 20, 21]. Interesting effects arise such as dynamical phase transitions [53, 54], the quantum Zeno effect (where increasing the strength of noise slows down quantum dynamics [55, 56]) and the quantum Mpemba effect (where non-equilibrium quantum systems relax faster when they are further from equilibrium [57, 58, 59, 60]). Thus, studying the effect of noise on quantum circuits is of crucial importance to understand the experiments performed on NISQ devices [61]. Here, we study a circuit that can be solved exactly in the presence of dephasing noise.

Models: To begin with, we construct an integrable quantum circuit that corresponds to a unitary operator

$$\begin{split}\mathbb{U}^{\eta}(\lambda,d)&=K_{0,1}^{\eta}(\lambda,d)\prod_{j=1}^{L/2-1}\check{R}_{2j,2j+1}^{\eta}(\lambda)\check{R}_{2j-1,2j}^{\eta}(\lambda),\end{split}$$ (1)

where $L$ is the number of qubits, $\check{R}_{j,j+1}$ given explicitly in [62] is the braid form of the trigonometric $R$-matrix of the six-vertex model [63, 36] acting on the $j$ and $j+1$ qubits, $K_{01}^{\eta}(\lambda,d)=\check{R}_{01}^{\eta}(-d)\check{R}^{\eta}_{01}(\lambda+d)$ is the boundary $K$-matrix that solves the reflection algebra [64, 65, 66], and $d$, $\lambda$ and $\eta$ are free parameters. In the Trotterziation limit [67, 68], choosing $\lambda=-i\delta t\sinh\eta$, this quantum circuit simulates the unitary evolution $\exp(-iHt)=\lim_{\delta t\to 0}\mathbb{U}(\lambda,d,\eta)^{\frac{t}{\delta t}}$  for an integrable Hamiltonian $H$ [25, 27, 69], where

$$\begin{split}H=&\sum_{j=1}^{L-1}\frac{1}{2}(\sigma_{j}^{x}\sigma_{j+1}^{x}+\sigma_{j}^{y}\sigma_{j+1}^{y}+\Delta\sigma_{j}^{z}\sigma_{j+1}^{z})\\ &+\frac{J}{2}\left(\sigma_{0}^{x}\sigma_{1}^{x}+\sigma_{0}^{y}\sigma_{1}^{y}+\Delta^{\prime}\sigma_{0}^{z}\sigma_{1}^{z}\right)\\ \end{split}$$ (2)

with $\Delta\mathord{=}\cosh(\eta)$, $J\mathord{=}\frac{\sinh^{2}(\eta)\cosh(d)}{\sinh^{2}(\eta)-\sinh^{2}(d)}$ and $\Delta^{\prime}\mathord{=}\frac{\cosh(\eta)}{\cosh(d)}$. This Hamiltonian describes the $XXZ$ spin chain with a boundary impurity. At the isotropic XXX ($\Delta\mathord{=}\Delta^{\prime}\mathord{=}1)$ and the XX ($\Delta\mathord{=}\Delta^{\prime}\mathord{=}0$) points, the model is integrable for arbitrary boundary coupling $J$. For generic $\Delta$, the boundary coupling ($J$) and the boundary anisotropy parameter ($\Delta^{\prime}$) are related and expressed in terms of a single free parameter $d$ to maintain the integrability [70, 71, 72, 73, 74, 75, 76]. In the ground state, the impurity is screened due to the Kondo effect when $\bar{d}\mathord{=}d/\eta$ is purely imaginary or $0<\bar{d}<\frac{1}{2}$ and it is screened by a bound mode exponenitially localized at the boundary when $1>\bar{d}>\frac{1}{2}$ [70].

To study the effect of dephasing noise on the integrable quantum circuit, we solve the continuous-time dynamics of the noisy integrable Hamiltonian using the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation, $\mathcal{L}[\rho(t)]=i\partial_{t}\rho(t)$ [20, 19] that governs the dynamics of the density matrix $\rho(t)$ with a Liouvillian [77, 78, 79, 80]

$$\mathcal{L}(\rho)=-\left[\rho,H\right]+i\sum_{j=1}^{L-1}\gamma\left(L_{j}\rho L_{j}^{\dagger}-\frac{1}{2}\left\{L^{\dagger}_{j}L_{j},\rho\right\}\right),$$ (3)

where the jump operators are $L_{j}\mathord{=}\sigma^{z}_{j}$, $\gamma$ is the dephasing rate, $H$ is the Hamiltonian in Eq. (2), $\{A,B\}$ and $[A,B]$ denote the anticommutator and the commutator, respectively. In the following, we will evolve the system from two initial states, Néel state and domain wall state expressed as: $\rho_{N/D}=(\prod_{z\in\mathcal{S}_{N/D}}\sigma^{+}_{z})\ket{0}\bra{0}(\prod_{z\in\mathcal{S}_{N/D}}\sigma^{-}_{z})$ where $\mathcal{S}_{N/D}$ denote the set of sites with $z\text{ }\mathord{=}\text{ }0,2,4,...,L$ for the Néel (N) state $\rho_{N}$, and $z\text{ }\mathord{=}\text{ }0,1,2,...,L/2$ for the domain wall (D) state $\rho_{D}$. Here, $\ket{0}$ refers to the state with all spins down (qubits $0$) state. In both cases, we compute the time-dependent expectation value of magnetization at site $j$, $S^{z}_{j}(t)\mathord{=}\mathrm{Tr}\{\sigma^{z}_{j}\rho(t)\}$.

XX chain: Let us start from the XX point as it allows for a simpler treatment with similar phenomenology [71]. It is convenient to work in the fermionic language via a Jordan-Wigner transformation $c^{\dagger}_{j}\mathord{=}\sigma^{+}_{j}\prod_{l<j}\sigma^{z}_{l}$. The magnetization is then related to the two-point correlator $G(x,y;t)=\mathrm{Tr}[c^{\dagger}_{x}c_{y}\rho(t)]$ as $S^{z}_{j}(t)=2G(j,j;t)-1$. A closed set of equations gives the dynamics of the two-point correlators because the Hamiltonian and the jump operators are both quadratic in fermions [81]. Therefore, the dynamics of the two-point correlation function in the noisy XX model can be mapped to a two-particle quantum mechanical problem [62]. The correlator can be vectorized as a state (see [62]) $|G(t))\in\mathbb{C}^{L^{2}}$ where $|G(t))=\sum_{x,y}G(x,y;t)|x,y)$ and $|x,y)$ forms an orthonormal basis in $\mathbb{C}^{L^{2}}$. Its evolution from an initial state $|G_{in})=\sum_{x,y}\mathrm{Tr}[G(x,y{;0})\rho_{N/D}]|x,y)$ reads $|G(t))=e^{-i\hat{h}t}|G_{in})$ with Liouvillian $\hat{h}=\sum_{x,y,x^{\prime},y^{\prime}}h^{x,y}_{x^{\prime},y^{\prime}}|x,y)(x^{\prime},y^{\prime}|$ and

$$h^{x,y}_{x^{\prime},y^{\prime}} =J_{x}(\delta_{x,x^{\prime}+1}+\delta_{x,x^{\prime}-1})-J_{y}(\delta_{y,y^{\prime}+1}+\delta_{y,y^{\prime}-1})$$ (4)

$$+4i\gamma(\delta_{x,x^{\prime}}\delta_{y,y^{\prime}}-1)+2i\gamma(\delta_{x^{\prime},0}+\delta_{y^{\prime},0}),$$

where $J_{x}=1+\delta_{x,0}(J-1)$. The Liouvillian right eigenmodes $|G_{E})=\sum_{x,y}G_{E}(x,y)|x,y)$ satisfy the stationary Schödinger equation $\hat{h}|G_{E})=E|G_{E})$ with eigenvalue $E$ in terms of which we expand the time-dependent state $|G(t))=\sum_{E}(\bar{G}_{E}|G_{in})e^{-iEt}|G_{E})$. Here $(\bar{G}_{E}|$ is the left eigenmode of the Liouvillian satisfying $(\bar{G}_{E}|\hat{h}^{\dagger}=E^{*}(\bar{G}_{E}|$ and it forms a bi-orthognal basis with $|G_{E})$ as $\sum_{E}|G_{E})(\bar{G}_{E}|=\mathbb{I}$. The time-dependent impurity magnetization (at site $j=0$) can be written as

$$S^{z}_{0}(t)=2\sum_{E}\alpha_{E}^{N/D}e^{-iEt}-1,$$ (5)

where the amplitude $\alpha_{E}^{N/D}=G_{E}(0,0)(\bar{G}_{E}|G_{in})$ depends on the initial state $\rho_{in}=\rho_{N/D}$ as follows $(\bar{G}_{E}|G_{in})=\sum_{z{\in\mathcal{S}_{N/D}}}\bar{G}_{E}(z,z)$. Thus, $S^{z}_{0}(t)$ requires the knowledge of the Liouvillian spectrum $E$ and of the amplitude $\alpha_{E}$ which encodes the initial conditions.

Figure 1: (a) The parameter space $(J,\gamma)$ is divided into two distinct regimes: the bound-mode regime, and the Kondo regime. The real part of the double-bound mode energy is computed for $L=30$ is given in the heatmap. For large $\gamma\mathord{\gg}1$, the boundary between two regimes is almost $J\mathord{=}\gamma/2$ as shown in the dashed line in (a). The Liouvillian spectrum and amplitude $\alpha_{E}=G_{E}(0,0)(\bar{G}_{E}|G_{in})$ for Néel state in the two-particle sector with dephasing rate $\gamma=0.25$ and impurity coupling $J=0.5$ (b) in the Kondo regime, $\gamma=2$, $J=0.5$ (c) in the Kondo regime, and $\gamma=2$, $J=3$ (d) in the bound mode regime due to the existence of the double-bound mode. The data points labeled in the legend with different shapes and numbers represent four different kinds (respectively correspond to the four kinds of solutions solved analytically in the main text). The scattering color plots are numerical results solved with system size $L=40$ where the heatmap shows absolute values of amplitude $\alpha_{E}^{N}$ for each eigenmode. The modes with larger $\alpha_{E}^{N}$ dominate the impurity dynamics: the double-bound mode in the bound mode regime (d), and the modes within the diffusive band in the Kondo regime (b, c). With a small noise $\gamma<1$, the diffusive band is mixed with the other bands (b), and when $\gamma>1$ the diffusive band is separated from the other bands (c,d).

Figure 2: The dynamics of the XX and XXX impurity models evolving from a Néel initial state under a Lindbladian evolution with and without dephasing noise. Different colors represent different dephasing rates $\gamma$, solid lines for $J=0.5$ in the Kondo regime, and dashed lines for $J=3$ in the bound-mode regime. The upper panel shows the data for the XXX chain for $L=10$ and the bottom panel for the XX chain for $L=40$. For both XX and XXX chains, in the bound mode phase, the signal decays faster as the dephasing rate increases. However, in the Kondo phase, the impurity decay is slower as the dephasing rate increases. The insets in each panel show the data on a log-linear scale displaying the same power law decay independent of $J$ in the long-time limit.

We now proceed to solve for the spectrum $E$ of the Liouvillian $\hat{h}$ defined in Eq. (4). In the absence of the impurity, this Liouvillian can be mapped to an integrable non-Hermitian Hubbard model after a unitary transformation [81, 82, 83, 84, 85], here we generalize it to include the impurity as discussed in [62]. The Liouvillian in Eq. (4) can thus be solved via the Bethe Ansatz (BA) [86, 87, 88]. It has eigenvalue $E(k,q)=2\cos k+2\cos q-4i\gamma$ with the eigenmodes wavefunction taking the BA form

$$G_{k,q}(x,y)=\frac{(-1)^{y}}{\mathcal{N}}\mathcal{S}\sum_{\sigma,\tau=\pm 1}F_{\sigma k,\tau q}(x,y)e^{i\sigma kx+i\tau qy},$$ (6)

where $F_{\sigma k,\tau q}(x,y)=A_{\sigma k,\tau q}\theta(y-x)+B_{\sigma k,\tau q}\theta(x-y)$, $\mathcal{N}$ is the normalization factor such that $|G_{{k,q}}|^{2}=1$ and $\mathcal{S}$ is the symmetrizer. The Liouvilian in Eq.(4) also has antisymmetrized eigenmodes [82], but only symmetrized eigenmodes participate in dynamics since the initial configuration is symmetric, i.e. $G(x,y;0)=G(y,x;0)$. Here the amplitudes are related to the boundary scattering matrix $K(k)$ and the bulk scattering matrix $S(k,q)$ as: $A_{\sigma k,\tau q}=S(\sigma k,\tau q)B_{\sigma k,\tau q}$ and $A_{+k,\tau q}=K(k)A_{-k,\tau q}$ [81, 82, 83, 84] with

$$\begin{split}&K(k)=-\frac{e^{ik}+e^{-ik}-2i\gamma-J^{2}e^{ik}}{e^{ik}+e^{-ik}-2i\gamma-J^{2}e^{-ik}},\\ &S(k,q)=\frac{\sin k-\sin q+2\gamma}{\sin k-\sin q-2\gamma}.\end{split}$$ (7)

The quantization condition for quasimomenta $k$ and $q$ are obtained from the Bethe Ansatz equations (BAE) [84]

$$\begin{split}e^{2ik(L+1)}&=K(k)S(k,q)S(-k,q),\\ e^{2iq(L+1)}&=K(q)S(q,k)S(-q,k).\end{split}$$ (8)

If $(k,q)$ is a solution of the BAE (8), then $(\pm k+2\pi,\pm q+2\pi)$ is also a solution. This allows us to restrict the Brillouin zone to $\Re k\in(0,\pi)$ with $\Re k$ ($\Im k$) as the real (imaginary) part of $k$. Furthermore, its complex conjugate $(k^{*},q^{*})$ also being a solution, indicates that if $E$ is a Liouvillian eigenvalue, then $-E^{*}$ is also an eigenvalue [82]. The quasimomenta $(k,q)$ satisfying the BAEs fall into one of the following four cases in the limit $L\to\infty$. Here, we briefly summarize the resuls and details can be found in [62].

1.Free band: When both S-matrices are regular $S,K\neq 0,\infty$, the quasimomenta $k$ and $q$ are complex, with imaginary parts vanishing as $\sim\frac{1}{L}$, i.e $k=\frac{\pi n_{1}}{L+1}+\mathcal{O}(L^{-1})$ and $q=\frac{\pi n_{2}}{L+1}+\mathcal{O}(L^{-1})$, for $n_{1},n_{2}=1,\ldots,L$. These solutions exist for all parameters.

2. Bound mode band: When $K\mathord{=}\infty$, quasimomenta are $k_{\pm}=i\log\left(\frac{\pm\sqrt{J^{2}-1-\gamma^{2}}-i\gamma}{J^{2}-1}\right)$ and $q=\frac{\pi n}{L+1}+\mathcal{O}(L^{-1})$ for $n=1,\ldots,L$. If $\Im k_{\pm}<0$, it is a bound mode; otherwise, it is not a valid solution. When both $\Im k_{\pm}<0$, another mode with $E=2\cos k_{+}+2\cos k_{-}-4i\gamma$ exists.

3. Diffusive band: When $S\mathord{=}\infty$, the eigenspectra are $E_{n}=4i\sqrt{\gamma^{2}-\sin\left(\frac{\pi n}{2L}\right)^{2}}-4i\gamma+\mathcal{O}(L^{-3})$ with $n=0,1,2,\ldots,L$. Gapless excitations with $E_{n}=-i\frac{\pi^{2}n^{2}}{2\gamma L^{2}}+\mathcal{O}(L^{-3})$ [81, 82], dominate the late time dynamics. These solutions with $\Re E_{n}=0$ exist for all parameters.

4. Double-bound mode: There are two solutions where the particles are boundary-localized. When $\gamma\mathord{<}1$, the eigenvalues can be solved perturbtively $E_{\pm}=4\cos k_{\pm}+4i\gamma\left[\sum_{x=1}^{L}\mathcal{N}^{2}e^{i4k_{\pm}x}+\frac{\mathcal{N}}{J^{2}}-1\right]+\mathcal{O}(\gamma^{2})$, where $k_{\pm}=\frac{i}{2}\log(J^{2}-1)+\pi/2\pm\pi/2$ and $\mathcal{N}=\frac{1}{J^{2}}+\frac{e^{-2ik_{\pm}}}{1-e^{-2ik_{\pm}}}$. When $\gamma>1$, the eigenvalues solved with an effective two-site Liouvillian are $E_{\pm}=\pm\sqrt{4J^{2}-\gamma^{2}}-i\gamma$.

Fig. 1 shows the spectrum of the Liouvillian in the Kondo and bound regimes. We present the cross-over diagram that delineates these two regimes in Fig. 1(a) through the value of the real part of the Double bound mode $\Re E_{DBM}$. A non-zero $\Re E_{DBM}$ signifies the bound mode regime, while a zero value indicates the Kondo regime. We further analyze the spectrum’s real and imaginary parts in the Kondo regime in Fig. 1(a) and (b) as well as the bound mode regime in Fig. 1 (c). We see the free band (spectrum 1; circles), the bound mode band (spectrum 2; crosses), and the diffusive band (spectrum 3; star) in the Kondo regime, where the majority of the weight in $\alpha_{E}^{N}$ lies in the diffusive band. In the bound mode regime, we also see the appearance of double bound modes (spectrum 4; diamonds) that have the largest weight in the coefficients $\alpha_{E}^{N}$ and hence they dominate the early-time impurity dynamics $S^{z}_{0}(t)$.

In [71], we showed that at the XX point without noise, there is a boundary phase transition at $J=\sqrt{2}$. When $J<\sqrt{2}$, the impurity is screened by many-body Kondo effect whereas when $J>\sqrt{2}$, the impurity is screened by single particle bound mode localized at the edge. Here, we turn to study the dynamical characterization of the boundary phases. In the bound mode phase ($J>\sqrt{2}$), the time-dependent magnetization shows oscillatory behavior $S^{z}_{0}(t)\sim\cos(\frac{J^{2}}{\sqrt{J^{2}-1}}t)$ when $t\to\infty$. However in the Kondo phase ($J<\sqrt{2}$), the impurity magnetization leaks to the bulk with the time scale proportional to the Kondo time $t_{K}=\sqrt{|J^{2}-1|}/J^{2}$. Detailed calculation showing this sharp phase transition in dynamical quantities is in [62].

When noise is introduced ($\gamma>0$), there is no longer a sharp phase transition. Instead, the impurity magnetization $S^{z}_{0}(t)$ decays similarly in both regimes at late times, exhibiting the Zeno effect, where the dynamics slows down as noise strength increases. This behavior is driven by the diffusive band as it has largest overlap $\alpha_{E}^{N}$. The modes in the diffusive band have $|\Im E(\gamma)|\propto\gamma^{-1}$, so the dynamics slows down as $\gamma$ increases and exhibits the quantum Zeno effect. In contrast, the early-time dynamics differ between the two regimes. In the bound mode regime, the impurity dynamics is dominated by the double-bound mode which has the largest overlap $\alpha_{E}^{N}$ and $|\Im E(\gamma)|\propto\gamma$. Thus, the impurity decays faster as $\gamma$ increases. However, in the Kondo regime, the early-time dynamics are also influenced by the diffusive band—similar to the long-time limit where the decay rate decreases with increasing noise, resulting the quantum Zeno effect. The model parameters $J$ and $\gamma$ corresponding to these two regimes are shown in Fig. 1(a).

Numerical results for $S^{z}_{0}(t)$ starting from the Néel state at the XX point are shown in the bottom panel of Fig. 2 with two situations when the system is in the Kondo regime $J=0.5$ and in the bound-mode regime $J=3$. As seen from the data, for the early time in the bound-mode regime, the noise adds a decaying envelope that is accelerated as the dephasing rate is increased. This is an expected result because the boundary mode is essentially a two-qubit state, and qubits subjected to dephasing noise decay faster as the decay rate increases. In the Kondo regime or after the double bound mode decays i.e. $t\sim\gamma^{-1}$, the noise slows the leakage of magnetization and leads to a quantum Zeno effect as the dephasing rate increases. In the long time limit, independent of the boundary coupling $J$, the impurity magnetization falls off as $S^{z}_{0}(t)\sim e^{-t/\tau}$ with $\tau\propto\frac{1}{\gamma L^{2}}$ in both regimes. This behavior of the impurity in the Kondo regime is consistent with the heating of the bulk induced by Hermitian Lindblad operators [89, 90, 91], driving the impurity to (asymptotically free) weak coupling behavior [92, 93, 94, 95].

Figure 3: Data collapse for $S^{z}(x,t)=S^{z}(\xi)$ with $\xi=xt^{-\frac{1}{2}}$ in the XX case (filled markers, with system size $L=40$) and XXX case (hollow markers, with system size $L=10$) starting from a domain wall initial state (as shown in the color plot in the right inset for XX case with $\gamma=2$) under various dephasing rates $\gamma$. Different marker shapes correspond to the spin configuration at different time slices, and different marker color refers to different dephasing rates. The left inset shows the dependence of the diffusion constant $D$ on $\gamma$ (filled markers for XX case and hollow for XXX case) and the dashed lines are the analytical result for the noisy XX chain $D=\frac{1}{2\gamma}$.

Although no noise is applied at the impurity site, the Zeno effect on impurity is present in the Kondo regime at all times and at late times ($t>\gamma^{-1}$) in the bound mode phase, because noise changes the bulk transport property from ballistic when the system is at the XX point to diffusive [96] with diffusion constant $D=\frac{1}{2\gamma}$ as shown in [82].

Heisenberg (XXX) chain: Next, we examine the Trotterized dynamics of the interacting circuit corresponding to the isotropic Heisenberg model ($\Delta=1$) with a boundary impurity. The equilibrium boundary phases have been studied in [72, 76, 97]. We show here that the bound mode and Kondo phases also exist dynamically, manifesting themselves in a similar way as in the XX case: the bound mode phase exists when $J>4/3$ where the impurity magnetization oscillates in time, and the Kondo phase when $J<4/3$ where it decays as shown in Fig. 2 in the noiseless case. The bulk transport properties in the noiseless case belong to the KPZ universality class as shown in [98, 99].

Turning on the noise we find that it impacts both the bulk and boundary dynamics similarly as in the XX case discussed above. To show the bulk dynamics is diffusive for XXX case, we prepare the system in a domain wall state where all sites on the left (right) half are up (down) and compute the magnetization $S^{z}_{j}(t)$ as time evolves. When the system is diffusive, in the continuum limit as lattice spacing $a\to 0$ and system size $L\to\infty$, the magnetization $S^{z}_{j}(t)\to S^{z}(x,t)$ satisfies the diffusion equation i.e. $(\partial_{t}-D\partial_{x}^{2})S^{z}(x,t)=0$, with initial condition $S^{z}(x,t\mathord{=}0)=\mathrm{Sign}(L/2-x)$. The numerical results for the dynamics of both XX and XXX cases are shown in Fig. 3, where all the data collapses on the curve $S^{z}(\xi)=-\mathrm{Erf}(\xi/\sqrt{D})$ with $\xi=x/\sqrt{t}$ which is the solution of diffusion equation from the domain wall initial configuration. For large values of $\gamma$, we find clear numerical evidence that $D=\frac{1}{2\gamma}$ as in the XX case, but for smaller values of $\gamma$, it deviates (as shown in the left inset in Fig. 3) due to finite size and time effects [62] and we therefore expect that the relation $D=\frac{1}{2\gamma}$ holds for each value of $\gamma$ considered. At the same time, the bulk diffusive dynamics reflects itself in the impurity dynamics as shown in the impurity magnetization dynamics from the initial Néel state with and without noise in the upper panel of Fig. 2. In the presence of noise, both Kondo and bound-mode regimes exist, and the impurity in each regime shows a similar behavior as in the free XX chain. In the bound mode regime, the impurity magnetization initially decays faster as the noise strength is increased, and after the bound modes vanish ($t>\gamma^{-1}$), the impurity experiences the quantum Zeno effect whereas in the Kondo regime, the impurity magnetization exhibits the Zeno effect in all time scales.

Conclusion: We found that dephasing noise affects impurity behavior differently in the Kondo and bound mode regimes. In the Kondo regime, increasing $\gamma$ slows the impurity magnetization decay, leading to a quantum Zeno effect. In the bound mode regime, magnetization initially decays faster but eventually exhibits the Zeno effect, similar to the Kondo regime. Additionally, the KPZ dynamics [98, 99] in the isotropic Heisenberg model become diffusive when the diffusion rate is introduced [100, 101].

Acknowledgment: We thank T. Giamarchi, V. Alba, F. Essler, and D. Huse for useful discussions. J.H.P. is partially supported by the Army Research Office Grant No. W911NF-23-1-0144 and the Office of Naval Research grant No. N00014-23-1-2357 (J.H.P.).

In the supplementary material, we first map the circuit dynamics to the Lindbladian evolution. Then at the XX point, we show that the corresponding Lindbladian can be mapped to an integrable Hubbard Hamiltonian with imaginary onsite interaction in and boundary hopping defect in doubled Liouville-Fock space. We then solve the dynamics of the noiseless model with the free-fermion technique. After that, we solve the non-Hermition Hubbard model in the two-particle sector to discuss the dynamics in the presence of noise. We also discuss the derivation of the Bethe-Ansatz equations for the complex onsite Hubbard model in an arbitrary particle sector. In the final section, we show the details of numerical results.

I Mapping from Quantum circuits to Liouvillian dynamics

Consider the trigonometric six vertex R-matrix [36, 63]

$$R(u)=\left(\begin{array}[]{cccc}1&0&0&0\\ 0&\frac{\sinh(u)}{\sinh(u+\eta)}&\frac{\sinh(\eta)}{\sinh(u+\eta)}&0\\ 0&\frac{\sinh(\eta)}{\sinh(u+\eta)}&\frac{\sinh(u)}{\sinh(u+\eta)}&0\\ 0&0&0&1\\ \end{array}\right)$$ (S.1)

which is a solution of the Yang-Baxter equation [86]

$$R_{12}(u-v)R_{13}(u)R_{23}(v)=R_{23}(v)R_{13}(u)R_{12}(u-v).$$ (S.2)

Notice that the Yang-Baxter equation remains satisfied if we shift $u_{i}\to u_{i}-\theta_{i}$ where $\theta_{i}$ are arbitrary inhomogeneous parameters.

Let us consider the following two single-row transfer matrices

$$T_{A}(u) =R_{A,L}(u-\theta_{L})R_{A,L-1}(u-\theta_{L-1})\cdots R_{A,2}(u-\theta_{2})R_{A,1}(u-\theta_{1})$$

$$\hat{T}_{A}(u) =R_{A,1}(u+\theta_{1})R_{A,2}(u+\theta_{2})\cdots R_{A,L-1}(u+\theta_{L-1})R_{A,L}(u+\theta_{L}),$$ (S.3)

where $1,2,\ldots,L$ are the labels for the physical space and $A$ is an auxiliary space.

Furthermore, denoting the physical space of impurity as $0$, we define

$$K_{A}(u)=R_{A,0}(u-\theta_{0}-d)R_{A,0}(u+\theta_{0}+d),$$ (S.4)

such that the $K$ matrices satisfy the reflection equations

$$R_{ij}(\lambda-\mu)K_{iA}(\lambda)R_{ji}(\lambda+\mu)K_{jA}(\mu)=K_{jA}(\mu)R_{ij}(\lambda+\mu)K_{iA}(\lambda)R_{ji}(\lambda-\mu).$$ (S.5)

Now, we define the monodromy matrix

$$\Xi(u)=T_{A}(u)K_{A}(u)\hat{T}_{A}(u).$$ (S.6)

The trace of the monodromy matrix over the auxiliary space is defined as the double-row transfer matrix [66]

$$t(u)=\operatorname{tr}_{A}\Xi(u).$$ (S.7)

We fix the inhomogeneity parameters [28]

$$\theta_{0}=\theta_{2}=\theta_{4}=\ldots=\theta_{L}=\frac{\lambda}{2}\quad\quad\theta_{1}=\theta_{3}=\theta_{5}=\ldots=\theta_{L-1}=-\frac{\lambda}{2}.$$ (S.8)

Such that when the transfer matrix Eq.(S.7) is evaluated at $u=\lambda$, we obtain

$$\mathbb{U}=t(k)=\mathcal{U}_{0,1}U_{2,3}U_{4,5}\ldots U_{L-1,L}U_{1,2}U_{3,4}U_{5,6}\ldots U_{L-2,L},$$ (S.9)

where

$$U_{j,j+1} =\check{R}_{j,j+1}(\lambda)=P_{j,j+1}R_{j,j+1}(\lambda),$$

$$\mathcal{U}_{0,1} =R_{0,1}(-d)R_{0,1}(\lambda+d)=\check{R}_{0,1}(-d)\check{R}_{0,1}(\lambda+d).$$ (S.10)

Here $P_{ij}=\frac{1}{2}(\mathbb{I}_{i,j}+\sigma^{x}_{i}\sigma^{x}_{j}+\sigma^{y}_{i}\sigma^{y}_{j}+\sigma^{z}_{i}\sigma^{z}_{j})$ is the permutation matrix and explicit form of the $\check{R}$ matrix reads

$$\check{R}^{\eta}_{i,j}(\lambda)=\frac{1}{2}\left(\mathbb{I}_{i,j}+\sigma^{z}_{i}\sigma^{z}_{j}\right)+\frac{\sinh(\lambda)}{2\sinh(\eta+\lambda)}\left(\sigma^{x}_{i}\sigma^{x}_{j}+\sigma^{y}_{i}\sigma^{y}_{j}\right)+\frac{\sinh(\eta)}{2\sinh(\lambda+\eta)}\left(\mathbb{I}_{i,j}-\sigma^{z}_{i}\sigma^{z}_{j}\right).$$ (S.11)

These local unitary can be written as

$$U_{j,j+1} =e^{-i\frac{J}{2}\delta t(\sigma^{x}_{j}\sigma^{x}_{j+1}+\sigma^{y}_{j}\sigma^{y}_{j+1}+\Delta(\sigma^{z}_{j}\sigma^{z}_{j+1}-1))},$$ (S.12)

$$\mathcal{U}_{0,1} =e^{-i\frac{J^{\prime}}{2}\delta t(\sigma^{x}_{j}\sigma^{x}_{j+1}+\sigma^{y}_{j}\sigma^{y}_{j+1}+\Delta^{\prime}(\sigma^{z}_{j}\sigma^{z}_{j+1}-1))}$$ (S.13)

upon identifying

$$e^{iJ\delta t(\Delta\mp 1)}=\frac{e^{\eta}\pm e^{\lambda}}{e^{\lambda+\eta}\pm 1}\quad\quad\text{and}\quad\quad e^{iJ^{\prime}\delta t(\Delta^{\prime}\mp 1)}=\frac{\left(e^{d+\eta}\pm 1\right)\left(e^{d+\lambda}\pm e^{\eta}\right)}{\left(e^{d}\pm e^{\eta}\right)\left(e^{d+\lambda+\eta}\pm 1\right)}.$$ (S.14)

Figure S.1: Illustrations of an integrable quantum circuit with spin-1/2 qubits in an initial classical Néel state (empty circles are qubits in the up states and black filled circles are in the down qubit state). Each rectangle represents a local unitary matrix; blue ones are the braid form of the trigonometric six-vertex R-matrix $\check{R}(\lambda)$ and red ones are the boundary reflection matrix $K(\lambda,d)=\check{R}(-d)\check{R}(\lambda+d)$. The dashed rectangle shows one layer of the time evolution operator $\mathbb{U}^{\eta}(\lambda,d)$ given by Eq.(1).

The circuit in Fig. S.1 corresponds to the unitary matrix of the following form

$$\begin{split}&\mathbb{U}(\delta t)=\mathcal{U}_{0,1}\prod_{j\in 2\mathbb{Z}}U_{j,j+1}\prod_{j\in 2\mathbb{Z}-1}U_{j,j+1}=1-iH\delta t+\mathcal{O}(\delta t^{2}),\end{split}$$ (S.15)

where in the Trotterization limit, we obtain an integrable Hamiltonian with boundary impurity of the form

$$H=\sum_{i=1}^{L-1}\frac{J}{2}(\sigma_{i}^{x}\sigma_{i+1}^{x}+\sigma_{i}^{y}\sigma_{i+1}^{y}+\Delta(\sigma_{i}^{z}\sigma_{i+1}^{z}-1))+\frac{J^{\prime}}{2}\left(\sigma_{0}^{x}\sigma_{1}^{x}+\sigma_{0}^{y}\sigma_{1}^{y}+\Delta^{\prime}(\sigma_{0}^{z}\sigma_{1}^{z}-1)\right),$$ (S.16)

where the parameters in two representations are related by $\Delta=\cosh(\eta),\Delta^{\prime}=\frac{\cosh(\eta)}{\cosh(d)}$ and $J^{\prime}=J\frac{\sinh^{2}(\eta)\cosh(d)}{\sinh^{2}(\eta)-\sinh^{2}(d)}$.

Notice that both $J^{\prime}$ and $\Delta^{\prime}$ are expressed through a single free parameter $d$ hence the model is not integrable for arbitrary $\Delta^{\prime}$ and $J^{\prime}$. However, at the isotropic $\eta=0$ [76, 72] and XX $\eta=i\frac{\pi}{2}$ [71] limit, the boundary and bulk anisotropy parameter becomes equal i.e. $\Delta^{\prime}=\Delta$ and $J^{\prime}$ becomes an arbitrary free parameter.

Taking into account the presence of noise while applying the gates, the density matrix of the qubits can be described by a discrete-time evolution as

$$\rho(t+\delta t)=\sum_{j}M^{\dagger}_{j}\rho(t)M_{j},$$ (S.17)

with $M_{0}=1-\Delta t(K+iH)$, for $j\neq 0$, $M_{j}=\sqrt{\delta t}L_{j}$ and $K=\frac{1}{2}\sum_{j}L^{\dagger}_{j}L_{j}\delta t$ and $L_{j}$ are the Lindbladian jump operators. Here, we will concentrate solely on the case where $L_{j}=\sigma^{z}_{j}$. Taking $\delta t\to 0$, we obtain the Lindbladian-master equation [19]

$$\frac{d}{dt}\rho(t)=\mathcal{L}(\rho),$$ (S.18)

where the Liouvillian operator can be written as

$$\mathcal{L}(\rho)=-i\left[H,\rho\right]+\gamma\sum_{i}\left[L^{\dagger}_{i}\rho L_{i}-\frac{1}{2}\left\{L^{\dagger}_{i}L_{i},\rho\right\}\right]$$ (S.19)

with $\{A,B\}=AB+BA$ and $[A,B]=AB-BA$ being the anticommutator and the commutator, respectively.

II XX chain point of the Hamiltonian

At the XX chain point i.e when $\eta=i\pi/2$, we end up with a Hamiltonian of a XX chain with an XX impurity [71]

$$H=\sum_{i=1}^{L-1}\frac{1}{2}\left(\sigma^{x}_{i}\sigma^{x}_{i+1}+\sigma^{y}_{i}\sigma^{y}_{i+1}\right)+\frac{J}{2}\left(\sigma^{x}_{0}\sigma^{x}_{1}+\sigma^{y}_{0}\sigma^{y}_{1}\right).$$ (S.20)

After performing the Jordan-Wigner transformation, one can write the model in terms of Fermionic operators i.e. the Hamiltonian becomes

$$H=\sum_{i=0}^{L-1}\left(c^{\dagger}_{i}c_{i+1}+c^{\dagger}_{i+1}c_{i}\right)+J\left(c^{\dagger}_{0}c_{1}+c^{\dagger}_{1}c_{0}\right),$$ (S.21)

and the Liouvillian operator takes the form

$$\mathcal{L}(\rho)=-i\left[H,\rho\right]+4\gamma\sum_{i=1}^{L-1}\left(c^{\dagger}_{i}c_{i}\rho c^{\dagger}_{i}c_{i}-\frac{1}{2}\{\rho,c^{\dagger}_{i}c_{i}\}\right).$$ (S.22)

The density matrix $\rho$ is an operator defined in the $2^{N}$ dimensional Hilbert space which can be written as $\rho=\sum_{m,n}\rho_{mn}\ket{m}\bra{n}$. We can rewrite it as a vector in a $4^{N}$ dimensional Hilbert space, so-called the Liouville-Fock space, as $\ket{\rho}=\sum_{m,n}\rho_{mn}\ket{n,m}\rangle$ without losing information [85]. Then the Liouvillian operator acts as an operator in the enlarged vector space in the following way

$$\mathcal{L}\ket{\rho}\rangle=\frac{d}{dt}\ket{\rho}\rangle,$$ (S.23)

with the Liouvillian as superoperator

$$\mathcal{L}=-iH+i\tilde{H}+4\gamma\sum_{i=1}^{L-1}\left(c^{\dagger}_{i}c_{i}\tilde{c}^{\dagger}_{i}\tilde{c}_{i}-\frac{1}{2}c^{\dagger}_{i}c_{i}-\frac{1}{2}\tilde{c}^{\dagger}_{i}\tilde{c}_{i}\right).$$ (S.24)

Here we introduce a new superoperator acting on the Liouville-Fock space in the following way:

$$\begin{split}c^{\dagger}_{i}\ket{n,m}\rangle&=c^{\dagger}_{i}\ket{n}\bra{m}\\ \tilde{c}^{\dagger}_{i}\ket{n,m}\rangle&=\ket{n}\bra{m}c_{i}\end{split}$$ (S.25)

If we rename $c^{\dagger}_{i}=c^{\dagger}_{i,\uparrow}$ and $\tilde{c}^{\dagger}_{i}=c^{\dagger}_{i,\downarrow}$, and perform the unitary transformation

$$U=\prod_{i=1}^{N/2}(1-2\tilde{c}^{\dagger}_{2i}\tilde{c}_{2i}),$$ (S.26)

the Liouvillian operator becomes the Hubbard model with imaginary coupling [81]

$$\mathcal{L}=-iJ(c^{\dagger}_{0,\sigma}c_{1,\sigma}+c^{\dagger}_{1,\sigma}c_{0,\sigma})-i\sum_{i=1}^{L-1}(c^{\dagger}_{i,\sigma}c_{i+1,\sigma}+c^{\dagger}_{i+1,\sigma}c_{i,\sigma})+\gamma\sum_{i=1}^{L}(2c^{\dagger}_{i,\uparrow}c_{i,\uparrow}-1)(2c^{\dagger}_{i,\downarrow}c_{i,\downarrow}-1)-\gamma(L-1),$$ (S.27)

where $\{c^{\dagger}_{i,\sigma},c_{j,\sigma^{\prime}}\}=\delta_{ij}\delta_{\sigma\sigma^{\prime}}$, and all other combinations anti-commute. In this formalism, the Lindbladian-master equation becomes an eigenvalue problem

$$\mathcal{L}\ket{\rho_{n}}\rangle=E_{n}\ket{\rho_{n}}\rangle.$$ (S.28)

One thing worthy to mention is that the dynamics of 2$k$-point correlation function $G^{x_{1},...,x_{k}}_{y_{1},...,y_{k}}(t)=\text{Tr}\left\{c^{\dagger}_{x_{1}}...c^{\dagger}_{x_{k}}c_{y_{1}}...c_{y_{k}}\rho(t)\right\}$, under this mapping, can be shown to be governed by the 2$k$-particle Schrodinger equation [81]. Particularly for the two-point correlation function we are interested in,

$$i\frac{d}{dt}G(x,y,t)=\left[J_{x}\Delta_{x}+J_{y}\Delta_{y}+4i\gamma\delta_{x,y}+2i\gamma(\delta_{x}+\delta_{y}-2)\right]G(x,y,t),$$ (S.29)

where we introduced a compact notation $J_{x}\mathord{=}J$ when $x\mathord{=}0$, $J_{x}\mathord{=}1$ otherwise and $\Delta f(x)=f(x-1)+f(x+1)$.

III Noiseless XX-impurity Quench

Let us first discuss the dynamics of the model in the absence of the noise. We consider the quench dynamics of the free model governed by a quantum mechanical Hamiltonian

$$h=\sum_{j=0}^{L-1}J_{j}\Delta_{j},$$ (S.30)

where $J_{x}\mathord{=}J$ when $x\mathord{=}0$, $J_{x}\mathord{=}1$ other wise and $\Delta_{j}f(j)=f(j+1)+f(j-1)$ is the discrete derivative. This first quantized Hamiltonian has eigenfunction

$$F_{k}(j)=\frac{1}{J_{j}}(Ae^{ikj}+Be^{-ikj}).$$ (S.31)

The boundary condition when $x=0$ gives  [71]

$$S_{b}(k)=\frac{B}{A}=-\frac{(J^{2}-1)e^{ik}-e^{-ik}}{(J^{2}-1)e^{-ik}-e^{ik}}$$ (S.32)

and $x=L$ gives the quantization condition of $k$

$$\frac{(J^{2}-1)e^{ik}-e^{-ik}}{(J^{2}-1)e^{-ik}-e^{ik}}=e^{i2(L+1)k}.$$ (S.33)

There exists a bound mode with $\Im k>0$, described by $B=0$ such that the wave function is normalizable in the thermodynamic limit, i.e. $(J^{2}-1)e^{ik}-e^{-ik}=0$ which gives $k_{BM}=\frac{i}{2}\log(J^{2}-1)$ with bound mode energy $E=2\cos k=\frac{J^{2}}{\sqrt{J^{2}-1}}$ and $\Im k>0$ gives $J>\sqrt{2}$. Likewise, $A=0$ gives $k_{BM}=-\frac{i}{2}\log(J^{2}-1)$. After studying the spectrum and equilibrium physics, now we can study the quench dynamics of this problem. Since the model is non-interacting, we need to solve only the single-mode evolution and take the algebraic sum of each mode at the end of the calculation. We prepare the initial state as $F(j,t=0)=F_{0}(j)=\delta_{j,j_{0}}$, and we solve for $F(j,t)$. The dynamics of a wave function is given by the following equation

$$\begin{split}F(x,t)&=\sum_{\ell,k}F^{*}_{k}(\ell)F_{0}(\ell)F_{k}(j)e^{i\epsilon_{k}t}\\ &=\sum_{k}F^{*}_{k}(j_{0})F_{k}(j)e^{i\epsilon_{k}t}.\end{split}$$ (S.34)

Let us look at the simplest case where we prepare the fermion (spin flip) at the impurity site initial $j_{0}=0$ and ask how the fermion number (magnetization) changes as time evolves, In the Kondo phase $J<\sqrt{2}$

$$n_{0}(t)=F(0,t)^{2}=\left(\int\frac{dk}{\pi}\frac{J^{2}\sin^{2}(k)}{J^{4}-4(J^{2}-1)\cos^{2}(k)}e^{-i2t\cos k}\right)^{2},$$ (S.35)

with $|F_{k}(0)|^{2}=\frac{1}{J^{2}}(A+B)^{2}=\frac{J^{2}\sin^{2}(k)}{J^{4}-4(J^{2}-1)\cos^{2}(k)}$. Changing the variable $\epsilon=2\cos k$, we obtain

$$F_{Kondo}(0,t)=\frac{1}{2J^{2}\pi}\int_{-2}^{2}d\epsilon\frac{\sqrt{4-\epsilon^{2}}}{1-\mathrm{Sign}(J^{2}-1)t_{K}^{2}\epsilon^{2}}e^{-i\epsilon t},$$ (S.36)

where the Kondo time $t_{K}=\frac{\sqrt{|J^{2}-1|}}{J^{2}}$ is the only scale in this integral. We verify this result numerically as shown in Fig.S.2. In the Bound Mode Phase, the same quantity can be expressed as the following

$$F_{BM}(0,t)=\frac{1}{2J^{2}\pi}\int_{-2}^{2}d\epsilon\frac{\sqrt{4-\epsilon^{2}}}{1-\mathrm{Sign}(J^{2}-1)t_{K}^{2}\epsilon^{2}}e^{-i\epsilon t}+\frac{J^{2}-2}{J^{2}-1}\cos(E_{b}t).$$ (S.37)

The last cosine term is from the bound mode with energy $E_{b}=\pm\frac{J^{2}}{\sqrt{J^{2}-1}}$.

Figure S.2: Dynamics of the spin at the impurity site with initial Néel state. The system is in bound mode phase $J=3$ (left panel) where the analytic curve of the time-dependent spin profile is given by Eq. (S.40) and numerical cure is calculated for $L=40$ with exact diagonalization. The system is in the Kondo phase with $J=0.5$ (right panel) where the analytic curve of time-dependent spin profile is given by Eq. (S.40) and numerical cure is calculated for $L=40$ with exact diagonalization.

In the Kondo phase, the wavefunction at the $\ell$ spin, if we flip $j$ initially, is given by

$$F_{K}(j,\ell,t)=\int_{0}^{\pi}\frac{dk}{2\pi}\frac{1}{J_{\ell}J_{j}}(e^{-ikj}+S_{b}(k)e^{ikj})(e^{ik\ell}+S_{b}(k)^{*}e^{-ik\ell})e^{-i2t\cos k}.$$ (S.38)

In the bound mode phase,

$$F_{B}(j,\ell,t)=F_{K}(j,\ell,t)+\frac{J^{2}}{J_{j}J_{\ell}}\frac{J^{2}-2}{J^{2}-1}\left[\frac{e^{-iE_{B}t}}{(\sqrt{J^{2}-1})^{j+\ell}}+\frac{e^{iE_{B}t}}{(-\sqrt{J^{2}-1})^{j+\ell}}\right].$$ (S.39)

The magnetization at site $\ell$ with $\{j_{x}\}$ spin initially flipped is given by

$$S^{z}(\ell,t)=2\sum_{j\in\{j_{x}\}}|F(j,\ell,t)|^{2}-1.$$ (S.40)

with $F=F_{K}$ or $F=F_{B}$ depending on which phase the system is in.

IV Liouvillian Spectrum

In the previous section we solved XX-impurity Model quenched from an initial state with only some spins flipped, i.e. $\ket{\psi_{0}}=\prod_{j_{x}}\sigma^{+}_{j_{x}}|\downarrow,\downarrow,...\downarrow\rangle$. Now we want to ask the question of how noise affects this model. Using Eq. (S.29), we need to study the two-particle sector of this model to solve for the dynamics of the two-point correlation function.

In the Liouville-Fock space, our initial state (initial condition of the correlation matrix) with spin-flip at site $j_{x}$ corresponds to the initial density matrix $\ket{\rho_{0}}\rangle=c^{\dagger}_{j_{x}\uparrow}c^{\dagger}_{j_{x}\downarrow}\ket{0}\rangle$. If we initially flip more spins, we just need to superpose them, i.e. $\ket{\rho_{0}}\rangle=\sum_{x_{j}}c^{\dagger}_{x_{j}\uparrow}c^{\dagger}_{x_{j}\downarrow}\ket{0}\rangle$. Due to the $U(1)\times U(1)$ symmetry of the Liouvillian operator, the number of fermions and total magnetization are conserved separately. Thus, we only need to diagonalize the eigenvalue problem in Eq. (S.28) in two particles $N=2$ and one spin-flip $M=1$ sector.

The time-independent Schrödinger equation read from Eq.(S.29) is

$$EG(x,y)=\left[J_{x}\Delta_{x}+J_{y}\Delta_{y}+4i\gamma\delta_{x,y}+2i\gamma(\delta_{x}+\delta_{y}-2)\right]G(x,y),$$ (S.41)

with eigenvector

$$G_{k,q}(x,y)=\frac{1}{\mathcal{N}}\mathcal{S}\sum_{\alpha_{1},\alpha_{2}=\pm 1}e^{i\alpha_{1}kx+i\alpha_{2}qy}\left[A(\alpha_{1}k,\alpha_{2}q)\theta(x-y)+B(\alpha_{2}q,\alpha_{1}k)\theta(y-x)\right],$$ (S.42)

where $\mathcal{N}$ is the normalization factor such that $1=\sum_{x,y}|G_{k,q}(x,y)|^{2}$ and $\mathcal{S}$ is the symmetrizer. Note if the transformation in Eq.(S.26) is not performed, then there is a factor $(-1)^{y}$ multiplying the wavefunction shown in the main text. The amplitudes are related as follows:

$$\begin{split}A(-k,q)&=K(k)A(k,q),\\ A(k,-q)&=K(q)A(k,q),\\ A(k,q)&=S(k,q)B(k,q),\end{split}$$ (S.43)

where $K(k)$ is given by the boundary condition $EF(0,x)=JF(1,x)+2i\gamma F(0,x)$ which gives

$$K(k)=-\frac{e^{ik}+e^{-ik}-2i\gamma-J^{2}e^{ik}}{e^{ik}+e^{-ik}-2i\gamma-J^{2}e^{-ik}},$$ (S.44)

and the bulk interaction scattering matrix is that of Hubbard models given by

$$S(k,q)=\frac{\sin k-\sin q+2\gamma}{\sin k-\sin q-2\gamma}.$$ (S.45)

If we impose the boundary condition on a finite chain, then $k$ and $q$ are eigenstates if they satisfy the Bethe Ansatz equations

$$e^{2ik(L+1)} =K(k)S(k,q)S(-k,q)$$ (S.46)

$$e^{2iq(L+1)} =K(q)S(q,k)S(-q,k),$$ (S.47)

which when written explicitly takes the form

$$\begin{split}-\frac{2\cos k-2i\gamma-J^{2}e^{ik}}{2\cos k-2i\gamma-J^{2}e^{-ik}}\left(\frac{\sin k-\sin q+2\gamma}{\sin k-\sin q-2\gamma}\right)\left(\frac{\sin k+\sin q-2\gamma}{\sin k+\sin q+2\gamma}\right)&=e^{i2(L+1)k}\\ -\frac{2\cos q-2i\gamma-J^{2}e^{iq}}{2\cos q-2i\gamma-J^{2}e^{-iq}}\left(\frac{\sin k-\sin q-2\gamma}{\sin k-\sin q+2\gamma}\right)\left(\frac{\sin k+\sin q-2\gamma}{\sin k+\sin q+2\gamma}\right)&=e^{i2(L+1)q}.\end{split}$$ (S.48)

The solution of $k$ and $q$ can be written as $k=k^{R}+ik^{I}$ and $q=q^{R}+iq^{I}$.

IV.1 1. The free band (Roots with vanishing Imaginary part)

By taking the norm of Eq. (S.48), it is easy to see that

$$k^{I}=-\frac{1}{2(L+1)}\ln\lvert S_{b}(k)S(k,q)S(-k,q)\rvert,$$ (S.49)

and the same for $k\to q$. One can see that as long as the scattering matrix is not singular, the imaginary part of the Bethe roots vanishes in the thermodynamic limit, and results in a real band of quasimomenta. In the thermodynamic limit, the imaginary part of the roots is given by replacing $k=\frac{\pi n_{1}}{L+1}+ik^{I}$ and $q=\frac{\pi n_{2}}{L+1}+iq^{I}$.

IV.2 Roots with Finite Imaginary Part

As mentioned earlier, when the scattering matrix is singular, the Bethe roots have a non-vanishing imaginary part in the thermodynamic limit. Now let us consider different situations with different parts of the S-matrix being singular.

Figure S.3: The upper two panels show the Liouvillian spectrum (top-left) and amplitude for the wavefunction in the double-bound mode (top-right) for $J=3$, $\gamma=0.5$ and $L=30$. The perturbative calculation matches the exact diagonalization result better since $\gamma$ is small. The lower two panels show the Liouvillian spectrum (bottom-left) and amplitude for the wavefunction in the double-bound mode (bottom-right) for $J=3$, $\gamma=2$ and $L=30$. The two-site Hamiltonian calculation matches the exact diagonalization result better since $\gamma$ is big. Also, wavefunctions for the both cases show that the particles are localized at the boundary more strongly and when $\gamma$ is bigger.

IV.3 2. Boundary Mode Band

In the case where the boundary $S-$matrix $K$ is singular, the bound mode exists and the singular term can be written as

$$2\cos k-2i\gamma-J^{2}e^{ik}=0,$$ (S.50)

which leads to the constraint $\Im k>0$. Changing the variable to $e^{ik}=x$, the equation becomes

$$(1-J^{2})x^{2}-2i\gamma x+1=0,$$ (S.51)

which has a solution

$$e^{ik}=\frac{\pm\sqrt{J^{2}-1-\gamma^{2}}-i\gamma}{J^{2}-1}:=e^{-\phi},$$ (S.52)

under the assumption $|e^{ik}|>1$. Plugging in $k$ into the other equation gives the solution of $q$ as

$$-\frac{2\cos q-2i\gamma-J^{2}e^{iq}}{2\cos q-2i\gamma-J^{2}e^{-iq}}\left(\frac{\sin k-\sin q-2\gamma}{\sin k-\sin q+2\gamma}\right)\left(\frac{\sin k+\sin q-2\gamma}{\sin k+\sin q+2\gamma}\right)=e^{i2(L+1)q}.$$ (S.53)

In the thermodynamic limit, $q=\frac{n_{2}\pi}{L+1}$ and $q=i\phi$.

IV.4 3. Diffusive Band

In the $L\to\infty$ limit, without lack of generality, we choose $\Im k_{1}<0$. Then the bulk scattering matrix being singular indicates the solutions satisfy the string hypothesis in the way $k_{2}=k_{1}^{*}+\pi$. These solutions form a band of solutions whose corresponding eigenvalues have a vanishing real part. Now we define $k=\Re k+i\Im k$, then the S-matrix being singular also gives the relationship between $\Re k$ and $\Im k$ as

$$\Re k=\pi-\arcsin\frac{\gamma}{\cosh\Im k}.$$ (S.54)

As a result, the eigenenergy is $E=2i\tanh\Im k-4i\gamma$, and the Bethe Ansatz equation can be written as the form of a Bethe-Gaudin-Takahashi (BGT) equation [82] as following

$$\left(\frac{2\cos\Re k\cosh\Im k-i2\sin\Re k\sinh\Im k-2i\gamma-J^{2}e^{i(\Re k-\Im k)}}{2\cos\Re k\cosh\Im k-i2\sin\Re k\sinh\Im k-2i\gamma-J^{2}e^{-i(\Re k-\Im k)}}\right)^{2}\left(\frac{2i\sinh\Im k\cos\Re k+2\gamma}{2i\sinh\Im k\cos\Re k-2\gamma}\right)^{2}=e^{-i4(L+1)\Re k}.$$ (S.55)

To solve this equation, it is easier to change the variable $\Im k=\mathrm{arccosh}\frac{\gamma}{z}$ and take the logarithm on both sides of the equation which gives

$$2\arcsin(z_{n})+\mathcal{O}(L^{-1})=-\pi\frac{n}{L}$$ (S.56)

with $n=0,1,2,...,L-2$. This equation which look extremely simple gives the spectrum $E_{n}=4i\sqrt{\gamma^{2}-\sin\left(\frac{n\pi}{2L}\right)^{2}}-4i\gamma+\mathcal{O}(L^{3})$, which match the numerical solution extremely well as shown in Fig. S.3. This figure shows that the analytical calculation fits the exact diagonalization reasonably well when $\gamma>1$ which the system is experiencing a true diffusive bulk. However, when $\gamma<1$, the diffusive band mixes with the free band and results in a worse but reasonable prediction because we are not taking care about the contribution from the impurity.

Noticing $\Im k<0$ gives the constraint $\Im E>-4\gamma$, and all solutions beyond this limit are not physical and should be dropped. Since $\cosh\Im k>1$, the largest imaginary energy is when $\Im k=0$, which gives $\Re k=\pi-\arcsin\gamma$ and thus resulting the energy $E_{max}=4i\sqrt{\gamma^{2}-1}-4i\gamma$. This is the lower bound of the dissipation band, and when $|E_{max}|>4\gamma$ which can only happen when $\gamma<1$ the dissipation band mixes with the free band.

IV.5 4. Double-Boundary Modes

As discussed in the main text, there is another solution when both the boundary S-matrix and the bulk S-matrix are singular. Such a mode physically describes two particles (in the imaginary Hubbard model) that are both exponentially localized at the of the chain and they experience a strong imaginary Hubbard interaction. The existence of these two states is discussed in the main text with a phase diagram in Fig.1(d). However, we are not able to separate these two singularities, so we treat them in two limits. When $\gamma<1$, the double-bound mode is perturbatively connected to the case when both particles form a noiseless bound mode, which means that

$$G_{k_{0},k_{0}}(x,y)=\mathcal{N}\frac{1}{J_{x}}e^{-ik_{0}x}\frac{1}{J_{y}}e^{-ik_{0}y}$$ (S.57)

with $k_{0}=\frac{i}{2}\log(J^{2}-1)+\pi$ or $k_{0}=\frac{i}{2}\log(J^{2}-1)$ and $\mathcal{N}=\frac{1}{J^{2}}+\frac{e^{-2ik_{0}}}{1-e^{-2ik_{0}}}$ being the normalization factor. Then the eigenvalue of this double-bound mode reads

$$E=4\cos k_{0}-4i\gamma+4i\gamma\left[\sum_{x=1}^{L}G_{k_{0},k_{0}}(x,x)^{2}+G_{k_{0},k_{0}}(0,0)\right]+\mathcal{O}(\gamma^{2}),$$ (S.58)

and when $\gamma>1$, the bound mode is more localized on the boundary (decays faster) as one can see from the single boundary bound mode. Then we can treat it by solving a two-site Hamiltonian ($L=2$) as

$$h=\begin{pmatrix}0&-J&-J&0\\ -J&-2i\gamma&0&-J\\ -J&0&-2i\gamma&-J\\ 0&-J&-J&0\end{pmatrix},$$ (S.59)

which has two boundary bound mode with energy $E_{\pm}=-i\gamma\pm\sqrt{4J^{2}-\gamma^{2}}$. It is a double-bound mode only when $\Re E\neq 0$ which gives the boundary between the Kondo phase and the bound mode phase at the large limit $\gamma$ as $J=\gamma/2$.

V N-particle Bethe Ansatz

As mentioned earlier, the XX chain with boundary impurity maps to the open boundary condition Hubbard model with boundary impurity given by

$$H=\sum_{j=1}^{L-1}\sum_{\sigma=\uparrow,\downarrow}(c^{\dagger}_{j,\sigma}c_{j+1,\sigma}+c^{\dagger}_{j+1,\sigma}c_{j,\sigma})+J(c^{\dagger}_{0,\sigma}c_{1,\sigma}+c^{\dagger}_{1,\sigma}c_{0,\sigma})+4i\gamma\sum_{j=1}^{L}n_{j,\uparrow}n_{j,\downarrow}-2i\gamma\sum_{j=1}^{L}\sum_{\sigma=\uparrow,\downarrow}n_{j,\sigma}.$$ (S.60)

As shown in [71], the boundary scattering matrix is

$$S_{b_{j}}(k)=\frac{-e^{2ik_{j}}J^{2}+e^{2ik_{j}}+1}{-J^{2}+\left(1+e^{2ik_{j}}\right)}=e^{2ik_{J}}\frac{2\cos(k_{j})-iJ^{2}\sin(k_{j})-J^{2}\cos(k_{j})}{2\cos(k_{j})+iJ^{2}\sin(k_{j})-J^{2}\cos(k_{j})},$$ (S.61)

and the Bulk S-matrix, as first calculated by Lieb and Wu for real values of U, is given by

$$S_{j,l}=\frac{\sin(k_{j})-\sin(k_{l})-2\gamma P_{j,l}}{\sin(k_{j})-\sin(k_{l})-2\gamma},$$ (S.62)

where $P_{j,l}$ is the permutation operator.

Now, the usual Bethe Ansatz routine can be employed. It is very important to notice that the boundary $S-$matrix does not depend on the spin variable but the bulk $S-$matrix depends on the spin variable (through the permutation operator). The Yang-Baxter equation is therefore trivially satisfied, thereby proving that the model is integrable.

Following [87, 66], we construct the N-particle Bethe Ansatz of the model using the functional Bethe Ansatz method. Let us consider the rational solution of the six-vertex model

$$R_{ij}(u)=u+\eta P_{ij},$$ (S.63)

where we will consider the crossing parameter $\eta=-2\gamma$.

Then for the open boundary condition, we consider the double-row transfer matrix

$$T_{0}(\lambda) =R_{0,N}(\lambda-\sin(k_{N}))R_{0,N-1}(\lambda-\sin(k_{N_{1}}))\cdots R_{0,2}(\lambda-\sin(k_{2}))R_{0,1}(\lambda-\sin(k_{1}))$$

$$\hat{T}_{0}(\lambda) =R_{0,1}(\lambda+\sin(k_{1}))R_{0,2}(\lambda+\sin(k_{2}))\cdots R_{0,N-1}(\lambda+\sin(k_{N_{1}}))R_{0,N}(\lambda+\sin(k_{N}))$$

Now, we define the monodromy matrix

$$\Xi(\lambda)=T_{0}(\lambda)\hat{T}(\lambda).$$ (S.64)

The trace of the monodromy matrix over the auxiliary space is defined as the double-row transfer matrix

$$t(\lambda)=\operatorname{tr}_{0}\Xi(\lambda).$$ (S.65)

It is quite easy to see that the transfer matrix forms a one-parameter family of commuting operators i.e.

$$[t(\lambda),t(\rho)]=0.$$ (S.66)

First notice that using coordinate Bethe Ansatz, one can immediately construct the eigenvalue problem for the transfer matrix

$$\tau_{j}= S_{j-1,j}\left(k_{j-1},k_{j}\right)\cdots S_{1,j}\left(k_{1},k_{j}\right){S_{b_{j}}}\left(k_{j}\right)S_{j,1}\left(-k_{j},k_{1}\right)S_{j,j-1}\left(-k_{j},k_{j-1}\right)$$ (S.67)

$$S_{j,j+1}\left(-k_{j},k_{j+1}\right)\cdots S_{j,N}\left(-k_{j},k_{N}\right)S_{N,j}\left(k_{N},k_{j}\right)\cdots S_{j+1,j}\left(k_{j+1},k_{j}\right).$$

Upon identification of the relation between the spectral parameter valued double-row transfer matrix $t(\lambda)$ and the transfer matrix obtained from the coordinate Bethe Ansatz $\tau_{j}$ via

$$\tau_{j}/S_{b_{j}}(k_{j})=\frac{t(-\sin(k_{j}))}{\eta(\sin(k_{j})-\eta)\prod_{l\neq j}^{N}(\sin(k_{j})-\sin(k_{l})-\eta)(\sin(k_{j})+\sin(k_{l})-\eta)}$$ (S.68)

Once again, the boundary $S-$matrix $S_{b_{j}}$ is independent of a spin variable; it is just a phase that trivially commutes with bulk $S-$matrices.

The eigenvalue of the transfer matrix satisfies the $T-Q$ relation

$$\Lambda(\lambda) =\prod_{j=1}^{N}\frac{2\lambda+2\eta}{2\lambda+\eta}\left(\lambda-\sin k_{j}+\eta\right)\left(\lambda+\sin k_{j}+\eta\right)\frac{Q(\lambda-\eta)}{Q(\lambda)}$$

$$\quad+\prod_{j=1}^{N}\frac{2\lambda}{2\lambda+\eta}\left(\lambda-\sin k_{j}\right)\left(\lambda+\sin k_{j}\right)\frac{Q(\lambda+\eta)}{Q(\lambda)},$$ (S.69)

where the $Q$ function is

$$Q(\lambda)=\prod_{l=1}^{M}(\lambda-\lambda_{l})(\lambda+\lambda_{l}+\eta).$$ (S.70)

Imposing the regularity condition of the $T-Q$ relation immediately gives the Bethe equation

$$\frac{\lambda_{j}+\eta}{\lambda_{j}}\prod_{j=1}^{N}\frac{\lambda_{j}-\sin(k_{j})+\eta}{\lambda_{j}-\sin(k_{j})}\frac{\lambda_{j}+\sin(k_{j})+\eta}{\lambda_{j}+\sin(k_{j})}=-\prod_{l=1}^{M}\frac{\lambda_{j}-\lambda_{l}+\eta}{\lambda_{j}-\lambda_{l}-\eta}\frac{\lambda_{j}+\lambda_{l}+2\eta}{\lambda_{j}+\lambda_{l}}.$$ (S.71)

To make the equations more symmetric, we change the variable $\lambda\to\lambda-\eta/2$, such that the Bethe equations becomes

$$\frac{\lambda_{j}+\frac{\eta}{2}}{\lambda_{j}-\frac{\eta}{2}}\prod_{j=1}^{N}\frac{\lambda_{j}-\sin(k_{j})+\frac{\eta}{2}}{\lambda_{j}-\sin(k_{j})-\frac{\eta}{2}}\frac{\lambda_{j}+\sin(k_{j})+\frac{\eta}{2}}{\lambda_{j}+\sin(k_{j})-\frac{\eta}{2}}=-\prod_{l=1}^{M}\frac{\lambda_{j}-\lambda_{l}+\eta}{\lambda_{j}-\lambda_{l}-\eta}\frac{\lambda_{j}+\lambda_{l}+\eta}{\lambda_{j}+\lambda_{l}-\eta}.$$ (S.72)

From Eq. (S.68), upon using the shift $\lambda\to\lambda-\eta/2$ we obtain

$$e^{-2ik_{j}L}/S_{b_{j}}(k_{j})=\prod_{l=1}^{M}\frac{\sin k_{j}-\lambda_{l}+\frac{\eta}{2}}{\sin k_{j}-\lambda_{l}-\frac{\eta}{2}}\frac{\sin k_{j}+\lambda_{l}+\frac{\eta}{2}}{\sin k_{j}+\lambda_{l}-\frac{\eta}{2}}.$$ (S.73)

Or,

$$e^{-2ik_{J}(L+1)}\frac{2\cos(k_{j})+iJ^{2}\sin(k_{j})-J^{2}\cos(k_{j})}{2\cos(k_{j})-iJ^{2}\sin(k_{j})-J^{2}\cos(k_{j})}=\prod_{l=1}^{M}\frac{\sin k_{j}-\lambda_{l}+\frac{\eta}{2}}{\sin k_{j}-\lambda_{l}-\frac{\eta}{2}}\frac{\sin k_{j}+\lambda_{l}+\frac{\eta}{2}}{\sin k_{j}+\lambda_{l}-\frac{\eta}{2}}.$$ (S.74)

Putting $\eta=-2\gamma$, we obtain

$$e^{2ik_{j}(L+1)}\frac{2\cos(k_{j})-iJ^{2}\sin(k_{j})-J^{2}\cos(k_{j})}{2\cos(k_{j})+iJ^{2}\sin(k_{j})-J^{2}\cos(k_{j})}=\prod_{l=1}^{M}\frac{\sin k_{j}-\lambda_{l}+\gamma}{\sin k_{j}-\lambda_{l}-\gamma}\frac{\sin k_{j}+\lambda_{l}+\gamma}{\sin k_{j}+\lambda_{l}-\gamma},$$ (S.75)

and the rapidities satisfy

$$\frac{\lambda_{j}+\gamma}{\lambda_{j}-\gamma}\prod_{j=1}^{N}\frac{\lambda_{j}-\sin(k_{j})+\gamma}{\lambda_{j}-\sin(k_{j})-\gamma}\frac{\lambda_{j}+\sin(k_{j})+\gamma}{\lambda_{j}+\sin(k_{j})-\gamma}=-\prod_{l=1}^{M}\frac{\lambda_{j}-\lambda_{l}+2\gamma}{\lambda_{j}-\lambda_{l}-2\gamma}\frac{\lambda_{j}+\lambda_{l}+2\gamma}{\lambda_{j}+\lambda_{l}-2\gamma}.$$ (S.76)

If one is interested in the many-point correlation function, one needs to solve the many-particle BAE given by Eq. (S.75) and Eq. (S.76). In the main text, we only considered the two-point correlation function, and we leave the general problem of the 2n-point correlation function for future work. The Bethe equations for the real Hubbard interaction were first obtained in  [84].

Figure S.4: The magnetization $S^{z}_{j}(t)$ for time $t$ and site $j$ in the XXX case with system size $L=10$ (right panel) and XX case with system size $L=40$ (left panel) for dephasing rate $\gamma=3$.

VI Numerical details

Here we will show the detailed calculation to get the diffusion constant $D$ numerically for the Heisenberg (XXX) and XX case. We start from a Doman wall configuration initially with $\rho_{D}=\sum_{z\in\mathcal{S}_{D}}\sigma^{+}_{z}\ket{0}\bra{0}\sigma^{+}_{z}$ and $\mathcal{S}_{D}=0,1,2,...,L/2$. Then we evolve the state $\rho(t)$ via the Lindbladian master equation (S.19) from the initial state $\rho(0)=\rho_{D}$. Then we calculated the magnetization profile $S^{z}_{j}(t)$ for each time $t$ and site $j$. The results for $\gamma=3$ in the XX case and XXX case are shown in Fig. S.4. For the XXX case, the domain wall hit the boundary when $t\sim 50$ and resulting in a finite size effect.

If the bulk is diffusive, then in the limit lattice spacing $a\to 0$ and system size $L\to\infty$ $S^{z}_{j}(t)\to S^{z}(x,t)$ should satisfy the diffusion equation $(\partial_{t}-D\partial_{x}^{2})S^{z}(x,t)=0$. The initial condition is a domain wall meaning $S^{z}(x,0)=\mathrm{Sign}(L/2-x)$ and the solution reads $S^{z}(x,t)=-\mathrm{Erf}(\frac{x}{\sqrt{Dt}})$. This is obvious that for different times and spaces, the function collapses on the curve with a normalized variable $\xi=x/\sqrt{t}$ in the way that $S^{z}(\xi)=-\mathrm{Erf}(\frac{\xi}{\sqrt{D}})$. The system is not perfectly diffusive for early-time dynamics and due to the finite size effect, we cannot evolve the system for an extremely long time. Therefore, we will fit the magnetization $S^{z}(\xi)$ to $-\mathrm{Erf}(\xi/\sqrt{D(t)})$ for a finite time $t$ to get the diffusion constant $D(t)$. To get the long-time behavior, we expand $D(t)=D(\infty)+\frac{D_{1}}{t}+\mathcal{O}(\frac{1}{t^{2}})$ and scale the diffusion constant when $t\to\infty$. The fitting and scaling are shown in Fig.S.5. We found for both XX and XXX case $D_{0}$ which is the diffusion constant for long-time dynamics is $D_{0}=\frac{1}{2\gamma}$ for big gamma. As for small $\gamma$, due to the finite size effect, we are not able to conclude for the XXX case. For the XX case, where we can go for larger system sizes, we can show $D_{0}=\frac{1}{2\gamma}$.

Figure S.5: Scaling of the diffusion constant $D$ with time $t$ and $\gamma$ for the XXX circuit. The left panel shows $D$ as a fit parameter obtained by fitting the spin profile to the $S^{z}(\xi)=\mathrm{Erf}(\xi/\sqrt{D(t)})$. The right panel shows $D(t)$ as a function for different times and we scale $D(\infty)$ by taking $1/t\to 0$. The dashed line is the fitting for $1/t\to 0$ and different colors for different $\gamma$.