Inhomogeneous MPA and exact steady states of boundary driven spin chains at large dissipation

TL;DR

This paper presents an exact solution for a boundary-driven XYZ spin-1/2 chain with strong dissipation, using novel site-dependent Lax operators and an inhomogeneous matrix product ansatz, revealing complex matrix structures.

Contribution

It introduces new site-dependent Lax operators and an inhomogeneous MPA for solving the XYZ spin chain in the Zeno limit, expanding understanding of dissipative quantum systems.

Findings

Exact steady state expressed via inhomogeneous MPA

Matrices cannot be simplified into a tridiagonal form

Results applicable to open XYZ and eight-vertex models

Abstract

We find novel site-dependent Lax operators in terms of which we demonstrate exact solvability of a dissipatively driven XYZ spin-1/2 chain in the Zeno limit of strong dissipation, with jump operators polarizing the boundary spins in arbitrary directions. We write the corresponding nonequilibrium steady state using an inhomogeneous MPA, where the constituent matrices satisfy a simple set of linear recurrence relations. Although these matrices can be embedded into an infinite-dimensional auxiliary space, we have verified that they cannot be simultaneously put into a tridiagonal form, not even in the case of axially symmetric (XXZ) bulk interactions and general nonlongitudinal boundary dissipation. We expect our results to have further fundamental applications for the construction of nonlocal integrals of motion for the open XYZ model with arbitrary boundary fields, or the eight-vertex model.

Equations

$$H=\sum _{n=1}^{N-1}{h}_{n,n+1},h_{n,n+1}={\vec{\sigma }}_n\cdot \mathrm{J}{\vec{\sigma }}_{n+1},$$

$$[H+{\vec{h}}_{\mathrm{l}}\cdot {\vec{\sigma }}_1+{\vec{h}}_{\mathrm{r}}\cdot {\vec{\sigma }}_N,R]=0,$$

$${\mathbf{L}}_n:={\vec{\sigma }}_n\cdot {\vec{L}}_n=\sum _{\alpha \in \mathcal{J}}{\sigma }_n^{\alpha }{L}_n^{\alpha }$$

$$[h_{n,n+1},{\mathbf{L}}_n{\mathbf{L}}_{n+1}]=i(I_n{\mathbf{L}}_{n+1}-{\mathbf{L}}_n{I}_{n+1}),$$

$$\sum _{\alpha ,\beta \in \mathcal{J}}[h_{n,n+1},{\sigma }_n^{\alpha }{\sigma }_{n+1}^{\beta }]L_n^{\alpha }{L}_{n+1}^{\beta }=$$

$$=i\sum _{\alpha \in \mathcal{J}}\left({\sigma }_{n+1}^{\alpha }{I}_n{L}_{n+1}^{\alpha }-{\sigma }_n^{\alpha }{L}_n^{\alpha }{I}_{n+1}\right).$$

$$\Omega =⟨0|{\mathbf{L}}_1\cdots {\mathbf{L}}_N|\psi ⟩,$$

$$\displaystyle\big{[}\sigma_{n}^{\alpha}A_{n}\big{]}\!\otimes\!\big{[}\sigma_{n}^{\beta}B_{n}\big{]}\!:=\!\sigma_{n}^{\alpha}\sigma_{n}^{\beta}A_{n}\otimes B_{n}$$

$${\mathbb{L}}_n={\mathbf{L}}_n\otimes {\mathbf{L}}_n^{\ast }:=\sum _{\alpha ,\beta \in \mathcal{J}}{\sigma }_n^{\alpha }{\sigma }_n^{\beta }{L}_n^{\alpha }\otimes (L_n^{\beta }{)}^{\ast },$$

$$R=⟨0,\bar{0}|{\mathbb{L}}_1\dots {\mathbb{L}}_N|\psi ,\bar{\psi }⟩.$$

$$\displaystyle\begin{gathered}\big{[}H+\vec{h}_{\rm l}\cdot\vec{\sigma}_{1}+\vec{h}_{\rm r}\cdot\vec{\sigma}_{N},R\big{]}=\bra{0,\bar{0}}\mathbb{F}_{1}\mathbb{L}_{2}\ldots\mathbb{L}_{N}\ket{\psi,\bar{\psi}}+\bra{0,\bar{0}}\mathbb{L}_{1}\ldots\mathbb{L}_{N-1}\mathbb{F}_{N}\ket{\psi,\bar{\psi}},\\ \mathbb{F}_{1}=[2\,\vec{h}_{\rm l}\cdot\vec{L}_{1}+i\,I_{1}]\otimes\mathbf{L}^{*}_{1}-\mathbf{L}_{1}\otimes[2\,\vec{h}_{\rm l}\cdot\vec{L}^{*}_{1}-i\,I_{1}],\quad\mathbb{F}_{N}=[2\,\vec{h}_{\rm r}\cdot\vec{L}_{N}-i\,I_{N}]\otimes\mathbf{L}^{*}_{N}-\mathbf{L}_{N}\otimes[2\,\vec{h}_{\rm r}\cdot\vec{L}^{*}_{N}+i\,I_{N}].\end{gathered}$$

$$⟨0|[2{\vec{h}}_{\mathrm{l}}\cdot {\vec{L}}_1+i{I}_1]=0,[2{\vec{h}}_{\mathrm{r}}\cdot {\vec{L}}_N-i{I}_N]|\psi ⟩=0.$$

$${\vec{L}}_n\times \mathrm{J}{\vec{L}}_{n+1}=\frac{1}{2}{\vec{L}}_n{I}_{n+1},\mathrm{J}{\vec{L}}_n\times {\vec{L}}_{n+1}=\frac{1}{2}{I}_n{\vec{L}}_{n+1}.$$

$$⟨k;{\mathcal{A}}_n|:=⟨k;{\mathcal{A}}_{n-1}|I_n,0\le k\le n-1$$

