Charge and spin current statistics of the open Hubbard model with weak coupling to the environment

TL;DR

This paper extends previous work to analyze charge and spin current statistics in the open Hubbard model with multiple baths, revealing universal behaviors in nonequilibrium quantum systems.

Contribution

It generalizes the calculation of current statistics to the Hubbard model with multiple baths, up to next-to-leading order in weak coupling, indicating potential universality across models.

Findings

Charge and spin current statistics computed for the open Hubbard model.

Universal features observed in nonequilibrium behavior.

Extension of methods from the XXZ model to Hubbard model.

Abstract

Based on generalization and extension of previous work [Phys. Rev. Lett. {\bf 112}, 067201 (2014)] to multiple independent markovian baths we will compute the charge and spin current statistics of the open Hubbard model with weak system-bath coupling up to next-to-leading order in the coupling parameter. The physical results are related to those for the $XXZ$ model in the analogous setup implying a certain universality which potentially holds in this class of nonequilibrium models.

Equations

$$\frac{\mathrm{d}\rho (t)}{\mathrm{d}t}=-\mathrm{i}[H,\rho (t)]+\epsilon \left({\hat{\mathcal{D}}}^{\mathrm{jump}}\rho (t)+{\hat{\mathcal{D}}}^{\mathrm{diss}}\rho (t)\right),$$

$${\hat{\mathcal{D}}}^{\mathrm{jump}}\rho :=\sum _{\mu }{L}_m\rho L_m^{†},{\hat{\mathcal{D}}}^{\mathrm{diss}}\rho :=\frac{1}{2}\sum _m\{L_m^{†}{L}_m,\rho \},$$

$${\hat{\mathcal{D}}}_{\vec{\chi }}^{\mathrm{jump}}$$

$${\hat{\mathcal{D}}}^{\mathrm{diss}}$$

$${\hat{\mathcal{D}}}_{d,m,+}^{\mathrm{jump}}\rho$$

$${\hat{\mathcal{D}}}_{d,m,-}^{\mathrm{jump}}\rho$$

$${\hat{\mathcal{D}}}_{d,m,+}^{\mathrm{diss}}\rho$$

$${\hat{\mathcal{D}}}_{d,m,-}^{\mathrm{diss}}\rho$$

$$[M_d,L_{m_{\mu }}^{d,\nu }]=\nu L_{m_{\mu }}^{d,\nu }$$

$$\frac{\mathrm{d}\rho (t)}{\mathrm{d}t}={\hat{\mathcal{L}}}_{\vec{\chi }}\rho (t):=-\mathrm{i}[H,\rho (t)]+\epsilon \left({\hat{\mathcal{D}}}_{\vec{\chi }}^{\mathrm{jump}}\rho (t)+{\hat{\mathcal{D}}}^{\mathrm{diss}}\rho (t)\right).$$

$$\left[-\mathrm{i}\mathrm{ad}H+\epsilon \left({\hat{\mathcal{D}}}_{\vec{\chi }}^{\mathrm{jump}}+{\hat{\mathcal{D}}}^{\mathrm{diss}}\right)\right]\rho (\vec{\chi })=\lambda (\vec{\chi })\rho (\vec{\chi }),$$

$${\langle I^{m}_{d}\rangle}_{c}:=\lim_{t\to\infty}\frac{1}{t}{\langle[N_{d}(t)]^{m}\rangle}_{c}=\frac{\partial^{m}\lambda(\vec{\chi})}{\partial({\rm i}\chi_{d})^{m}}\Big{|}_{\vec{\chi}\to 0}.$$

$$\rho (\vec{\chi })=\sum _{p=0}^{\infty }(\mathrm{i}\epsilon {)}^p{\rho }^{(p)},\lambda (\vec{\chi })=\sum _{p=1}^{\infty }{\epsilon }^{2p-1}{\lambda }^{(2p-1)}.$$

$$(\mathrm{ad}H){\rho }^{(0)}$$

$$(\mathrm{ad}H){\rho }^{(1)}+{\hat{\mathcal{D}}}_{\vec{\chi }}{\rho }^{(0)}$$

$${\rho }^{(0)}=\sum _k{\alpha }_k{Q}_k.$$

$$-{\sigma }_{m_{+}}^{d,\mathrm{z}}+{\sigma }_{m_{-}}^{d,\mathrm{z}}\in \mathrm{Im}\mathrm{ad}H,$$

$$P{\sigma }_{m_{+}}^{d,\mathrm{z}}={\sigma }_{m_{-}}^{d,\mathrm{z}}P,\mathrm{or}P{\sigma }_{m_{+,-}}^{d,\mathrm{z}}=-{\sigma }_{m_{+,-}}^{d,\mathrm{z}}P,$$

