Spin-1/2 XXZ chain coupled to two Lindblad baths: Constructing nonequilibrium steady states from equilibrium correlation functions

TL;DR

This paper shows that for the spin-1/2 XXZ chain, the nonequilibrium steady state under weak driving can be constructed from equilibrium correlation functions, bridging closed and open system approaches.

Contribution

It introduces a method to derive nonequilibrium steady states from equilibrium correlations, enabling analysis of larger systems and clarifying discrepancies in transport coefficient extraction.

Findings

Nonequilibrium steady states can be constructed from equilibrium correlations.

The approach allows studying larger systems than traditional methods.

Potential pitfalls exist in extracting transport coefficients from finite systems.

Abstract

State-of-the-art approaches to extract transport coefficients of many-body quantum systems broadly fall into two categories: (i) they target the linear-response regime in terms of equilibrium correlation functions of the closed system; or (ii) they consider an open-system situation typically modeled by a Lindblad equation, where a nonequilibrium steady state emerges from driving the system at its boundaries. While quantitative agreement between (i) and (ii) has been found for selected model and parameter choices, also disagreement has been pointed out in the literature. Studying magnetization transport in the spin-1/2 XXZ chain, we here demonstrate that at weak driving, the nonequilibrium steady state in an open system, including its buildup in time, can remarkably be constructed just on the basis of correlation functions in the closed system. We numerically illustrate this direct correspondence of closed-system and open-system dynamics, and show that it allows the treatment of comparatively large open systems, usually only accessible to matrix product state simulations. We also point out potential pitfalls when extracting transport coefficients from nonequilibrium steady states in finite systems.

Equations

$$H=J\sum _{r=1}^N(S_r^x{S}_{r+1}^x+S_r^y{S}_{r+1}^y+\Delta S_r^z{S}_{r+1}^z),$$

$$⟨S_r^z(t)S_{r^{\prime }}^z(0){⟩}_{\text{eq}}=\frac{\text{tr}[e^{iHt}{S}_r^z{e}^{-iHt}{S}_{r^{\prime }}^z]}{2^N}.$$

$$\dot{\rho }(t)=\mathcal{L}\rho (t)=i[\rho (t),H]+\mathcal{D}\rho (t),$$

$${\cal D}\,\rho(t)=\sum_{j}\alpha_{j}\Big{(}L_{j}\rho(t)L^{\dagger}_{j}-\frac{1}{2}\{\rho(t),L_{j}^{\dagger}L_{j}\}\Big{)}\ ,$$

$$L_1=S_{{\text{B}}_1}^{+},$$

$$L_2=L_1^{†}=S_{{\text{B}}_1}^{-},$$

$$L_3=S_{{\text{B}}_2}^{+},$$

$$L_4=L_3^{†}=S_{{\text{B}}_2}^{-},$$

$$⟨S_r^z(t)⟩=\text{tr}[\rho (t)S_r^z],$$

$$H_{\text{eff}}=H-i\gamma +i\gamma \mu (n_{{\text{B}}_1}-n_{{\text{B}}_2}),$$

$$H_{\text{eff}}\approx H-i\gamma$$

$$|\psi (t)⟩\approx e^{-\gamma t}{e}^{-iHt}|\psi (0)⟩,$$

$$|{\psi }^{\prime }(t)⟩=\frac{L_j|\psi (t)⟩}{\Vert L_j|\psi (t)⟩\Vert },$$

$$p_j=\frac{{\alpha }_j{\Vert L_j|\psi (t)⟩\Vert }^2}{{\sum }_j{\alpha }_j{\Vert L_j|\psi (t)⟩\Vert }^2}.$$

$$⟨S_r^z(t)⟩\approx \frac{1}{{\text{T}}_{\text{max}}}\sum _{\text{T}=1}^{{\text{T}}_{\text{max}}}\frac{⟨{\psi }_{\text{T}}(t)|S_r^z|{\psi }_{\text{T}}(t)⟩}{{\Vert |{\psi }_{\text{T}}(t)⟩\Vert }^2},$$

$$|\psi (0)⟩\propto \sum _j{c}_j|{\varphi }_j⟩,$$

$$d_{r,\text{T}}(t)=2\mu \sum _j$$

$$A_j$$

$$\text{with}{a}_j$$

$$⟨S_r^z(t)⟩\approx \frac{1}{{\text{T}}_{\text{max}}}\sum _{\text{T}=1}^{{\text{T}}_{\text{max}}}{d}_{r,\text{T}}(t),$$

$$|\psi (0)⟩\to |{\psi }^{\prime }⟩\propto L_1|\psi (0)⟩.$$

$$d_r(t)\equiv \frac{⟨{\psi }^{\prime }(t)|S_r^z|{\psi }^{\prime }(t)⟩}{{\Vert |{\psi }^{\prime }(t)⟩\Vert }^2}\approx ⟨{\psi }^{\prime }|e^{iHt}{S}_r^z{e}^{-iHt}|{\psi }^{\prime }⟩,$$