$${\mathcal{A}}_0=\mathbb{C}⟨0;{\mathcal{A}}_0|,{\mathcal{A}}_n={\mathcal{A}}_{n-1}{I}_n\oplus \mathbb{C}⟨n;{\mathcal{A}}_n|.$$

$${\mathcal{A}}_0=\mathbb{C}⟨0|,{\mathcal{A}}_n={\mathcal{A}}_{n-1}\oplus \mathbb{C}⟨n|.$$

$$L_n^{\alpha }=\sum _{k=0}^{n-1}\sum _{l=0}^n{L}_{n;k,l}^{\alpha }|k⟩⟨l|,I_n=\sum _{k=0}^{n-1}|k⟩⟨k|,$$

$$L_1^x$$

$$L_1^y$$

$$L_1^z$$

$$r=({\xi }^2+{\eta }^2)({\omega }_{xy}{\eta }^2-1)({\omega }_{xy}{\xi }^2+1)({\omega }_{xz}{\xi }^2+{\omega }_{yz}{\eta }^2+1).$$

$$L_n^x=\frac{1}{2}(L_n^{+}+L_n^{-}),L_n^y=\frac{1}{2i}(L_n^{-}-L_n^{+}),$$

$$\begin{array}{c}{L}_n^z=\sum _{k=0}^{n-1}|k⟩⟨k+1|,\\ L_n^{\pm }=\pm {\eta }^{\mp 1}\sum _{k=0}^{n-1}\sum _{l=0}^n{\left(\frac{\pm i}{2\cos \gamma }\right)}^{k-l+1}{M}_{n;k,l}|k⟩⟨l|,\end{array}$$

$$\begin{array}{c}{M}_{n;k,l}=\frac{(\xi -{\xi }^{-1})P_{n,k+1,l}(\cos \gamma )}{(\xi +{\xi }^{-1})\sin \gamma }-\frac{2P_{n,k+1,l-1}(\cos \gamma )}{(\xi +{\xi }^{-1})\cos \gamma },\\ P_{n,k,l}(x)=\sum _{m=0}^l(-1{)}^m(\genfrac{}{}{0ex}{}{n-k}{m}\left)\right(\genfrac{}{}{0ex}{}{n-m-1}{l-m})x^{n-2m}.\end{array}$$

$$\begin{array}{rl}{K}_n^{+}{L}_{n+1}^{\pm }& =q^{\pm 1}{L}_n^{\pm }{K}_{n+1}^{+},\\ K_n^{-}{L}_{n+1}^{\pm }& =q^{\mp 1}{L}_n^{\pm }{K}_{n+1}^{-},\\ L_n^{+}{L}_{n+1}^{-}-L_n^{-}{L}_{n+1}^{+}& =\frac{K_n^{+}{K}_{n+1}^{+}-K_n^{-}{K}_{n+1}^{-}}{q-q^{-1}},\\ K_n^{+}{K}_{n+1}^{-}& =K_n^{-}{K}_{n+1}^{+}.\end{array}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\rho (t)=-i[H^{\prime },\rho (t)]+\Gamma {\mathcal{D}}_{\mathrm{l}}[\rho (t)]+\Gamma {\mathcal{D}}_{\mathrm{r}}[\rho (t)],$$

$${\mathcal{D}}_{\mu }[\rho ]=2k_{\mu }\rho k_{\mu }^{†}-\{k_{\mu }^{†}{k}_{\mu },\rho \},$$

$$k_{\mathrm{l},\mathrm{r}}=({\vec{n}}_{\mathrm{l},\mathrm{r}}^{\prime }+i{\vec{n}}_{\mathrm{l},\mathrm{r}}^{\prime \prime })\cdot {\vec{\sigma }}_{0,N+1}$$

$$\vec{n}(\theta ,\varphi )=(\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,\cos \theta ).$$

$${\mathcal{D}}_{\mu }[{\rho }_{\mu }]=0,\mathrm{tr}[{\rho }_{\mu }\vec{\sigma }]={\vec{n}}_{\mu },$$

$$H^{\prime }=H+h_{0,1}+h_{N,N+1}.$$

$${\rho }^{(0)}={\rho }_{\mathrm{l}}\otimes R\otimes {\rho }_{\mathrm{r}},$$

Full text

I Introduction

Matrix product ansatz (MPA) is arguably one of the most useful theoretical concepts in statistical and quantum physics of one-dimensional locally interacting systems. It appears in a diverse variety of contexts, ranging from exact form of the ground state for a certain type of non-integrable spin-1 chains (the so-called AKLT model AKLT1 ; AKLT2 of valence bond solids) to exact description of the non-equilibrium steady states of both classical interacting Markov chains (e.g. simple exclussion processes Derrida ; Evans and driven cellular automata Carlos ), as well as Lindblad equation in quantum integrable systems TP2011 ; JPAreview ; in special cases it can even describe the full time evolution Katja . Moreover, it enters a general description of the so-called finitely-correlated-states Werner , as well as the variational ansatz for a classical simulation of equilibrium and time-dependent quantum states (aka DMRG-related methods Vidal ; Schollwoeck ). In all known theoretical applications of MPA, the constituent matrices of the ansatz are position independent, and satisfy certain bulk cancellation condition, related to a particular matrix representation of either Yang-Baxter or Zamolodchikov-Faddeev algebra.

In the context of boundary-driven open quantum systems, i.e. integrable spin chains with Lindblad jump operators that act only on the boundary sites, it is particularly challenging to understand the maximal set of dissipative boundary processes for which the non-equilibrium steady state density matrix can be written exactly. So far, this has only been possible (for bulk integrable models such as XXZ spin-$1/2$ chain or Fermi-Hubbard model) for the so-called pure-source/pure-sink boundaries, or boundaries which target opposite longitudinal directions JPAreview . Note, though, that a global SU(2) symmetry allows for a more general solvable boundary processes in the isotropic XXX model KPS2013 . It has, however, remained an open question if and how these exact steady state solutions fit into the general framework of integrablity. For example, except in the special case of dissipatively driven noninteracting models fabian , the solvable dissipatively driven boundaries cannot be generated using the solutions of the ubiquitous reflection equations sklyanin , which constitute the standard framework for generating integrable boundaries in the coherent (nondissipative, Hamiltonian) setting or in the quantum Hamiltonian formulation of classical Markov processes.

