Convergence radius of perturbative Lindblad driven non-equilibrium steady states

TL;DR

This paper investigates the convergence properties of perturbative expansions for non-equilibrium steady states in Lindblad-driven spin chains, revealing how convergence radius varies with system size and parameters.

Contribution

Develops a systematic formal approach to compute perturbative expansions of non-equilibrium steady states in open quantum spin chains, analyzing convergence behavior in different regimes.

Findings

Convergence radius shrinks with increasing system size in coupling expansion.

Convergence radius exceeds 1 in driving parameter expansion.

Radius approaches 1 as anisotropy shifts from XY to Ising regime.

Abstract

We address the problem of analyzing the radius of convergence of perturbative expansion of non-equilibrium steady states of Lindblad driven spin chains. A simple formal approach is developed for systematically computing the perturbative expansion of small driven systems. We consider the paradigmatic model of an open $XXZ$ spin 1/2 chain with boundary supported ultralocal Lindblad dissipators and treat two different perturbative cases: (i) expansion in system-bath coupling parameter and (ii) expansion in driving (bias) parameter. In the first case (i) we find that the radius of convergence quickly shrinks with increasing the system size, while in the second case (ii) we find that the convergence radius is always larger than $1$, and in particular it approaches $1$ from above as we change the anisotropy from easy plane ($XY$) to easy axis (Ising) regime.

Tables

Table 1. Table 1: N 𝑁 N dependence of the critical index n 0 subscript 𝑛 0 n_{0} values for different scenarios: (A) perturbation in ε 𝜀 \varepsilon parameter, with anisotropy Δ > 0 Δ 0 \Delta>0 , (B) perturbation in ε 𝜀 \varepsilon parameter, with anisotropy Δ = 0 Δ 0 \Delta=0 , and (C) perturbation in μ 𝜇 \mu parameter, with anisotropy Δ > 0 Δ 0 \Delta>0 .

$N$	2	3	4	5	6
$dim\mathscr{H}=2^{2N}$	16	64	256	1024	4096
(A): $\epsilon$-perturbation, $\Delta >0$	2	6	26	98	N/A
(B): $\epsilon$-perturbation, $\Delta =0$	2	4	6	8	10
(C): $\mu$-perturbation, $\Delta >0$	2	4	12	36	N/A

Equations

$$H=\sum _{j=1}^{N-1}{\sigma }_j^x{\sigma }_{j+1}^x+{\sigma }_j^y{\sigma }_{j+1}^y+\Delta {\sigma }_j^z{\sigma }_{j+1}^z,$$

$$\frac{d\rho (t)}{dt}=-i\left[H,\rho (t)\right]+\sum _{k}2L_k\rho (t)L_k^{†}-\left\{L_k^{†}{L}_k,\rho (t)\right\},$$

$$-i\left(\mathrm{ad}H\right){\rho }_{\infty }+\epsilon \hat{\mathcal{D}}{\rho }_{\infty }=0,$$

$$\hat{\mathcal{D}}:=\frac{1}{2}\left(1+\mu \right){\hat{\mathcal{D}}}_{+}+\frac{1}{2}\left(1-\mu \right){\hat{\mathcal{D}}}_{-},$$

$${\hat{\mathcal{D}}}_{\pm }\rho :=2{\sigma }_1^{\pm }\rho {\sigma }_1^{\mp }-\left\{{\sigma }_1^{\mp }{\sigma }_1^{\pm },\rho \right\}+2{\sigma }_N^{\mp }\rho {\sigma }_N^{\pm }-\left\{{\sigma }_N^{\pm }{\sigma }_N^{\mp },\rho \right\}.$$

$${\rho }_{\infty }=\sum _{n=0}^{\infty }(i\epsilon {)}^n{\rho }^{(n)}={\rho }^{(0)}+\sum _{n=1}^{\infty }(i\epsilon {)}^n{\rho }^{(n)},$$