$$\frac{d_r(t)}{2}\approx ⟨S_r^z(t)S_{{\text{B}}_1}^z(0){⟩}_{\text{eq}}$$

$$\frac{{\bar{d}}_r(t)}{2}=$$

$$+(p_3-p_4)⟨S_r^z(t)S_{{\text{B}}_2}^z(0){⟩}_{\text{eq}}$$

$${\bar{d}}_r(t)=\mu ⟨S_r^z(t)S_{{\text{B}}_1}^z(0){⟩}_{\text{eq}}-\mu ⟨S_r^z(t)S_{{\text{B}}_2}^z(0){⟩}_{\text{eq}}$$

$${\bar{d}}_{{\text{B}}_1}(\tau -0^{+})=-{\bar{d}}_{{\text{B}}_2}(\tau -0^{+})$$

$$y_1$$

$$y_2$$

$$y_3$$

$$y_4$$

$$p_1=p_4=\frac{1}{2}\frac{(1+\mu )y_1}{(1+\mu )y_1+(1-\mu )y_2}$$

$$p_2=p_3=\frac{1}{2}\frac{(1-\mu )y_1}{(1+\mu )y_1+(1-\mu )y_2}$$

$$p_1-p_2=p_4-p_3=\frac{\mu -2{\bar{d}}_{{\text{B}}_1}(\tau -0^{+})}{2-4\mu {\bar{d}}_{{\text{B}}_1}(\tau -0^{+})}.$$

$$A_{\tau }=\frac{(p_1-p_2)-{\bar{d}}_{{\text{B}}_1}(\tau -0^{+})}{\mu },$$

Full text

Introduction. Our understanding of the properties of many-body quantum systems out of equilibrium has seen remarkable advances in the last decades thanks to various experimental and theoretical breakthroughs [1, 2, 3, 4, 5]. Central questions are concerned with the emergence of particular (thermal or nonthermal) steady states in the long-time limit, but also with the (universal) properties of the actual nonequilibrium process towards such states in the course of time [2, 3, 4, 5]. Broadly speaking, these and related questions are usually studied in two different scenarios: (i) the system of interest is perfectly isolated from its environment and evolves unitarily in time; (ii) the system’s time evolution is nonunitary due to an explicit coupling to an external bath which can affect the dynamics (see, e.g., Ref. [6, 7, 8]).

