Diffusion and entanglement in open quantum systems

Recently, transport properties of nonequilibrium quantum systems have been extensively studied [1–3]. Of particular interest are dissipative quantum systems [4–12], which in many cases lead to diffusive transport —ubiquitous in both quantum and classical physics.

The macroscopic fluctuation theory (MFT) [13], a universal coarse grained theory, has been rigorously shown to completely encapsulate the dynamics of classical diffusive systems [14–16]. The MFT allows to calculate nonequilibrium long-range correlations [17], large deviations [18], fluctuation-induced forces [19], predict dynamical phase transitions [20–24], their associated critical exponents [16,25,26] and many more [27,28]. The present study is motivated by a salient question: Does a universal description of quantum diffusive dynamics, comparable to the MFT, exist?

Let us briefly recap the MFT. Rather than describing a diffusive process at the mesoscopic scale through a master equation, the MFT takes a coarse grained approach where, away from criticality, only conserved fields survive. To make this more concrete, consider a 1D diffusive system of interacting particles with particle density n. Conservation of particles in the bulk implies the continuity equation $\partial_t n = -\partial_x j$. Then, the diffusion equation results from the combination of the continuity equation and Fick's law of the mean current $j = -D \partial_x n + \sigma E$, where $D,\sigma$ are the diffusivity and conductivity encapsulating the dynamics and E is an external field. To infer what the correlations in the system are, one must go beyond the mean description and include fluctuations. The MFT asserts that the fluctuations of the process at the coarse grained level are entirely described by supplementing the current in Fick's law with a white noise term of magnitude proportional to the conductivity. Therefore, regardless of the microscopic description, the dynamics of the coarse grained mean behaviour as well as the fluctuations are entirely encapsulated by $D,\sigma$. This is the universal structure of the MFT.

A universal coarse grained theory of quantum diffusion is required to go beyond the classical MFT, e.g., to exhibit entanglement. First, it is necessary to provide a microscopic model that exhibits entanglement and diffusive dynamics at the coarse grained level [29], expanding the MFT's mean evolution structure. Then, one can verify whether the structure of the fluctuations remains universal.

In this letter we introduce a microscopic model described by a quantum master equation that demonstrates at the coarse grained level both diffusive dynamics and entanglement at the steady state. The mean evolution is a set of non-local diffusion equations. This example demonstrates that a general theory describing diffusive open quantum system goes beyond the MFT. Furthermore, it provides us with a working example that will allow in turn to develop such a diffusive theory of open quantum systems.

Since simple diffusive systems are Markovian and dissipative, a convenient initial starting point would be to consider Markovian open quantum systems, given by the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation [30]. Consider a quantum system, coupled to a large environment with fast relaxation times. The back-action of the system on the environment can thus be neglected on the relevant time scale of the system. The evolution of the system's density matrix follows $\partial_t \rho = \mathcal{L}(\rho)$, where

$$\begin{eqnarray} \begin{array}{rcl} \mathcal{L}(\rho) &= - i [H,\rho] + \displaystyle\sum_k L_{M_k} (\rho),\\ L_M (\rho) &= M \rho M^\dagger - \displaystyle\frac{1}{2} \{ M^\dagger M, \rho\}. \end{array} \end{eqnarray} \tag{ 1 }$$