$$\left(\mathrm{ad}H\right){\rho }^{(n)}=-\hat{\mathcal{D}}{\rho }^{(n-1)},\forall n\ge 1,$$

$$\hat{\mathcal{D}}\rho^{(n_{0})}=\sum_{j=1}^{n_{0}}c_{j}\Big{(}\hat{\mathcal{D}}\rho^{(j-1)}\Big{)}\ .$$

$${\rho }^{(n_0+1)}=\sum _{j=1}^{n_0}{c}_j{\rho }^{(j)}.$$

$$(\mathrm{ad}H){\rho }^{(n_0+1)}$$

$${\rho }^{(n_0+2)}=\sum _{j=1}^{n_0}{c}_j{\rho }^{(j+1)}=\sum _{j=2}^{n_0+1}{c}_{j-1}{\rho }^{(j)},$$

$${\rho }^{(n_0+k)}=\sum _{j=1}^{n_0}{c}_j{\rho }^{(j+k-1)}=:\sum _{j=1}^{n_0}{R}_j^{(n_0+k)}{\rho }^{(j)},$$

$$R_j^{(n_0+2)}=R_{j-1}^{(n_0+1)}+c_j{R}_{n_0}^{(n_0+1)},\forall 1\le j\le n_0,$$

$$R_j^{(n_0+k+1)}=R_{j-1}^{(n_0+k)}+c_j{R}_{n_0}^{(n_0+k)},\forall 1\le j\le n_0,\forall k\ge 1,$$

$$M=\left(\begin{array}{ccccc}0& & & & c_1\\ 1& 0& & & c_2\\ & \ddots & \ddots & & ⋮\\ & & 1& 0& c_{n_0-1}\\ & & & 1& c_{n_0}\end{array}\right),$$

$${\mathbf{R}}^{(n_0+k)}=M^k{\mathbf{R}}^{(n_0)},\forall k\ge 1.$$

$${\mathbf{R}}^{(n)}=M^{n-1}{\mathbf{R}}^{(1)},\forall n\ge 1.$$

$${\mathbf{R}}^{(\infty )}=\sum _{n=1}^{\infty }(i\epsilon {)}^n{\mathbf{R}}^{(n)}=\left(\sum _{n=1}^{\infty }(i\epsilon {)}^n{M}^{n-1}\right){\mathbf{R}}^{(1)}=i\epsilon (\text{mathbbm}{1}_{n_0}-i\epsilon M{)}^{-1}{\mathbf{R}}^{(1)},$$

$${\rho }_{\infty }=\rho (0)+\sum _{j=1}^{n_0}{\mathbf{R}}_j^{(\infty )}{\rho }^{(j)}.$$

$$\lambda ={\left(\underset{1\le j\le n_0}{\max }\{|{\lambda }_j|\}\right)}^{-1},$$

$$\epsilon \hat{\mathcal{D}}=\frac{\epsilon }{2}\left({\hat{\mathcal{D}}}_{+}+{\hat{\mathcal{D}}}_{-}\right)+\mu \frac{\epsilon }{2}\left({\hat{\mathcal{D}}}_{+}-{\hat{\mathcal{D}}}_{-}\right)=:{\hat{\mathcal{D}}}_0+\mu {\hat{\mathcal{D}}}_{\mu },$$

$$-i\Big{[}\left(\mathrm{ad}H\right)+i\hat{\mathcal{D}}_{0}\Big{]}\rho_{\infty}+\mu\hat{\mathcal{D}}_{\mu}\rho_{\infty}=0\ ,$$

$${\rho }_{\infty }=\sum _{n=0}^{\infty }(i\mu {)}^n{\rho }_{\mu }^{(n)},$$

$$T_{\mu }{\rho }_{\mu }^{(n)}=-{\hat{\mathcal{D}}}_{\mu }{\rho }_{\mu }^{(n-1)},\forall n\ge 1,$$

$${\lambda }_{\mu }={\left(\underset{1\le j\le n_0}{\max }\{|({\lambda }_{\mu }{)}_j|\}\right)}^{-1},$$