In systems with a global conservation law, a fundamental role is played by transport processes [9]. Quantum transport is also a prime example of a research question that is explored both from a closed-system and an open-system perspective. In closed systems, a widely used approach is linear response theory, where the Kubo formula allows for the extraction of transport coefficients from equilibrium correlation functions, which can be studied in the time or frequency domain and in real or momentum space [9]. While nonintegrable systems are expected to exhibit normal diffusion [10, 11, 12], the concrete calculation of diffusion constants for specific models turns out to be a hard task in practice. This difficulty has been one of the motivations for the development of sophisticated numerical methods [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. Moreover, some classes of models can generically feature anomalous subdiffusion or superdiffusion in certain parameter regimes [28, 29, 30, 2, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41].

In contrast, when studying transport in an open-system setting, the model of interest is often coupled at its edges to two reservoirs, e.g., at different temperatures or chemical potentials, leading to a nonequilibrium steady state in the long-time limit. Then, the profile and currsnt of this steady state yield information on the transport behavior [42, 43, 44, 45]. A popular description of such an open system is provided by the Lindblad quantum master equation [6], not least since it allows for efficient numerical simulations based on matrix product states, giving access to comparatively large system sizes [46, 47, 48, 49, 50, 51, 39]. While quantitative agreement of transport coefficients according to the Lindblad description with those from closed-system approaches has been found for selected models and parameter regimes [52, 53, 54], also disagreement has been pointed out in the literature [55], and there is no proof that both approaches have to agree [56, 55, 57, 58, 9].

From a physical perspective, computed transport coefficients for a given system should of course be independent of the method employed. In fact, some of us have recently shown that the dynamics of closed and open systems can be connected with each other in a certain simple setting. Specifically, Ref. [54] considered an initially homogeneous system coupled locally to a single Lindblad bath, which induces a net magnetization into the system. Remarkably, it was shown that if the Lindblad driving is weak, the flow of the magnetization in the open system, i.e., the broadening of the nonequilibrium density profile, can be described by an appropriate superposition of equilibrium correlation functions in the closed system. Building on this result, we here go beyond Ref. [54] in a crucial point and explore the more common situation of two Lindblad baths inducing a nonequilibrium steady state. Considering magnetization transport in the paradigmatic spin-$1/2$ XXZ chain as an example, we demonstrate that the steady state in the open system can be constructed on the basis of correlation functions in the closed system. We support our analytical results by large-scale numerical simulations and show that our scheme enables efficient unravelings of Lindblad equations for systems with up to $36$ sites, which are usually only accessible with matrix product state techniques.

Closed System. We consider the one-dimensional XXZ model, which is described by the Hamiltonian

[TABLE]

where $S_{r}^{j}$ (${j=x,y,z}$) are spin-1/2 operators at site $r$, ${J>0}$ is the antiferromagnetic coupling constant, and $\Delta$ denotes the anisotropy in the $z$ direction. Moreover, $N$ is the number of sites and we employ periodic boundary conditions, ${S_{N+1}^{j}\equiv S_{1}^{j}}$. The XXZ chain conserves the global magnetization, $[H,\sum_{r}S_{r}^{z}]=0$, and we will particularly focus on the regime ${\Delta>1}$, where it is well-established that spin transport is diffusive [9]. This diffusive transport behavior can, for instance, be seen in the Gaussian shape of the infinite-temperature spin-spin correlation function at $\Delta=1.5$ [59], see Fig. 1(a),

[TABLE]

The root-mean-squared displacement of the above grows as ${\Sigma(t)\propto\sqrt{t}}$, see Fig. 1(b), where ${\Sigma^{2}(t)=\sum_{r}(r-r^{\prime})^{2}C_{rr^{\prime}}(t)-[\sum_{r}(r-r^{\prime})C_{rr^{\prime}}(t)]^{2}}$ and ${C_{rr^{\prime}}(t)=4\langle S_{r}^{z}(t)S_{r^{\prime}}^{z}\rangle_{\text{eq}}}$. Moreover, a diffusion coefficient can be defined as ${2D(t)=\tfrac{d}{dt}\Sigma^{2}(t)}$ [60]. As shown in Fig. 1(b), $D(t)$ takes on a constant value ${D/J\approx 0.6}$ for ${tJ\lesssim 10}$, which is approximately independent of time (and system size [61, 59]) and consistent with other results in the literature [62, 63, 64, 51].

In the following, we will show that the equilibrium correlation function $\langle S_{r}^{z}(t)S_{r^{\prime}}^{z}(0)\rangle_{\text{eq}}$ in Eq. (2) is not only central to transport in the closed system, but can remarkably be used to predict the buildup of a nonequilibrium steady state in an open-system situation where the spin chain is weakly driven by two Lindblad baths. While we focus on the integrable XXZ chain as a concrete example due to its interesting transport properties, we expect our conceptual findings to apply to a wider range of models. In particular, while our derivation [54, 65] is largely model-independent, it implicitely assumes sufficiently fast local equilibration, which should be even better fulfilled in nonintegrable chaotic systems.

Open System. Let us consider a scenario, where the XXZ chain is explicitly coupled to an environment. We describe this setting with a Lindblad equation,

[TABLE]

which consists of a coherent time evolution of the density matrix $\rho$ with respect to $H$ and an incoherent damping term,

[TABLE]

with non-negative rates $\alpha_{j}$, Lindblad operators $L_{j}$, and the anticommutator $\{\bullet,\bullet\}$. While the derivation of this equation can be a subtle task for a given microscopic model [43, 66], it is the most general form of a time-local quantum master equation, which maps any density matrix to a density matrix, i.e., which preserves trace, hermiticity, and positivity [6]. Here, we choose [9]

[TABLE]

where $\gamma$ is the system-bath coupling and $\mu$ is the driving strength. $L_{1}$ and $L_{2}$ are local Lindblad operators at site $\text{B}_{1}$ and flip a spin up and down, respectively. $L_{3}$ and $L_{4}$ act similarly on another site $\text{B}_{2}$. In the following, we set ${\text{B}_{1}=1}$ and ${\text{B}_{2}=N/2+1}$. Note that we still consider periodic boundary conditions. However, our approach can also generally be applied to open boundaries with the two baths at the system’s edges ${\text{B}_{1}=1}$ and ${\text{B}_{2}=N}$, and we present results for this setting in [65]. For ${\mu>0}$, the first (second) bath induces a net polarization of $\mu/2$ ($-\mu/2$), leading to a steady state in the long-time limit with a characteristic density profile and a constant current. Note that, while the Lindblad modeling (3) - (8) is standard in the context of transport in quantum lattice models [9], there exist other approaches to open-system dynamics which can also address potential non-Markovian effects [67].

In addition to the long-time limit, we are interested in the temporal buildup of the steady state. Thus, we study the time evolution of local densities

[TABLE]

which depends on the parameters of the system $H$, but also on the bath parameters $\gamma$ and $\mu$. As an initial state, we here consider a homogeneous situation with ${\rho(0)\propto\mathbb{1}}$ being the infinite-temperature ensemble.

Quantum-trajectory approach. One possibility to solve the Lindblad equation is given by the concept of stochastic unraveling, which relies on pure states $|\psi\rangle$ rather than density matrices [68, 69]. It consists of an alternating sequence of stochastic jumps with one of the Lindblad operators and deterministic evolutions governed by an effective Hamiltonian ${H_{\text{eff}}=H-\tfrac{i}{2}\sum_{j}\alpha_{j}\,L_{j}^{\dagger}L_{j}}$. For our choice of Lindblad operators,

[TABLE]

with ${n_{r}=S_{r}^{+}S_{r}^{-}=S_{r}^{z}+\mathbb{1}/2}$. For weak driving ${\mu\ll 1}$, the time scale on which the last term in Eq. 10 affects the dynamics is much longer than the typical time scale between jumps. Thus, the effective Hamiltonian can be approximated as

[TABLE]

and the time evolution of a pure state reads

[TABLE]

i.e., apart from the scalar damping term, the dynamics is generated by the closed system $H$ only. The approximation in Eq. (12) is one of the main ingredients to establish a correspondence between the dynamics of the isolated and the weakly-driven XXZ chain below. For larger values of $\mu$, the effective Hamiltonian generating the dynamics of $|\psi(t)\rangle$ also involves the two operators $n_{\text{B}_{1}}$ and $n_{\text{B}_{2}}$, cf. Eq. (10).

Naturally, since $H_{\text{eff}}$ is a non-Hermitian operator, the norm of a pure state is not conserved as a function of time. As a consequence, for a given $\varepsilon$ drawn at random from a uniform distribution ${]0,1]}$, there is a time, where the condition ${\left\lVert|\psi(t)\rangle\right\rVert^{2}>\varepsilon}$ is first violated. At this time, a jump with one of the Lindblad operators occurs and the new and normalized pure state reads