$${\rho }^{(0)}=\frac{1}{(2d{)}^n}\text{mathbbm}1,$$

$${\lambda }^{(1)}=\sum _{d}2(-1+\cos ({\chi }_d)-\mathrm{i}{\mu }_d\sin ({\chi }_d)),$$

$$H=-2\sum _{s,x}(c_{s,j}^{†}{c}_{s,j+1}+c_{s,j+1}^{†}{c}_{s,j})+u\sum _j(2n_{\uparrow ,j}-1)(2n_{\downarrow ,j}-1)$$

$$c_{\uparrow ,j}=P_{j-1}^{({\sigma }^1)}{\sigma }_j^{1,-},c_{\downarrow ,j}=P_n^{({\sigma }^1)}{P}_{j-1}^{({\sigma }^2)}{\sigma }_j^{2,-}$$

$$H=2\sum _{j=1}^{n-1}\sum _{d=1}^2\left({\sigma }_j^{d,+}{\sigma }_{j+1}^{d,-}+{\sigma }_j^{d,-}{\sigma }_{j+1}^{d,+}\right)+\frac{u}{2}\sum _{j=1}^n{\sigma }_j^{1,\mathrm{z}}{\sigma }_j^{2,\mathrm{z}},$$

$$h_{j,j+1}=2\sum _{d=1}^2\left({\sigma }_j^{d,+}{\sigma }_{j+1}^{d,-}+{\sigma }_j^{d,-}{\sigma }_{j+1}^{d,+}\right)$$

$$+\frac{u}{2}\left({\sigma }_j^{1,\mathrm{z}}{\sigma }_j^{2,\mathrm{z}}+{\sigma }_{j+1}^{1,\mathrm{z}}{\sigma }_{j+1}^{2,\mathrm{z}}\right)$$

$$h_{j,j+1}^{{\sigma }^1}:=2{\sigma }_j^{1,+}{\sigma }_{j+1}^{1,-}+2{\sigma }_j^{1,-}{\sigma }_{j+1}^{1,+},$$

$$h_{j,j+1}^{{\sigma }^2}:=2{\sigma }_j^{2,+}{\sigma }_{2,j+1}^{-}+2{\sigma }_j^{2,-}{\sigma }_{j+1}^{2,+}$$

$${\mathbf{S}}^{+}=\sqrt{2}\sum _{k=0}^{\infty }\left(|k^{+}⟩⟨k+{\frac{1}{2}}^{+}|+|k+{\frac{1}{2}}^{-}⟩⟨k+1^{-}|\right),$$

$${\mathbf{S}}^{-}=\sqrt{2}\sum _{k=0}^{\infty }(-1{)}^k\left(|k+{\frac{1}{2}}^{+}⟩⟨k^{+}|+|k+1^{-}⟩⟨k+{\frac{1}{2}}^{-}|\right),$$

$$\displaystyle{\mathbf{S}}^{0}=\sum_{k=0}^{\infty}\bigl{(}{|2k^{+}\rangle}{\langle 2k^{+}|}+{|2k\!+\!\frac{1}{2}^{+}\rangle}{\langle 2k\!+\!\frac{1}{2}^{+}|}$$

$$\displaystyle\qquad+{|2k\!+\!1^{-}\rangle}{\langle 2k\!+\!1^{-}|}+{|2k\!+\!\frac{1}{2}^{-}\rangle}{\langle 2k\!+\!\frac{1}{2}^{-}|}\bigr{)}$$

$$+\varphi \sum _{k=1}^{\infty }\left(|2k-{\frac{1}{2}}^{+}⟩⟨2k-{\frac{1}{2}}^{+}|+|2k^{-}⟩⟨2k^{-}|\right),$$

$$\displaystyle{\mathbf{S}}^{\rm z}=\sum_{k=1}^{\infty}\bigl{(}{|2k\!-\!1^{+}\rangle}{\langle 2k\!-\!1^{+}|}+{|2k\!-\!\frac{1}{2}^{+}\rangle}{\langle 2k\!-\!\frac{1}{2}^{+}|}$$

$$\displaystyle\qquad+{|2k^{-}\rangle}{\langle 2k^{-}|}+{|2k\!+\!\frac{1}{2}^{-}\rangle}{\langle 2k\!+\!\frac{1}{2}^{-}|}\bigr{)}$$