Full text

I Introduction

Open quantum systems approach Breuer and Petruccione (2002) has generated a great deal of interest in recent years, in particular since the theory is not only able to efficiently describe subsystems of large quantum systems in thermal equilibrium, but also captures non-equilibrium physics of systems driven out of and even far from equilibrium. For example, one can efficiently model the coherent quantum transport problem in terms of a Lindblad equation with incoherent quantum jump processes limited to the ends of the chain Wichterich et al. (2007). It has been shown that in the case when the bulk dynamics is completely integrable, one can often write exact solutions for the many-body density matrix of the boundary driven steady state Prosen (2015) and these, in turn, generate interesting new non-equilibrium physics, such as quasilocal conservation laws Ilievski et al. (2016), persistent currents, etc. Remarkably, this rather formal approach resolved a long lasting debate (starting with the work of Zotos and collaborators. Castella et al. (1995)) on existence of anomalous (ballistic) spin transport at high-temperatures Prosen (2011a), a problem of not only theoretical but also experimental interest Heidrich-Meisner et al. (2007); Hlubek et al. (2010).

However, exact analytic solutions for the non-equilibrium steady state are only possible for a rather limited set of Lindblad dissipators, such as for example, pure magnetic source on one end, and pure magnetic sink on the other end Prosen (2015). In cases of more general and generic boundaries, one may attempt to consider formal perturbative expansions of the steady state in a parameter which breaks the integrability (solvability).

In this paper, we address the problem of calculating the radius of convergence of such perturbative expansions for a general finite open quantum system. We consider a specific and widely studied case of open $XXZ$ spin $1/2$ chain with magnetic pump boundaries which depend on two parameters, (i) the system-bath coupling parameter, and (ii) the driving (bias) parameter which generates the magnetic current in the non-equilibrium steady state (NESS). We demonstrate that, when expanding in the coupling parameter, the case in which the first two orders are explicitly and analytically known Prosen (2011a), the convergence radius of the series quickly, probably exponentially, shrinks with increasing the chain length, so the solution may have little relevance for the correctly scaled thermodynamic limit. On the other hand, if expanding in the driving parameter, around the so-called linear response regime, explicit computations strongly suggest that the radius of convergence remains finite and bounded by $1$ from below, so this perturbative expansion should be well relevant even in the thermodynamic limit.

The paper is structured as follows. In section II we develop a systematic approach for calculating the radius of convergence for perturbative expansion of NESS for finite open quantum systems. In sections III and IV we present numerical results on the example of boundary driven $XXZ$ chain, for coupling and driving expansions, respectively, and conjecture the asymptotic behaviors.

II Evaluating the radius of convergence

We study an open anisotropic Heisenberg $XXZ$ spin 1/2-chain with a constant nearest-neighbor interaction, given by the Hamiltonian

[TABLE]

where $\sigma_{j}^{x,y,z}$, for $j=1,\ldots,N$ are Pauli operators acting on a Hilbert tensor product space $\mathcal{H}=(\mathbb{C}^{2})^{\otimes N}$. We also use the standard notation $\sigma^{0}=\mathbbm{1}_{2}$. The dynamics of the system is given by a Markovian master equation in the Lindblad form

[TABLE]

and in this case we have Lindblad operators mimicking magnetic reservoirs acting only at the boundaries of the chain, namely $L_{1,2}=\sqrt{\frac{1}{2}\varepsilon(1\pm\mu)}\ \sigma_{1}^{\pm}$ and $L_{3,4}=\sqrt{\frac{1}{2}\varepsilon(1\mp\mu)}\ \sigma_{N}^{\mp}$. We stress here the different roles of the two parameters in the Lindblad operators: $\varepsilon$ is the strength of coupling of the spin baths at the boundaries of the chain, while $|\mu|\leq 1$ tells us how strong is the non-equilibrium driving force acting at the edges of the system. If $\mu=0$, spin baths acting on $j=1$ and $j=N$ are symmetric, and the system will reach an equilibrium state after a sufficient long time (in our case, this will be an infinite temperature Gibbs state – completely mixed state). On the other hand, if $\mu=\pm 1$, we have maximal non-equilibrium driving force, and the analytical solution for the stationary state is well known in this case Prosen (2011b); Karevski et al. (2013).