[TABLE]

where the specific jump is chosen with probability

[TABLE]

After this jump, the next deterministic evolution takes place. This sequence of stochastic jumps and deterministic evolutions leads to a particular trajectory $|\psi_{\text{T}}(t)\rangle$. The time-dependent density matrix according to the Lindblad equation can eventually be approximated by the average over different trajectories $\mathrm{T}$. Thus, expectation values read

[TABLE]

where ${\text{T}_{\text{max}}}$ is the number of trajectories.

In order to mimic the homogeneous state ${\rho(0)\propto\mathbb{1}}$, we use random pure states as initial condition for the stochastic unraveling,

[TABLE]

where the real and imaginary parts of the coefficients $c_{j}$ in some given basis $|\phi_{j}\rangle$ are drawn at random according to a Gaussian probability distribution with zero mean. Crucially, by exploiting the concept of quantum typicality [70, 71, 72, 73, 74, 59, 75], expectation values ${\bra{\psi}\bullet\ket{\psi}}$ of local observables evaluated within such random states can be related to infinite-temperature averages $\text{tr}[\bullet]/2^{N}$. This is used in the following to connect the equilibrium correlation functions $\langle S_{r}^{z}(t)S_{r^{\prime}}^{z}(0)\rangle_{\text{eq}}$ [Eq. (2)] to the dynamics $\langle S_{r}^{z}(t)\rangle$ in the open system [Eq. (9)].

Constructing steady states from correlation functions. In Ref. [54], it was demonstrated that individual quantum trajectories of the open system can be described by closed-system equilibrium correlation functions if the driving by the Lindblad bath is weak. We here build on this result and apply it to the case of two Lindblad baths leading to a nonequilibrium steady state. While we relegate details of the derivation to the supplemental material [65], we find that for small coupling $\gamma$ and weak driving $\mu$, the local magnetization dynamics within a single trajectory T can be approximated as ${d_{r,\text{T}}(t)\approx\bra{\psi_{\text{T}}(t)}S_{r}^{z}\ket{\psi_{\text{T}}(t)}/\left\lVert|\psi_{\text{T}}(t)\rangle\right\rVert^{2}}$, where

[TABLE]

with ${\mathcal{C}_{r}(t)\equiv\langle S_{r}^{z}(t)S_{\text{B}_{1}}^{z}\!(0)\rangle_{\text{eq}}-\langle S_{r}^{z}(t)S_{\text{B}_{2}}^{z}\!(0)\rangle_{\text{eq}}}$. Here, ${\langle\bullet\rangle_{\text{eq}}=\text{tr}[\bullet]/2^{N}}$ denotes the infinite-temperature ensemble, $\Theta(t)$ is the Heavyside function, and the sum runs over the jump times $\tau_{j}$ of the particular trajectory T. Moreover, the amplitudes $A_{j}$ in Eq. (17) read

[TABLE]