$$+\varphi \sum _{k=0}^{\infty }\left(|2k+{\frac{1}{2}}^{+}⟩⟨2k+{\frac{1}{2}}^{+}|+|2k+1^{-}⟩⟨2k+1^{-}|\right),$$

$$\mathbf{X}$$

$$X_k(\varphi ,w)=\left(\begin{array}{cc}-(w+ku)w& 1-(w+ku)w(1-{\varphi }^2)\\ -kuw& 1-kuw(1-{\varphi }^2),\end{array}\right).$$

Full text

I Introduction

One of the key challenges of modern physics is understanding non-equilibrium behavior of both quantum and classical system EisertFriesdorfGogolin . In contrast to ordinary equilibrium statistical physics many approaches for understanding non-equilibrium physics exist, as there are many ways to drive a system from equilibrium. From quenches to open quantum systems one is interested often in similar concepts, such as transport of a certain quantity, that is the behavior of current of this quantity. The standard approach is to compute the expectation value of an observable (typically the current in the case of transport) with respect to some state of the system. In open quantum systems openbook this is usually the expectation value in the non-equilibrium steady state (NESS) to which the system evolves can evolve to in the infinite time limit. However, this approach is not the whole story and one may gain significantly more physical insights by studying the whole probability distribution of the observable. This is usually more difficult than computing just the expectation value.

This approach is additionally complicated in the quantum setting as measurements of an observable (such as the current) influence the time evolution of the system. One possible resolution to this issue is to look at a framework in which the non-equilibrium driving (via coupling to the baths) effectively performs the measurement as well. The framework of Markovian quantum master equations (Lindblad equation) is particularly well-suited for setting such a protocol. The study of transport in this context has attracted a lot of attention recently. For some very recent examples of such studies see Refs. Arenas1 ; Ilievski ; EverestLesanovsky ; ZnidScarVar ; ManzanoChuangCao ; KPS2 ; IlievskiZunkovic ; Prosen2 ; WolffSheikhanKollath and a recent review Prosenreview in the context of exact solutions.

The Lindblad master equation Lindblad for the system’s density operator $\rho(t)$ can be written as

[TABLE]

where $H$ is the Hamiltonian of the system generating the unitary part of the time evolution, and the dissipator $\hat{\cal D}$ implements the action of environment on the system. The latter is split into the quantum jumps and the dissipation

[TABLE]

where, as in the previous work Prosen1 ; countstatweak , the system-bath coupling parameter $\varepsilon$ is assumed to be small enough to allow the use of perturbation theory. The applicability of this method, namely the convergence radius of the perturbation series in this parameter is discussed in LemosProsen , where the authors point out issues of exponentially (in system size) decreasing size of the radius of convergence when one does perturbation for small $\varepsilon$. Thus, the methods relying on perturbation expansion in $\varepsilon$ cannot be applied in the thermodynamic limit, but can still be useful for finite system sizes.

Full counting statistics is the method which allows us to formulate in an elegant and efficient way the computation of cumulants of a current of a conserved quantity scully ; levitov ; rev , provided that certain conditions, elaborated in the next section hold. For some very recent results, as well as the related large deviation theory, in the open quantum framework in this area see, e.g., Refs. Garrahan ; ZnidaricFCS ; ZnidaricFCS2 ; ZnidaricFCS3 ; Monthus ; PigeonXuereb ; ManzanoHurtado ; MK ; flucthydro .

In this paper we will study the multi-bath generalization of the results obtained in Ref. countstatweak . Namely, we will show how a simple generalization of the full counting statistics for spin-1/2 systems weakly coupled to multiple baths can be achieved in the leading order in $\varepsilon$.

Afterwards, we will give a solution to the full counting statistics in next-to-leading (third) order for the open 1D (fermi) Hubbard model. The Hubbard model is a paradigm of many-body quantum physics. It describes numerous fundamental transport phenomena, in particular its two dimensional version is believed to capture superconductivity in cuprates. The 1D Hubbard model is integrable and exactly solvable through the Bethe ansatz BaxterBook ; EsslerBook ; LiebWu . By mapping the fermionic Hubbard model to a spin ladder we will find the full counting statistics of both the spin and charge currents in the third order of system-bath coupling.

The solution is based on the construction ProsenHubbard of the exact steady state density operator for the open maximally driven Hubbard model, which was later put into a Lax form in Ref. PopkovProsen . We show how to generate the solution to our counting problem using this infinite dimensional Lax operator. The solution, much like the previously discussed one for the $XXZ$ spin chain, is given by almost conserved charges which break the $\mathbb{Z}_{2}$ symmetry of the model (in this case the particle-hole symmetry) ProsenHubbardCharge .

