# Dissipative Dynamics in Open Fermionic Chains

A. I. Karanikas1 and G. E. Pavlou2111Present address: Department of Physics,
University of Crete, 71003, Heraklion, Greece 1 National and Kapodistrian
University of Athens, Physics Department, Nuclear & Particle Physics Section
Panepistimiopolis, Ilissia 15771 Athens, Greece  2 Institute of Applied and
Computational Mathematics (IACM), Foundation for Research and Technology
Hellas (FORTH), IACM/FORTH, GR-71110 Heraklion, Greece  $^1$
<EMAIL_ADDRESS>$^2$<EMAIL_ADDRESS>
###### Abstract

By merging the Feynman-Vernon’s approach with the out-of-equilibrium Keldysh-
Schwinger formalism, we construct the reduced generating functional through
which all the time-dependent correlation functions of an open fermionic system
can be directly derived by applying the appropriate functional derivatives. As
a concrete example, we investigate the transverse Ising model, we derive the
covariance matrix at the steady state of the system and we investigate its
critical behavior.

##  1 Introduction

The unavoidable influence of the environment on open quantum systems is of
great importance for understanding their physical properties and for handling
their practical applications. When a system is embedded in a controlled or
uncontrolled environment, induced decoherence is one of the major issues for
storing and processing quantum information. Among these systems, for reasons
of both technical and theoretical origin, fermionic and spin chains with
short-range interactions, play a central role as their dynamics reveal a quite
interesting complexity even at their simple 1D versions. Isolated systems, the
Hamiltonian of which is quadratic in fermionic or bosonic degrees of freedom,
have been extensively studied as they are exactly solvable and the structure
of their ground state offers the explanation basis for highly non-trivial
phenomena as quantum phase transitions of topological character [1, 2].
However, the embedding of such a system into a bath induces major changes in
the behavior of the system’s correlations in a way that is not easy to
calculate even if the isolated system is completely solvable. When a system is
isolated, the structure of its ground state as well as the energy gap between
this state and the excited ones has a central role for the system’s properties
[3]. In the case of an open system, the role of the ground state is played by
the so-called steady state that is, the state at which the reduced density
matrix relaxes at the infinite time limit [4, 5]. The study of the steady
state’s quantum properties as well as the rate at which this it is approached,
is of great importance for both theoretical and practical reasons. A driven
approach to the steady state has been used to deal with quantum information
processing [6, 7, 8, 9]. The issue is standardly probed in the Lindbland
master equation framework in which the steady state can be defined as the
(right) eigenstate of the Lindblad super-operator with zero eigenvalues [10,
11].

In a recent paper [12], we presented a scheme that combines the Feynman-
Vernon’s influence functional technique with the out-of-equilibrium Keldysh
formalism. This construction enables the direct calculation of the
environment’s impact on the correlation functions of an open bosonic or
fermionic quantum system. In this approach, the focus is not on the reduced
density matrix but on the correlation functions per se. This permits the
calculation of functions that contain more than a single time variable [13].

In the present work, we apply this technique to probe the dynamics of
quadratic 1D fermionic chains that when isolated, can be mapped on a spin
chain model. The focus of the calculation is on the subsystem’s covariance
matrix [5], the quantity that encaptures the properties of reduced density
matrix. For the quadratic case, we find the steady state at which the system
relaxes as well as the rate at which this state is approached.

The structure of the paper is the follows: In Section 2 we briefly present the
basic ingredients of our formalism which is characterized by the the
introduction of the reduced generating functional. The latter is written as a
coherent state path integral over paths parametrized in terms of the Keldysh-
Schwinger complex time variable. The environment is simulated by a collection
of fermionic harmonic oscillators and, by assuming that it interacts linearly
with the system, its degrees of freedom are integrated out. The resulting
quantity, called influence functional, fully expresses the impact of the
environment on the dynamics of the system. In this way, we construct the
reduced generating functional through which the calculation of the covariance
matrix, as well the calculation of any reduced correlation function, is
immediate via functional differentiations. The results of this Section are
quite general not depending on the specific system under consideration.

In Section 3, we apply the aforementioned formalism for a system that is
described by a Hamiltonian which is quadratic when written in terms of
Majoranas variables and we derive the general form of the generating
functional through which the all-important covariance matrix can be
straightforwardly calculated. We also discuss the Markovian limit, on which
the analytical calculations of the next Section are based.

Section 4 refers to a fermionic system which, when isolated, is characterized
by a quantum phase transition. In the framework of the present formalism, we
examine its properties when it is part of a compound system. The site-
translational symmetry it possesses, facilitates the diagonalization of the
corresponding Hamiltonian providing, thus, a concrete example of the
calculations presented in Section 3. We derive the covariance matrix and we
examine its analytic properties at the steady state limit. The main result of
this Section is the confirmation that the non-analyticity presented in the
ground state correlations of the isolated system, remains in the steady state
of the open system indicating the possible persistence of quantum
correlations.

Finally, in Section 5 we present the conclusions of the current investigation.
In A, we present some technical details that were left out in the main text.

##  2 The Reduced Generating Functional

We consider a compound fermionic system consisting of two parts: The
(sub)system $S$ and its environment $E$. The dynamics of the isolated $S+E$
system is controlled by a Hamiltonian of the form:

$\hat{H}=\hat{H}_{S}\otimes\hat{I}_{E}+\hat{I}_{S}\otimes\hat{H}_{E}+\hat{H}_{I}$
(1)

The Hamiltonian operators in the last expression are defined in terms of
fermionic creation and annihilation operators:

$\hat{H}_{S}=\hat{H}\left(\hat{a}_{S}^{\dagger},\hat{a}_{S}\right),~{}~{}~{}\hat{H}_{E}=\hat{H}\left(\hat{a}_{E}^{\dagger},\hat{a}_{E}\right),~{}~{}~{}\hat{H}_{I}=\hat{H}\left(\hat{a}_{S}^{\dagger},\hat{a}_{S};\hat{a}_{E}^{\dagger},\hat{a}_{E}\right)$
(2)

The physical picture we adopt is the following: Up to a moment $t_{in.}$ the
environment and the system are independent of each other
$\left(\hat{H}_{I}\left(t\leq t_{in.}\right)=0\right)$ and the density
operator of the compound system can be expressed in the product form
$\hat{\rho}\left(t_{in.}\right)=\hat{\rho}_{E}\left(t_{in.}\right)\otimes\hat{\rho}_{S}\left(t_{in.}\right)$.
Before the moment $t_{in.}$ the environment is in equilibrium at temperature
$T=\beta^{-1}$ meaning that
$\hat{\rho}_{E}=\left(Z_{E}\left(\beta\right)\right)^{-1}\exp\left(-\beta\hat{H}_{E}\right)$.
After the initialization of the interaction, the parts of the compound system
entangle, and time evolution becomes non-unitary.

As probes for the system’s dynamics, one usually considers correlations of the
form:

$G_{jk}\left(\rho;t_{2},t_{1}\right)=Tr\left[\hat{\rho}\left(t_{in.}\right)\hat{T}\left(\hat{a}_{S,j}^{\dagger}\left(t_{2}\right)\hat{a}_{S,k}\left(t_{1}\right)\right)\right].$
(3)

In the last expression, $\hat{a}_{S,j}\left(t\right)$ is a fermionic
Heisenberg operator that refers to the system, the subscript $j$ is a site or
space index, $\hat{\rho}\left(t_{in.}\right)$ is the initial density operator
of the compound system, $\hat{T}$ is the time ordering operator and the trace
operation refers to both the environment’s and the system’s degrees of
freedom.

Needless to say, the calculation of 3 and of any reduced correlation function,
is not a trivial task. All the efforts, analytical or numerical, for
confronting the issue are restricted in the framework of Lindblad’s master
equation [10, 11]. In the present study we shall adopt a more general
formalism for the study of the system’s dynamics, based on the introduction of
a functional that generates, via the application of the appropriate
derivatives, the reduced functions we are interested in.

To set the stage, we begin by noting that correlators of Heisenberg field
operators are characterized by a forward-backward time structure (see Eq. 3).
As long as the system is in equilibrium and isolated this structure is easily
taken into account by assuming an adiabatic evolution of the non-interacting
vacuum. However, when the system is out of equilibrium or open, the adiabatic
hypothesis is not valid. In that case, the dynamics has been confronted, a
long ago, via the Keldysh-Schwinger formalism [14, *Keldysh1965]. In the
present case, we examine a quantum system in contact with its surroundings or,
put in other words, we focus on a (presumably small) part of an isolated
compound system. As it is obvious, the situation shares a lot with the out of
equilibrium dynamics. Taking this into account, we have presented in [12] a
formalism based on a reduced version of the Keldysh-Schwinger correlation
functions. It is based on an extension of Feynman-Vernon’s technique along the
Keldysh time contour, an extension that focuses on the calculation of
correlation functions bypassing the determination of the reduced density
matrix per se.

The scheme begins with the interpretation of the trace operation in Eq. 3 in
the (over complete) coherent state basis
$\left|\mathbf{z}_{S}\right>\equiv\otimes_{j\in S}\left|z_{j}\right>$ for the
system and
$\left|\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upzeta}_{E}\right>\equiv\otimes_{\mu\in
E}\left|\zeta_{\mu}\right>$ for the environment, followed by the construction
of the generating functional for system’s correlators.

Apart from the details about the peculiarities of defining path integration
over coherent states, the key ingredient in the aforementioned construction is
that the paths entering the generating functional are parametrized along the
so-called Keldysh time contour [14, *Keldysh1965].

The Keldysh time contour is defined on the complex plane along a closed
contour $P$ that encircles the real $t$ axis running from $t_{in.+}\equiv
t_{in.}+i0$ to $t_{in.-}\equiv t_{in.}-i0$. The contour consists of two
straight lines. The first one, denoted as $L_{+}$, joins the point $t_{in.+}$
to an arbitrary time instance $T_{+}=T+i0$. Along this line the time variable
is denoted as $t_{+}$. The second line joins $t_{-}=T-i0$ and along this line,
defined as $L_{-}$, time is denoted as $t_{-}$. In the case of thermal initial
states, the contour is extended by a complex time line, running parallel to
the imaginary axis, from $t_{in.-}$ to $t_{in.}-i\beta$, where $\beta^{-1}$ is
the temperature of the corresponding thermal state. This third line is denoted
as $L_{\beta}$, while the extended contour is denoted as $C$.

Figure 1: The Keldysh Contour

A natural ordering arises along this configuration since times along $L_{+}$
are considered to be before times along $L_{-}$, which are also considered to
be before times along $L_{\beta}$ (see Fig. 1). Taking as granted the
technical details presented in Refs. [12, 13]. we quote here the result for
the reduced generating functional:

$\eqalign{Z_{S}\left[\bar{\mathbf{{J}}},\mathbf{J}\right]&=\smallint
d^{2}\mathbf{z}_{S}d^{2}\mathbf{z}^{\prime}_{S}\left<-\mathbf{z}_{S}\right|\hat{\rho}_{S}\left(t_{in.+}\right)\left|\mathbf{z}^{\prime}_{S}\right>\cr&\times\int\limits_{\mathbf{z}\left(t_{in.+}\right)=\mathbf{z}_{S}}^{\bar{\mathbf{{z}}}\left(t_{in.-}\right)=\bar{\mathbf{{z}}}^{\prime}_{S}}{D^{2}\mathbf{z}\left(t_{P}\right)e^{iS_{S}\left[\bar{\mathbf{{z}}},\mathbf{z}\right]-E\left(\bar{\mathbf{{z}}}_{S},\mathbf{z}_{S}\right)+i\int\limits_{P}{dt_{P}\left(\mathbf{\bar{J}}\cdot\mathbf{\dot{z}}+\mathbf{\bar{z}}\cdot\mathbf{\dot{J}}\right)}}I_{E}\left[\bar{\mathbf{{z}}},\mathbf{z}\right]}.}$
(4)

At this point explanations are needed for the last expression to become
transparent: The Grassmann fields
$\mathbf{z}\left(t_{P}\right)\equiv\mathbf{z}_{P}\left(t\right)$ entering in
the path integral are defined along the Keldysh contour $P$ (assuming that the
system is at zero temperature):

$\mathbf{z}\left(t_{P}\right)=\left\\{\begin{array}[]{c}\mathbf{z}\left(t_{+}\right)\equiv\mathbf{z}_{+}~{}~{}\mathrm{
along}~{}~{}L_{+}\\\
\mathbf{z}\left(t_{-}\right)\equiv\mathbf{z}_{-}~{}~{}\mathrm{along}~{}~{}L_{-}\\\
\end{array}\right..$ (5)

These variables are considered as independent and they are integrated
separately. The action $S$ refers to the system and assumes the following
form:

$S_{S}=\int\limits_{P}{dt_{P}}\left[\frac{i}{2}\left(\mathbf{\bar{z}}\cdot\mathbf{\dot{z}}-\mathbf{\dot{\bar{z}}}\cdot\mathbf{z}\right)-H_{S}\left(\mathbf{\bar{z}},\mathbf{z}\right)\right].$
(6)

In Eq. 4 the action is accompanied by the surface term

$E\left(\mathbf{\bar{z}}^{\prime}_{S},\mathbf{z}_{S}\right)=\frac{1}{2}\left[\left|\mathbf{\bar{z}}^{\prime}_{S}\right|^{2}+\left|\mathbf{z}_{S}\right|^{2}\right]-\frac{1}{2}\left[\mathbf{\bar{z}}^{\prime}_{S}\cdot\mathbf{z}\left(t_{in.-}\right)+\mathbf{\bar{z}}\left(t_{in.+}\right)\cdot\mathbf{z}_{S}\right].$
(7)

The generating functional depends on the initial state of the system through
the factor
$\left<-\mathbf{z}_{S}\right|\hat{\rho}_{S}\left(t_{in.+}\right)\left|\mathbf{z}^{\prime}_{S}\right>$.
When the system is initially in its ground state,
$\hat{\rho}_{S}\left(t_{in.+}\right)=|GS\rangle\langle GS|$, this factor
together with the surface term are integrated out and the generating
functional reduces to the simpler form:

$Z_{S}^{\left(0\right)}\left[\mathbf{\bar{J}},\mathbf{J}\right]=\int{D^{2}\mathbf{z}\left(t_{P}\right)e^{iS_{S}\left[\mathbf{\bar{z}},\mathbf{z}\right]+i\int\limits_{P}{dt_{P}\left(\mathbf{\bar{J}}\cdot\mathbf{z}+\mathbf{\bar{z}}\cdot\mathbf{J}\right)}}I_{E}\left[\mathbf{\bar{z}},\mathbf{z}\right]}.$
(8)

Sources have been added in Eq. 4 to produce the system’s correlation functions
through functional differentiation. For example, the correlation function 3
can be derived by using the formula:

$G_{jk}\left(\rho;t_{2},t_{1}\right)=\left.\frac{\delta^{2}\ln
Z_{S}^{\left(0\right)}}{\delta
J_{j}\left(t_{1+}\right)\delta\bar{J}_{k}\left(t_{2+}\right)}\right|_{J=0}.$
(9)

The functional $I_{E}$, appearing in Eq. 4 is produced after tracing out the
environmental degrees of freedom. It is the reason we consider $Z_{S}$ as
reduced. It incorporates the influence of the environment on the system’s
dynamics and has the same significance as the Feynman-Vernon’s influence
functional for the calculation of the reduced density matrix. It is defined as

$I_{E}\left[\mathbf{\bar{z}},\mathbf{z}\right]=\frac{1}{Z_{E}\left(\beta\right)}\int\limits_{AP}{D\bm{\zeta}\left(t_{C}\right)e^{i\int\limits_{C}{dt_{C}\left(\mathbf{\bar{\bm{\zeta}}}\cdot\mathbf{\dot{\bm{\zeta}}}-H_{E}\left(\mathbf{\bar{\bm{\zeta}}},\bm{\zeta}\right)\right)-i\int\limits_{C}{dt_{C}H_{I}\left(\mathbf{\bar{\bm{\zeta}}},\bm{\zeta};\mathbf{\bar{z}},\mathbf{z}\right)}}}}.$
(10)

Here, the paths are parametrized along the full Keldysh contour $C=P\cup
L_{\beta}$ with the agreement that interactions are absent along the thermal
branch $L_{\beta}$. The antiperiodic boundary conditions are induced by the
fermionic trace operation. As expected, the calculation of the influence
functional is not a trivial task. In the present work, we adopt the
simplifying assumption that the environment can be simulated by a collection
of (fermionic) harmonic oscillators which interact linearly with the system:

$\eqalign{\hat{H}_{E}=\sum_{\mu}{E_{\mu}}\left(\hat{a}_{E,\mu}^{\dagger}\hat{a}_{E,\mu}-\frac{1}{2}\right),\cr\hat{H}_{I}=\sum_{j\in
S,\mu\in
E}{\left(\lambda_{j\mu}\hat{a}_{E,\mu}^{\dagger}\hat{a}_{S,j}+\bar{\lambda}_{j\mu}\hat{a}_{S,j}^{\dagger}\hat{a}_{E,\mu}\right)}}.$
(11)

This interaction results in fermion number non-conservation and is responsible
for the dissipation effects we are interested in. In this case, the influence
functional can be exactly calculated [13]:

$I_{E}=C_{E}^{-1}\exp\left[-\int\limits_{P}{dt_{P}\int\limits_{P^{\prime}}{dt^{\prime}_{P^{\prime}}\bar{\mathbf{z}}\left(t_{P}\right)\mathbf{\Delta}\left(t_{P}-t^{\prime}_{P^{\prime}}\right)\mathbf{z}\left(t^{\prime}_{P^{\prime}}\right)}}\right].$
(12)

The different indexing in time integration indicates the fact that the time
parameters in the last expression can run along different branches of the
Keldysh contour. The constant appearing in Eq. 12, $C_{E}^{-1}=\prod_{\mu\in
E}{2\cosh\left(\beta E{{{}_{\mu}}/{2}}\right)}$ , is of no significance for
the calculation of correlation functions. The crucial quantities are the
matrix elements $\left(\mathbf{\Delta}\right)_{jk}=\Delta_{jk}$ which assume
the form [13]:

$\eqalign{\Delta_{jk}\left(t_{P}-t^{\prime}_{P^{\prime}}\right)&=\sum_{\mu\in
E}{\bar{\lambda}_{j\mu}\lambda_{k\mu}}\left[\Theta\left(t_{P}-t^{\prime}_{P^{\prime}}\right)-\frac{1}{1+e^{\beta
E_{\mu}}}\right]e^{-i\left(t_{P}-t^{\prime}_{P^{\prime}}\right)E_{\mu}}=\cr&=\int\limits_{0}^{\infty}{dED_{jk}\left(E\right)}\left[\Theta\left(t_{P}-t^{\prime}_{P^{\prime}}\right)-\frac{1}{1+e^{\beta
E}}\right]e^{-i\left(t_{P}-t^{\prime}_{P^{\prime}}\right)E}.}$ (13)

Here we defined:

$D_{jk}\left(E\right)=\sum_{\mu\in
E}{\bar{\lambda}_{j\mu}}\lambda_{k\mu}\delta\left(E-E_{\mu}\right).$ (14)

In the following we shall assume that we can write:

$\sum_{\mu\in E}{\bar{\lambda}_{j\mu}\lambda_{k\mu}}=g_{jk}\sum_{\mu\in
E}{\left|\gamma_{\mu}\right|^{2}},~{}~{}~{}g_{jk}=\bar{g}_{kj}.$ (15)

As it is obvious the coefficients $g_{jk}$ characterize the strength of the
interaction between the sites $j,k$ of the system, induced by the environment.
The two extreme choices $g_{jk}=g\delta_{jk}$ and $g_{jk}=g~{}~{}\forall j,k$
correspond, the first one, to the case where each site of the system interacts
independently with the environment, and the second one to the case in which
all the sites of the system interact with each other with the same strength,
irrespective of their distance.

Thus, the influence of the environment, as it is encoded in Eq. 12, yields a
contribution to the generating functional which, in general, is nonlocal both
in space and time indices. Inserting Eq. 12 into Eq. 8 we find the expression
that we indent to use for the calculation in the current investigation:

$Z_{S}^{\left(0\right)}\left[\mathbf{\bar{J}},\mathbf{J}\right]=C_{E}^{-1}\int{D^{2}\mathbf{z}\left(t_{P}\right)e^{iS_{S}\left[\mathbf{\bar{z}},\mathbf{z}\right]+iS_{INF.}\left[\mathbf{\bar{z}},\mathbf{z}\right]+i\int\limits_{P}{dt_{P}\left(\mathbf{\bar{J}}\cdot\mathbf{z}+\mathbf{\bar{z}}\cdot\mathbf{J}\right)}}}$
(16)

with

$S_{INF.}=i\int\limits_{P}{dt_{P}\int\limits_{P^{\mathrm{{}^{\prime}}}}{dt^{\prime}_{P^{\prime}}}}\mathbf{\bar{z}}\left(t_{P}\right)\mathbf{\Delta}\left(t_{P}-t^{\prime}_{P^{\prime}}\right)\mathbf{z}\left(t^{\prime}_{P^{\prime}}\right).$
(17)

The generating functional in Eq. 16 can, in principle at least, be exactly
calculated if the system’s Hamiltonian is quadratic. In such a case the
covariance matrix can also be exactly calculated. We shall confront the issue
in the next chapter. At this point, is worth note that through the (reduced)
generating functional correlations of the form
$\left<\hat{O}_{1,S}\left(t_{1}\right)\hat{O}_{2,S}\left(t_{2}\right)...\hat{O}_{n,S}\left(t_{n}\right)\right>_{\rho}$
can be calculated. This must be contrasted to the case of calculations based
on Lindblad’s equation where one considers correlations that contain a single
time variable [16].

##  3 Quadratic Hamiltonians

In this chapter we are interested in systems the dynamics of which are
described by Hamiltonians that, when written in terms of Majorana variables,
assume the quadratic form:

$H_{S}=\frac{i}{4}\sum_{j,k=1}^{N}{\sum_{u,v=0,1}{\gamma_{j,u}}}A_{ju,kv}\gamma_{k,v}=\frac{i}{4}\sum_{\alpha,\beta=1}^{2N}{\gamma_{\alpha}}A_{\alpha\beta}\gamma_{\beta}\equiv\frac{i}{4}\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}\mathbf{A}\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}.$
(18)