where ${A_{j}\to 1/2}$ for ${d_{\text{B}_{1},\text{T}}(\tau_{j}-0^{+})\to 0}$. Note that, due to the symmetry $d_{\text{B}_{1},\text{T}}(\tau_{j}-0^{+})=-d_{\text{B}_{2},\text{T}}(\tau_{j}-0^{+})$, only $\text{B}_{1}$ enters the above expressions. Equation (17) is the main result of this Letter. It predicts the magnetization dynamics in the open system by suitably superimposing equilibrium correlation functions of the closed system involving the two bath sites $\text{B}_{1}$ and $\text{B}_{2}$. In particular, from Eq. (17), the trajectory-averaged magnetization dynamics follows as

[TABLE]

where each $d_{r,\text{T}}(t)$ is evaluated for a different sequence ${(\tau_{1},\tau_{2},\ldots)}$ of $\tau_{j}$. Given the exponential damping in Eq. (12), the $\tau_{j}$ can be generated as ${\tau_{j+1}=\tau_{j}-\ln\varepsilon_{j+1}/2\gamma}$, where $\varepsilon_{j+1}$ are random numbers drawn from a box distribution ${]0,1]}$. If the correlation functions ${\langle S_{r}^{z}(t)S_{\text{B}_{1}}^{z}\!(0)\rangle_{\text{eq}}}$ and ${\langle S_{r}^{z}(t)S_{\text{B}_{2}}^{z}\!(0)\rangle_{\text{eq}}}$ are known, it is thus straightforward to evaluate Eq. (20) for a large number of sequences.

Numerical Illustration. We now test our theoretical prediction and its accuracy for a specific example, namely the spin-$1/2$ XXZ chain with ${\Delta=1.5}$, ${N=20}$, and periodic boundary conditions. The baths are located at $B_{1}=1$ and $B_{2}=11$ and we focus on small coupling ${\gamma/J=0.1}$ and weak driving ${\mu=0.1}$. Additional data for other values of $\Delta$, $\gamma$, and $\mu$, as well as for open boundary conditions can be found in [65].

Our theoretical prediction (20) is carried out numerically for ${{\cal O}(10^{4}-10^{5})}$ different sequences of jump times, which turns out to be sufficient to obtain negligibly small statistical errors. For comparison, we simulate the exact dynamics of the open system by performing a stochastic unraveling of the Lindblad equation. We stress that while (20) is derived in the limit of weak driving, cf. Eq. (12), the stochastic unraveling is here performed for the full $H_{\text{eff}}$ in Eq. (10).

In Fig. 2, we depict the outcome of the comparison. In Fig. 2(a), we show the time evolution of the local magnetization ${\langle S_{r}^{z}(t)\rangle}$ for different sites $r$. The site dependence of the steady-state profile is depicted in Fig. 2(b) and is well described by a linear function, except for the sites located exactly at the bath contacts. Importantly, we observe a remarkably good agreement between our prediction (20) and the exact open-system dynamics for all times up to ${tJ=50}$, where the steady-state profile is already established. This confirms our main result (17). We also note that for larger values of $\gamma$ and $\mu$, deviations are expected to become more pronounced, see [65].

For ${N\gg 20}$, stochastic unraveling cannot be carried out, since the required average over many trajectories becomes unfeasible. In contrast, our theoretical prediction (20) can be evaluated for larger system sizes, since only the equilibrium correlation functions are needed in Eq. (17). In particular, by relying on quantum typicality [26, 24], we simulate these correlation functions for up to $N=36$ lattice sites on Jülich’s “JUWELS” supercomputer. As shown in Figs. 3(a) and 3(b), we are thus able to describe the buildup of a nonequilibrium steady state in a $N=36$ XXZ chain weakly driven by Lindblad baths at sites $\text{B}_{1}=1$ and $\text{B}_{2}=19$. Open-system simulations for such system sizes are typically only accessible with matrix product state techniques, which are in turn usually restricted to open boundary conditions.

On the extraction of transport coefficients. In the nonequilibrium steady state, the diffusion constant can be calculated as ${D=-\langle j_{r}\rangle/\nabla\langle S_{r}\rangle}$ for some site $r$ in the bulk away from the bath sites. Here, $j_{r}$ is the local spin-current operator. Its expectation value can be expressed as ${\langle j_{r}\rangle=\frac{\text{d}}{\text{d}t}\langle\delta S^{z}_{\text{B}_{1}}\!(t)\rangle/2}$, where ${\langle\delta S^{z}_{\text{B}_{1}}\!(t)\rangle}$ is the magnetization injected by the first bath (see [65] for more details), and the factor $1/2$ takes into account that magnetization can flow to the left and to the right of this bath, due to periodic boundary conditions. As shown in Fig. 3(c), ${\langle\delta S^{z}_{\text{B}_{1}}\!(t)\rangle}$ grows linearly in time (i.e., a constant $\langle j_{r}\rangle$) and we can thus evaluate the diffusion constant in the steady state as the ratio of the slopes in Figs. 3(b) and 3(c). In this way, we obtain a value ${D/J\approx 0.99}$ which differs notably from what we found earlier in the context of Fig. 1.