Here $[ \cdot, \cdot]$ and $\{ \cdot, \cdot \}$ are correspondingly the commutation and anti-commutation relations. H is a Hermitian operator and may contain both the system's Hamiltonian and its interaction with the environment. The set of M_k operators, orthonormal under the trace norm $\textrm{Tr} M^\dagger_k M_k' =\delta_{k,k'}$, are related to the interaction with the environments. One can also follow the evolution of operators in the Heisenberg picture, given by the adjoint equation $\partial_t O = \mathcal{L}^* (O)$. If all the M_k are Hermitian, $\mathcal{L}^* (\bullet) = +i [H,\bullet] + \sum_k L_{M_k}(\bullet)$.

To explore quantum diffusion, we introduce the "selective dephasing" model. Consider a periodic spin chain of $\Omega>2$ sites, with the Pauli matrices $\sigma^{x,y,z,\pm}_k$ operating at site k. The selective dephasing model has the GKSL form (1). We set H to be the XX Hamiltonian $H_{\textrm{XX}} = \varepsilon \sum^\Omega_{k=1 } \sigma^x_k \sigma^x_{k+1} + \sigma^y_k \sigma^y_{k+1}$. The operators $M_k = \sqrt{\eta \nu_f} \sigma^z_k$ correspond to dephasing on the sites $k\in \mathcal{A}$, where $\nu_f$ has units of inverse time and η is a dimensionless parameter that controls the strength of the dephasing. Three choices of $\mathcal{A}$, which select the dephasing sites, are considered:

$$\begin{eqnarray} \mathcal{A}_\Omega &\equiv \{ 1,2,\ldots,\Omega \}, \quad \mathcal{A}_{\{ \xi \}} \equiv \{ j\in \mathcal{A}_\Omega | j\neq \xi,\xi+2 \}, \nonumber\\ \mathcal{A}_{\textrm{odd}} &\equiv \{ 1,3,5,\ldots,\Omega-1 \}, \end{eqnarray} \tag{ 2 }$$

where for $\mathcal{A}_{\textrm{odd}}$, Ω is assumed to be even. Note that the exact value of ξ is irrelevant due to the periodic boundary conditions we consider. This GKSL model can be derived using the singular coupling limit for spins coupled to a set of thermalized bosonic baths [30] and for fermions in an optical lattice [31].

The $\mathcal{A}_{\Omega}$ model has been exhaustively studied [9,32–34] and it even admits a Bethe ansatz solution [10]. First, we recap the derivation of the diffusive dynamics of the $\mathcal{A}_\Omega$ model. The dephasing leaves no quantum fingerprints at long times. Then, we show that protecting some of the sites from dephasing as in the $\mathcal{A}_{\{ \xi \} },\mathcal{A}_{\textrm{odd}}$ models may maintain the diffusive dynamics and exhibits quantum fingerprints at long times.

For the analysis, we define the operators

$$\begin{eqnarray} \begin{array}{rcl} P^\pm _k &=& \sigma^\pm _k \sigma^\mp _k , \quad \mathbb{P}_k = P^+ _k \displaystyle\prod_{l\in \mathcal{A}_\Omega/\{ k \} } P^- _l ,\\ \mathbb{S}^m _k &=& (\sigma^+ _k \sigma^- _{k+m} +\sigma^- _k \sigma^+ _{k+m} ) \displaystyle\prod_{l\in \mathcal{A}_\Omega/\{ k, k+m \} } P^- _l . \end{array} \end{eqnarray} \tag{ 3 }$$

In the following, the long time dynamics of the selective dephasing models will be shown to be governed by these operators. Furthermore, the $\mathbb{S}^m_k$ will be used as witnesses of entanglement.

The operator evolution for $P^+_i$ is given by the discrete conservation equation $\partial_t P^+_i = -J_i + J_{i-1}$, where $J_i = 2i \varepsilon ( \sigma^+_i\sigma^-_{i+1} - \sigma^-_i\sigma^+_{i+1} )$. Therefore, on a ring $Q = \sum_{k\in \mathcal{A}_\Omega} \textrm{Tr} P^+_k \rho(s)$ is conserved for any η. This is true for any $\mathcal{A}$ dephasing model. Trying to evaluate $\partial_t J_i$ reveals that the equations do not close [34], so this naive approach fails. Fortunately, the dynamics simplifies as states not in the kernel of $\sum_{k \in \mathcal{A}_\Omega} L_{M_k}$ are quickly filtered [9,34] (see the Supplementary Material Supplementarymaterial.pdf (SM)). States in the kernel are called pointer states. The resulting long time dynamics, quickly attained for large η values, is governed by transitions between these pointer states. Hence, we derive a perturbative expression for the pointer states dynamics at large η for a general GKSL equation of the form

$$\begin{eqnarray} \mathcal{L}(\rho) &= \mathcal{L}_S (\rho) + \eta \mathcal{L}_b (\rho),\\ \mathcal{L}_S (\rho) &= - i [H,\rho], \quad \mathcal{L}_b (\rho) = \nu_f \displaystyle\sum_{k\in \mathcal{A}} L_{M_k}(\rho).\nonumber\\ \end{eqnarray} \tag{ 4 }$$

For large η, the density matrix quickly converges into a mixture of pointer states —the states in the kernel of $\mathcal{L}_b$. Then, the density matrix can only evolve on a long time scale [35–37]. A slow dynamics, where the density matrix spreads to other pointer states emerges if $\Pi_0 \mathcal{L}_S \Pi_0 =0$ and $\mathcal{L}_b = \sum_{\nu\leq0} \nu \Pi_\nu$. Here $\Pi_\nu$ are projectors and specifically $\Pi_0$ is the projector into the kernel. To leading order $\partial_s\rho=\mathfrak{u} \rho$, where [9,34] (see the SM)

$$\begin{equation} \mathfrak{u} = - \Pi_0 \mathcal{L}_S (\mathcal{L}^\perp_b)^{-1} \mathcal{L}_S, \end{equation} \tag{ 5 }$$

with the rescaled time $s=t/\eta$ at the limit of $t,\eta \rightarrow \infty$. $(\mathcal{L}^\perp_b)^{-1} = \sum_{\nu<0} \nu^{-1} \Pi_\nu$ is the inverse to the restriction of $\mathcal{L}_b$ outside of the kernel. The above assumptions are valid in all the variants of the selective dephasing model (see the SM). Finding the projector to the kernel $\Pi_0$ may be easier than finding a basis for the pointer states as one needs to worry about positivity and unit trace. For Hermitian M_k, the operator evolution is $\partial_s O = \mathfrak{u} O$. Hence, dealing with pointer operators, i.e., $O = \Pi_0 O$ may be easier than dealing with pointer states.

To analyze the effective dynamics, it is essential to identify the kernel of $\mathcal{L}_b$. For $\mathcal{A}_{\Omega}$, this kernel is spanned by the $2^\Omega$ pointer states $\prod_k P^{\pm}_k$. One can give a classical interpretation to the effective dynamics of (5) in this case. Interpreting the operators $P^\pm$ as projectors into an occupied (empty) lattice site state implies that the pointer states span all the configurations of a lattice gas. Namely, each site is either occupied by a particle or not. The dynamics corresponds to the simple symmetric exclusion process (SSEP), where particles can jump to empty neighboring lattice sites with rate $D= \nicefrac{2\varepsilon^2}{\nu_f}$ [38,39]. The SSEP is diffusive and classical. So, we have gained no insight on the behaviour of diffusive open quantum systems.

By protecting some of the sites from the dephasing noise, one expects a larger kernel, which in turn could lead to richer dynamics. This is the reasoning behind the $\mathcal{A}_{\{ \xi \}}, \mathcal{A}_{\textrm{odd}}$ models. We restrict the study here to the already interesting case Q = 1. For the analysis of the qualitatively similar $Q\neq1$ see the SM.

The kernel becomes larger than in the $\mathcal{A}_\Omega$ case as it contains, e.g., the pointer operator $\mathbb{S}^2_\xi$ (see (3)), as well as the states $\mathbb{P}_k$. The effective dynamics for the $\mathcal{A}_{\{ \xi \}}$ dephasing provides us with a closed set of pointer operator equations,

$$\begin{eqnarray} \begin{array}{rcl} \partial_s \mathbb{P}_k &=& D \Delta \mathbb{P}_k, \quad k \neq \xi-1,\xi,\xi+1,\xi+2,\\ \partial_s \mathbb{P}_k &=& 2D \Delta \mathbb{P}_k + D \mathbb{S}^2 _\xi, \quad k = \xi,\xi+2,\\ \partial_s \mathbb{P}_{\xi+1} &=& 2D \Delta \mathbb{P}_{\xi+1} - 2D\mathbb{S}^2 _{\xi},\\ \partial_s \mathbb{P}_{\xi-1} &=& D \Delta \mathbb{P}_{\xi-1} +D(\mathbb{P}_{\xi}-\mathbb{P}_{\xi-1}),\\ \partial_s \mathbb{S}^2 _\xi &=& -2D \Delta \mathbb{P}_{\xi+1} -4D\mathbb{S}^2 _\xi, \end{array} \end{eqnarray} \tag{ 6 }$$

where Δ is the discrete Laplacian operator.

Note that in the slow dynamics regime, operators outside of the kernel have already decayed in the Heisenberg picture, i.e., $\partial_s O =0$ if $\Pi_0 O =0$. Therefore, more information on the slow dynamics can be recovered from the evolution equations of the (non-pointer operators) $\mathbb{S}^2_{k\neq\xi}$. Since $\Pi_0\mathbb{S}^2_{k\neq\xi} =0$, we obtain a set of identities

$$\begin{eqnarray} \begin{array}{rcl} \Delta \mathbb{P}_{k+1}&=&0, \quad |k-\xi|\geq 3,\\ \Delta \mathbb{P}_{\xi+3} &=& \mathbb{P}_{\xi+3} - \mathbb{P}_{\xi+2},\\ \Delta \mathbb{P}_{\xi-1} &=& \mathbb{P}_{\xi-1} - \mathbb{P}_{\xi},\\ \mathbb{S}^2 _\xi &=& \Delta \mathbb{P}_\xi = \Delta\mathbb{P}_{\xi+2}. \end{array} \end{eqnarray} \tag{ 7 }$$

The slow dynamics identities (7) are at the operator level and hence do not depend on initial conditions. We denote by $S_{\xi}(s) = \mathrm{Tr}\mathbb{S}^2_\xi \rho(s)$ and $q_k (s) = \mathrm{Tr}\mathbb{P}_k \rho(s)$ the expectation values associated with the pointer states. From the $\Omega-1$ identities (7) and from the conserved $Q=1= \sum_k q_k$, we recover a single solution $q_{k\neq \xi+1} = \frac{1}{\Omega} (1-S_\xi)$ and $q_{\xi+1} = \frac{1}{\Omega}( 1 + (\Omega-1)S_\xi)$. The density matrix of the effective dynamics is a (positive semidefinite, Hermitian and trace 1) combination of the operators $\mathbb{P}_k$ and $\mathbb{S}^2_\xi$. The Peres-Horodecki criterion asserts that a two two-level system is entangled if and only if the partial transpose to the density matrix is not positive definite [40] (see the SM). We use the Peres-Horodecki criterion to probe for bipartite entanglement in our system. In the effective dynamics bipartite entanglement can occur in the $\mathcal{A}_{\{ \xi \} }$ model only between the spins in $\xi,\xi+2$ and only if $S_\xi \neq 0$ (see the SM). Hence, $S_\xi$ serves as a witness of entanglement in the system.

We evaluate the values $S_\xi$ and q_k using the evolution equations (6). They correspond to the trivial solution $S_\xi=0$ and $q_k =1/\Omega$ in the effective dynamics. Hence, entanglement quickly decays. This suggests, as expected, that a microscopic protection from the environment is not enough to retain an interesting hydrodynamic quantum behaviour (see the SM for numerical verification).

As before, we study the effective dynamics at large η to obtain the evolution within the kernel. Notice that $\mathbb{P}_k$ for $k\in \mathcal{A}_\Omega$ and the operators $\mathbb{S}^{2m}_{2k}$ are in the kernel for $k,m=1,2,\ldots,\Omega/2$. For Q = 1, we obtain a set of evolution equations and identities corresponding to the operators in the kernel. Here, we keep track also of $\mathbb{S}^{2m}_k$ for odd k, even though they are not in the kernel and thus have vanishing expectation values in the effective dynamics. This highlights the emerging diffusive picture, but does not change any expectation value (see fig. 1 and the SM). The (diffusive) evolution equations are

$$\begin{eqnarray} \begin{array}{rcl} \partial_s \mathbb{P}_{k} &=& 2D \Delta \mathbb{P}_{k} - D \Delta \mathbb{S}^2 _{k-1},\hspace{8pc}\\ \partial_s \mathbb{S}^2 _{k} &=& -2D\Delta \mathbb{P}_{k+1} + 2D\Delta \mathbb{S}^2 _{k}-D\Delta \mathbb{S}^4 _{k-1},\\ \partial_s \mathbb{S}^{2+2m} _{k} &=& - D\Delta \mathbb{S}^{2m} _{k+1} +2D\Delta \mathbb{S}^{2m+2} _{k}- D\Delta \mathbb{S}^{4+2m} _{k-1}, \end{array} \end{eqnarray} \tag{ 8 }$$

for $2\leq 2+2m \leq \Omega/2$. The periodicity of the system helps to truncate eq. (8) as $\mathbb{S}^{\Omega/2 +2 }_k = \mathbb{S}^{\Omega/2 - 2}_{k+ (\Omega/2 +2)}$ for $\Omega/2$ even, and $\mathbb{S}^{\Omega/2 +1 }_k = \mathbb{S}^{\Omega/2 - 1}_{k+ (\Omega/2 +1)}$ for $\Omega/2$ odd.

Zoom In Zoom Out Reset image size

To find the steady state solution, we use the translational invariance of the system (in discrete jumps of 2). Define $q_{\rm o} = \mathrm{Tr} \rho(s) \mathbb{P}_{2k+1}$, $q_{\rm e} = \mathrm{Tr} \rho(s) \mathbb{P}_{2k}$ and $S^{2m} =\mathrm{Tr} \rho(s) \mathbb{S}^{2m}_{2k}$. From (8) and for $\Omega/2$ even, we find the steady state values $S^{2m} = -S^{2m+2}$ as well as $q_{\rm o}-q_{\rm e} = \frac{1}{2} S^2$. This is numerically confirmed in the SM.

Similarly to the $A_{\{\xi\}}$ dephasing, the expectation value of S^2m is a witness of entanglement in the steady state. Note that the steady state density matrix is constructed using the pointer operators $\mathbb{P}_k, \mathbb{S}^{2m}_{2k}$. Then, using the Peres-Horodecki criterion shows that any two even sites $2k',2k''$ are bipartite entangled if the expectation of $\mathbb{S}^{2(k''-k')}_{2k'}$ is non-vanishing (see the SM).

For $\Omega/2$ odd, periodicity implies that $S^{\Omega/2+1} = S^{\Omega/2-1}$. Then, the steady state equations result in $S^{2m}=0, q_{\rm o}=q_{\rm e} = 1/\Omega$ and thus the system will not be entangled at the steady state. In what follows, we only study $\Omega/2$ even, where the entanglement can survive at the steady state. The effective evolution equations (8) suggest that the expectation values of $\sum_{k \in \mathcal{A}_\Omega}\mathbb{S}^{2m}_{k}$ are conserved quantities in the effective dynamics. However, this is not the case for the t time evolution. It implies that only beyond a transient period $t \gtrsim \eta \nu_f/ \varepsilon^2$, $\sum_{k\in \mathcal{A}_\Omega}\mathbb{S}^{2m}_{k}$ is conserved.

Using the effective conservation, we define a witness of entanglement in the system $B(t) = \sum_{k\in\mathcal{A}_\Omega}\mathrm{Tr} \rho(t) \mathbb{S}^2_k$. Using this witness, fig. 1 and the SM numerically demonstrate that entanglement survives the long time limit for generic initial conditions, provided the initial conditions have some nonzero entanglement value (other initial conditions were tested with similar results). Notice that the entanglement survives on all length scales of the system. That is, the bipartite entanglement spreads to all length scales even if it is initially restricted to a local region. Hence, we have successfully found a model where long-range entanglement can be observed for diffusive transport.

Contrary to the classical diffusion equation description, the hydrodynamic equations for the quantum diffusion of the $\mathcal{A}_{\textrm{odd}}$ dephasing model can be expected to be non-local due to the long-ranged entanglement. For $\Omega/2$ even, we define the coarse grained length $x=j/\Omega$, and time $\tau = s D / \Omega^2$ scales as well as the coarse grained operators $\mathbb{P}(x,t),\mathbb{S}^{2m}(x,\tau)$. Using the vector representation $v = ( \mathbb{P}, \mathbb{S}^2, \mathbb{S}^4, \ldots \mathbb{S}^{\Omega/2} )$, the non-local hydrodynamic diffusion equation stemming from (8) is

$$\begin{equation} \partial_\tau v(x) = \partial_{xx} \mathcal{D}_0 v(x) + \partial_{xx} \mathcal{D}_1 v(x+\Omega/2), \end{equation} \tag{ 9 }$$

where

$$\begin{eqnarray*} \hspace{-22pt}\mathcal{D}_1 &= \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 0 & \ldots & & \\ \vdots & \ddots & &\\ & & 0 &0 \\ & & -1 &0 \end{array}\right),\nonumber \end{eqnarray*}$$