The model is exactly solvable using a Bethe ansatz Bethe (1931); Korepin et al. (1997), but this closed system solution no longer applies when it is driven far out of equilibrium with a dynamics given by a Lindblad equation. In two papers Prosen (2011a, b), one of us made a progress and calculated analytically some physical quantities for the driven $XXZ$ model. In the present paper, we initially address to Prosen (2011a), where the author formally expanded the solution as a perturbative series in $\varepsilon$, related to the strength of the coupling between the baths and the boundaries of the chain, but as it has been already noted there, the “convergent properties of perturbation series are unknown”. For sake of understanding, we recall the main steps of such construction. First, the NESS density operator is given by $\rho_{\infty}=\lim_{t\to\infty}\rho(t)$, so it is a fixed point for the dynamics (2)

[TABLE]

where $\left(\mathrm{ad}H\right)\rho:=\left[H,\rho\right]$, and the dissipator here is defined as

[TABLE]

with

[TABLE]

The NESS solution for this regime of weak coupling is then formally expanded as a series in $\varepsilon$

[TABLE]

and naïvely replacing the series above into Eq. (3) we obtain a recurrence relation

[TABLE]

which allows us to evaluate the terms of the sequence $\{\rho^{(n)}\}$ given the initial condition $\rho^{(0)}=2^{-N}\mathbbm{1}_{2^{N}}$.

In Prosen (2011a), the author evaluated the first and second order terms, but that is not the point here. He raises three questions about some possible mathematical problems in constructing this sequence: (i) it is not clear if recurrence relation (6) should have unique solutions for each $n$. (ii) It is not also clear if it is always possible, for each $n$, to find out the next term $\rho^{(n+1)}$. In order to do so, we must prove that $\hat{\mathcal{D}}\rho^{(n)}\in\mathrm{Im}\ (\mathrm{ad}H)$, for any $n$. (iii) We do not know the radius of convergence for (5). We stress now that we no longer concern about (i) and (ii), since the uniqueness of NESS is guaranteed by Evans theorem Evans (1977) – the application of Evans theorem to the present model is well discussed in Ref. Prosen (2015). So the only open problem left is (iii).

We start constructing the sequence term by term using the recurrence relation (6), up to $\rho^{(n_{0})}$, where $n_{0}$ is the unknown index such that, for the first time, the vector $\hat{\mathcal{D}}\rho^{(n_{0})}$ can be written as a linear combination of the previous ones, i.e.

[TABLE]

Just to make it clearer, we emphasize now the crucial definition of $n_{0}$: starting from the first term of the series (5), namely the vector $\rho^{(0)}$, we use recurrence relation (6) to find out the next one, $\rho^{(1)}$. Once it was found, we check if vectors $\hat{\mathcal{D}}\rho^{(0)}$ and $\hat{\mathcal{D}}\rho^{(1)}$ are linearly dependent: if so, we have found that $n_{0}=1$ Ker ; otherwise we keep using Eq. (6) again to find next terms of the series, one by one, until we find a $\rho^{(n_{0})}$ such that, as we did on (7), the vector $\hat{\mathcal{D}}\rho^{(n_{0})}$ may be written as a linear combination of $\hat{\mathcal{D}}\rho^{(0)},\ldots,\hat{\mathcal{D}}\rho^{(n_{0}-1)}$, which are linearly independent by definition. Of course, once fixed the system size $N$, the Hilbert space is finite with $\dim\mathcal{H}=2^{2N}$, so we always will have such index $n_{0}\leq 2^{2N}-1$ – this minus 1 refers to the first term of the series $\rho^{(0)}$, which is orthogonal to all other vectors $\rho^{(n)}$ Ker . But $n_{0}$ is in practice (much) smaller than the dimension of the Hilbert space, we list some values in section III, see table 1. If we manage to find the coefficients $c_{j}$ of this linear combination above, then we can obtain any further term of the sequence. We claim that the next term is given by