II Full counting statistics of open systems with multiple baths

We aim at computing the full spin or charge exchange (and current) statistics in the limit of weak system-bath coupling when the Markovian approximation for the evolution of system’s density matrix $\rho(t)$ is appropriate openbook . In contrast to previous work countstatweak we will consider multiple, $D$ species of spins-1/2 labeled by a superscript $d=1,2\ldots D$ associated with the physical site $j$, i.e., $\sigma^{d,\pm}_{j}$, $\sigma^{d,{\rm z}}_{j}$, where these are the standard Pauli operators acting on a tensor product space $(\mathbb{C}^{2})^{\otimes Dn}$. Also, assume that our system is coupled to several sets of baths, i.e., in contrast to the previous work studied in countstatweak we allow for the system to be coupled to more than two baths. Let $H$ conserve one or more extensive observables with integer spectrum $M_{d}$, e.g. magnetization, charge, particle number etc, $[H,M_{d}]=0$. The subscript $d$ in $M_{d}$ denotes that $M_{d}$ is composed solely of spin operators of species $d$.

The Lindblad jump operators composing the dissipator, which are assumed to act ultra-locally (on one spin site only), represent the action of each bath on the system. For every spin species $d$ there exists one pair of baths corresponding to and driving the current of that spin species. We will then consider the flow (current) of observable $M_{d}$ from a subset of baths ${\cal I}$ (with the jump operator composed of one element in the pair of spin operators of species $d$ for every $d$, say $\sigma^{d,+}$) to its complement ${\cal O}$ (with the jump operator composed of the other element in the pair of spin operators of species $d$ for every $d$, say $\sigma^{d,-}$). We write $\tau_{\mu}:=1$ for $\mu\in{\cal I}$, and $\tau_{\mu}:=-1$ for $\mu\in{\cal O}$. Let us denote the amount of quantity $M_{d}$, say $\sum_{j}\sigma^{d,{\rm z}}_{j}$, transported in time $t$ from ${\cal I}$ to ${\cal O}$ by $N_{d}(t)$. In the open system’s framework, $N_{d}(t)$ is exactly the sum of $\delta_{\mu}$ for all the jumps for $\mu\in{\cal I}$ minus the sum of $\delta_{\mu}$ for all the jumps of $\mu\in{\cal O}$. In other words, we let the jump operators change $M_{d}$ by $\delta_{\mu}\in\mathbb{Z}$, and if they drive current of $M_{d}$ from ${\cal I}$ to ${\cal O}$ we say that the current is driven in the positive direction ($\tau_{\mu}=1$) and if the drive current of $M_{d}$ in the other direction then we say that the current is driven in the negative direction ($\tau_{\mu}=-1$).

In order to apply the counting field formalism we introduce now a counting vector field $\vec{\chi}=\{\chi_{1},\chi_{2}\,\ldots,\chi_{d},\ldots\}$ where ${\chi}_{d}$ counts, for spin species $d$, how much of charge $M_{d}$ was transferred during the time evolution of the system. For simplicity we will take $D$ symmetric current driving dissipators acting on ”complementary” sites of spin species $d$, each pair of sites being labeled by $m_{\nu}$ and each site in the pair by $\nu=\pm$. Then we may split the dissipator in the following manner,

[TABLE]

where,

[TABLE]

and,

[TABLE]

where the index $m$ corresponds to the pair of sites on which the dissipator acts, the sub-subscript, $\mu=\pm$ in $m_{\mu}$, denotes on which site in the pair of sites $m$ a Lindblad jump operator, $L_{m_{\mu}}^{d,\nu}=\sigma^{d,\nu}_{m_{\mu}}$, acts and $\nu=\pm$ in $L_{m_{\mu}}^{d,\nu}=\sigma^{d,\nu}_{m_{\mu}}$, denotes either a raising or lowering operator

[TABLE]

and $d$ is the species of spin on which the Lindblad jump operator acts. The sign of the product $\mu\nu=\pm$ determines whether a Lindblad operator $L_{m_{\mu}}^{d,\nu}$ drives the current of $M_{d}$ in the positive or negative direction, i.e., $\tau_{\mu}=1$ or $\tau_{\mu}=-1$, respectively.

With this modification the Lindblad equation (1) is deformed in the following way,

[TABLE]

The eigenvalue problem for the leading (with maximal real part) eigenvalue $\lambda(\vec{\chi})$ of this deformed Liouvillian ${\hat{\cal L}}_{\vec{\chi}}$ is given as,

[TABLE]