This discrepancy may be explained by finite-size effects. Specifically, as also apparent in Fig. 1(a), the equilibrium correlation function $\langle S_{r}^{z}(t)S_{r^{\prime}}^{z}(0)\rangle_{\text{eq}}$ is affected by finite-size effects already at ${tJ\sim 20}$, where the broadening of the density profile has explored the full system. These effects then likely translate into the steady state in the weakly-driven open system and its finite-$N$ estimate of the diffusion constant. Such finite-size effects demonstrate that care must be taken when extracting transport properties both in closed and open systems. Importantly, we stress that the main conceptual result of our work, i.e., establishing a connection between weakly-driven Lindblad dynamics and closed quantum systems, remains unabated. In the supplemental material [65], we provide more details on this issue: Specifically, one can assume an ideal situation where the closed system behaves perfectly diffusive without finite-size corrections (in contrast to Fig. 1), in which case the equilibrium correlation functions $\langle S_{r}^{z}(t)S_{r^{\prime}}^{z}(0)\rangle_{\text{eq}}$ follow analytically as damped modified Bessel functions [65]. Using this idealized Ansatz, we find that the nonequilibrium steady state indeed yields the same diffusion constant as the closed system.

We note that one can extract a diffusion coefficient also from the finite-time dynamics of the open system, even before the steady state is established, via ${D(t)=\partial_{t}\langle S_{r}^{z}(t)\rangle/\nabla^{2}\langle S_{r}^{z}(t)\rangle}$, where ${\nabla^{2}\langle S_{r}^{z}(t)\rangle=\langle S_{r-1}^{z}(t)\rangle-2\langle S_{r}^{z}(t)\rangle+\langle S_{r+1}^{z}(t)\rangle}$. We are able to find a $D(t)$ in Fig. 4 that exhibits an approximately constant plateau ${D/J\approx 0.6}$ for ${tJ\lesssim 25}$ (while at longer $t$ the behavior becomes uncontrolled due to dividing two small numbers), consistent with our analysis of the closed system in Fig. 1.

Conclusion. Considering the example of magnetization transport in the spin-$1/2$ XXZ chain, we have connected linear response theory to the dynamics in an open quantum system driven by two Lindblad baths. Specifically, building on Ref. [54], we have shown that, at weak driving, the nonequilibrium steady state and its buildup in time can be constructed by suitably superimposing equilibrium correlation functions of the closed system.

Conceptually, our results for a specific model might reflect the natural expectation that transport coefficients obtained from closed-system and open-system approaches should agree with each other, at least if the driving is sufficiently weak. While we have presented data for systems with periodic boundary conditions, we provide additional results in [65], where we consider the more common case of open boundaries with Lindblad driving at the edge spins. In particular, we find that our main result (20) works convincingly also in this case and is in good agreement with state-of-the-art simulations based on time-evolving block decimation [47, 48]. From a practical perspective, our results enable the treatment of quite large open systems, which are usually not accessible by full stochastic unraveling. It would be an interesting attempt to generalize our setting to other jump operators, e.g., dephasing noise with $L_{j}=S_{j}^{z}$, and other questions beyond quantum transport.

Acknowledgments. We sincerely thank J. Wang for fruitful discussions. Our research has been funded by the Deutsche Forschungsgemeinschaft (DFG), projects 397107022 (GE 1657/3-2), 397300368 (MI 1772/4-2), and 397067869 (STE 2243/3-2), within DFG Research Unit FOR 2692, grant no. 355031190. J. R. acknowledges funding from the European Union’s Horizon Europe research and innovation programme, Marie Skłodowska-Curie grant no. 101060162, and the Packard Foundation through a Packard Fellowship in Science and Engineering. We gratefully acknowledge the Gauss Centre for Supercomputing e.V. for funding this project by providing computing time on the GCS Supercomputer JUWELS [77] at Jülich Supercomputing Centre (JSC). Z. L. and S. N. acknowledge support by the projects J1-2463 and P1-0044 program of the Slovenian Research Agency, EU via QuantERA grant T-NiSQ, and also computing time for the TEBD calculations at the supercomputer Vega at the Institute of Information Science (IZUM) in Maribor, Slovenia. We also acknowledge computing time at the HPC3 at University Osnabrück, which has been funded by the DFG, grant no. 456666331.

Details on the derivation of the theoretical prediction

In this section, we are going to sketch the derivation of our theoretical prediction given in Eq. (17) of the main text, which is an extension of the derivation in Ref. [54], where only a single Lindblad bath was considered, instead of the two baths treated here.

To this end, let us for the moment consider a simple scenario featuring a jump with the Lindblad operator $L_{1}$ immediately at some time $t$, say ${t=0}$. Then,

[TABLE]