In this paper we propose a new direction for a general construction of integrable incoherent boundaries of interacting quantum chains. We construct local Lax operators that at first sight appear not to be related to the standard solutions of the Yang-Baxter equation for the eight-vertex model takhtajan ; sklyaninxyz and have auxiliary dimension which, unlike in the usual analytic scenario, differs from site to site, i.e. auxiliary spaces at different physical sites are manifestly non-isomorphic. These site-dependent Lax operators generate a conserved transfer matrix with an inhomogeneous matrix product structure.

As a straightforward application of our result we use this mechanism to solve the problem of a boundary driven anisotropic XYZ spin-$1/2$ chain in the limit of strong dissipation (the so-called Zeno regime), where the driving mechanism polarises the boundary-localised degrees of freedom in a fixed direction of arbitrary choice.

The paper is organized in two parts. In the first part we introduce the model and the novel Lax operators and show, how they can be used to construct operators that commute with the model’s Hamiltonian. Our main technical tool is to show the validity of a generalized divergence condition for the Lax operators that guarantees cancellation of unwanted terms in the bulk. The divergence condition appears as an infinite set of recurrence relations that can be solved once the initial seed is provided. Complete analytical ansatz is established rigorously for a special case of the XXZ spin-$1/2$ chain. In the more general case of XYZ model, we only explicitly provide the seed for the recurrence. The complexity of these equations currently only allows us to treat this second case as a numerical scheme. The second part of the letter deals with applications. Here, we introduce the dissipative boundary processes with arbitrary polarization and, in the limit of strong dissipation, treat them using solutions of the recurrence relations. In the XXZ case we provide the explicit inhomogeneous matrix product form of the non-equilibrium steady state, while in the XYZ case we provide a carefully empirically verified (conjectured) computational recipe for its construction.

II Inhomogeneous Matrix Product Ansatz

We consider a quantum spin-$1/2$ chain on a $N$-site one-dimensional lattice. Each spin is acted upon by Pauli matrices $\sigma^{\alpha}\in{\rm End}(\mathbb{C}^{2})$, where $\alpha\in{\cal J}=\{x,y,z\}$. For each $n\in\{1,2,\ldots N\}$ we denote the local one-site operators by $\sigma^{\alpha}_{n}=\mathbb{1}_{2^{n-1}}\otimes\sigma^{\alpha}\otimes\mathbb{1}_{2^{N-n}}$, where $\mathbb{1}_{d}$ is a $d\times d$ identity matrix. The dynamics that we will consider in this exposition is generated by the anisotropic Heisenberg Hamiltonian (also known as XYZ model)

[TABLE]

which acts over the total Hilbert space ${\cal H}=(\mathbb{C}^{2})^{\otimes N}$. Here and in the following ${\rm J}={\rm diag}(J_{x},J_{y},J_{z})$ denotes the diagonal tensor of spin coupling constants (i.e. the anisotropy tensor) and $\vec{\sigma}_{n}=(\sigma^{x}_{n},\sigma^{y}_{n},\sigma^{z}_{n})$. For ${\rm J}=\mathbb{1}_{3}$ the model describes the isotropic Heisenberg magnet (XXX model), while ${\rm J}={\rm diag}(1,1,\cos\gamma)$ is a parametrization of the axially symmetric XXZ model.

We will present a new MPA, by means of which one can construct a positive, semi-definite operator $R=\Omega\Omega^{\dagger}\in{\rm End}({\cal H})$ that satisfies

[TABLE]

where $\vec{h}_{\rm l},\vec{h}_{\rm r}\in\mathbb{R}^{3}$ are arbitrary magnetic field polarizations on the left and the right hand side of the spin chain, respectively. The aim is to construct conserved quantities, unrelated to the off-diagonal Bethe ansatz ODBA1 ; ODBA2 that reproduces the integrable hierarchy of the Hamiltonian (1) with arbitrary boundary magnetic fields. The main physical application that we shall discuss here is in the context of boundary driven spin chains, also described in the parallel work inPrep19 , but one may, in addition, use it for encoding generic (inhomogeneous) conservation laws with a finite spatial correlation structure. We start by describing the site-dependent MPA for $\Omega$ and the mechanism responsible for the validity of Eq. (2).

II.1 Inhomogeneous cancellation mechanism

Consider a sequence $\{{\cal A}_{n}\}_{n=0}^{N}$ of auxiliary vector spaces with dimensions ${\rm dim}({\cal A}_{n})=n+1$ and let $L^{\alpha}_{n},\,I_{n}\in{\rm Lin}({\cal A}_{n-1},{\cal A}_{n})$ for $\alpha\in{\cal J}$ be linear maps between them. Denoting $\vec{L}_{n}=(L_{n}^{x},L_{n}^{y},L_{n}^{z})$ for each $n\in\{1,2,\ldots N\}$, we define Lax operator

[TABLE]

as an element of ${\rm Lin}({\cal H}\otimes{\cal A}_{n-1},{\cal H}\otimes{\cal A}_{n})$. A sequence of such site-dependent Lax operators will constitute the MPA for the factor $\Omega$ of the operator $R$ satisfying (2). The key to solving this equation is the inhomogeneous Sutherland equation

[TABLE]

in which the commutator on the left hand side concerns only the operators acting on ${\cal H}$; the equation should be read as

[TABLE]

For the ansatz we now set

[TABLE]

where $\bra{0}\in{\cal A}_{0}$ and $\ket{\psi}=\sum_{n=0}^{N}\psi_{n}\ket{n}$ are boundary vectors that we will identify later. The right boundary vector is inferred from $\bra{\psi}=\sum_{n=0}^{N}\psi_{n}^{*}\bra{n}\in{\cal A}_{N}$ by invoking the duality relation: for $\bra{k}\in{\cal A}_{N}$, the dual vector is defined through $\braket{k}{l}=\delta_{k,l}$.