$$\begin{eqnarray} \mathcal{D}_0 &= \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} 2 & -1 & 0 & 0 & 0 & \ldots& 0 \\ -2 & 2 & -1 & 0 & 0 & \ldots & 0 \\ 0 & -1 & 2 & -1 & 0 & \ldots & 0 \\ 0 & 0 & -1 & 2 & -1 & \ldots & 0 \\ \vdots& & \ddots & \ddots& & & \\ 0 & 0 & \ldots& & 0 & -1 & 2 \end{array}\right).\nonumber \end{eqnarray}$$

One can show, using Fourier analysis and the Perron-Frobenius theorem, that the system of equations is stable.

In this letter we have introduced and studied the dynamics of the selective dephasing model with three selective dephasing schemes. Accentuating the dephasing part of the dynamics allows analytical treatment of the long time dynamics.

A local protection from dephasing, as in the $\mathcal{A}_{\{ \xi \}}$ dephasing model, is not sufficient to keep quantum fingerprints at the steady state. The dephasing only at odd sites together with the Hamiltonian dynamics results in effectively conserved numbers. These new conserved numbers are witnesses of entanglement at the steady state at any length scale. Furthermore, the coarse grained dynamics is encapsulated in (9), a hydrodynamic non-local diffusion equation.

Note that the effective dynamics cannot predict the steady state values of the conserved numbers from the initial state. This is due to the microscopic time $t\lesssim \eta \nu_f/\varepsilon^2$, where there is no conservation. However, we stress that the effectively conserved numbers generally do not vanish if some bipartite entanglement is introduced in the initial condition (fig. 1 and the SM). We further note that the perturbation theory in (5) may break the quantum dynamical semigroup structure [41]. For Q = 1, we verify that effective evolution of the density matrix in the pointer state subspace satisfies a GKSL equation (1) according to $\partial_s \Pi_0 \rho = \frac{1}{\nu_f} \Pi_0 L_{H_{XX}} (\rho)$. Further evidence for the validity of the perturbation theory is provided in the SM.