The Majorana variables entering the last expression are defined according to
the rule:

$\gamma_{2j-1}=\bar{z}_{S,j}+z_{S,j}\equiv\gamma_{j,0},~{}~{}~{}\gamma_{2j}=i\left(\bar{z}_{S,j}-z_{S,j}\right)\equiv\gamma_{j,1}.$
(19)

In Eq. 18 we adopted the compact notation of Ref. [5]
$\alpha=\left(j,u\right)$ in which $j=1,...,N$ and $u=0,1$. In accordance, in
the present Section, vectors and matrices are defined in a $2N$ dimensional
space:

$\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}\equiv\left(\gamma_{1},\gamma_{2},...,\gamma_{2N-1},\gamma_{2N}\right)=\left(\gamma_{1,0},\gamma_{1,1},...,\gamma_{N,0},\gamma_{N,1}\right).$
(20)

The $2N\times 2N$ antisymmetric matrix $\mathbf{A}\doteq
A_{\alpha\beta}=A_{jk,uv}=-A_{kj,vu}=-A_{\beta\alpha}$ describes the
interaction between the sites of the system. In the quadratic case, the
generating functional 16 can be exactly calculated by minimizing the action:

$S=S_{S}+S_{INF.}+S_{J}.$ (21)

Due to form 18 of the Hamiltonian, it is convenient to express $S$ in terms of
Majorana variables:

$\eqalign{S=&\frac{i}{4}\int\limits_{P}{dt_{P}}\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}\left(t_{P}\right)\left(\mathbf{I}\partial_{t_{P}}-\mathbf{A}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}\left(t_{P}\right)\cr&+\frac{i}{4}\int\limits_{P}{dt_{P}}\int\limits_{P^{\prime}}{dt^{\prime}_{P^{\prime}}}\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}\left(t_{P}\right)\mathbf{D}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}\left(t^{\prime}_{P^{\prime}}\right)\cr&+\frac{1}{2}\int\limits_{P}{dt_{P}}\mathbf{f}\left(t_{P}\right)\cdot\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}\left(t_{P}\right).}$
(22)

In the last expression we introduced the sources:

$f_{j,0}=\bar{J}_{j}-J_{j},~{}~{}~{}f_{j,1}=i\left(\bar{J}_{j}+J_{j}\right)$
(23)

and we defined

$\mathbf{D}\doteq
D_{\alpha\beta}=D_{ju,kv}=\Delta_{jk}\left(1-\sigma^{y}\right)_{uv}.$ (24)

By minimizing Eq. 22 we get the following “classical” equation:

$\left(\mathbf{I}\partial_{t_{P}}-\mathbf{A}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}^{cl.}\left(t_{P}\right)+\int\limits_{P^{\prime}}{dt^{\prime}_{P^{\prime}}}\tilde{D}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}^{cl.}\left(t^{\prime}_{P^{\prime}}\right)=-i\mathbf{f}\left(t_{P}\right)$
(25)

in which

$\tilde{D}_{\alpha\beta}\left(t_{P},t^{\prime}_{P^{\prime}}\right)=\frac{1}{2}\left(D_{\alpha\beta}\left(t_{P},t^{\prime}_{P^{\prime}}\right)-D_{\beta\alpha}\left(t^{\prime}_{P^{\prime}},t_{P}\right)\right)=-\tilde{D}_{\beta\alpha}\left(t_{P},t^{\prime}_{P^{\prime}}\right).$
(26)

By introducing Green’s function $\Gamma_{\alpha\beta}$ as the causal solution
of Green’s equation

$\eqalign{&\int\limits_{P^{\prime\prime}}{dt^{\prime\prime}_{P^{\prime\prime}}}\left[\left(\mathbf{I}\partial_{t_{P}}-\mathbf{A}\right)\delta\left(t_{P}-t^{\prime\prime}_{P^{\prime\prime}}\right)+\mathbf{\tilde{D}}\left(t_{P},t^{\prime\prime}_{P^{\prime\prime}}\right)\right]\mathbf{\Gamma}\left(t^{\prime\prime}_{P^{\prime\prime}},t^{\prime}_{P^{\prime}}\right)=\cr&=-\mathbf{I}i\delta\left(t_{P}-t^{\prime}_{P^{\prime}}\right).}$
(27)

the generating functional 16 assumes the form:

$Z_{S}^{\left(0\right)}\left[\mathbf{f}\right]=C_{E}^{-1}\exp\left[\frac{i}{4}\int\limits_{P^{\prime}}{dt^{\prime}_{P^{\prime}}\int\limits_{P}{dt_{P}\mathbf{f}\left(t_{P}\right)\mathbf{\Gamma}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\mathbf{f}\left(t^{\prime}_{P^{\prime}}\right)}}\right].$
(28)

As it is obvious, having in hand the solution of Eq. 27 and the generating
functional 28, one can calculate, via a unitary transformation of the sources,
all the correlation functions of the system.

In the present study, we focus on the real and antisymmetric covariance matrix
that is defined as [5]:

$C_{jk,uv}\left(t\right)=\mathrm{Tr}\left\\{\hat{\rho}\left(t\right)\frac{i}{2}\left[\hat{\gamma}_{j,u},\hat{\gamma}_{k,v}\right]\right\\}=\mathrm{Tr}_{\mathrm{S}}\left\\{\hat{\rho}_{S}^{R}\left(t\right)\frac{i}{2}\left[\hat{\gamma}_{j,u},\hat{\gamma}_{k,v}\right]\right\\}.$
(29)

The key ingredient in this expression is the reduced density matrix

$\hat{\rho}_{S}^{R}\left(t\right)=Tr_{E}\left\\{e^{-i\left(t-t_{in.}\right)\hat{H}}\hat{\rho}\left(t_{in.}\right)e^{i\left(t-t_{in.}\right)\hat{H}}\right\\}$
(30)

which, in the present case, is a Gaussian state and, as such, all its
properties are encoded into the covariance matrix 29.

Knowledge of the covariance matrix permits physically important issues to be
addressed. One of the most challenging refers to the rate at which
correlations diminish due to decoherence. Closely related to it is the
approach of the reduced density matrix to the, so-called, steady state,
$\hat{\rho}_{S}^{R}\left(t\right)\stackrel{{\scriptstyle[}}{{t\rightarrow\infty}}]{}{\rightarrow}\hat{\rho}_{S,0}^{R}$
and the consequent approach of the covariance matrix to the corresponding
steady state covariance:
$C_{jk,uv}\left(t\right)\stackrel{{\scriptstyle[}}{{t\rightarrow\infty}}]{}{\rightarrow}C_{jk,uv}^{\left(0\right)}$.
At the thermodynamic limit $N\rightarrow\infty$, the knowledge of the analytic
structure of these functions is of great importance as it can reveal possible
quantum phase transitions.

The covariance matrix can be straightforwardly calculated through the
generating functional 28 by applying functional derivatives as follows:

$C_{jk,uv}\left(t\right)=-4i\left.\frac{\delta^{2}\ln
Z_{S}^{\left(0\right)}}{\delta f_{ju}\left(t_{+}\right)\delta
f_{kv}\left(t_{+}\right)}\right|_{J=0}.$ (31)

In the present approach, the key function for achieving the final result is
the solution of Green’s equation 27. However, this equation, in its general
form, is nonlocal both in space and time indices, a fact that makes its
analytical solution a very difficult task. Even so, it can be written in a way
more convenient for calculations as the real antisymmetric matrix
$A_{\alpha\beta}$ can always be diagonalized through an orthogonal
transformation:

$\mathbf{Q}^{\bot}\mathbf{AQ}=\oplus_{j=1}^{N}\left(\begin{array}[]{cc}0&\varepsilon_{j}\\\
-\varepsilon_{j}&0\\\ \end{array}\right),~{}~{}\varepsilon_{j}\geq
0,~{}~{}\mathbf{Q}\in O\left(2N\right).$ (32)

Thus, by using the unitary transformation

$\mathbf{V}=\mathbf{Q}\oplus^{N}\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&1\\\
i&-i\\\ \end{array}\right)$ (33)

we can write

$\mathbf{V}^{\dagger}\mathbf{AV}=\oplus^{N}\left(\begin{array}[]{cc}i\varepsilon_{j}&0\\\
0&-i\varepsilon_{j}\\\ \end{array}\right)\equiv i\mathbf{E}.$ (34)

By writing

$\eqalign{\sum_{j,u}{\gamma_{j,u}}V_{jk,uv}=\psi_{k,v}\rightarrow\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}\cdot\mathbf{V}=\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi},~{}~{}\mathbf{K}=\mathbf{V}^{\dagger}\mathbf{\tilde{D}V},\cr\mathbf{j}=\mathbf{f}\cdot\mathbf{V}\rightarrow
j_{k,v}=\sum_{j,u}{f_{j,u}}V_{jk,uv}}$ (35)

the action 22 can be written as follows:

$\eqalign{S&=\frac{i}{4}\int\limits_{P}{dt_{P}}\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}\left(t_{P}\right)\left(\mathbf{I}\partial_{t_{P}}-i\mathbf{E}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}\left(t_{P}\right)\cr&+\frac{i}{4}\int\limits_{P}{dt_{P}}\int\limits_{P^{\prime}}{dt^{\prime}_{P^{\prime}}}\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}\left(t_{P}\right)\mathbf{K}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}\left(t^{\prime}_{P^{\prime}}\right)\cr&-\frac{1}{4}\int\limits_{P}{dt_{P}}\mathbf{\bar{j}}\left(t_{P}\right)\cdot\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}\left(t_{P}\right)-\frac{1}{4}\int\limits_{P}{dt_{P}}\bar{\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}}\left(t_{P}\right)\cdot\mathbf{j}\left(t_{P}\right)}.$
(36)