where $\rho(\vec{\chi})$ is the corresponding right eigenvectors, which will be assumed to be unique (as per the conditions of the Evans theorem Evans ). We have also defined the superoperator ${\,{\rm ad}\,}Hx:=[H,x]$ which acts on the space of operators. Note that $\rho(\vec{\chi}=0)$ corresponds to the non-equilibrium steady state of the undeformed Lindblad equation (1). We remark that the steady state can be also alternatively degenerate and symmetry reducible sym1 ; sym2 ; sym3 . Such a case is studied in Ref. ManzanoHurtado .

The leading eigenvalue $\lambda(\vec{\chi})$ in (9) will be the cumulant generating function for the currents $I_{d}$ of $M_{d}$,

[TABLE]

To give an intuitive picture of why this method works we will now give a short elaboration for only one spin species $D=1$, writing $\vec{\chi}\equiv\chi$. The method of full counting statistics can be understood most easily by observing that since $H$ commutes with the quantity of interest, and the Lindblad operators change the value of this quantity by well-defined (classical) amounts one may project the density matrix of the system $\rho(t)$ on a subspace of $N$ spin transfers between the pair of baths in time $t$ which we denote as $\rho_{N}(t)$. Thus, the probability of $N$ spin transfers in time $t$ is given as $P_{N}(t)={{\rm tr}}\rho_{N}(t)$. We can then perform a Fourier transform in $N$ and then find that the dissipator in the original Lindblad equation (1) is modified to include the counting field in the form written in Eq. (5). Then in the long-time limit where $\rho(\chi,t)\approx e^{\lambda(\chi)t}\rho(\chi,t=0)$, the leading (with largest real part) eigenvalue $\lambda(\chi)$ is the cumulant generating function for the current (10).

The perturbation expansion is given by

[TABLE]

The first two orders obtained by plugging in (11) in are,

[TABLE]

where we also defined the total dissipator as $\hat{\cal D}_{\vec{\chi}}:=\hat{\cal D}^{\rm jump}_{\vec{\chi}}+\hat{\cal D}^{\rm diss}$. In general any linear combination of operators which commute with $H$ satisfy the zeroth-order equation (12). So we may write in general this solution using a linearly independent combination of $Q_{k}$ such that $[H,Q_{k}]=0$,

[TABLE]

Note that since that we have assumed that the corresponding right eigenvector is unique not every possible solution to the zeroth-order equation (12) can be the correct leading order in the perturbation expansion (11). In order to find the correct one, we also demand Prosen1 ; countstatweak that the solution to the first-order equation (13) exists. To do this we will demand a generalization of the simple algebraic condition on the image space of ${\,{\rm ad}\,}H$ (see countstatweak ) to hold,

[TABLE]

for every spin species $d$ and for every pair of sites $m$ on which the Lindblad jump operators act. This requirement will be fulfilled if a generalization of the previous parity symmetry condition studied in countstatweak holds. That is, an operator $P$ ($P^{2}=\mathbbm{1})$ should exist, which satisfies either

[TABLE]

and commutes with the Hamiltonian $[P,H]=0$ together with all the conserved operators $[P,Q_{k}]=0$. The proof of this statement is the same as in Ref. countstatweak : Since ${\,{\rm ad}\,}H$ is Hermitian its image is orthogonal w.r.t. Hilbert-Schmidt inner-product to all $Q_{k}$. Thus the condition (15) is equivalent to: ${{\rm tr}}(-\sigma^{d,{\rm z}}_{1}+\sigma^{d,{\rm z}}_{n})Q_{k}=0,\forall Q_{k}$. From this we may directly see that ${{\rm tr}}(-\sigma^{d,{\rm z}}_{m_{+}}+\sigma^{d,{\rm z}}_{m_{-}})Q_{k}={{\rm tr}}P(-\sigma^{d,{\rm z}}_{m_{+}}+\sigma^{d,{\rm z}}_{m_{-}})Q_{k}P=-{{\rm tr}}(-\sigma^{d,{\rm z}}_{m_{+}}+\sigma^{d,{\rm z}}_{m_{-}})Q_{k}$. From which follows that ${{\rm tr}}(-\sigma^{d,{\rm z}}_{m_{+}}+\sigma^{d,{\rm z}}_{m_{-}})Q_{k}=0$.

It is then trivial to show that given,

[TABLE]

and,

[TABLE]

a solution to the first-order equation (13) exists. This can be done by plugging the expressions (17) and (18) into (13) and showing that (15) holds.

By performing a suitable transformation $\vec{\chi}\to-{\rm i}\vec{\chi}$ we switch to (provided there are no issues with analyticity) what is called large-deviation theory Oono ; Touchette .