Let $\otimes$ denote a (partial) tensor product of two copies of the auxiliary space that acts as an ordinary matrix multiplication over the physical (quantum) space ${\cal H}$. For $A_{n},B_{n}\in{\rm Lin}({\cal A}_{n-1},{\cal A}_{n})$ and arbitrary $\alpha,\beta\in{\cal J}$ it is defined as

[TABLE]

and then extended by linearity. Introducing a two-point Lax operator

[TABLE]

an element of ${\rm Lin}({\cal H}\otimes{\cal A}_{n-1}^{\otimes 2},{\cal H}\otimes{\cal A}_{n}^{\otimes 2})$, we can now write the MPA for the whole operator $R$:

[TABLE]

Here $\ket{\psi,\bar{\psi}}:=\ket{\psi}\otimes(\ket{\psi})^{*}$, while $(\bullet)^{*}$ denotes complex conjugation over the auxiliary space and hermitian conjugation over the physical space: $\mathbf{L}_{n}^{*}=\sum_{\alpha}(\sigma_{n}^{\alpha})^{\dagger}(L^{\alpha}_{n})^{*}=\sum_{\alpha}\sigma_{n}^{\alpha}(L^{\alpha}_{n})^{*}$. As shown in Appendix A, utilisation of the Sutherland equation (4) now results in expression

[TABLE]

We now see that the operator $R=\Omega\Omega^{\dagger}$, where $\Omega$ is given in (6), commutes with the Hamiltonian with boundary magnetic fields, i.e. satisfies equation (2), provided that

[TABLE]

In what follows, we will (i) provide unique solutions to the inhomogeneous Sutherland equation (4), thus specifying the MPA, and (ii) show an example of an interesting and nontrivial physical application where the boundary equations (13) can be solved.

II.2 Solution of the cancellation mechanism

A straightforward calculation shows, that Sutherland equation (4) is component-wise equivalent to a pair of discrete Landau-Lifshitz equations

[TABLE]

Fixing $\vec{L}_{n}$ and operators $I_{n},I_{n+1}$, this is an overdetermined set of linear equations for $\vec{L}_{n+1}$. In this paper we will demonstrate (partly prove) and use the following:

Proposition**.**

For a fixed initial datum $\vec{L}_{1}$ (seed), which depends on two free complex parameters, there exist a solution to recurrence relations (14), unique up to a choice of basis in each ${\cal A}_{n}$.

II.2.1 Specifying the bases of the auxiliary spaces

To specify the basis in each auxiliary space, we assume that operators $I_{n}$ are non-degenerate, choose basis $\{\bra{k;{\cal A}_{n-1}}\}_{k=0}^{n-1}$ of ${\cal A}_{n-1}$, and then define

[TABLE]

as the first $n$ elements of the basis of ${\cal A}_{n}$. Next, assuming non-degeneracy of some other Lax component, say $L_{n}^{z}$, we define an additional basis vector $\bra{n;{\cal A}_{n}}:=\bra{n-1;{\cal A}_{n-1}}L^{z}_{n}$, to get the sequence

[TABLE]

If we interpret the operator $I_{n}$ as an inclusion map $I_{n}:{\cal A}_{n-1}\hookrightarrow{\cal A}_{n}$ we can denote all basis elements by $\bra{k}$, irrespective of the auxiliary space. The entire auxiliary sequence is then embedded in an infinite-dimensional linear space ${\cal A}_{\infty}={\rm lsp}\{\bra{k}\}_{k=0}^{\infty}$ and corresponds to

[TABLE]

Introducing the dual basis $\{\ket{k}\}$ through the orthogonality relation $\braket{k}{l}=\delta_{k,l}$, the operators $L^{\alpha}_{n}$ and $I_{n}$ can now be represented as rectangular $n\times(n+1)$ matrices

[TABLE]

where $L^{z}_{n;n-1,l}=\delta_{n,l}$. Having specified the basis, the boundary vectors of the MPA (6) are now identified. On the left hand side we simply take the state $\bra{0}$ that, by subsequent action of $L_{1}^{z},L_{2}^{z},\ldots$, generates the entire basis in the embedded space ${\cal A}_{\infty}$, while on the right hand side $\ket{\psi}=\sum_{n=0}^{N}\psi_{n}\ket{n}$, where the complex parameters $\psi_{n}$ will be fixed by the second of the boundary equations (13).

II.2.2 Solutions of the recurrence and connection to integrability

Solving the nonlinear coupled equation (14) for $\vec{L}_{1}$, at $n=1$ and for arbitrary spin coupling constants $J_{\alpha}$ we obtain – up to either trivial or equivalent solutions – the following two-parametric solution for the seed ($\xi,\eta\in\mathbb{C}$):

[TABLE]

where $\omega_{\alpha\beta}:=4\,(J_{\alpha}^{2}-J_{\beta}^{2})$ and we have denoted

[TABLE]

Using a symbolic computer algebra we have checked, that the overdetermined linear equations (14) now generate unique $\vec{L}_{n}$ for $n=2,3\ldots N$. Each matrix element of any auxiliary Lax component $L_{n}^{\alpha}$ is of the form $p(\xi,\eta)+q(\xi,\eta)\sqrt{r}$, where $p,q$ are some rational functions. Unfortunately the complexity of the solution quickly increases with $n$ and we were unable to determine its explicit analytic structure. Hence, for $n>5$ and arbitrary $J_{\alpha}$ one can only efficiently solve the recurrence equations (14) numerically.

Nevertheless, for the special case of XXZ model, where $J_{x}=J_{y}=1,J_{z}=\cos\gamma$, $\gamma\in\mathbb{R}$ (or $i\mathbb{R}$), the recurrence (14) can in fact be explicitly analytically solved (see Appendix B). After changing the basis by writting

[TABLE]

the solution reads