In this case the classical equation of motion assumes the form:

$\left(\mathbf{I}\partial_{t_{P}}-i\mathbf{E}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}^{\mathbf{cl}.}\left(t_{P}\right)+\int\limits_{P^{\prime}}{{dt^{\prime}}_{P^{\prime}}}\mathbf{K}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}^{\mathbf{cl}.}\left(t^{\prime}_{P^{\prime}}\right)=-i\mathbf{j}\left(t_{P}\right).$
(37)

Thus, by defining the function
$\mathbf{L}=\mathbf{V}^{\mathbf{\dagger}}\mathbf{\Gamma V}$, the Green’s
equation 27 is recasted as:

$\left(\mathbf{I}\partial_{t_{P}}-i\mathbf{E}\right)\mathbf{L}\left(t_{P},t^{\prime}_{P^{\prime}}\right)+\int\limits_{P^{\prime\prime}}{{dt^{\prime\prime}}_{P^{\prime\prime}}}\mathbf{K}\left(t_{P},{t^{\prime\prime}}_{P^{\prime\prime}}\right)\mathbf{L}\left({t^{\prime\prime}}_{P^{\prime\prime}},t^{\prime}_{P^{\prime}}\right)=-\mathbf{I}i\delta\left(t_{P}-t^{\prime}_{P^{\prime}}\right)$
(38)

The last expression is a set of equations coupled due to the presence of the
kernel $\mathbf{K}$, the dissipation kernel. In the current paper we adopt an
approximation that can lead to the exact solution of the Green’s Eqs. 38. It
is the same approximation on which the Lindbland approach is founded namely
the, so-called, Markovian approximation [16]. It is based on the hypothesis
that the environment acts as a memoryless bath and, consequently, that the
influence functional can be considered as local in time. The validity of such
an approximation is based on the existence of two characteristic time scales
[17]. The first one, $\tau_{E}$, refers to the decay of environmental
correlations, while the second one, $\tau_{S}$, characterizes the frequencies
$\varepsilon_{S}~{}\tau_{S}^{-1}$ after which system’s dynamics are screened
out. If the last scale is much larger than the first one, the function 30
behaves much like as a delta function:

$\Delta_{jk}\left(t_{P}-t^{\prime}_{p^{\prime}}\right)\rightarrow
g_{jk}\Delta_{PP^{\prime}}\left(\varepsilon_{S}\right)\delta\left(t-t^{\prime}\right)$
(39)

The matrix elements $\Delta_{PP^{\prime}}$ are defined as [13]:

$\eqalign{\Delta_{++}=-i\delta E+\Gamma\left(\frac{1}{2}-b\right),\\\
\Delta_{--}=i\delta E+\Gamma\left(\frac{1}{2}-b\right),\mathrm{ }\\\
\Delta_{+-}=-\Gamma b,\\\ \Delta_{-+}=\Gamma\left(1-b\right)\\\ }$ (40)

with:

$\Gamma=2\pi D\left(\varepsilon_{S}\right),~{}~{}\delta
E=\mathrm{Pr}.\int\limits_{0}^{\infty}{dE}\frac{D\left(E\right)}{E-\varepsilon_{S}},~{}~{}b=\frac{1}{1+e^{\beta\varepsilon_{S}}}$
(41)

The Markovian limit makes the dissipation kernel diagonal in the time
variables. Even so, expression 38 is a set of $4\times\left(2N\right)^{2}$
coupled equations the solution of which is a rather complicated task. However,
when the Hamiltonian is translationally invariant Eqs. 38 can be exactly
solved. This is the task of the next Section.

##  4 Kitaev’ s chains

In the present Section we examine a simple -yet interesting- fermionic model
in one spatial dimension, namely the $1D$ superconducting chain [1, 18, 19]:

$\eqalign{\hat{H}_{S}=\cr\sum_{j=1}^{N}{\left[-w\left(\hat{a}_{j}^{\dagger}\hat{a}_{j+1}+\hat{a}_{j+1}^{\dagger}\hat{a}_{j}\right)+\Delta\left(\hat{a}_{j}\hat{a}_{j+1}+\hat{a}_{j+1}^{\dagger}\hat{a}_{j}^{\dagger}\right)-\mu\left(\hat{a}_{j}^{\dagger}\hat{a}_{j}-\frac{1}{2}\right)\right]}.}$
(42)

Although simple and completely solvable when isolated, this model is quite
rich due to the underlying quantum phase transition. After the Jordan-Wigner
transformation

$\hat{a}_{j}=\prod_{k=1}^{j-1}{\left(-\sigma_{k}^{z}\right)}\sigma_{j}^{-},\hat{a}_{j}^{\dagger}=\prod_{k=1}^{j-1}{\left(-\sigma_{k}^{z}\right)}\sigma_{j}^{+},\sigma_{j}^{\pm}=\frac{1}{2}\left(\sigma_{j}^{x}\pm
i\sigma_{j}^{y}\right)$ (43)

the fermionic system is mapped onto the anisotropic $XY$ spin system:

$\hat{H}_{S}=\sum_{j}{\left(\frac{\Delta-w}{2}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+\frac{\Delta+w}{2}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{y}-\frac{\mu}{2}\hat{\sigma}_{j}^{z}\right)}.$
(44)

For reasons of simplicity, we choose $w=-\Delta=1$ and we write $\mu=2h$. In
this case, the spin system 44 reduces to the transverse Ising model:

$\hat{H}_{S}=-\sum_{j}{\left(\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+h\hat{\sigma}_{j}^{z}\right)}$
(45)

while its fermionic ancestor reads:

$\hat{H}_{S}=-\sum_{j=1}^{N}{\left[\left(\hat{a}_{j}^{\dagger}-\hat{a}_{j}\right)\left(\hat{a}_{j+1}^{\dagger}+\hat{a}_{j+1}\right)+2h\left(\hat{a}_{j}^{\dagger}\hat{a}_{j}-\frac{1}{2}\right)\right]}.$
(46)

In the presence of the environment, the dynamics of these systems become quite
complicated studied mainly in the framework of the Lindbland master equation.
From this point of view the model 42, and all its integrable variants,
provides the ideal stage for presenting the analytic formalism introduced in
Sections 2 and 3. However, a note is needed at this point: When the system we
are interested in is isolated, the Hamiltonians 42 and 44 are mathematically
equivalent meaning that the calculation of the relevant correlation functions
can be based either on the first or the second one. When the system is part of
a larger compound system, it is the mutual interaction that defines the
physical degrees of freedom. In the case under study the environment is a
fermionic bath meaning that the physical degrees of freedom for the
description of the system are fermionic and the relevant quantum Hamiltonian
is 44.

The classical Hamiltonian representing the fermionic system 46 in the
integral7, can be written in the form:

$\eqalign{H_{S}&=\sum_{j}{\left[\left(z_{j}z_{j+1}-\bar{z}_{j}\bar{z}_{j+1}\right)+\left(z_{j}\bar{z}_{j+1}-\bar{z}_{j}z_{j+1}\right)-2h\bar{z}_{j}z_{j}\right]}=\cr&=\frac{i}{4}\sum_{jk,uv}{\gamma_{j,u}A_{jk,uv}\gamma_{k,v}}}$
(47)

The matrix $A$ entering the Majorana expression of the Hamiltonian is defined
as:

$\eqalign{A_{j,j+1;0,0}=-A_{j+1,j;0,0}=2,A_{j,j+1;0,1}=-A_{j+1,j;1,0}=-2,\cr
A_{j,j;0,1}=-A_{j,j;1,0}=2h.}$ (48)

The first step for obtaining the solution of the classical equation 27 is the
diagonalization of system’s Hamiltonian. For the isolated case, this task is
accomplished via the implementation of the discrete Fourier transform [14,
*Keldysh1965, 19]:

$z_{j}=\frac{1}{\sqrt{N}}\sum_{m=0}^{N-1}{e^{i\varphi_{m}j}}c_{m}.$ (49)

The exact form of the phase $\varphi_{m}=\frac{2\pi}{N}\left(m+\kappa\right)$
depends on the number of fermions in the system. When this number is even,
$\kappa=1/2$, while $\kappa=0$ when odd. The first case is connected with the
antiperiodic boundary condition $z_{j+N}=-z_{j}$ and the latter with the
periodic one $z_{j+N}=z_{j}$ [19]. When dissipation is present, the fermion
number is not conserved meaning that, strictly speaking, the transformation 49
should not be used. However, at the thermodynamic limit $N\rightarrow\infty$
which we are interested in, the summation in Eq. 47 is extended over
$\mathbb{Z}$ and the discrete transform 49 can be replaced by its continuum
version:

$z_{j}=\frac{1}{\sqrt{2\pi}}\int\limits_{-\pi}^{\pi}{d\varphi}e^{i\varphi
j}c_{\varphi},~{}~{}~{}c_{\varphi}=\frac{1}{\sqrt{2\pi}}\sum_{j=-\infty}^{\infty}{e^{-i\varphi
j}}z_{j}.$ (50)

These transformations are going to be helpful only if the system-environment
interaction does not destroy the translational invariance presented in the
isolated system. To preserve this symmetry, we shall assume that
$g_{jk}=g_{j-k}=g_{k-j}$.

For the Majorana variables, the Fourier transforms read as follows:

$\gamma_{j,u}=\frac{1}{\sqrt{2\pi}}\int\limits_{-\pi}^{\pi}{d\varphi}e^{i\varphi
j}\gamma_{\varphi,u},~{}~{}~{}\gamma_{\varphi,0}=\bar{c}_{-\varphi}+c_{\varphi},~{}~{}~{}\gamma_{\varphi,1}=i\left(\bar{c}_{-\varphi}-c_{\varphi}\right).$
(51)

Inserting expressions 51 into Eq. 47, the Hamiltonian can be rewritten in
terms of the conjugate fields in the following form:

$H_{S}=\int\limits_{0}^{\pi}{d\varphi}\bar{\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}}_{\varphi}\mathbf{A}_{\varphi}\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}_{\varphi}.$
(52)

To arrive at the last expression, we used the translational invariance of the
Hamiltonian to write $A_{jk,uv}=A_{j-k,uv}$, we defined the vectors
$\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}_{\varphi}=\left(\gamma_{\varphi,0},\gamma_{\varphi,1}\right)$
and the matrices:

$\left(\mathbf{A}_{\varphi}\right)_{uv}\equiv
A_{\varphi,uv}=\sum_{j}{e^{-i\varphi
j}}A_{j,uv}=-A_{-\varphi,vu},A_{j,uv}=\frac{1}{2\pi}\int\limits_{-\pi}^{\pi}{d\varphi
e^{i\varphi j}}A_{\varphi,uv}.$ (53)