III Charge and spin statistics of the open Hubbard model at weak system-bath coupling

The first-order result of the previous section was universal and only depended on the local nature of jump operators and existence of the parity symmetry. In this section we show that one can compute higher order in $\varepsilon$ for a multispecies integrable model, namely the 1D Hubbard model. The fermionic Hubbard model with $n$ sites is given by the following Hamiltonian,

[TABLE]

where $c_{s,j}$ are the canonical Fermi operators on site $j$ and $j\in\{1\ldots n\}$, $s\in\{\uparrow,\downarrow\}$ is the spin of each fermion, $u$ is a non-dimensional parameter determining the interaction strength, and $n_{s,j}=c^{\dagger}_{s,j}c_{s,j}$ is the local charge density.

We now apply the Wigner-Jordan transformation,

[TABLE]

where $P^{(\sigma^{1})}_{j}:=\sigma_{1}^{1,\rm z}\sigma_{2}^{1,{\rm z}}\cdots\sigma_{j}^{1,\rm z}$, $P^{(\sigma^{2})}_{j}:=\sigma_{1}^{2,\rm z}\sigma_{2}^{2,{\rm z}}\cdots\sigma_{j}^{2,\rm z}$, and of course $\sigma^{1,0}_{j}\equiv\sigma^{2,0}_{j}\equiv\mathbbm{1}$. Thus we transform the fermionic Hubbard model into a spin ladder,

[TABLE]

We will also likewise define two local ”kinetic energy” densities for the two species of spin,

[TABLE]

Note that the Hubbard Hamiltonian conserves the numbers of spin-up and spin-down electrons, $M_{\sigma^{1}}=\sum_{j=1}^{n}\frac{1}{2}(\sigma^{1,{\rm z}}_{j}+\mathbbm{1})$, $M_{\sigma^{2}}=\sum_{j=1}^{n}\frac{1}{2}(\sigma^{2,{\rm z}}_{x}+\mathbbm{1})$, $[H,N_{\sigma^{1,2}}]=0$.

We now continue for sake of being self-contain by quickly reviewing the results of Refs. ProsenHubbard and PopkovProsen on the open Hubbard model where it has been shown that there exists an infinite dimensional Lax operator seemingly unrelated to the known one discovered by Shastry Shastry .

All the following definitions are given as in Ref. PopkovProsen . Define ${\mathbf{S}},\acute{{\mathbf{S}}},\grave{{\mathbf{S}}},{\mathbf{T}},\acute{{\mathbf{T}}},\grave{{\mathbf{T}}}\in{\rm End}({\cal H}_{\rm a}\otimes{\cal H}_{\rm p})$, and ${\mathbf{X}},{\mathbf{Y}}\in{\rm End}({\cal H}_{\rm a})$ (scalars over ${\cal H}_{\rm p}$), where ${\cal H}_{\rm a}$ is an infinite dimensional auxiliary space. We use labelling the basis vectors of the auxiliary space as ${\cal V}=\{0^{+},\frac{1}{2}^{+},\frac{1}{2}^{-},1^{-},1^{+},\frac{3}{2}^{+},\frac{3}{2}^{-},2^{-},2^{+}\ldots\}$ so thatfootnote ${\cal H}_{\rm a}={\rm lsp}\{{|v\rangle};v\in{\cal V}\}$. Introduce another definition, namely the components ${\mathbf{S}}^{s},{\mathbf{T}}^{s}\in{\rm End}({\cal H}_{\rm a})$, ${\mathbf{S}}=\sum_{s\in\{+,-0,{\rm z}\}}{\mathbf{S}}^{s}\sigma^{1,s}$, ${\mathbf{T}}=\sum_{s\in\{+,-0,{\rm z}\}}{\mathbf{T}}^{s}\sigma^{2,s}$ and likewise for $\acute{{\mathbf{S}}},\grave{{\mathbf{S}}},\acute{{\mathbf{T}}},\grave{{\mathbf{T}}}$.

Then ${\mathbf{S}}$ is given by the following matrix representations,

[TABLE]

where $\phi\in\mathbb{C}$ is a free parameter. Defining an operator ${\mathbf{G}}$ which interchanges spin species in physical space, i.e., ${\mathbf{G}}\sigma^{1,s}{\mathbf{G}}=\sigma^{2,s}$, ${\mathbf{G}}\sigma^{2,s}{\mathbf{G}}=\sigma^{1,s}$ and in the auxiliary space operates as ${\mathbf{G}}{|k^{\pm}\rangle}:={|k^{\pm}\rangle}$, ${\mathbf{G}}{|k\!+\!\frac{1}{2}^{\pm}\rangle}:={|k\!+\!\frac{1}{2}^{\mp}\rangle}$, $k\in\mathbb{Z}^{+}$, ${\mathbf{T}}$ is given by ${\mathbf{T}}={\mathbf{G}}{\mathbf{S}}{\mathbf{G}}$.