[TABLE]

where we have introduced the following symbols:

[TABLE]

Here, the free variables $\xi,\eta\in\mathbb{C}$ provide a different parametrization than those in Eqs. (19), where the spin coupling constants $J_{\alpha}$ are arbitrary. For general $\xi,\eta$ we have checked that this solution of the Sutherland equation (4) [or, equivalently, the recurrence relations (14)] cannot be reduced to any known solution of the Yang-Baxter equation by means of local twists in the auxiliary spaces ${\cal A}_{n}$. Instead, recombining the Lax component $L^{z}_{n}$ and the inclusion operator $I_{n}$ as $K_{n}^{\pm}=\tfrac{i}{2}(I_{n}\pm 2\sin\gamma L_{n}^{z})$ and denoting $q=e^{i\gamma}$, we notice that Eq. (14) is equivalent to the inhomogeneous quantum group relations

[TABLE]

These relations define an inhomogeneous analogue of the $q$-deformed spin algebra ${\cal U}_{q}(sl_{2})$ and thus provide a yet-unexplored manifestation of the XXZ spin-$1/2$ chain integrability structure, albeit now with an extensive (in system size) number of generators. Note that, when site-independent, relations (26) describe the $sl_{2}$ spin algebra as $q\to 1$ ($XXX$ model). On the other hand, while the matrices $L^{\alpha}_{n}$ of our ansatz (24) do become elementwise site-independent in this limit, they do not reduce to the standard spin ladder operators. In the following we illustrate the facility of the recurrence scheme (14) for solving a boundary driven Lindblad equation in the regime of strong dissipation.

III Application: Quantum Zeno limit of the boundary driven XYZ chain

We wish to use the ansatz, described in the preceding section, to construct the nonequilibrium steady state (NESS) of the Lindblad equation

[TABLE]

at large dissipation strength $\Gamma$, where ${\cal D}_{\mu}[\rho]$, $\mu\in\{{\rm l},{\rm r}\}$, denote the dissipators at the left and right ends of the chain of $N+2$ sites, which we label by [math] and $N+1$, respectively. They are of the form

[TABLE]

with the two jump operators

[TABLE]

targeting polarizations $\vec{n}_{\mu}=\vec{n}(\theta_{\mu},\phi_{\mu})$, where

[TABLE]

Real vectors $\vec{n}^{\prime}_{\mu}=\vec{n}(\frac{\pi}{2}-\theta_{\mu},\pi+\phi_{\mu})$ and $\vec{n}^{\prime\prime}_{\mu}=\vec{n}(\frac{\pi}{2},\phi_{\mu}-\frac{\pi}{2})$, together with $\vec{n}_{\mu}$, form an orthonormal basis of $\mathbb{R}^{3}$. The targeted states of the dissipators are single-site pure states $\rho_{\mu}=|\psi_{\mu}\rangle\langle\psi_{\mu}|$, such that

[TABLE]

i.e. $k_{\mu}\ket{\psi_{\mu}}=0$ ($\ket{\psi_{\mu}}$ should not be confused with the right boundary vector $\ket{\psi}$ of the MPA). The Hamiltonian is now provided by Eq. (1) extended by two sites:

[TABLE]

The problem of constructing the NESS in the limit of strong dissipation (Zeno limit) has been rigorously examined, but not solved in 2018PopkovZenoDynamics . In the limit $\Gamma\to\infty$, when the unitary part of the dynamical equation (27) can be neglected, the NESS should obviously be of the form

[TABLE]

where $R$ is some operator acting on the Hilbert space ${\cal H}$ of the internal degrees of freedom, labeled with $1,2,\ldots N$. For large but finite $\Gamma$, we can proceed perturbatively by expanding $\rho_{\infty}=\sum_{k\geq 0}\Gamma^{-k}\rho^{(k)}$. Plugging the expansion into the Lindblad equation (27), demanding stationarity ${\rm d}\rho_{\infty}/{{\rm d}t}=0$ and comparing the orders of $\Gamma^{-1}$, we get

[TABLE]

which is automatically satisfied, as well as a sequence of equations

[TABLE]

Equations (35) in particular state that $[H^{\prime},\rho^{(k)}]$ belongs to the image space of the dissipator (superoperator) ${\cal D}_{\rm l}+{\cal D}_{\rm r}$, which only acts on the boundary degrees of freedom, labeled by [math] and $N+1$. Since normalization of the density matrix is conserved by the Lindblad equation (27) and separately by its unitary part, respectively ${\rm d}(\operatorname{tr}\!{[}\rho(t)])/{{\rm d}t}=0$, and $\operatorname{tr}\!{[}H^{\prime},\rho(t)]=0$, this yields additional condition $\mathop{\mathrm{tr}_{0,N+1}}\limits[H^{\prime},\rho^{(k)}]=0$ NJPdecoherence2015 . For $k=0$ it explicitly reads

[TABLE]

where $H_{\cal D}$ is the dissipation-projected Hamiltonian that acts on sites $1,2,\ldots N$ and takes the following form:

[TABLE]

Although this condition seems rather insignificant at first, it remarkably constitutes the core of our solution to the problem of the strongly boundary driven spin chain.

Indeed, we have arrived at the problem defined in the first section, which can be solved by our ansatz (6), if equations (13) are satisfied for the boundary field orientations $\vec{h}_{\rm l}={\rm J}\vec{n}_{\rm l}$ and $\vec{h}_{\rm r}={\rm J}\vec{n}_{\rm r}$. The left boundary equation $\bra{0}[2\,({\rm J}\vec{n}_{\rm l})\cdot\vec{L}_{1}+i\,I_{1}]=0$, in reality a set of two equations for two variables, completely fixes the parameters $\eta$ and $\xi$ in the components $L^{\alpha}_{1}$. In the XXZ case, the solution reads

[TABLE]