[TABLE]

Indeed, if we replace it into the recurrence relation (6), we get

[TABLE]

so it is the solution for the recurrence relation in this step. Moreover, we can easily prove that $\hat{\mathcal{D}}\rho^{(n_{0}+1)}\in\mathrm{Im}\ (\mathrm{ad}H)$, since the Range set is a subspace of $\mathcal{H}$, so we are able to keep constructing the sequence term by term. The next one is found using the same coefficients,

[TABLE]

and we can easily show that this is the solution following the same steps that we have just used in Eq. (8).

We could do this repeatedly from now on and get all the terms of (5) up to any desired order $n>n_{0}$, but the critical step now is to rewrite them as a linear combination of the elements of the set $\mathcal{B}_{n_{0}}=\left\{\rho^{(1)},\ldots,\rho^{(n_{0})}\right\}$. In other words, although $\mathcal{B}_{n_{0}}$ is not a basis for the whole Hilbert space $\mathcal{H}$, the NESS solution is essentially in the subspace $\mathcal{H}_{n_{0}}$ spanned by this set, except for its zero-th order term $\rho^{(0)}$. Precisely, we have \big{(}\rho_{\infty}-\rho^{(0)}\big{)}\in\mathcal{H}_{n_{0}}. We now aim to rewrite each term as a linear combination of these vectors

[TABLE]

for any $k\geq 1$, where in the last equality in Eq. (11) above we have just defined the coefficients $R_{j}^{(n_{0}+k)}$ for our new $n_{0}$-dimensional vector $\mathbf{R}^{(n_{0}+k)}$. For $k=1$ is trivial to see that $R_{j}^{(n_{0}+1)}=c_{j}$. For $k=2$, we can match equations (10) and (8) to show that

[TABLE]

where by convention we have defined $R_{0}^{(n_{0}+1)}=0$. Using an analogous procedure, we can straightforward show that the relation above holds for any $k\in\mathbb{N}$, i.e.

[TABLE]

where once more we have defined, by convention, $R_{0}^{(n_{0}+k)}=0$.

We can now compactly rewrite it for $n\geq n_{0}$ as $\mathbf{R}^{(n+1)}=M\mathbf{R}^{(n)}$, if we define $M$ as the $n_{0}\times n_{0}$ square matrix below

[TABLE]

and by induction we can easily get

[TABLE]

But we can do even more if we apply it backwards: it is obvious to see that in this basis we have $\mathbf{R}^{(1)}=(1,0,\ldots,0)^{T}$, and also that $M\mathbf{R}^{(1)}=(0,1,0,\ldots,0)^{T}=\mathbf{R}^{(2)}$, and so on. Thus one can rewrite

[TABLE]

Now we define $\mathbf{R}^{(\infty)}$ as the operator \big{(}\rho_{\infty}-\rho^{(0)}\big{)} spanned in our basis $\mathcal{B}_{n_{0}}$, and from this we can formally get

[TABLE]

where the last equality is true if the series in $\varepsilon$ converges. We have now a very compact and elegant form to express the NESS, and please note that $\mathbf{R}^{(\infty)}$ gives the components of NESS in the basis $\mathcal{B}_{n_{0}}$, i.e.

[TABLE]

Although on right hand side of (15) we only have $\mathbf{R}^{(1)}$, which is essentially the operator $\rho^{(1)}$ written in the basis $\mathcal{B}_{n_{0}}$, we remind that one needs to construct the sequence $\{\rho^{(1)},\ldots,\rho^{(n_{0})}\}$ in order to find the coefficients $c_{1},\ldots,c_{n_{0}}$ defined in Eq. (7) to construct the matrix $M$, and then evaluate the resolvent $(\mathbbm{1}_{n_{0}}-i\varepsilon M)^{-1}$. But we recap that we have started with a formal series (5) which has been rewritten in a compact form in the last equality of (15). Therefore the original series converges if the series for the resolvent converges, and this will happen if $\varepsilon<\lambda$, with