The operator ${\mathbf{X}}$ is given by,

[TABLE]

where $X_{k}=\{X^{\nu,\nu^{\prime}}_{k}\}_{\nu,\nu^{\prime}\in\{-,+\}}$ are $2\times 2$ matrices

[TABLE]

and $w\in\mathbb{C}$ is another free parameter. The operator ${\mathbf{Y}}$ is defined as,

[TABLE]

Demanding that ${\mathbf{X}}$ is invertible, that is $w\neq 0$, $\det X_{k}\neq 0$, we may define implicitly,

[TABLE]

and finally

[TABLE]

With these definitions (i.e., Eqs. 24 -29) it was shown in PopkovProsen that the following relations hold (recall the definition in (23)),

[TABLE]

where the subscripts in ${\mathbf{S}}_{j}$ (${\mathbf{T}}_{j}$) indicate that they act locally on spin site $j$ of spin species 1 (2).

Define a Lax operator ${\mathbf{L}}$ and another operator which will play the role of its ”derivative” $\widetilde{{\mathbf{L}}}$,

[TABLE]

Then using the above relations, Eqs. 30 - 34 and the definition in (23), it can be shown that the so-called Sutherland-Shastry relation holds,

[TABLE]

where $h_{j,j+1}$ is the local energy density operator given by (21).

Now define a central operator,

[TABLE]

Through this operator we may calculate a quadratically extensive ProsenHubbardCharge almost conserved charge,

[TABLE]

By deriving the Sutherland-Shastry relation (37) and using definition (38) it can be shown that this operator $Z$ is indeed almost conserved,

[TABLE]

Define an operator $S$ which is a product of spin species 2 - flip $P=\prod_{j=1}^{n}\sigma^{2,\rm{x}}_{j}$ and the transformation $u\to-u$, $S^{2}=\mathbbm{1}$. It is obvious that $[H,S]=0$. Let,

[TABLE]

and then by multiplying (40) from the left and the right by $S$ we find,

[TABLE]

This completes the review of the work PopkovProsen we need to continue and now we can finally state and solve our problem. We will study driving of the Hubbard with a total of eight Lindblad operators,

[TABLE]

where $d=1,2$ designates either spin species 1 or 2 with appropriate driving parameters $\mu_{1,2}$. The rest of the setup is defined in the previous section and the perturbation expansion is the same as before (11).

The Hubbard Hamiltonian fulfils the parity symmetry requirement (16) with $P=\prod_{d=1}^{2}\prod_{j=1}^{n}\sigma^{d,\rm{x}}_{j}$ being a global spin flip operator in the $\rm{z}$-basis.

Clearly, the zeroth order solution is $\rho^{(0)}=4^{-n}\mathbbm{1}$ and the first order solution is given through the conserved charges in (39). By applying the total dissipator $\hat{\cal D}_{\vec{\chi}}\rho^{(0)}$ and using the algebraic condition (15) we find that the solution to the first order equation (13) is,

[TABLE]

where,

[TABLE]

Moving on to the second order equation,

[TABLE]

we recall that like for the $XXZ$ spin chain Prosen1 $Z_{1,2}$ contain no terms with factors $\sigma_{j}^{d,{\rm z}}$ ($d=1,2$). We then only need to check (47) for $Z_{d}\in\{\mathbbm{1},\sigma^{d,+}_{1,n},\sigma^{d,-}_{1,n}\}$. We take a quadratic ansatz for the second order solution,

[TABLE]

The calculation for the second-order coefficients is now straightforward, but tedious, and we can find the expressions for the coefficients $c^{(2)}_{j,k},c^{(2)}_{j^{\dagger},k},c^{(2)}_{j,k^{\dagger}},c^{(2)}_{j^{\dagger},k^{\dagger}}.$ They are quite long and not directly relevant for the rest of the discussion so we omit them here and instead will write the explicit expression in Appendix  V.

Like before this determines $\rho^{(2)}$ up to the addition of conserved quantities $Q_{k}$, where $[H,Q_{k}]=0$, that is, $\rho^{(2)^{\prime}}=\rho^{(2)}+\sum_{k}\alpha_{k}Q_{k}$. However, like in countstatweak they are irrelevant because via the third order equation,

[TABLE]

we find that,