In the general XYZ model, the solution to the left boundary equation (13) exists as well and is unique for our choice of bases in ${\cal A}_{0}$ and ${\cal A}_{1}$. Alternatively, we can choose a gauge, different than in Eqs. (19), (20) and (21), in which the seed that generates the solution to the recurrence (14) becomes explicitly dependent on the left-edge polarisation axes $\vec{n}_{\rm l}$,

[TABLE]

and satisfies the left boundary equation by construction.

Having specified the parameters, thus fixing the ansatz in the bulk of the system we now turn to the right boundary equation in (13), i.e. $[2\,({\rm J}\vec{n}_{\rm r})\cdot\vec{L}_{N}-i\,I_{N}]\ket{\psi}=0$, which determines $\ket{\psi}$. Writing $\ket{\psi}=\sum_{n=0}^{N}\psi_{n}\ket{n}$, with $\psi_{0}=1$, this is a set of $N$ linear equations for $N$ unknowns $\psi_{n}$. One solution always exists and seems to be unique for generic values of the boundary angles $\theta_{\mu}$ and $\phi_{\mu}$. In particular cases, for example, for XXZ chain with $\theta_{\rm r}=\phi_{\rm r}=0$, it can easily be computed analytically:

[TABLE]

In general, we compute it numerically.

When unique, the resulting operator $\rho^{(0)}=\rho_{\rm l}\otimes R\otimes\rho_{\rm r}$ indeed reproduces the NESS of the Lindblad equation (27) in the Zeno limit: (i) In special cases, where the latter is known analytically 2012XYtwist , we find it in complete agreement with our ansatz. (ii) In generic cases, we resort to comparison with numerically exact NESS, computed via a method proposed in 2018PopkovZenoDynamics , which yields equivalence up to the preset numerical precision. (iii) For finite values of the dissipation strength $\Gamma$, the ansatz $\rho^{(0)}$ converges towards the NESS of the finite-$\Gamma$ Lindblad equation (27), i.e. towards the solution of

[TABLE]

as shown in Fig. 1, again indicating that the ansatz is correct. The right-hand-side plot on Fig. 1 also shows that operators $R$ and $H_{\cal D}$ are functionally independent, i.e. $R\neq f(H_{\rm D})$ for at least a piece-wise smooth function $f$, in turn implying nontriviality of our ansatz.

Note that there are also cases, in which the right boundary vector $\ket{\psi}$ of the MPA is not unique. We hypothesize this to happen in measure-zero subset of the parameter space. Even in this case, however, we find that the Zeno NESS is correctly reproduced by our ansatz for a specific choice of the right boundary vector. Resolving this issue analytically requires considering higher orders $\rho^{(k)}$ of the perturbative expansion, which is out of our present scope.

The MPA expression for $\rho^{(0)}$ allows for an efficient computation of local observables, such as magnetization profiles and spin current, for previously inaccessible system sizes; see Fig. 2. In Fig. 3 we plot the phase diagram of the spin current exhibiting high sensitivity with resonance spiking as a function of anisotropy parameter. For a detailed analysis of the problem we refer the reader to Ref. inPrep19 .

III.1 Zeno limit with asymmetric dissipation rates

Apart from the symmetric dissipative action on both boundaries, we can as well consider Lindblad problem with infinitely large, but different dissipation rates at the left and the right edge. To this end, we renormalize the Lindblad jump operators as

[TABLE]

where $0<\kappa$ is the measure of the left-right asymmetry. The ratio of the effective dissipative rates is then fixed to $\kappa$. We can now study the Zeno limit of the Lindblad master equation as a function of $\kappa$, for the general XYZ model.

Heuristically we observe the following remarkable fact: nonequilibrium steady state for the XYZ model does not depend on the asymmetry $\kappa$ of the dissipation rates, as long as both dissipation rates go to infinitity, i.e. in the Zeno limit $\Gamma\to\infty$ (for any finite dissipation $\Gamma$ the steady state of course depends on $\kappa$). This property is rather exceptional and is related to a subtle property of the Zeno effective dynamics of the XYZ model, described below.

It has been shown in 2018PopkovZenoDynamics that the effective Zeno dynamics of a quantum system is governed by (a) the dissipation projected Hamiltonian $H_{\cal D}$ (37), (b) by a classical Markov process with rates $w_{\alpha,\beta}$, calculated using the eigenstates $\ket{\alpha}$ of $H_{\cal D}$ and some auxiliary operators $g_{\mu}$, calculated from the dissipator,

[TABLE]

For instance, under generic assumption of a non-degenerate spectrum, the bulk of the NESS density matrix is diagonal in the basis $\ket{\alpha}$, $R=\sum_{\alpha}p_{\alpha}\ket{\alpha}\bra{\alpha}$, where the ‘probability vector’ $\{p_{\alpha}\}$ is given as an invariant state (fixed point) of a classical Markov process

[TABLE]

If two systems have the same dissipation-projected Hamiltonians, and the solutions of Eq. (42) are also the same, then the Zeno limit of NESS is the same as well. This is exactly the situation we have: firstly, the dissipation projected Hamiltonian is just determined by the kernel of the dissipator and therefore does not depend on the asymmetry $\kappa$. It is given in Eq. (37). Secondly, for our choice of the dissipators, the sum over $\mu$ in (43) consists of two terms, each one associated with a separate boundary; the respective $g_{\rm l}(\kappa),g_{\rm r}(\kappa)$ can be calculated using the method developed in 2018PopkovZenoDynamics and read

[TABLE]

where the unit vectors $\vec{n}_{\mu}^{\prime},\vec{n}_{\mu}^{\prime\prime}$, for $\mu\in\{{\rm l},{\rm r}\}$, have been introduced after Eq. (30).

For a general “MPA-integrable” case we now numerically observe

[TABLE]

The properties (46) and (47) are very special; they hold if the bulk is described by a homogeneous XYZ (integrable) Hamiltonian. We checked, for example, that if we take a nonintegrable Hamiltonian (i.e., switch on integrability breaking terms), they are no longer satisfied.