[TABLE]

where $\lambda_{j}$ are the eigenvalues of the matrix $M$. Now we have an approach to evaluate the radius of convergence $\lambda$ for the series (5).

III Weak system-bath coupling regime

The theoretical procedure to find the exact analytic expression (15), described in the previous section, is quite simple, or better saying, straightforward for small lengths of the chain $N$. We used Mathematica to perform explicit computations for $N=2,3,4,5$ which give us information about the radius of convergence $\lambda$ defined by Eq. (16). We emphasize here that all calculations are done in exact arithmetic after setting some values for parameters $\mu$ and $\Delta$, but are becoming practically unfeasible for $N\geq 6$. Anyway, from the data obtained we could conjecture the behavior for $\lambda$ dependence on $N$ and $\Delta$.

First of all, the software must compute the critical index $n_{0}$ defined in equation (7). We list in Table 1 some values for $n_{0}$ in three different perturbation regimes, and the case (A) corresponds to our first expansion in $\varepsilon$ parameter with anisotropy $\Delta>0$.

We have also checked that the critical index $n_{0}$ does not depend on parameter $\Delta$, except for a sharp change at $\Delta=0$ – see case (B) on table 1. We remind that $\dim\mathcal{H}=2^{2N}$, so this first step shows us that we only need a small number of vectors on $\mathcal{H}$ to express the NESS solution $\rho_{\infty}$, and this allowed us to optimize computational resources. Even so, as previously said, it was not possible to find $n_{0}$ for a spin chain with size $N\geq 6$, but anyway we can see some pattern showing up when $N$ varies from 2 to 5 and conjecture how the radius of convergence $\lambda$ behaves as $N$ grows.

Once the indices $n_{0}$ are now known, we can – by means of Mathematica – find the sequence $\left\{\rho^{(1)},\ldots,\rho^{(n_{0})}\right\}$ using recurrence relation (6), and then evaluate the coefficients $c_{j}$ as in equation (7). We fixed $\mu=1/2$ for all the evaluations, and ran the system size $N$ for $2\leq N\leq 5$, and anisotropy $\Delta>0$ from $10^{-3}$ to $10^{3}$. As described in section II, after finding the coefficients $c_{j}$, we construct the matrix $M$ and calculate its spectral radius $\lambda$. In figure 1, we have chosen different fixed values for the anisotropy $\Delta$, and then we plot a graph of $\log\lambda$ vs $N$, and the results suggest that $\lambda$ decays with the system size $N$, and apparently faster than exponentially. Moreover, for $\Delta>1$, the radius of convergence $\lambda$ decays faster when we increase the anisotropy $\Delta$, as we can see in Figure 1. When $0<\Delta<1$, $\lambda$ also decays with $N$, as exemplified by data for $\mu=1/2$ plotted in Figure 1. We note that for maximum driving $\mu=1$, the NESS is known analytically to be a polynomial in $\varepsilon$ of order $2N-2$ Prosen (2011b), so the $\varepsilon$-expansion has trivially infinite radius of convergence there for any $\Delta$.