[TABLE]

and due to the form of the dissipator only terms of the form $O\in\{\mathbbm{1},(-\sigma^{d,{\rm z}}_{1}+\sigma^{d,{\rm z}}_{n})\}$ in $\rho^{(2)^{\prime}}$ can contribute to (50). Now since $(-\sigma^{d,{\rm z}}_{1}+\sigma^{d,{\rm z}}_{n})\in{\rm Im}{\,{\rm ad}\,}H$ for $d=1,2$ due to the existence of a solution to the first order equation (13) merely $\rho^{(2)}$ is sufficient to determine the third order correction,

[TABLE]

The full form of $\lambda^{(3)}$ is long and is given as Eq. (69) in the appendix. Now let us finally turn to computing the third-order correction to the current cumulants. Note that there are two physically relevant currents we can study, first the charge current (corresponding to setting \chi=\chi_{1}=\chi_{2}\) and secondly the spin current (corresponding to setting \chi=\chi_{1}=-\chi_{2}\). In the first case we have for charge current for odd and even $n$,

[TABLE]

The spin current cumulants $\left<I^{n}_{(3),\rm{sp}}\right>_{\rm{c}}$ can be obtained from the charge current ones simply by setting $\mu_{1}\to-\mu_{1}$ in (52). We plot them in Fig.  1.

It is interesting to note that if we set $\mu=\mu_{1}=\mu_{2}$ the charge current cumulants become quite similar to the ones for the $XXZ$ spin chain countstatweak . They only differ up to a constant replacing $u$ in (52). This reflects potentially a deep symmetry connection between the two models. Indeed, since the derivation only depended on the fact that $Z_{1,2}$ contain only terms of the form $Z_{d}\in\{\mathbbm{1},\sigma^{d,+}_{1,n},\sigma^{d,-}_{1,n}\}$ on sites $1$ and $n$ this could hint at a possible universal result for this class of integrable models.

We can also calculate the charge and spin current connected correlations using the explicit expression for $\lambda^{(1)}$ and $\lambda^{(3)}$ (eqs. (18), and (69) in the appendix, respectively) by using a generalization of (10) for multiple counting fields,

[TABLE]

where in the case of two counting fields $d_{k}=1,2$ and we will again transform into the spin-charge picture using $\chi_{1}=(\chi_{\rm{ch}}+\chi_{\rm{sp}})/2$ and $\chi_{2}=(\chi_{\rm{ch}}-\chi_{\rm{sp}})/2$. The first few are,

[TABLE]

and,

[TABLE]

The general expressions for arbitrary order correlations between two currents ate given in appendix. One can also calculate higher order connected correlations from the expression for $\lambda^{(1)}$ and $\lambda^{(3)}$ (eqs. (18) and (69) in the appendix, respectively, give the two-point connected correlation functions).

IV Conclusion

In this paper we generalized the results of Ref. countstatweak for the full counting statistics of parity symmetric spin systems in the weak-system bath coupling regime to multiple baths and spin species in the leading order. Then later we calculated the the next-to-leading order for the charge and spin current statistics for the case of the Hubbard model in the next-to-leading order, as well as the charge-spin correlations. The form of the charge current cumulants are similar to those in the $XXZ$ case countstatweak indicating a possible universal feature of these class of integrable models. Based on recent work LemosProsen likely do not hold in the thermodynamic limit due to issues with this perturbation expansion, however they should be applicable to finite size system. One would thus desire a fully non-perturbative extension for the $XXZ$ and the Hubbard model, which currently remain elusive.

V Appendix: The second-order correction to the full counting statistics of the open Hubbard model

V.1 Explicit coefficients of the second order correction to the counting eigenvector

Here we will write out the coefficient in the second-order correction to the eigenvector corresponding to the leading eigenvalue of the deformed Liouvillian of the open Hubbard model Eq. (48) and the full form of the third-order correction to the current cumulant generating function.

The coefficients multiplying operators which are not conjugate transposed are,

[TABLE]

The mixed ones are,

[TABLE]

and,

[TABLE]

Continuing,

[TABLE]

while the rest of coefficients are zero. Note that several pairs are similar in form, though not exactly the same.

The full form of the third-order correction to the cumulant generating function is given as,

[TABLE]

V.2 General spin-charge current correlations in the driven Hubbard model

We finally give the full form of the spin-charge current correlations. First the odd-odd ones,

[TABLE]

then the even-odd ones,

[TABLE]

the odd-even ones can be obtained by setting $\mu_{2}\to-\mu_{2}$ and interchanging $m$ and $k$, and the even-even ones are,

[TABLE]