Due to the Kolmogorov criterion (47), the steady state probabilities $p_{\alpha}$ satisfy the detailed balance condition $p_{\alpha}/p_{\beta}=w_{\beta,\alpha}/w_{\alpha,\beta}$. Then, due to Eq. (46), we have

[TABLE]

Since the probabilities $p_{\alpha}(\kappa)$ are normalized, $\sum_{\alpha}p_{\alpha}(\kappa)=1$, the steady state of the associated Markov process for the asymmetric boundary driving is the same as in the case of symmetric driving: $p_{\alpha}(\kappa)=p_{\alpha}(1)$. Consequently, the Zeno limit of NESS will remain the same, i.e., it will be independent of the asymmetry $\kappa$ in the dissipation at the left and the right boundary.

IV Discussion

The Sutherland equation – divergence condition (4) and the boundary equations (13) are two crucial ingredients in the construction of conservation laws and nonequilibrium steady states of boundary driven spin chains. Here we have proposed a generalized, inhomogeneous Sutherland equation, in which the Lax matrices of the MPA explicitly depend on the lattice site. We have demonstrated the applicability of the resulting MPA by generating the nonequilibrium steady state of a boundary driven XYZ spin-$1/2$ chain with strong dissipative spin-polarizing boundary baths. Generically, our ansatz (6) can be also used as a tool to construct nontrivial conservation laws for the open spin chain with arbitrary nondiagonal boundary fields (2).

The structure of constituent matrices of our ansatz (6) is very different from that of previously treated Lax operators, which satisfy the celebrated Yang-Baxter equation. Besides having a site-dependent auxiliary structure, our Lax operators cannot be put into a tridiagonal form, even after all of the nonisomorphic local auxiliary spaces ${\cal A}_{n}$ are embedded into a joint infinite-dimensional auxiliary vector space. For example, it can be checked that our explicit representation (24,25) cannot be reduced to the highest weight representation of the ${\cal U}_{q}(sl_{2})$ quantum group symmetry of the XXZ model, which has been used to solve Lindblad equation for the longitudinal TP2011 or transverse KPS2013 dissipative boundaries. In other words, the Lax structure proposed here, seems to correspond to a new representation of the underlying symmetry algebra, in which the auxiliary space is not fixed to some ${\cal U}_{q}(sl_{2})$ module, but rather corresponds to a ladder of linear vector spaces, transitions between which are represented by matrices of our ansatz. Similarly, we expect that for the anisotropic XYZ model our inhomogeneous Lax operators and nonequilibrium dissipative solutions go beyond the off-diagonal Bethe ansatz which diagonalizes the closed Hamiltonian ODBA1 ; ODBA2 . It is left as an open future problem to find explicit analytic expression for the inhomogeneous Lax operators in the general XYZ case, presumably in terms of Jacobi elliptic functions.

We acknowledge discussions with M. Petkovšek and V. Romanovsky. The work has been supported by European Research Council (ERC) through the advanced grant 694544 – OMNES and the grant P1-0402 of Slovenian Research Agency (ARRS). V.P. also acknowledges support by the DFG grant KL 645/20-1.

Appendix A Cancellation mechanism and the boundary equations

In this appendix we elaborate on the boundary equations, that need to be satisfied in order for the commutation relation

[TABLE]

to hold for an operator $R=\Omega\Omega^{\dagger}$, with $\Omega=\bra{0}\mathbf{L}_{1}\mathbf{L}_{2}\ldots\mathbf{L}_{N}\ket{\psi}$. The Hamiltonian is given by Eq. (1), while the inhomogeneus Lax operators $\mathbf{L}_{n}=\vec{\sigma}_{n}\cdot\vec{L}_{n}$ satisfy the so-called divergence condition, given in Eq. (4). A straightforward application of the latter yields

[TABLE]

where $(\bullet)^{*}$ denotes complex conjugation over the auxiliary space and hermitian conjugation over the physical space and $\ket{\bar{\psi}}:=(\ket{\psi})^{*}$. For example, for $\alpha\in{\cal J}$ we have $(\sigma_{n}^{\alpha}L^{\alpha}_{n})^{*}:=(\sigma_{n}^{\alpha})^{\dagger}(L_{n}^{\alpha})^{*}=\sigma_{n}^{\alpha}(L_{n}^{\alpha})^{*}$. Using $\mathbb{L}=\mathbf{L}\otimes\mathbf{L}^{*}=\sum_{\alpha,\beta\in{\cal J}}\sigma_{n}^{\alpha}\sigma_{n}^{\beta}\,L_{n}^{\alpha}\otimes(L_{n}^{\beta})^{*}$, where $\otimes$ denotes the tensor product over auxiliary spaces and ordinary matrix multiplication over the physical space ${\cal H}$, we can rewrite this as

[TABLE]

On the other hand we have

[TABLE]

where the commutators can be explicitly rewritten as

[TABLE]

and similarly $[\vec{h}_{\rm r}\cdot\vec{\sigma}_{N},\mathbb{L}_{N}]=2\,(\vec{h}_{\rm r}\cdot\vec{L}_{N})\otimes\mathbf{L}^{*}_{N}-2\,\mathbf{L}_{N}\otimes(\vec{h}_{\rm r}\cdot\vec{L}^{*}_{N})$. Putting everything together, we get

[TABLE]

If the boundary equations $\bra{0}[2\,\vec{h}_{\rm l}\cdot\vec{L}_{1}+iI_{1}]=0$ and $[2\,\vec{h}_{\rm r}\cdot\vec{L}_{N}-iI_{N}]\ket{\psi}=0$ are satisfied, the operator $R$ commutes with the Hamiltonian $H+\vec{h}_{\rm l}\cdot\vec{\sigma}_{1}+\vec{h}_{\rm r}\cdot\vec{\sigma}_{N}$.

Appendix B Proof of the ansatz in the XXZ case