To understand what happens to our approach when $\Delta$ gets closer to zero, we used the procedure described in section II to study the solutions when we have exactly $\Delta=0$ in Eq. (1), namely the $XX$ model. Repeating all the steps for this simpler model, we are now able to find the $n_{0}$ indices for system size $N\leq 6$, and they are listed on case (B) of table 1. Although the solution for $XX$ model is well known Prosen (2011b), we use our approach here to understand what happens if we keep using smaller positive values for the anisotropy in our perturbation, trying to get closer to the limit $\Delta\downarrow 0$. For $\Delta=0$, the values for $n_{0}$ change drastically, as we can see comparing cases (A) and (B) in table 1. Except for $N=2$, we can see that $n_{0}$ is (much) smaller when $\Delta=0$. By the way, these first values indicate that the critical index increases linearly with the system size, $n_{0}=2(N-1)$. Anyway, fixing some $N$, our results suggest that the square matrix $M$ is much larger for any $\Delta>0$, no matter how small the anisotropy parameter is, so $M$ has to dramatically shrink when we solve the $XX$ model using this approach. We could observe that this affects the behavior of the inverse of the largest eigenvalue when we increase $N$, this may explain the behavior of $\lambda$ as we approach $\Delta\downarrow 0$ for $N=5$. Moreover, we were able to find out that for $N$ from 2 up to 6, the matrix $M$ has only two different eigenvalues, $\pm 1/2$, each one with $(N-1)$ geometric multiplicity. Therefore, when we increase $N$, the radius of convergence is constant $\lambda_{0}=2$ for $XX$ model.

IV Linear Response Regime

In section II, we have described a procedure to study the radius of convergence $\lambda$ for a perturbation series (5) for bath-system coupling parameter $\varepsilon$. We can easily adopt it to study the solutions for the linear response regime, where we use another parameter as the perturbative one. Here, again we start from the solution in the equilibrium regime, when we have symmetric magnetic baths ($\mu=0$) coupled to the boundaries. So we would expect that, at least for small values of the driving force $|\mu|$, we could study the NESS solution perturbatively. In other words, we are changing the perturbation parameter from $\varepsilon$ to $\mu$.

We start rewriting the fixed point equation (3). From (4), after some easy manipulations, we have

[TABLE]

where in the last equality we defined two new operators, $\hat{\mathcal{D}}_{0}$ and $\hat{\mathcal{D}}_{\mu}$. One can easily check that $\hat{\mathcal{D}}_{0}^{\dagger}=\hat{\mathcal{D}}_{0}$. We get

[TABLE]

and we define the linear operator $T_{\mu}:=\left(\mathrm{ad}H\right)+i\hat{\mathcal{D}}_{0}$, just to make expressions simpler. In a completely analogous way, we can try a formal series on $\mu$ as solution for NESS

[TABLE]

where the terms of the sequence are now labeled with a $\mu$ index just to distinguish that we are now performing an expansion on this perturbative parameter. Once again, if we naïvely substitute it on fixed point equation (18), we find a similar recurrence relation

[TABLE]

and again we used the same initial condition $\rho^{(0)}=2^{-N}\mathbbm{1}_{2^{N}}$, which refers to the equilibrium solution for $\mu=0$. But at this point we make some comments to show that we have a more comfortable situation now than in the previous $\varepsilon$-perturbation: one can easily see that even if $T_{\mu}$ is not a self-adjoint operator, its real part, $\mathrm{ad}H$, as its imaginary part, $\hat{\mathcal{D}}_{0}$, are both self-adjoint operators. As consequence, we can prove that $\mathrm{Ker}\ (T_{\mu})$ is a one-dimensional subspace of $\mathcal{H}$ spanned by the identity vector $\mathbbm{1}_{2^{N}}=2^{N}\rho^{(0)}$. In other words, we no longer have to worry about degeneracy constructing the elements for the sequence one by one. We also may prove that we can find any element for the sequence Eva , so we are only concerned about convergence properties for the series (19). In this sense, the perturbation is even simpler for $\mu$ parameter. The rest of the argument follows exactly the same line as in section II.