It’s not difficult to find that:

$\mathbf{A}_{\varphi}=-\frac{i}{2}\left(\begin{array}[]{cc}0&h+e^{-i\varphi}\\\
-\left(h+e^{i\varphi}\right)&0\\\ \end{array}\right).$ (54)

By using the unitary transformation

$\mathbf{V}_{\varphi}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&1\\\
-i&i\\\ \end{array}\right)e^{-i\theta_{\varphi}\sigma^{x}}$ (55)

with

$\cos 2\theta_{\varphi}=\frac{h+\cos\varphi}{\varepsilon_{\varphi}},~{}~{}\sin
2\theta_{\varphi}=\frac{\sin\varphi}{\varepsilon_{\varphi}},~{}~{}\varepsilon_{\varphi}=\sqrt{\left(h+\cos\varphi\right)^{2}+\sin^{2}\varphi}$
(56)

the matrix 54 can be put into diagonal form:

$\mathbf{A}_{\varphi}=-\frac{1}{2}\mathbf{V}_{\varphi}\left(\begin{array}[]{cc}\varepsilon_{\varphi}&0\\\
0&-\varepsilon_{\varphi}\\\ \end{array}\right)\mathbf{V}_{\varphi}^{\dagger}.$
(57)

By defining

$\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}_{\varphi}={\mathbf{V}^{\mathbf{\dagger}}}_{\varphi}\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}_{\varphi},~{}~{}~{}\bar{\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}}_{\varphi}=\mathbf{\bar{\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\upgamma}}}_{\varphi}\mathbf{V}_{\varphi}$
(58)

the Hamiltonian 52 simplifies as follows:

$H_{S}=-\frac{1}{2}\int\limits_{0}^{\pi}{d\varphi\varepsilon_{\varphi}}{\bar{\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}}}_{\varphi}\sigma^{z}\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}_{\varphi}.$
(59)

Introducing the function

$\tilde{g}_{\varphi}=\frac{1}{\sqrt{2\pi}}\sum_{j}{e^{i\varphi
j}}g_{j}=\tilde{g}_{-\varphi}$ (60)

the influence of the environment is written as

$S_{INF.}=\frac{i}{2}\int\limits_{0}^{\pi}{d\varphi}\int\limits_{P}{dt_{P}}\int\limits_{P^{\mathrm{{}^{\prime}}}}{dt^{\prime}_{P^{\prime}}}{\bar{\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}}}_{\varphi}\left(t_{P}\right)\mathbf{K}_{\varphi}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}_{\varphi}\left(t^{\prime}_{P^{\prime}}\right)$
(61)

where

$\mathbf{K}_{\varphi}=\tilde{g}_{\varphi}\mathbf{V}_{\varphi}^{\dagger}\left(\begin{array}[]{cc}\Delta^{\left(-\right)}&i\Delta^{\left(+\right)}\\\
-i\Delta^{\left(+\right)}&\Delta^{\left(-\right)}\\\
\end{array}\right)\mathbf{V}_{\varphi},~{}~{}\Delta^{\left(\pm\right)}=\Delta\left(t_{P},t^{\prime}_{P^{\prime}}\right)\pm\Delta\left(t^{\prime}_{P^{\prime}},t_{P}\right).$
(62)

Thus, the action 21 written in terms of the fields 53 gets the form:

$\eqalign{S=&\frac{1}{2}\int\limits_{0}^{\pi}{d\varphi}\int\limits_{P}{dt_{P}}\bar{\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}}_{\varphi}\left(i\partial_{t_{P}}+\varepsilon_{\varphi}\sigma^{z}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}_{\varphi}\cr&+\frac{i}{2}\int\limits_{0}^{\pi}{d\varphi}\int\limits_{P}{dt_{P}}\int\limits_{P^{\prime}}{dt^{\prime}_{P^{\prime}}}\bar{\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}}_{\varphi}\left(t_{P}\right)\mathbf{K}_{\varphi}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}_{\varphi}\left(t^{\prime}_{P^{\prime}}\right)\cr&+\frac{1}{2}\int\limits_{0}^{\pi}{d\varphi}\int\limits_{P}{dt_{P}\left(\mathbf{\bar{j}}_{\varphi}\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}_{\varphi}+\bar{\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}}_{\varphi}\mathbf{{j}}_{\varphi}\right)}.}$
(63)

For the source term we followed the notation (see Eq. 23):

$\mathbf{\bar{j}}_{\varphi}=\mathbf{f}_{\varphi}\mathbf{V}_{\varphi},~{}~{}~{}\mathbf{f}_{\varphi}=\left(f_{\varphi,0},f_{\varphi,1}\right),~{}~{}~{}f_{\varphi,u}=\frac{1}{\sqrt{2\pi}}\sum_{j}{e^{i\varphi
j}f_{j,u}}$ (64)

By minimizing Eq. 62 we get the classical equation:

$\left(i\partial_{t_{P}}+\varepsilon_{\varphi}\sigma^{z}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}_{\varphi}^{cl.}\left(t_{P}\right)+i\int\limits_{P^{\prime}}{dt^{\prime}_{P^{\prime}}}\mathbf{K}_{\varphi}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\uppsi}_{\varphi}^{cl.}\left(t^{\prime}_{P^{\prime}}\right)=-\mathbf{j}_{\varphi}\left(t_{P}\right).$
(65)

The corresponding Green’s equation reads:

$\left(i\partial_{t_{P}}+\varepsilon_{\varphi}\sigma^{z}\right)\mathbf{L}_{\varphi}\left(t_{P},t^{\prime}_{P^{\prime}}\right)+i\int\limits_{P^{\prime\prime}}{dt^{\prime\prime}_{P^{\prime\prime}}}\mathbf{K}_{\varphi}\left(t_{P},t^{\prime\prime}_{P^{\prime\prime}}\right)\mathbf{L}_{\varphi}\left(t^{\prime\prime}_{P^{\prime\prime}},t^{\prime}_{P^{\prime}}\right)=\mathbf{I}\delta\left(t_{P}-t^{\prime}_{P^{\prime}}\right).$
(66)

Having in hand the solution of the last equation and using Eq. 58 the
generating functional gets the form

$Z_{S}^{\left(0\right)}\left[\mathbf{j}\right]=C_{E}^{-1}\exp\left[-\frac{i}{2}\int\limits_{0}^{\pi}{d\varphi}\int\limits_{P^{\mathrm{{}^{\prime}}}}{dt^{\prime}_{P^{\prime}}\int\limits_{P}{dt_{P}\mathbf{\bar{j}}_{\varphi}\left(t_{P}\right)\mathbf{L}_{\varphi}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\mathbf{j}\left(t^{\prime}_{P^{\prime}}\right)}}\right].$
(67)

Taking into account 58 the generating functional can be rewritten in the form
28 with

$\Gamma_{jk,uv}\left(t_{P},t^{\prime}_{P^{\prime}}\right)=\frac{1}{4\pi}\int\limits_{0}^{\pi}{d\varphi
e^{i\varphi\left(j-k\right)}\left[\mathbf{V}_{\varphi}\mathbf{L}_{\varphi}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\mathbf{V}_{\varphi}^{\dagger}\right]_{uv}}.$
(68)

The Green’s equation 66 with causal boundary conditions can be exactly solved.
In A we present the details of the calculation. Here it is enough to quote the
result for $t_{+}=t^{\prime}_{+}=t$:

$\eqalign{\mathbf{L}_{\varphi}=-i\left[\frac{1}{2}-\sin^{2}\theta_{\varphi}\left(1-e^{-\left|\tilde{g}_{\varphi}\right|\Gamma\left(t-t_{in.}\right)}\right)\right]\sigma^{z}\\\
~{}~{}~{}~{}~{}~{}~{}-\frac{i}{2}\frac{\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}{\varepsilon_{\varphi}^{2}+\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}\sin^{2}\left(2\theta_{\varphi}\right)\left(1-e^{-\left|\tilde{g}_{\varphi}\right|\Gamma\left(t-t_{in.}\right)}\right)\left(\begin{array}[]{cc}\cos^{2}\theta_{\varphi}&0\\\
0&\sin^{2}\theta_{\varphi}\\\ \end{array}\right)}$ (69)

By applying Eq. 31 it is straightforward to find the covariance matrix:

$\eqalign{C_{jk,uv}\left(t\right)=\frac{2}{\pi}\int\limits_{0}^{\pi}{d\varphi}\left(\begin{array}[]{cc}0&\cos\left(\varphi\left(j-k\right)-2\theta_{\varphi}\right)\\\
-\cos\left(\varphi\left(j-k\right)+2\theta_{\varphi}\right)&0\\\
\end{array}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\left[\frac{1}{2}-\sin^{2}\theta_{\varphi}\left(1-e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left(t-t_{in.}\right)}\right)\right]\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathrm{
}+\frac{1}{2\pi}\int\limits_{0}^{\pi}{d\varphi}\frac{\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}{\varepsilon_{\varphi}^{2}+\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}\sin^{2}\left(2\theta_{\varphi}\right)\left(1-e^{-\left|\tilde{g}_{\varphi}\right|\Gamma\left(t-t_{in.}\right)}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\left(\begin{array}[]{cc}\sin\varphi\left(j-k\right)&\cos\left(\varphi\left(j-k\right)-2\theta_{\varphi}\right)\cos\left(2\theta_{\varphi}\right)\\\
-\cos\left(\varphi\left(j-k\right)+2\theta_{\varphi}\right)\cos\left(2\theta_{\varphi}\right)&\sin\varphi\left(j-k\right)\\\
\end{array}\right)}$ (70)

At the limit $t\rightarrow\infty$, we get the covariance at the steady state:

$\eqalign{C_{jk,uv}^{\left(0\right)}=\frac{1}{\pi}\int\limits_{0}^{\pi}{d\varphi}cos\left(2\theta_{\varphi}\right)\left(\begin{array}[]{cc}0&\cos\left(\varphi\left(j-k\right)-2\theta_{\varphi}\right)\\\
-\cos\left(\varphi\left(j-k\right)+2\theta_{\varphi}\right)&0\\\
\end{array}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}+\frac{1}{2\pi}\int\limits_{0}^{\pi}{d\varphi}\frac{\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}{\varepsilon_{\varphi}^{2}+\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}\sin^{2}\left(2\theta_{\varphi}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}\times\left(\begin{array}[]{cc}\sin\varphi\left(j-k\right)&\cos\left(\varphi\left(j-k\right)-2\theta_{\varphi}\right)\cos\left(2\theta_{\varphi}\right)\\\
-\cos\left(\varphi\left(j-k\right)+2\theta_{\varphi}\right)\cos\left(2\theta_{\varphi}\right)&\sin\varphi\left(j-k\right)\\\
\end{array}\right)}$ (71)

The existence of the second term in the rhs of Eqs. 69 \- 71 is due to the
fact that the dissipation kernel in Eq. 62 has non-zero off-diagonal entries.
An important observation is that both the functions 70 and 71 are non-analytic
at the point $\left|h\right|=1$. Firstly, we neglect the off-diagonal
contribution of the dissipation kernel by assuming that
$\left(g\Gamma/2\right)\ll 1$. This assumption simplifies the results without
qualitatively changing them. In this case the steady-state covariance reads:

$C_{jk,uv}^{\left(0\right)}\simeq\frac{1}{\pi}\int\limits_{0}^{\pi}{d\varphi}\cos\left(2\theta_{\varphi}\right)\left(\begin{array}[]{cc}0&\cos\left(\varphi\left(j-k\right)-2\theta_{\varphi}\right)\\\
-\cos\left(\varphi\left(j-k\right)+2\theta_{\varphi}\right)&0\\\
\end{array}\right)$ (72)

If the environment were absent ($\Gamma=0$) the corresponding result would
have the form (see Eq. 70):

$C_{jk,uv}^{\left(0\right)}\simeq\frac{1}{\pi}\int\limits_{0}^{\pi}{d\varphi}\left(\begin{array}[]{cc}0&\cos\left(\varphi\left(j-k\right)-2\theta_{\varphi}\right)\\\
-\cos\left(\varphi\left(j-k\right)+2\theta_{\varphi}\right)&0\\\
\end{array}\right).$ (73)

This function represents the covariance at the ground state of the isolated
system and it is non-analytic. The integrals involved, are connected with the
complete elliptic integral of the second kind that has a branch cut at
$\left|h\right|=1$. This non-analytic behavior reflects the underlying quantum
phase transition in the closed system. When the system is open the situation
changes but non-analyticity persists. This is immediately confirmed by
considering 72 for $j=k$:

$C_{jj,uv}^{\left(0\right)}\simeq\frac{1}{\pi}\int\limits_{0}^{\pi}{d\varphi}\cos^{2}\left(2\theta_{\varphi}\right)\left(\begin{array}[]{cc}0&1\\\
-1&0\\\
\end{array}\right)=\frac{1}{\pi}\int\limits_{0}^{\pi}{d\varphi}\frac{\left(h+\cos\varphi\right)^{2}}{h^{2}+2h\cos\varphi+1}\left(\begin{array}[]{cc}0&1\\\
-1&0\\\ \end{array}\right)$ (74)

A simple calculation yields the result:

$C_{jj,uv}^{\left(0\right)}\left(h\right)\simeq\left[\frac{1}{2}\theta\left(1-\left|h\right|\right)+\left(1-\frac{1}{2h^{2}}\right)\theta\left(\left|h\right|-1\right)\right]\left(\begin{array}[]{cc}0&1\\\
-1&0\\\ \end{array}\right)$ (75)

and

$\frac{d}{d\left|h\right|}C_{jj,uv}^{\left(0\right)}\left(h\right)=\left(\begin{array}[]{cc}0&1\\\
-1&0\\\
\end{array}\right)\times\left\\{\begin{array}[]{cc}0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left|h\right|<1\\\
{{1}/{\left|h\right|^{3}}},~{}~{}~{}\left|h\right|>1\\\ \end{array}\right.$
(76)

Thus, in the presence of dissipation the steady state’s covariance matrix
possesses a finite discontinuity at the point $\left|h\right|=1$. This point
of non-analyticity coincides with the critical point of the, well-known,
quantum phase transition that characterizes the corresponding isolated system
and its equivalent, the transverse Ising model.

This result indicates that, although differently, non-analyticity remains in
the steady state of the open system despite the impact of the environment.

For $j-k=L\rightarrow+\infty$ a straightforward calculation yields the result:

$C_{jk,uv}^{\left(0\right)}\left(h\right)\underset{L\rightarrow+\infty}{\simeq}\frac{1-h^{2}}{2h^{2}}\left(\begin{array}[]{cc}0&h^{L}\theta\left(1-\left|h\right|\right)\\\
\frac{1}{h^{L}}\theta\left(\left|h\right|-1\right)&0\\\ \end{array}\right)$
(77)

and

$C_{jk,uv}^{\left(0\right)}\left(h\right)\underset{L\rightarrow-\infty}{\simeq}\frac{h^{2}-1}{2h^{2}}\left(\begin{array}[]{cc}0&\frac{1}{h^{\left|L\right|}}\theta\left(\left|h\right|-1\right)\\\
h^{\left|L\right|}\theta\left(1-\left|h\right|\right)&0\\\
\end{array}\right).$ (78)

Combining these expressions, we can infer that for
$\left|L\right|\rightarrow\infty$ and
$\left|\left|h\right|-1\right|\rightarrow 0$ the steady state covariance
behaves as [4, 20]:

$C_{jk,uv}^{\left(0\right)}\left(h\right)\stackrel{{\scriptstyle[}}{{\left}}|L\right|\rightarrow\infty;\left|h\right|\rightarrow
1]{}{\thicksim}\left|\left|h\right|-1\right|\exp\left(-\left|L\right|\left|\left|h\right|-1\right|\right).$
(79)

The last expression indicates the fact that the function
$\frac{d}{d\left|h\right|}C_{jk,uv}^{\left(0\right)}\left(h\right)$ is
characterized, at the limit $\left|h\right|\rightarrow 1$, by long range
correlations the length of which has the form
$\xi\sim\frac{1}{\left|\left|h\right|-1\right|}\rightarrow\infty$.

##  5 Conclusions

We developed a framework for the construction of the time-dependent
correlation functions of an open fermionic system. It is a framework, that is
based on the merge of the Feynman-Vernon’s technique with the Keldysh-
Schwinger formalism on the ground of the coherent state path integrals. By
integrating out the environmental degrees of freedom, we defined the
generating functional through which every system’s correlation function can be
derived by applying the appropriate functional derivatives. We presented the
details of this construction in the case of a system that is quadratic when
written in terms of Majorana variables. In the last Section of the work, we
applied the proposed formalism to examine a model which is completely solvable
when isolated namely, the transverse Ising model. We calculate the exact form
of the covariance matrix at the steady state. We confirm the it is non-
analytic at the same point at which the corresponding isolated model develops
a quantum phase transition and we examine the specific form in which this non-
analyticity is realized.

## Appendix A

In this Appendix we solve the Green’s equation 66 of the main text:

$\left(i\partial_{t_{P}}+\varepsilon_{\varphi}\sigma^{z}\right)\mathbf{L}_{\varphi}\left(t_{P},t^{\prime}_{P^{\prime}}\right)+i\int\limits_{P^{\prime\prime}}{dt^{\prime\prime}_{P^{\prime\prime}}}\mathbf{K}_{\varphi}\left(t_{P},t^{\prime\prime}_{P^{\prime\prime}}\right)\mathbf{L}_{\varphi}\left(t^{\prime\prime}_{P^{\prime\prime}},t^{\prime}_{P^{\prime}}\right)=\mathbf{I}\delta\left(t_{P}-t^{\prime}_{P^{\prime}}\right).$
(80)

The solution we are interested for, must obey the boundary conditions:

$\mathbf{L}\left(t_{P}=t_{in.+},t^{\prime}_{P^{\prime}}\right)=\mathbf{L}\left(t_{P},t^{\prime}_{P^{\prime}}=t_{in.-}\right)=0,~{}~{}~{}\mathbf{L}\left(t_{P},t^{\prime}_{P^{\prime}}\right)\stackrel{{\scriptstyle[}}{{\left|t-t^{\prime}\right|\rightarrow+\infty}}]{}{\rightarrow}0.$
(81)

The first demand expresses the causal character of propagation. The second one
is connected with the decoherence encoded into the third term of Eq. 80, the
term that brings the influence of the environment into the propagation.

The dissipation kernel appearing in last equation has the form:

$\mathbf{K}_{\varphi}=\tilde{g}_{\varphi}\left(\begin{array}[]{cc}\Delta\left(t_{P},t^{\prime}_{P^{\prime}}\right)\cos^{2}\theta_{\varphi}-\Delta\left(t^{\prime}_{P^{\prime}},t_{P}\right)\sin^{2}\theta_{\varphi}&-\frac{i}{2}\left(\Delta\left(t_{P},t^{\prime}_{P^{\prime}}\right)+\Delta\left(t^{\prime}_{P^{\prime}},t_{P}\right)\right)\sin\left(2\theta_{\varphi}\right)\\\
\frac{i}{2}\left(\Delta\left(t_{P},t^{\prime}_{P^{\prime}}\right)+\Delta\left(t^{\prime}_{P^{\prime}},t_{P}\right)\right)\sin\left(2\theta_{\varphi}\right)&\Delta\left(t_{P},t^{\prime}_{P^{\prime}}\right)\sin^{2}\theta_{\varphi}-\Delta\left(t^{\prime}_{P^{\prime}},t_{P}\right)\cos^{2}\theta_{\varphi}\\\
\end{array}\right)$ (82)

The solution of the Green’s equation is a matrix of Green’s functions:

$\mathbf{L}_{\varphi}=\left(\begin{array}[]{cc}L_{\varphi,00}&L_{\varphi,01}\\\
L_{\varphi,10}&L_{\varphi,11}\\\ \end{array}\right)$ (83)

The dissipation kernel splits the set of equations 80 to two groups of coupled
equations:

$\eqalign{\partial_{t_{P}}L_{\varphi,00}\left(t_{P},t^{\prime}_{P}\right)-i\varepsilon_{\varphi}L_{\varphi,00}\left(t_{P},t^{\prime}_{P}\right)+\int\limits_{P^{\prime\prime}}{{dt^{\prime\prime}}_{P}}K_{\varphi,00}\left(t_{P},{t^{\prime\prime}}_{P}\right)L_{\varphi,00}\left({t^{\prime\prime}}_{P},t^{\prime}_{P}\right)\\\
~{}~{}~{}~{}~{}~{}+\int\limits_{P^{\prime\prime}}{{dt^{\prime\prime}}_{P}}K_{\varphi,01}\left(t_{P},{t^{\prime\prime}}_{P}\right)L_{\varphi,10}\left({t^{\prime\prime}}_{P},t^{\prime}_{P}\right)=-i\delta\left(t_{p}-t^{\prime}_{P}\right)\\\
\partial_{t_{P}}L_{\varphi,10}\left(t_{P},t^{\prime}_{P}\right)+i\varepsilon_{\varphi}L_{\varphi,10}\left(t_{P},t^{\prime}_{P}\right)+\int\limits_{P^{\prime\prime}}{{dt^{\prime\prime}}_{P}}K_{\varphi,10}\left(t_{P},{t^{{}^{\prime\prime}}}_{P}\right)L_{\varphi,00}\left({t^{\prime\prime}}_{P},t^{\prime}_{P}\right)\\\
~{}~{}~{}~{}~{}~{}+\int\limits_{P^{\prime\prime}}{{dt^{\prime\prime}}_{P}}K_{\varphi,11}\left(t_{P},{t^{\prime\prime}}_{P}\right)L_{\varphi,10}\left({t^{\prime\prime}}_{P},t^{\prime}_{P}\right)=0}$
(84)

and

$\eqalign{\partial_{t_{P}}L_{\varphi,11}\left(t_{P},t^{\prime}_{P}\right)+i\varepsilon_{\varphi}L_{\varphi,11}\left(t_{P},t^{\prime}_{P}\right)+\int\limits_{P^{\prime\prime}}{{dt^{\prime\prime}}_{P}}K_{\varphi,10}\left(t_{P},{t^{\prime\prime}}_{P}\right)L_{\varphi,01}\left({t^{\prime\prime}}_{P},t^{\prime}_{P}\right)\\\
~{}~{}~{}~{}~{}~{}+\int\limits_{P^{\prime\prime}}{{dt^{\prime\prime}}_{P}}K_{\varphi,11}\left(t_{P},{t^{\prime\prime}}_{P}\right)L_{\varphi,11}\left({t^{\prime\prime}}_{P},t^{\prime}_{P}\right)=-i\delta\left(t_{p}-t^{\prime}_{P}\right)\\\
\partial_{t_{P}}L_{\varphi,01}\left(t_{P},t^{\prime}_{P}\right)-i\varepsilon_{\varphi}L_{\varphi,01}\left(t_{P},t^{\prime}_{P}\right)+\int\limits_{P^{\prime\prime}}{{dt^{\prime\prime}}_{P}}K_{\varphi,00}\left(t_{P},{t^{\prime\prime}}_{P}\right)L_{\varphi,01}\left({t^{\prime\prime}}_{P},t^{\prime}_{P}\right)\\\
~{}~{}~{}~{}~{}~{}+\int\limits_{P^{\prime\prime}}{{dt^{\prime\prime}}_{P}}K_{\varphi,01}\left(t_{P},{t^{\prime\prime}}_{P}\right)L_{\varphi,11}\left({t^{\prime\prime}}_{P},t^{\prime}_{P}\right)=0.}$
(85)

Restoring the explicit form of the Keldysh variables each of the groups 84
and85 splits into two distinct groups of four coupled equations. As the
calculations for each one of the $16$ groups have the same structure it is
enough to present here only those groups that yield the result for
$\mathbf{L}_{\varphi}\left(t_{+},t^{\prime}_{+}\right)$.

In the following, and for reasons of simplicity, we have set $b=0$ and $\delta
E=0$ in expressions 41 and we have used the abbreviations:

$A_{\varphi}=\tilde{g}_{\varphi}\frac{\Gamma}{2}\cos\left(2\theta_{\varphi}\right),~{}~{}~{}~{}B_{\varphi}=\tilde{g}_{\varphi}\frac{\Gamma}{2}\sin\left(2\theta_{\varphi}\right).$
(86)

After these, the first tetrad of equations results from Eq. 84 and reads as
follows:

$\eqalign{\left(\partial_{t}-i\varepsilon_{\varphi}+A_{\varphi}\right)L_{\varphi,00}\left(t_{+},t_{+}\right)+\Gamma\tilde{g}_{\varphi}\sin^{2}\theta_{\varphi}L_{\varphi,00}\left(t_{-},t_{+}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-iB_{\varphi}\left(L_{\varphi,10}\left(t_{+},t_{+}\right)-L_{\varphi,10}\left(t_{-},t_{+}\right)\right)=-i\delta\left(t-t^{\prime}\right),\\\
\Gamma\tilde{g}_{\varphi}\cos^{2}\theta_{\varphi}L_{\varphi,00}\left(t_{+},t_{+}\right)+\left(\partial_{t}-i\varepsilon_{\varphi}-A_{\varphi}\right)L_{\varphi,00}\left(t_{-},t_{+}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-iB_{\varphi}\left(L_{\varphi,10}\left(t_{+},t_{+}\right)-L_{\varphi,10}\left(t_{-},t_{+}\right)\right)=0,\\\
iB_{\varphi}\left(L_{\varphi,00}\left(t_{+},t_{+}\right)-L_{\varphi,00}\left(t_{-},t_{+}\right)\right)+\left(\partial_{t}+i\varepsilon_{\varphi}-A_{\varphi}\right)L_{\varphi,10}\left(t_{+},t_{+}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\Gamma\tilde{g}_{\varphi}\cos^{2}\theta_{\varphi}L_{\varphi,10}\left(t_{-},t_{+}\right)=0,\\\
iB_{\varphi}\left(L_{\varphi,00}\left(t_{+},t_{+}\right)-L_{\varphi,00}\left(t_{-},t_{+}\right)\right)+\Gamma\tilde{g}_{\varphi}\sin^{2}\theta_{\varphi}L_{\varphi,10}\left(t_{+},t_{+}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(\partial_{t}+i\varepsilon_{\varphi}+A_{\varphi}\right)L_{\varphi,10}\left(t_{-},t_{+}\right)=0.}$
(87)

The set of Eqs. 87 can be exactly solved. Without any boundary restriction it
is a straightforward exercise to find:

$\eqalign{L_{\varphi,00}\left(t_{+},t^{\prime}_{+}\right)=-i\cos^{2}\theta_{\varphi}e^{i\varepsilon_{\varphi}\left(t-t^{\prime}\right)-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left(t-t^{\prime}\right)}\theta\left(t-t^{\prime}\right)+i\sin^{2}\theta_{\varphi}e^{i\varepsilon_{\varphi}\left(t-t^{\prime}\right)-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left(t^{\prime}-t\right)}\theta\left(t^{\prime}-t\right)\\\
~{}~{}~{}~{}~{}~{}-\frac{i}{2}\frac{\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}{\varepsilon_{\varphi}^{2}+\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}\sin^{2}\left(2\theta_{\varphi}\right)\cos^{2}\theta_{\varphi}e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left|t-t^{\prime}\right|}\left(\cos\left(\varepsilon_{\varphi}\left|t-t^{\prime}\right|\right)+\tilde{g}_{\varphi}\frac{\Gamma}{2}\frac{\sin\left(\varepsilon_{\varphi}\left|t-t^{\prime}\right|\right)}{\varepsilon_{\varphi}}\right)}$
(88)

and

$L_{\varphi,10}\left(t_{+},t^{\prime}_{+}\right)=-\frac{1}{2}e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left(t-t^{\prime}\right)}\tilde{g}_{\varphi}\frac{\Gamma}{2}\frac{\sin\varepsilon_{\varphi}\left(t-t^{\prime}\right)}{\varepsilon_{\varphi}}\theta\left(t-t^{\prime}\right).$
(89)

The second group we are interested for is connected with Eq. 85 and assumes
the form:

$\eqalign{\left(\partial_{t}+i\varepsilon_{\varphi}-A_{\varphi}\right)L_{\varphi,11}\left(t_{+},t^{\prime}_{+}\right)+\Gamma\tilde{g}_{\varphi}\cos^{2}\theta_{\varphi}L_{\varphi,11}\left(t_{-},t^{\prime}_{+}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+iB_{\varphi}\left(L_{\varphi,01}\left(t_{+},t^{\prime}_{+}\right)-L_{\varphi,01}\left(t_{-},t^{\prime}_{+}\right)\right)=-i\delta\left(t-t^{\prime}\right),\\\
\Gamma\tilde{g}_{\varphi}\sin^{2}\theta_{\varphi}L_{\varphi,11}\left(t_{+},t^{\prime}_{+}\right)+\left(\partial_{t}+i\varepsilon_{\varphi}+A_{\varphi}\right)L_{\varphi,11}\left(t_{-},t^{\prime}_{+}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+iB_{\varphi}\left(L_{\varphi,01}\left(t_{+},t^{\prime}_{+}\right)-L_{\varphi,01}\left(t_{-},t^{\prime}_{+}\right)\right)=0,\\\
-iB_{\varphi}\left(L_{\varphi,11}\left(t_{+},t^{\prime}_{+}\right)-L_{\varphi,11}\left(t_{-},t^{\prime}_{+}\right)\right)+\left(\partial_{t}-i\varepsilon_{\varphi}+A_{\varphi}\right)L_{\varphi,10}\left(t_{+},t^{\prime}_{+}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\Gamma\tilde{g}_{\varphi}\sin^{2}\theta_{\varphi}L_{\varphi,01}\left(t_{-},t^{\prime}_{+}\right)=0,\\\
-iB_{\varphi}\left(L_{\varphi,11}\left(t_{+},t^{\prime}_{+}\right)-L_{\varphi,11}\left(t_{-},t^{\prime}_{+}\right)\right)+\Gamma\tilde{g}_{\varphi}\cos^{2}\theta_{\varphi}L_{\varphi,01}\left(t_{+},t^{\prime}_{+}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(\partial_{t}-i\varepsilon_{\varphi}-A_{\varphi}\right)L_{\varphi,01}\left(t_{-},t_{+}\right)=0.}$
(90)

The solution of this system can be deduced from the solution of the system 87
by observing that the coefficients in 90 can be produced from the coefficients
in (A.10) by making the replacement
$\theta_{\varphi}\rightarrow\theta_{\varphi}\pm{{\pi}/{2}}$. This makes easy
to find:

$\eqalign{L_{\varphi,11}\left(t_{+},t^{\prime}_{+}\right)=-i\sin^{2}\theta_{\varphi}e^{i\varepsilon_{\varphi}\left(t-t^{\prime}\right)-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left(t-t^{\prime}\right)}\theta\left(t-t^{\prime}\right)+i\cos^{2}\theta_{\varphi}e^{i\varepsilon_{\varphi}\left(t-t^{\prime}\right)-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left(t^{\prime}-t\right)}\theta\left(t^{\prime}-t\right)\\\
~{}~{}~{}~{}~{}~{}-\frac{i}{2}\frac{\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}{\varepsilon_{\varphi}^{2}+\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}\sin^{2}\left(2\theta_{\varphi}\right)\sin^{2}\theta_{\varphi}e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left|t-t^{\prime}\right|}\left(\cos\left(\varepsilon_{\varphi}\left|t-t^{\prime}\right|\right)+\tilde{g}_{\varphi}\frac{\Gamma}{2}\frac{\sin\left(\varepsilon_{\varphi}\left|t-t^{\prime}\right|\right)}{\varepsilon_{\varphi}}\right)}$
(91)

and

$L_{\varphi,01}\left(t_{+},t^{\prime}_{+}\right)=L_{\varphi,10}\left(t_{+},t^{\prime}_{+}\right).$
(92)

As discussed in the main text (see Eqs. 64 and 65) the Green’s matrix
$\mathbf{L}$ propagates a two-component fermion field:

$\psi_{\varphi,u}^{cl.}\left(t_{P}\right)=\sum_{\nu}{\int\limits_{P^{\prime}}{dt^{\prime}_{P^{\prime}}}}L_{\varphi,u\nu}\left(t_{P},t^{\prime}_{P^{\prime}}\right)j_{\varphi,\nu}\left(t^{\prime}_{P^{\prime}}\right).$
(93)

One of the components is propagating forward in time, while the second one,
backwards. Based on this observation, we construct, the appropriate for the
present problem, Green’s functions by imposing the following boundary
conditions:

$\eqalign{L_{\varphi,00}\left(t_{+}=t_{in.+},t^{\prime}_{+}\right)=0,~{}~{}~{}L_{\varphi,01}\left(t_{+}=t_{in.+},t^{\prime}_{+}\right)=0\\\
L_{\varphi,10}\left(t_{+},t^{\prime}_{+}=t_{in.+}\right)=0,~{}~{}~{}L_{\varphi,11}\left(t_{+},t^{\prime}_{+}=t_{in.+}\right)=0.}$
(94)

To obtain these functions we can add at the Green’s functions we found, the
general solution of the corresponding homogenous equation

$L_{\varphi,u\nu}\rightarrow
L_{\varphi,u\nu}+\sum_{j=1}^{4}{a_{\varphi,u\nu}^{\left(j\right)}}e^{ip_{j}t},~{}~{}~{}p_{j}=\pm\varepsilon_{\varphi}\pm
i\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma$ (95)

and determine the coefficients to satisfy 94 together with
$L_{\varphi,u\nu}\stackrel{{\scriptstyle[}}{{\left|t-t^{\prime}\right|\rightarrow+\infty}}]{}{\rightarrow}0$.
In this way we find:

$\eqalign{L_{\varphi,00}\left(t_{+},t^{\prime}_{+}\right)=-ie^{i\varepsilon_{\varphi}\left(t-t^{\prime}\right)-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left|t-t^{\prime}\right|}\theta\left(t-t^{\prime}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}+iF_{\varphi,00}\left(t-t^{\prime}\right)\left[e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left|t-t^{\prime}\right|}-e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left(t+t^{\prime}-2t_{in.}\right)}\right]}$
(96)

and

$\eqalign{L_{\varphi,11}\left(t_{+},t^{\prime}_{+}\right)=ie^{i\varepsilon_{\varphi}\left(t-t^{\prime}\right)-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left|t-t^{\prime}\right|}\theta\left(t^{\prime}-t\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}-iF_{\varphi,11}\left(t-t^{\prime}\right)\left[e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left|t-t^{\prime}\right|}-e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left(t+t^{\prime}-2t_{in.}\right)}\right]}$
(97)

The abbreviations in these expressions read as follows:

$\eqalign{F_{\varphi,00}=\sin^{2}\theta_{\varphi}e^{i\varepsilon_{\varphi}\left(t-t^{\prime}\right)}-\frac{1}{2}\frac{\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}{\varepsilon_{\varphi}^{2}+\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}\sin^{2}\left(2\theta_{\varphi}\right)\cos^{2}\theta_{\varphi}\\\
~{}~{}~{}~{}~{}~{}~{}~{}\times\left(\cos\left(\varepsilon_{\varphi}\left|t-t^{\prime}\right|\right)+\tilde{g}_{\varphi}\frac{\Gamma}{2}\frac{\sin\left(\varepsilon_{\varphi}\left|t-t^{\prime}\right|\right)}{\varepsilon_{\varphi}}\right)}$
(98)

and

$\eqalign{F_{\varphi,11}=\sin^{2}\theta_{\varphi}e^{i\varepsilon_{\varphi}\left(t-t^{\prime}\right)}+\frac{1}{2}\frac{\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}{\varepsilon_{\varphi}^{2}+\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}\sin^{2}\left(2\theta_{\varphi}\right)\sin^{2}\theta_{\varphi}\\\
~{}~{}~{}~{}~{}~{}~{}~{}\times\left(\cos\left(\varepsilon_{\varphi}\left|t-t^{\prime}\right|\right)+\tilde{g}_{\varphi}\frac{\Gamma}{2}\frac{\sin\left(\varepsilon_{\varphi}\left|t-t^{\prime}\right|\right)}{\varepsilon_{\varphi}}\right).}$
(99)

In the same way we find

$L_{\varphi,01}\left(t_{+},t^{\prime}_{+}\right)=-\frac{1}{2}e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\varGamma\left(t-t^{\prime}\right)}\tilde{g}_{\varphi}\frac{\Gamma}{2}\frac{\sin\varepsilon_{\varphi}\left(t-t^{\prime}\right)}{\varepsilon_{\varphi}}\theta\left(t-t^{\prime}\right)$
(100)

and

$\eqalign{L_{\varphi,10}\left(t_{+},t^{\prime}_{+}\right)=\frac{1}{2}e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left|t-t^{\prime}\right|}\tilde{g}_{\varphi}\frac{\Gamma}{2}\frac{\sin\varepsilon_{\varphi}\left|t-t^{\prime}\right|}{\varepsilon_{\varphi}}\theta\left(t^{\prime}-t\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}-e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left|t-t^{\prime}\right|}\tilde{g}_{\varphi}\frac{\Gamma}{2}\frac{\sin\varepsilon_{\varphi}\left|t-t^{\prime}\right|}{\varepsilon_{\varphi}}\left[e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left|t-t^{\prime}\right|}-e^{-\frac{1}{2}\left|\tilde{g}_{\varphi}\right|\Gamma\left(t+t^{\prime}-2t_{in.}\right)}\right].}$
(101)