In the XXZ case, the tensor of anisotropic spin-spin interactions becomes ${\rm J}={\rm diag}(1,1,\cos\gamma)$. Writing $L^{x}_{n}=\frac{1}{2}(L^{+}_{n}+L^{-}_{n})$ and $L^{y}_{n}=\frac{1}{2i}(L^{-}_{n}-L^{+}_{n})$, the discrete spatial Landau-Lifshitz equations given by Eq. (14) hold, if

[TABLE]

Our goal in this appendix is, to show that the ansatz

[TABLE]

satisfies algebraic relations (55). This will be done in two parts. Firstly, we will discuss three lemmas which will facilitate the proof of the relations themselves. The latter will be presented in the second part.

B.1 Lemmas

Lemma 1**.**

Polynomials given in (56), satisfy the following recurrence relations

[TABLE]

Proof.

This is a simple consequence of the Pascal rule for the binomial coefficients. ∎

Remark*.*

Note, that the recurrence relations hold irrespective of what integer $l$ is. For example, we have $P_{n,k,0}(x)=x^{n}$ and $P_{n,k,l}=0$ for $l<0$, which is consistent with the relations.

Lemma 2**.**

Binomial coefficients satisfy relations

[TABLE]

for $0\leq t\leq n-1$, $t^{\prime}\in\mathbb{Z}$ and

[TABLE]

for $-1\leq t\leq n-1$, $t^{\prime}\in\mathbb{Z}$.

Proof.

The first relation is

[TABLE]

In the last equality we have used one of the standard binomial sum identities: $\sum_{s=m}^{n}(-1)^{n-s}{n\choose s}{s\choose m}=\delta_{n,m}$. To do so, we have truncated the sums, $\sum_{s^{\prime}=0}^{n-1}\to\sum_{s^{\prime}=t^{\prime}-1}^{n-t-1}$ and $\sum_{s^{\prime}=0}^{n-1}\to\sum_{s^{\prime}=t^{\prime}}^{n-t-1}$, respectively. This is possible even for $t^{\prime}\leq 0$. In this case, the first sum will vanish, since ${a\choose b}=0$ if $b<0$. On the other hand, we can only change the upper bound from $n-1$ to $n-t-1$ if $t\geq 0$.

The second relation is

[TABLE]

Again, we have used the identity $\sum_{s=m}^{n}(-1)^{n-s}{n\choose s}{s\choose m}=\delta_{n,m}$, after truncating the sum according to $\sum_{s^{\prime}=0}^{n}\to\sum_{s^{\prime}=t^{\prime}}^{n-t-1}$. This is possible even for $t=-1$ and $t^{\prime}\leq 0$. ∎

Lemma 3**.**

For $0\leq k\leq n-1$ polynomials given by (56), satisfy the relations

[TABLE]

Proof.

We start by proving the first relation. We write out the left hand side:

[TABLE]

Since ${a\choose b}=0$ for $a<b$ or $b<0$, we can truncate the sum over $s$ at $n-1$ and extend sums over $t$ and $t^{\prime}$ up to $n-1$ and $n$, respectively. We get

[TABLE]

which, after using Lemma 2, becomes

[TABLE]

To produce the Kronecker deltas via identity $\sum_{s=m}^{n}(-1)^{n-s}{n\choose s}{s\choose m}=\delta_{n,m}$, we had to truncate the sum over $t$ at $n-k-1$. This is allowed by the assumption $k\geq 0$.

The second relation is even simpler to prove, again starting by writing out the left hand side:

[TABLE]

Since ${a\choose b}=0$ for $b>a$ or $b<0$, we can extend the sum over $t$ up to the maximum $s-1=n-1$. The sum over $t^{\prime}$ can be extended up to $n$. We then get

[TABLE]

and using Lemma 2,

[TABLE]

This completes the proof of the three lemmas. ∎

Remark*.*

The polynomial relations from Lemma 3 are trivially satisfied even for $l<0$, since then $P_{n,k,l}=0$.

B.2 Proof of the algebraic relations

We will now finally prove, that the ansatz (56) satisfies the algebraic relations (55), using the Lemmas 1 and 3 from the previous subsection.

Corollary 1**.**

Relations $L_{n}^{+}L_{n+1}^{-}-L_{n}^{-}L_{n+1}^{+}=i\,L_{n}^{z}I_{n+1}$ and $L_{n}^{+}L_{n+1}^{-}-L_{n}^{-}L_{n+1}^{+}=i\,I_{n}L_{n+1}^{z}$ are satisfied.

Proof.

Since they are equivalent, we will only prove the first one. Explicitly, the first relation reads

[TABLE]

where

[TABLE]

Because of the prefactor $(-1)^{k}-(-1)^{l}$, only the cases where $k-l$ is an odd integer need to be checked. Since there is no $\xi$-dependence on the right hand side of the relation (68), the second sum in $A_{k,l}$ should be zero. Note, that we can use Lemma 3 in all of the terms of the matrix element $A_{k,l}$. This gives

[TABLE]

For odd $k-l$ it becomes $A_{k,l}=(-1)^{k+1}\cos\gamma\,\delta_{k,l-1}$. Using this result we see that (68) is indeed satisfied. ∎

Corollary 2**.**

Relations $L_{n}^{z}L_{n+1}^{\pm}-\cos\gamma\,L_{n}^{\pm}L_{n+1}^{z}=\pm\frac{i}{2}\,L_{n}^{\pm}I_{n+1}$ are satisfied.

Proof.

Explicitly, they are both equivalent to

[TABLE]

where we note $M_{n;k,l-1}=0$, for $l=0$. They are obviously satisfied by courtesy of the first polynomial recurrence in Lemma 1. ∎

Corollary 3**.**

Relations $\cos\gamma\,L_{n}^{z}L_{n+1}^{\pm}-L_{n}^{\pm}L_{n+1}^{z}=\pm\frac{i}{2}\,I_{n}L_{n+1}^{\pm}$ are satisfied.

Proof.

Explicitly, they read

[TABLE]

Again, note $M_{n;k,l-1}=0$, for $l=0$. These relations are satisfied due to the second polynomial recurrence in Lemma 1. ∎