For a random initial state $|\psi(0)\rangle$, as given in Eq. (16), this jump results in a random superposition over a subset of pure states with a spin-up at site $\text{B}_{1}$, which mimics ${\rho(0)\propto\mathbb{1}+S^{z}_{\text{B}_{1}}}$.

At weak driving $\mu\ll 1$, the subsequent deterministic evolution before the next jump reads

[TABLE]

cf. Eq. (12) of the main text. Now, using the concept of typicality, Eq. (S2) can be rewritten as

[TABLE]

with ${S_{r}^{z}(t)=e^{iHt}S_{r}^{z}e^{-iHt}}$ and ${\langle\bullet\rangle_{\text{eq}}=\text{tr}[\bullet]/2^{N}}$ denoting the infinite-temperature ensemble [59]. Analogously, one can obtain such a relation for the other possible jumps with the Lindblad operators $L_{j}$, which then involve either ${\langle S_{r}^{z}(t)S_{\text{B}_{1}}^{z}\!(0)\rangle_{\text{eq}}}$ or ${\langle S_{r}^{z}(t)S_{\text{B}_{2}}^{z}\!(0)\rangle_{\text{eq}}}$. Note that in the derivation of Eq. (S3), we used the facts that $S_{r}^{z}=n_{r}-\mathbb{1}/2$, $(n_{r})^{2}=n_{r}$, and $\text{tr}[S_{r}^{z}(t)]=0$, see e.g., Ref. [24] for more details.

For the above homogeneous initial state $|\psi(0)\rangle$, the jump probabilities according to Eq. (14) are simply given by $p_{j}=\alpha_{j}/4\gamma$. Consequently, averaging over all 4 jumps possibilities,

[TABLE]

yields the theoretical prediction

[TABLE]

for the time evolution after the first and before the second jump.

Let us consider a second jump at a later time $\tau$. The corresponding jump probabilities $p_{j}$ can then be derived based on typicality arguments. To this end, we assume a random pure state ${|\psi(\tau-0^{+})\rangle}$ right before the jump with ${\bar{d}_{\text{B}_{1}}\!(\tau-0^{+})\neq 0}$ and ${\bar{d}_{\text{B}_{2}}\!(\tau-0^{+})\neq 0}$. Then, due to the symmetry

[TABLE]

we have

[TABLE]

with ${y_{1}+y_{2}+y_{3}+y_{4}=1}$. Thus, the jump probabilities read

[TABLE]

with $p_{1}+p_{2}+p_{3}+p_{4}=1$. A straightforward calculation yields

[TABLE]

Using this expression, we can define the amplitude

[TABLE]

to incorporate both the probabilities for the next jump as well as the fact that some magnetization is already induced at the bath site. With this, we can eventually formulate a theoretical prediction for the deterministic evolution after the second jump, in analogy to the case of a single Lindblad bath [54]. This prediction reads

[TABLE]

with the Heavyside function $\Theta(t)$ and the abbreviation

[TABLE]

Reiterating this procedure finally yields a generalization of Eq. (S15) to a sequence of jump times $\tau_{j}$,

[TABLE]

i.e., Eq. (17) in the main text.

As discussed in Ref. [54], the central assumption within the above derivation is that the system has sufficient time to equilibrate between two jumps or, in other words, the magnetization injected at the contact sites has to spread over some region of the system. This requirement means that, in addition to a weak driving $\mu$, one has to choose a small coupling $\gamma$. Still, it might happen that even for a small coupling the equilibration process is hampered, as it is the case for open boundary conditions in certain models and parameter regimes, see the discussion below for more details.

Perfect Diffusion

As mentioned in the main text, we attribute differences between diffusion constants in open and closed systems to finite-size effects at long times, where the steady state is established. To support this, it is instructive to consider the idealized assumption of perfect diffusion in the closed system. In this case, which has no finite-size effects at any time, the equilibrium correlation functions take on the simple form [9]

[TABLE]

where ${\cal I}_{r}(t)$ is the modified Bessel function of the first kind and of the order $r$. Choosing ${D_{\text{closed}}/J=0.6}$, we show our prediction for the open system in Fig. S1. Indeed, using the relationship ${D=-\langle j_{r}\rangle/\nabla\langle S_{r}\rangle}$, we find a corresponding diffusion constant ${D/J\approx 0.60}$ in the steady state [80]. Thus, for perfectly diffusive behavior without finite-size effects, the steady state yields the same transport coefficient as the equilibrium correlation function (in contrast to the realistic case discussed in the context of Fig. 3 in the main text).

Furthermore, as already discussed in the main text, it is possible to calculate $D(t)$ already at finite times before the steady state is established, via

[TABLE]

with

[TABLE]

which is just the diffusion equation for a lattice in one spatial dimension. Evaluating this expression for, e.g., site ${r=9}$ and time ${tJ=25}$ for the perfectly diffusive data in Fig. S1, we obtain the value

[TABLE]

which is again in good agreement with $D_{\text{closed}}$.