Combining these results, we find that for $t_{+}=t^{\prime}_{+}=t$:

$\eqalign{\mathbf{L}_{\varphi}=-i\left[\frac{1}{2}-\sin^{2}\theta_{\varphi}\left(1-e^{-\left|\tilde{g}_{\varphi}\right|\Gamma\left(t-t_{in.}\right)}\right)\right]\sigma^{z}\\\
~{}~{}~{}~{}~{}~{}~{}-\frac{i}{2}\frac{\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}{\varepsilon_{\varphi}^{2}+\left(\tilde{g}_{\varphi}\frac{\Gamma}{2}\right)^{2}}\sin^{2}\left(2\theta_{\varphi}\right)\left(1-e^{-\left|\tilde{g}_{\varphi}\right|\Gamma\left(t-t_{in.}\right)}\right)\left(\begin{array}[]{cc}\cos^{2}\theta_{\varphi}&0\\\
0&\sin^{2}\theta_{\varphi}\\\ \end{array}\right)}$ (102)

## References

## References

  * [1] T. J. Osborne, M. A. Nielsen, Entanglement in a simple quantum phase transition, Physical Review A 66 (3) (2002) 032110.  doi:10.1103/physreva.66.032110.
  * [2] G. Vidal, J. I. Latorre, E. Rico, A. Kitaev, Entanglement in quantum critical phenomena, Physical Review Letters 90 (22) (2003) 227902.  doi:10.1103/physrevlett.90.227902.
  * [3] S. Sachdev, Quantum Phase Transitions, 2nd Edition, Cambridge University Press, 2011\.  doi:10.1017/CBO9780511973765.
  * [4] M. Höning, M. Moos, M. Fleischhauer, Critical exponents of steady-state phase transitions in fermionic lattice models, Physical Review A 86 (1) (2012) 013606\.  doi:10.1103/physreva.86.013606.
  * [5] B. Horstmann, J. I. Cirac, G. Giedke, Noise-driven dynamics and phase transitions in fermionic systems, Physical Review A 87 (1) (2013) 012108.  doi:10.1103/physreva.87.012108.
  * [6] B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, P. Zoller, Preparation of entangled states by quantum markov processes, Physical Review A 78 (4) (2008) 042307.  doi:10.1103/physreva.78.042307.
  * [7] F. Verstraete, M. M. Wolf, J. I. Cirac, Quantum computation and quantum-state engineering driven by dissipation, Nature Physics 5 (9) (2009) 633–636.  doi:10.1038/nphys1342.
  * [8] F. Pastawski, L. Clemente, J. I. Cirac, Quantum memories based on engineered dissipation, Physical Review A 83 (1) (2011) 012304.  doi:10.1103/physreva.83.012304.
  * [9] C. A. Muschik, E. S. Polzik, J. I. Cirac, Dissipatively driven entanglement of two macroscopic atomic ensembles, Physical Review A 83 (5) (2011) 052312.  doi:10.1103/physreva.83.052312.
  * [10] T. Prosen, Third quantization: a general method to solve master equations for quadratic open fermi systems, New Journal of Physics 10 (4) (2008) 043026.  doi:10.1088/1367-2630/10/4/043026.
  * [11] T. Prosen, B. Žunkovič, Exact solution of markovian master equations for quadratic fermi systems: thermal baths, open XY spin chains and non-equilibrium phase transition, New Journal of Physics 12 (2) (2010) 025016\.  doi:10.1088/1367-2630/12/2/025016.
  * [12] I. Lyris, P. Lykourgias, A. I. Karanikas, Reduced dynamics in open bosonic and fermionic systems, Open Systems & Information Dynamics 28 (03) (sep 2021).  doi:10.1142/s1230161221500116.
  * [13] I. Lyris, P. Lykourgias, A. Karanikas, Fermionic path integrals and correlation dynamics in a 1d XY system, Annals of Physics 421 (2020) 168286.  doi:10.1016/j.aop.2020.168286.
  * [14] L. Keldysh, Zh. Eksp. Teor. Fiz. 47 (1964) 1515–1527.
  * [15] L. Keldysh, Sov. Phys. JETP 20 (1965) 2018.
  * [16] H.-P. Breuer, F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press, USA, 2007.
  * [17] J. Preskill, Quantum Information and Computation, CreateSpace Independent Publishing Platform, 2015.
  * [18] A. Kitaev, C. Laumann, Topological phases and quantum computation (Apr. 2009).  arXiv:0904.2771, doi:10.48550/ARXIV.0904.2771.
  * [19] T.-C. Wei, S. Vishveshwara, P. M. Goldbart, Global geometric entanglement in transverse-field XY spin chains: finite and infinite systems, Quantum Information and Computation 11 (3&4) (2011) 326–354.  doi:10.26421/qic11.3-4-10.
  * [20] J. Eisert, T. Prosen, Noise-driven quantum criticality (Dec. 2010).  arXiv:1012.5013, doi:10.48550/ARXIV.1012.5013.