From now on, the approach follows as we did on previous section: First, by means of Mathematica code, we use Eq. (7) to obtain the indices $n_{0}$ for this perturbation, which are listed in the case (C) of table 1. Again the software could not evaluate it for $N\geq 6$. Then we fixed $\varepsilon=1$ and changed the system size $N$ from 2 to 5, and anisotropy $\Delta$ from $2^{-14}$ to $2^{14}$. For each case, we ran the code to evaluate the radius of convergence $\lambda_{\mu}$. Just for the sake of better understanding, we remind the main steps to the reader: since the index $n_{0}$ is known for each $N$, we use recurrence relation (20) to construct the basis \Big{\{}\rho^{(1)}_{\mu},\ldots,\rho^{(n_{0})}_{\mu}\Big{\}}, reminding that the vector $\hat{\mathcal{D}}_{\mu}\rho^{(n_{0})}_{\mu}$ can be written as a linear combination of vectors $\hat{\mathcal{D}}_{\mu}\rho^{(0)}_{\mu},\ldots,\hat{\mathcal{D}}_{\mu}\rho^{(n_{0}-1)}_{\mu}$, in a completely similar way as we defined it on Eq. (7). Now we have found coefficients $(c_{\mu})_{j}$, so we can construct the matrix $M_{\mu}$, exactly as in (13), but with different size $n_{0}$ and with its last column given by $(c_{\mu})_{1},\ldots,(c_{\mu})_{n_{0}}$. Now we can obtain $\mathbf{R}^{(\infty)}_{\mu}$ as in Eq. (15), and to evaluate the radius of convergence for $\mu$-perturbation as

[TABLE]

where $(\lambda_{\mu})_{j}$ are the eigenvalues of the matrix $M_{\mu}$.

Here, in the linear response regime, our numerical results allow us to conjecture a very interesting behavior for $\lambda_{\mu}$. In Figure 2, for a fixed system size $N$ up to 5, we study how $\lambda_{\mu}$ depends on $\Delta$ in logarithmic scale. We can clearly see a transition between two behaviors as the black dashed guideline in Figure 2 indicates: for $0<\Delta<1$, we can see that $\lambda_{\mu}$ decreases as $\Delta$ increases. However, when we look to $\Delta>1$ we see that the radius of convergence still decreases as anisotropy increases, but very slightly. Moreover, from our results we can infer a lower bound $\lambda_{\min}=1$ for the radius of convergence in the linear response regime.

We can see that we always have $\lambda_{\mu}>1$, and it approaches the lower bound when $\Delta$ increases. These results may indicate that the radius of convergence for this perturbation series does not decay to zero with the system size, in contrast with the $\varepsilon$ perturbation. Actually, we conjecture that $\lambda_{\mu}\downarrow 1$ as $\Delta\to\infty$. If this statement is true, the perturbative solution in this linear response regime is reliable at least for $\lambda_{\min}=1$, independent of the anisotropy $\Delta$ and the system size $N$. So, in opposition to the $\varepsilon$-parameter expansion, where the radius of convergence decays at least exponentially with the system size $N$, here, in linear response regime our results indicates that in the $\mu$-parameter expansions the radius $\lambda_{\mu}$ remains larger than 1 as $N$ increases, so the expansion remains relevant, and maybe it is still suitable in thermodynamic limit. Moreover, the results indicate that a perturbative solution (19) can always reach the maximal driving solution $\mu=\pm 1$ regime, for example.

V Conclusions

We have elaborated on a formal and numerical analysis of radius of convergence for the perturbative solutions of NESS in boundary driven quantum spin chains. We have in particular expanded around the integrable points, where driven Lindblad equation allows for exact solutions. Even though we could only do exact numerical computations for relatively short chains (up to 5 sites), our results allow as to draw some general conclusions. For example, when expanding in the system-bath coupling strength parameter, the radius of convergence generally shrinks to zero very quickly by increasing the system size. On the other hand, when expanding in the driving (bias) parameter (the first order being just the linear response physics), then the radius of convergence appear to be uniformly lower bounded by $1$.

Acknowledgements.

This work was supported by grant number 249011/2013-1 of CNPq (Brazil), grants P1-0044, N1-0025 of Slovenian Research Agency, and ERC grant OMNES. HCFL thanks Alexandre C.L. Almeida for his help and comments.