The study of the effective dynamics of the selective dephasing model was restricted to its simplest case Q = 1. Onerous treatment is required for $Q\neq 1$, but similar results, i.e., diffusion and entanglement, are obtained (see the SM).

It is now clear that a coarse grained theory of diffusive open quantum systems must allow entanglement on all scales resulting in non-local evolution equations. Thus, such a theory will differ from the classical MFT, which is considered a local theory.

The present study, in particular the odd dephasing model, sets the stage for studying whether dynamics in diffusive quantum systems is indeed universal as its classical counterpart. Here, we considered the GKSL equation, which is an average equation. Following the degrees of freedom of the environment allows to inspect quantum trajectories, e.g., through quantum Langevin equations [42,43]. Finding whether the coarse grained description of the quantum trajectories is universal will be an important step towards a nonequilibrium theory of diffusive open quantum systems.

Recall that the selective dephasing model can be derived from a Hamiltonian dynamics through the singular coupling limit. Therefore, it is of interest to implement the selective dephasing model $\mathcal{A}_{\textrm{odd}}$ using cutting-edge experimental techniques [44–47]. How to take advantage of the entanglement propagation to all scales of the selective dephasing model is yet another open challenge.

Denis Bernard, Federico Carollo, Yaroslav Don, Juan P. Garrahan, Rob Jack, Toni Jin, Igor Lesanovsky, Katarzyna Macieszczak and Naftali Smith are acknowledged for fruitful discussions. The support of grant ANR-14-CE25-0003 is appreciated.