Open Boundaries and comparison with TEBD simulations

So far, we have focused on systems with periodic boundary conditions, which are the natural choice for closed systems, whereas state-of-the-art matrix product state approaches to open systems commonly rely on open boundary conditions with Lindblad driving at the systems’ edges.

In Fig. S2, we show a comparison between our theoretical prediction (20) and full stochastic unraveling, similar to Fig. 2 in the main text, but now for a XXZ chain with open boundaries, where the two baths are placed at the ends of the chain, ${\text{B}_{1}=1}$ and ${\text{B}_{2}=N}$. As shown in Fig. S2(a), the time evolution of the local magnetization ${\langle S_{r}^{z}(t)\rangle}$ is well captured by the prediction (20), though deviations start to appear at times ${tJ\gtrsim 20}$. These slight deviations might be caused by the fact that the equilibrium correlation functions can exhibit unusual behavior in the case of open boundary conditions. In particular, $\langle S_{\text{B}_{1}}^{z}\!(t)S_{\text{B}_{1}}^{z}\!(0)\rangle$ and $\langle S_{\text{B}_{2}}^{z}\!(t)S_{\text{B}_{2}}^{z}\!(0)\rangle$ do not fully decay for any $\Delta>1$ in the case of open boundary conditions due to the presence of a strong zero mode, where the edge spins retain memory of their initial conditions for very long times [78, 79].

In addition to our prediction (20) and the stochastic unraveling, we include in Fig. S2 numerical data obtained by a state-of-the-art matrix product state implementation based on time-evolving block decimation (TEBD) with time step $\mathrm{d}tJ=0.05$ and bond dimension $\chi=128$. Importantly, we find that this TEBD data is in perfect agreement with the results from stochastic unraveling.

The convincing agreement between TEBD, stochastic unraveling, and our theoretical prediction is further highlighted in Fig. S2(b), where the site dependence of the profile ${\langle S_{r}^{z}(t)\rangle}$ is shown for time ${tJ=50}$. Moreover, we find that this profile is well described by a linear function in the bulk, far away from the bath contacts. The injected magnetization ${\langle\delta S^{z}_{\text{B}_{1}}\!(t)\rangle}$ is shown in Fig. S2(c). The diffusion constant in the open system is again given by the ratio of the slopes in Figs. S2(b) and S2(c). This way, we find the value ${D/J\approx 0.63}$ in very good agreement with the value ${D/J\approx 0.6}$ in the closed system.

Injected magnetization

The first Lindblad bath injects magnetization at the corresponding contact site $r=\text{B}_{1}$. This magnetization can also be predicted and written as

[TABLE]

with

[TABLE]

cf. Eq. (17). This magnetization can then be related to the local currents, which are all the same in the steady state, i.e.,

[TABLE]

Thus, it is sufficient to know $\langle j_{\text{B}_{1}}\rangle$, which follows from the injected magnetization via

[TABLE]

where the factor $1/2$ takes into account that the injected magnetization can flow to the left and to the right of this bath, due to periodic boundary conditions. By the use of this expression, we find that, for the case discussed in Fig. 2, $\langle j_{r}\rangle/J\approx 0.0028$.

For comparison, we can calculate the local currents for the same model within the stochastic unraveling as

[TABLE]

for $\mathcal{O}(10^{4})$ trajectories. For long times $tJ\gtrsim 30$, where the steady state is established, we find a corresponding value of $\langle j_{r}\rangle/J\approx 0.003$, which is close to the predicted one above.

Large Coupling / Strong Driving

In the main text, we have focused on the case of small coupling ${\gamma/J=0.1}$ and weak driving ${\mu=0.1}$, where we have found a convincing agreement between the derived prediction and exact numerics in Fig. 2. To illustrate that deviations occur for larger coupling or stronger driving, we depict a corresponding comparison for (a) ${\gamma/J=1}$ (${\mu=0.1}$) and (b) ${\mu=1}$ (${\gamma/J=0.1}$) in Fig. S3. In both cases (a) and (b), deviations are visible already at finite times before the steady state is reached. Interestingly, in the case (a) of strong coupling, the overall agreement is still satisfactory and the profile in the steady state is predicted accurately.

Other Anisotropies

In Fig. 2 of the main text, we have provided a comparison of the magnetization dynamics for anisotropy ${\Delta=1.5}$. Complementarily, we show a comparison for anisotropies (a) ${\Delta=1.0}$ and (b) ${\Delta=0.5}$ in Fig. S4. For the case of ${\Delta=1.0}$, we again find a convincing agreement. While for the case of ${\Delta=0.5}$ the overall behavior of prediction and numerics is still similar, deviations are visible at long time scales. These deviations reflect that the prediction is not exact in a mathematical sense but involves physical assumptions (such as, e.g., on the equilibration properties of the involved equilibrium correlation functions) which may not hold perfectly in any given situation.

\close@column@grid