Analyzing the spectral density of a perturbed analog quantum simulator using Keldysh formalism

TL;DR

This paper investigates how decoherence affects the spectral density in analog quantum simulators by mapping them onto fermionic systems and applying Keldysh formalism, revealing the limits of spectral feature resolution.

Contribution

It introduces a diagrammatic Keldysh approach to analyze the spectral density of noisy quantum simulators, bridging quantum simulation and many-body physics.

Findings

Features stronger than the decoherence rate are resolvable.

Weaker spectral features are washed out by noise.

Keldysh approach agrees with master-equation calculations for small systems.

Abstract

Simulation of interacting electron systems is one of the great challenges of modern quantum chemistry and solid state physics. Controllable quantum systems offer the opportunity to create artificial structures which mimic the system of interest. An interesting quantity to extract from these quantum simulations is the spectral function. We map a noisy quantum simulator onto a fermionic system and investigate the influence of decoherence on the simulation of the spectral density using a diagrammatic approach on Keldysh contour. We show that features stronger than the single-qubit decoherence rate can be resolved, while weaker features wash out. For small systems, we compare our Keldysh approach to master-equation calculations.

Equations

$$G(\omega )=(G_0^{-1}(\omega )-\Sigma (\omega ){)}^{-1}.$$

$$\{c_i,c_j^{†}\}={\delta }_{ij},\{c_i,c_j\}=\{c_i^{†},c_j^{†}\}=0,$$

$$c_i=\prod _{j<i}\left(-{\sigma }_j^z\right){\sigma }_i^{-}=e^{i{\varphi }_i}{\sigma }_i^{-}.$$

$$H_q=\sum _{ij}{\hat{\Phi }}_{ij}{t}_{ij}{\sigma }_i^{+}{\sigma }_j^{-},$$

$$H_{q,\alpha }=\sum _i{\sigma }_i^{\alpha }{\hat{X}}_i^{\alpha },$$

$$⟨{\hat{X}}_i^{\alpha }(t){\hat{X}}_j^{\beta }(t^{\prime }){⟩}_0={\delta }_{ij}{\delta }_{\alpha \beta }⟨{\hat{X}}_i^{\alpha }(t){\hat{X}}_i^{\alpha }(t^{\prime }){⟩}_0$$

$$H_{f,z}=\sum _i{c}_i^{†}{c}_i^{}{\hat{X}}_i^z.$$

$$H_{f,x}=\sum _i\prod _{j<i}(1-2c_j^{†}{c}_j^{})(c_i^{†}+c_i^{}){\hat{X}}_i^x,$$

$${\mathcal{H}}_f=H_f+H_{f,z}+H_{f,x}+H_B.$$

$$H_f=\sum _{ij}{t}_{ij}{c}_i^{†}{c}_j^{},$$

$$\dot{\rho }=-$$

$$+\frac{{\Gamma }_{1,n}}{2}(2{\sigma }_n^{-}\rho {\sigma }_n^{+}-[{\sigma }_n^{+}{\sigma }_n^{-},\rho {]}_{+})=\mathcal{L}\rho ,$$

$$⟨\hat{A}(t+\tau )\hat{B}(t)⟩$$

$$⟨\hat{A}(t)\hat{B}(t+\tau )⟩$$

$$G$$

$$=-i⟨\prod _{k<n}(-{\sigma }_k^z){\sigma }_n^{-}(t+\tau )\prod _{l<n^{\prime }}(-{\sigma }_l^z){\sigma }_{n^{\prime }}^{-}(t)⟩.$$

$$G(t,t^{\prime })=\left(\begin{array}{cc}{G}^{+}(t,t^{\prime })& G^K(t,t^{\prime })\\ 0& G^{-}(t,t^{\prime })\end{array}\right),$$

$$G_{ij}^{\pm }(t,t^{\prime })$$

$$G_{ij}^K$$

$$G_{ii}^K(t,t^{\prime })=-i⟨({\sigma }_i^{-}(t){\sigma }_i^{+}(t^{\prime })+{\sigma }_i^{+}(t^{\prime }){\sigma }_i^{-}(t)⟩.$$

$$i{G}_{ij}^C(t,t^{\prime })$$

$$=\mathrm{Tr}\left[T_C{U}_C{\tilde{c}}_i(t){\tilde{c}}_j^{†}(t^{\prime })\rho (z_0)\right],$$

$${\rho }_S(z_0)=Z_S^{-1}{e}^{-\beta {\sum }_{ij}{a}_{ij}{c}_i^{†}{c}_j^{}}=Z_S^{-1}{e}^{-\beta H_{\mathrm{ini}}},$$

$$\tilde{O}(t)=e^{i(H_f+H_B)t}O{e}^{-i(H_f+H_B)t}.$$

$$U_C=T_C\exp \left\{-i{\int }_{C}d{t}^{\prime }[{\tilde{H}}_{f,x}(t^{\prime })+{\tilde{H}}_{f,z}(t^{\prime })]\right\},$$

$${\hat{G}}_C(z_1,$$

$$+{\int }_{C}d\bar{z}d{\bar{z}}^{\prime }{\hat{G}}_{C,0}(z_1,\bar{z})\Sigma (\bar{z},{\bar{z}}^{\prime }){\hat{G}}_C({\bar{z}}^{\prime },z_2).$$

$$\rho (t_0)=e^{-\beta {\sum }_{ij}{a}_{ij}{c}_i^{†}{c}_j}=e^{-\beta H_{\mathrm{ini}}}.$$

$$\left(i\frac{d}{dt}-{\hat{H}}_f(t)\right)G^{\pm }(t,t^{\prime })$$

$$\left(i\frac{d}{dt}-{\hat{H}}_f(t)\right)G^K(t,t^{\prime })$$

$$I^K(t,t^{\prime })=[{\Sigma }^{+}\circ G^K$$

$$G^K(0,0)=-i[1-2f({\hat{H}}_{\mathrm{ini}})],$$

$$G^{+}(t,t^{\prime })$$

$$G^K(t,t^{\prime })$$

$$G^{\pm }(\omega )$$

Full text

I Introduction

Using quantum mechanics we can describe many interesting systems in a wide range of fields, covering biology, chemistry, and (solid state) physics. The interacting many-body Schrödinger equation provides a microscopic description of atoms, molecules, and solids. However, due to the exponential growth of the Hilbert space with the system size, the exact solution of the full interacting Schrödinger equation becomes numerically demanding even for systems of moderate size.D. Poulin et al. (2015) This problem is inherent to the simulation of quantum systems on conventional computers, limiting their application to approximations of the full Schrödinger equation. For example, density functional theory (DFT)Hohenberg and Kohn (1964) has been used with great success to calculate electronic properties of solids and molecules on classical computers for the past decades. However, DFT, as other approximation schemes, is limited to weakly interacting systems and many interesting systems such as high-T${}_{\text{C}}$ superconductivityBednorz and Müller (1986) or transition metals,Jiang et al. (2012); Cramer and Truhlar (2009) are beyond the scope of DFT. The idea to use well-controllable quantum systems to overcome these difficultiesFeynman (1982) promises an exponential speedup compared to classical numerics.Lloyd (1996); Zalka (1998) Quantum simulationJohnson et al. (2014); Buluta and Nori (2009); Cirac and Zoller (2012) on different physical systems, such as trapped ions,Blatt and Roos (2012) ultra cold gases,Jaksch and Zoller (2005); Bloch et al. (2012) or superconducting circuits,Makhlin et al. (2001); Devoret and Schoelkopf (2013); Houck et al. (2012); Barends et al. (2015); Heras et al. (2014) have been proposed to solve problems in different fields, such as quantum chemistry Whitfield et al. (2011); Kassal et al. (2011) and strongly-correlated electrons.Barends et al. (2015); Wecker et al. (2015); Bauer et al. (2016) The idea behind quantum simulation is to map the Hamiltonian onto a controllable quantum system – the quantum simulator – that mimics the simulated system.

One fundamental roadblock on the way to feasible quantum simulation is the limited coherence time of qubits.Ithier et al. (2005); Devoret and Schoelkopf (2013) Universal digital quantum simulation by discretization of the time evolution relies on the application of gates.Barends et al. (2015); Heras et al. (2014); Lanyon et al. (2011); Barends et al. (2016) While quantum error correction is possible,Fowler et al. (2012); Knill and Laflamme (1997) the number of coherently applied gates needed to calculate molecules of relevant size exceeds the capability of quantum hardware available in the near future.Wecker et al. (2014) Analog quantum simulation is a promising alternative for first applications for quantum simulation.Manousakis (2002); Sarovar et al. (2017) Many proposals for analog quantum simulation of spin systems,Porras and Cirac (2004); Kim et al. (2010); Greif et al. (2013); Friedenauer et al. (2008); Simon et al. (2011), HubbardReiner et al. (2016); Hart et al. (2015) and HolsteinMei et al. (2013) models, exciton transport,Mostame et al. (2012) and many other models exist or have been realized on different plattforms already.Wenz et al. (2013); Trotzky et al. (2012); Zohar et al. (2013)

At present it remains unclear how to analyze the effect of decoherence in large systems of coupled qubits forming an analog quantum simulator.Hauke et al. (2012); Schwenk et al. (2016) In this article, we demonstrate the application of established many-body physics methods on systems of many qubits to analyze the coherence behavior of a quantum simulator consisting of many fault-prone qubits. An important quantity of interest we would like to extract from a quantum simulator is the spectral function, which can also be formulated in terms of the Green’s function $G(\omega)$. We expect intuitively that decoherence rates of single qubits determine the spectral resolution of a quantum simulator. Using Wick’s theorem and many-body perturbation theory, we formulate a connection between the ideal Green’s function $G_{0}$ and the simulated, decoherence afflicted perturbed Green’s function $G$:

[TABLE]

The self-energy $\Sigma$ describes the influence of the environment on the quantum simulator and is defined as the sum of all irreducible diagrams. Using this connection we investigate how single-qubit decoherence affects the spectral function of several multi-qubit systems, in and out of equilibrium.

We study systems of qubits which can be mapped to a fermionic system via a Jordan-Wigner transformation.Jordan and Wigner (1928) In Sec. II, we discuss both directions of this mapping. We first introduce the Jordan-Wigner transformation to map a fermionic system on a quantum simulator. After this we show how to map the coupling of qubits to the environment onto fermions using the same transformation. With this transformation we can map the noisy quantum simulator to a system of fermions coupled to a bath, similar to electron-phonon coupling. This system can be tackled with standard many-body perturbation theory, introduced in Sec. III. There, we use the transformation to fermionic systems to obtain a diagrammatic expansion in the coupling to the environment on Keldysh contour. In particular, we derive Eq. (1). In Sec. IV, we apply this method to a simple chain of qubits with dephasing due to a bosonic bath and decay due to two-level systems. We compare our results for the steady-state properties with results obtained by master-equation methods. We show that we can qualitatively understand the decoherence of large systems with our method. In Sec. V, we show that our method can also be used to model transient evolution of a chain of qubits after an initialization into a non-thermal state. Additionally, we compare our results with corresponding master-equation simulations. Conclusions are given in Sec. VI.

II General Discussion

In this section, we first describe the mapping of a fermionic system to a quantum simulator consisting of spin-$\frac{1}{2}$ qubits. Afterwards, we use this method to map the Hamiltonian of a noisy quantum simulator back to a fermionic system which we need to be able to apply many-body perturbation theory.

II.1 Fermionic system mapped to qubits

In second quantization, the Hamiltonian $H_{f}$ of an electronic system depends on the creation (annihilation) operator $c_{i}^{\dagger}$ ($c_{i}$) of an electron in state $i$, where $i$ contains information about site, orbital, and spin. To simulate such a Hamiltonian on an analog quantum simulator consisting of coupled qubits, we have to map the fermionic operators to a spin-$\frac{1}{2}$ system that describes the physical qubits. This physical difference reflects into different commutation relations for spin-$\frac{1}{2}$ Pauli operators $\hat{\sigma}_{\pm}$ and fermionic operators.Whitfield et al. (2011) While the latter always anti-commute

[TABLE]

spin operators obey mixed commutation relations. On different sites, operators commute while on-site operators anti-commute. There exist different mapping methods between fermions and spin operators, such as the Jordan-Wigner transformation,Jordan and Wigner (1928) Majorana representation,Schad et al. (2015) or the Bravyi-Kitaev transformation.Bravyi and Kitaev (2002) Here, we use the most common method, the one dimensional Jordan-Wigner transformation (JWT)Jordan and Wigner (1928)

[TABLE]

The JWT reproduces correct commutation relations and preserves the dimensionality of the Hilbert space. This transformation requires ordering of the operators and, unfortunately, introduces non-local interactions between all the qubits of the quantum simulator. For example, the one-particle term $\sum_{ij}t_{ij}c_{i}^{\dagger}c_{j}$ transforms into

[TABLE]

with the Jordan-Wigner string $\hat{\Phi}_{ij}=\prod_{k=\mathrm{min}(i,j)}^{\mathrm{max}(i,j)}(-\sigma_{k}^{z})$ that accounts for the correct parity of the system. The transformed Hamiltonian $H_{q}$ can can either be implemented on a digital quantum computer or an analog quantum simulatorManousakis (2002). Either way, from these simulations we can in principle extract the systems Green’s functions, for example $G^{<}_{0,ij}(t,t^{\prime})\propto\langle c_{i}^{\phantom{\dagger}}(t)c_{j}^{\dagger}(t^{\prime})\rangle_{0}$,Tian et al. (2016); Dallaire-Demers and Wilhelm (2016) where the subscript 0 refers to the noise-free result. Unfortunately, a realistic quantum simulator suffers from decoherence. Qubits couple to uncontrolled environmental degrees of freedom, such as two-level systems,Schriefl et al. (2006) and information is lost to these environmental degrees of freedom.

Decoherence of single qubits and small systems of qubits has been analyzed with quantum master equations (QME). Carmichael (2002) The QME approach allows for excellent quantitative results for small systems. Here, we want to qualitatively understand the behavior of large systems. Instead of a master-equation approach we will map the faulty quantum simulator back to a fermionic system and expand the disturbed Green’s function $G(\omega)$, obtained from the simulation on the faulty quantum simulator, in powers of the coupling to the noise. This is identical to an expansion in terms of the ideal Green’s function of the unperturbed quantum simulator $G_{0}(\omega)$ which we would like to obtain from the simulation, and establishes the connection between perturbed and ideal Green’s function, Eq. (1).

II.2 Decoherence as a fermionic problem

The qubits of a quantum simulator couple to and thus dissipate information into the noisy environment. The interaction with the environment adds a perturbation of the form

[TABLE]

to the quantum simulator Hamiltonian, where $\alpha=x,\,z$ and $\hat{X}_{i}^{\alpha}$ is a bath operator that describes the coupling between qubit $i$ and its respective bath. Longitudinal coupling $\propto\sigma^{z}_{i}$ describes (pure) dephasing, while transverse coupling $\propto\sigma_{i}^{x}$ is responsible for relaxation. Their dynamics are governed by free bath Hamiltonians $H_{B,\alpha}$. We assume that each qubit couples to individual, uncorrelated baths for decay and dephasing, such that

[TABLE]

holds for all free bath correlation functions. Since Wick’s theorem does not apply for spin operators, we map noise and quantum simulator Hamiltonians back to fermionic operators. For dephasing the Jordan-Wigner transformation, Eq. (3), simply yields an on-site term of the form

[TABLE]

This mapping is quite general and can be used for any noise inducing dephasing. While longitudinal coupling is easy to handle after the transformation, the transformation of transversal coupling induces non-local electron interactions of the form

[TABLE]

due to Jordan-Wigner strings. Additionally to non-locality, the mapping of transversal coupling of qubits to the bath with a Jordan-Wigner transformation produces a Hamiltonian that is odd in fermionic creation and annihilation operators. The linear term $\sim(c_{i}^{\phantom{\dagger}}+c_{i}^{\dagger})$ in the coupling stems from the violation of fermion number conservation due to the exchange of excitations between qubit and bath. This linearity induces another fundamental problem: One cannot replace the Dyson time-ordering that appears in the time evolution operator with a Wick time-ordering necessary for Wick’s theorem. This problem can be avoided if the bath operator $\hat{X}_{i}^{x}$ is odd in fermionic operators. In this case, the Hamiltonian becomes even in creation/ annihilation operators and we can use the standard perturbation theory with Wick time-ordering. In this paper, we discuss a system where this is the case. With dephasing, decay, and free bath Hamiltonian, the effective fermionic Hamiltonian that describes a noisy quantum simulator with fermionic operators takes the form

[TABLE]

We use this Hamiltonian to describe a noisy quantum simulator with the Keldysh contour technique.

III Decoherence on Keldysh contour

In this section, we calculate the non-equilibrium electronic Green’s functions of the noisy quantum simulator. The general idea is to use a non-equilibrium field theory to expand the Green’s functions in the coupling between system (quantum simulator) and bath to obtain the Dyson equation shown in Eq. (1). To check the validity of the mapping from a noisy quantum simulator to a many-body fermionic system, Eqs. (7) and (8), we compare the field theoretical results with master-equation calculations for small systems. We choose a non-interacting fermion system,

[TABLE]

in order to focus on effects of decoherence. Due to the assumption given in Eq. (6), we can treat dephasing and decay separately. In the first part of the section, we briefly show how to calculate the Green’s functions in a typical system-bath approach with the help of a master equation. On a classical computer, this method is feasible for small systems with a size of up to $~{}20-30$ qubits. In the second part, we use the mapping of the noise to fermions and obtain the Green’s functions in the presence of longitudinal coupling to a bath of harmonic oscillators with Ohmic spectral density. In the last part of the section, we show how to calculate decay at zero temperature due to a bath of two level systems, one of the major sources of decoherence for superconducting qubits.Martinis et al. (2005); Müller et al. (2015); Barends et al. (2014)

III.1 Master equation

For small systems the time evolution of the noisy quantum simulator with Hamiltonian $\mathcal{H}_{q}=H_{q}+H_{q,x}+H_{q,z}+H_{B}$ can be simulated approximately with quantum master equations. In the Lindblad form, the density matrix of the quantum simulator obeys the equation of motion

[TABLE]

with the dephasing rate $\Gamma_{2^{*},n}=S_{z,n}(0)$, and the decay rate $\Gamma_{1,n}=S_{x,n}(\epsilon_{n})$. The bath spectral density $S_{n}(\omega)$ will be defined below, in Eq. (38). We can formally integrate the master equation (11) and obtain $\rho(t)=e^{\mathcal{L}(t-t^{\prime})}\rho(t^{\prime})$, where $\mathcal{L}$ is the Lindblad super-operator. From here we use the quantum regression theoremCarmichael (2002) to obtain two-time correlation functions. For $\tau\geq 0$ the theorem states

[TABLE]

We implement Eq. (11) with the QuTip python packageJohansson et al. (2013) to calculate the disturbed fermionic Green’s functions, for example,

[TABLE]

Due to the presence of noise we expect differences to the ideal result. The method can be applied straightforwardly to time dependent problems.

To compare to the calculations obtained in the steady-state approach of the fermionic theory (see below for details), we assume that in the distant past, at $t=0$, our system was decoupled from the bath and in thermal equilibrium $\rho(0)=Z^{-1}e^{-\beta H_{q}}$, where $\beta$ is the inverse temperature. Subsequently, we allow the system to evolve under the Hamiltonian $\mathcal{H}_{q}$ for a time $t\gg\Gamma_{0}^{-1}$, with $\Gamma_{0}$ the smallest decoherence rate. During this evolution correlations between bath and qubit system build up and the entire system reaches a steady state at time $t$ with density matrix $\rho_{0}\equiv\rho(t)$. Because the system has reached a stationary state at time $t$, the Green’s function depends only on the time difference $\tau$, while the exact starting time $t\gg\Gamma_{2^{*}}^{-1}$ is insignificant. Thus, we can calculate the Fourier transform with respect to $\tau$ and compare the results with the Green’s function obtained with the fermionic field theory.

III.2 Fermionic perturbation theory

The results of a noisy quantum simulation are encoded in the perturbed non-equilibrium Green’s function (matrix in the Keldysh space)

[TABLE]

where $G^{\pm}$ are the retarded (+) and advanced (-), and $G^{K}$ the Keldysh Green’s function, defined as

[TABLE]

Due to the mapping between quantum simulator and fermionic system the fermionic Green’s function corresponds to certain qubit correlators that can be measured in an experiment. For example,

[TABLE]

To obtain Eq. (1), which connects the set of perturbed Green’s functions, Eq. (15), to the ideal ones, we follow Refs. [Rammer and Smith, 1986,Stefanucci and van Leeuwen, 2013] and first introduce the contour-ordered Green’s function

[TABLE]

which can be reduced to the real-time Green’s functions Eq. (15). Here, $T_{C}$ is the time-ordering operator along the contour and $\rho(z_{0})$ is the density matrix of quantum simulator and environment at the initial contour time $z_{0}$. In what follows we require that the density matrix at $z_{0}$ factorizes into bath and system $\rho(z_{0})=\rho_{S}(z_{0})\otimes\rho_{B}(z_{0})$ and the system is non-interacting, e.g.,

[TABLE]

while the bath is in thermal equilibrium. This condition ensures the validity of Wick’s theorem, which we will use in the next step. The exact choice of the contour depends on initial conditions of the problem. The most general contour is depicted in Fig. 1(a). We comment on the different choices for the contour in the following subsections. To obtain the second equality in Eq. (19) we changed to the interaction picture with respect to the coupling between simulator (system) and environment (bath) where operators in the interaction picture are defined as

[TABLE]

We further introduced the time evolution along the contour in the interaction picture,

[TABLE]

where the time integration runs along the contour $C$. We expand $U_{C}$ in powers of the interaction between bath and system, apply Wick’s theorem to both bath and system operators. We obtain a diagrammatic expansion of the full contour-ordered Green’s function in terms of the ideal contour-ordered Green’s function $G_{C.0}$ and bath correlation functions $D_{C}$. Defining the self-energy $\Sigma_{C}$ as the sum of all irreducible diagrams we find the contour Dyson equation

[TABLE]

An analytical continuation from contour to real times yields Dyson equations for the non-equilibrium Green’s functions of the noisy quantum simulator introduced in Eq. (15).

III.2.1 Interacting initial state

The full contour depicted in Fig. 1(a) covers the most general situation: At the start of the simulation, $z=t_{0}$, we prepare the system in an initial density matrix $\rho(t_{0})$ of the simulator. Because the initial density matrix on the contour has to be non-interacting, we have to include the time evolution along the vertical path from $z_{0}=t_{0}-i\beta$ with some initialization Hamiltonian $H_{M}(z)$ that evolves the non-interacting density matrix $\rho(z_{0})$ into the interacting matrix at $t_{0}$. Here, $H_{M}$ is a Hamiltonian defined in such a way that it evolves the system to the correct initial state along the vertical branch.Stefanucci and van Leeuwen (2013) We will not use this contour because we find that we can treat many experimentally relevant situations with simpler contours.

III.2.2 Non-interacting initial state

We assume that we prepare the simulator in an initial state that can be cast into the non-interacting form

[TABLE]

With this kind of initial state we can cover a lot of interesting situations.Stefanucci and van Leeuwen (2013) For example with $a_{ij}=\delta_{ij}(1-2\delta_{i1})$ and $\beta\to\infty$ we can implement a situation where qubit $i=1$ is in the up state while all remaining qubits are in the down state. Without initial correlations and non-interacting initial density matrices of the form introduced above, we can omit the vertical part of the contour (a) and evolve only along the real time forward and backward branch of the contour. For this case, we can cast the contour Dyson equation (23) into differential equations for the real-time Green’s functions (15):

[TABLE]

where we introduced the convolution $[A\circ B](t,t^{\prime})=\int_{t_{0}}^{T}d\bar{t}\,A(t,\bar{t})B(\bar{t},t^{\prime})$. Due to the structure of the contour, it suffices to set the upper limit of the convolution integrals to $T=\rm max(t,t^{\prime})$. This set of differential equations has to be solved with boundary conditions

[TABLE]

where $f(\hat{H}_{\rm ini})=[1+\exp(\beta\hat{H}_{\rm ini})]^{-1}$. Together with the identities

[TABLE]

we find a closed system of differential equations which can be solved either numerically or analytically. These equations are known as the Kadanoff-Baym equations (KBE). Except for special situations, the Green’s functions will depend on the initial time $t_{0}$ and on both time arguments $t$ and $t^{\prime}$ separately. This approach corresponds to time-dependent simulations on an analog quantum simulator.

III.2.3 Adiabatic switching

If we are not interested in time-dependent simulations but only on equilibrium properties such as the spectral function or decoherence rates, we can further simplify the contour. Since we are not interested in the evolution from some initial time $t_{0}$ we let $t_{0}\to-\infty$ and evolve the system on the contour 1(b).Kamenev (2011) In the distant past, system and bath were decoupled and in thermal equilibrium at temperature $T$. We switch the coupling between the bath and the system adiabatically from zero to the full value, i.e. $(H_{f,x}+H_{f,z})\to(H_{f,x}+H_{f,z})e^{-\eta|t|}$.

For a time-independent Hamiltonian $H_{f}$ this approach describes a steady state of the system. We assume that the Green’s function depends on the time difference only $G(t,t^{\prime})=G(t-t^{\prime})$, while components that depend on the center of mass time $\sim t+t^{\prime}\gg t-t^{\prime}$ have decayed. This corresponds to neglecting initial correlations of the system. In this case we can introduce the Fourier transformed function $G(\omega)=\int dte^{-i\omega t}G(t)$. The Dyson equation for the Fourier transformed Green’s functions simplifies to Eq. (1):Kamenev (2011)

[TABLE]

From retarded and advanced Green’s function follows the spectral weight

[TABLE]

a quantity we would like to extract in an experimental simulation. We use relation (29) between retarded and advanced self-energy and express $\Sigma^{\pm}=\hat{\Delta}(\omega)\pm i\Gamma(\omega)/2$ through the rate function $\hat{\Gamma}(\omega)=2\operatorname{Im}[\Sigma^{+}(\omega)]$ and the energy shift $\hat{\Delta}(\omega)=\operatorname{Re}[\Sigma^{+}(\omega)]$. The rate corresponds to the lifetime of excitations and determines the width of peaks in the spectral function. It compares to decoherence rates used in master equation calculations. We note that for a time independent Hamiltonian $\hat{H}_{f}$ the retarded and advanced Green’s function always depend only on the time difference $G^{\pm}(t,t^{\prime})=G^{\pm}(t-t^{\prime})$. This statement is independent from initial state preparation. Thus, for any time independent Hamiltonian we can define the Fourier transformed retarded and advanced functions. The Keldysh component $G^{K}$ carries information about the occupation of states and is connected to the distribution function $F$ of the system as Kamenev (2011)

[TABLE]

Thus, from $G^{K}(t,t)$ follows the distribution $F(t)$ at a given time $t$. We conclude this section by introducing the ideal Green’s functions of the fermionic system

[TABLE]

where $F(\omega)=1-2f(\omega)$ is the fermion distribution function, $f(\omega)=[1+\exp(\beta\omega)]^{-1}$ is the Fermi function and $\hat{H}_{f}$ is the matrix representation of the free Hamiltonian $H_{f}$ in the local qubit basis.

IV Adiabatic approach and steady-state analysis

In this section, we analyze systems in the adiabatic limit. We focus on decoherence and thus consider a non-interacting fermionic system. We calculate dephasing rates due to bosonic environment, analyze relaxation due to a bath of two-level systems, and demonstrate the influence of single qubit decoherence on the spectral density of large simulated systems. For small systems we compare the results obtained with our method with master-equation calculations.

IV.1 Dephasing due to a bosonic bath

Here, we calculate the Green’s functions for a quantum simulator subject to pure dephasing due to a bath of harmonic oscillators, where the qubits couple linearly to the displacement of the oscillators

[TABLE]

This model describes, e.g., coupling of a qubit to a resistive environment.Leggett et al. (1987) The bath of oscillators is characterized by its power spectral densityIthier et al. (2005)

[TABLE]

where $J_{i}(\omega)$ is the spectral function and $\langle\cdots\rangle_{0}$ is an average with respect to $\rho_{B}(t_{0}\to-\infty)$.

The interaction term in Eq. (7) with the noise operators in Eq. (37) is identical to a local interaction between electrons and phonons. For this type of system the lowest-order self-energy in Born approximation for the interaction with the bath [see Fig. 1(b)] is given byRammer and Smith (1986)

[TABLE]

Because bath operators of different qubits are uncorrelated, the self-energy is diagonal in the local qubit basis $\Sigma_{ij}\sim\delta_{ij}\Sigma_{ii}$. With the bath correlation functions

[TABLE]

we find the diagonal components of the self-energy

[TABLE]

From the self-energy we extract the lifetime

[TABLE]

with the spectral weight $\mathcal{A}_{0}(\omega)$ of the unperturbed simulator. In the following we calculate the self-energy and rate function for a single qubit. For a flat spectral density we expect to find the golden rule dephasing rate $\Gamma_{2}\propto S(0)$. After this we analyze a system of coupled qubits and compare the results with master-equation calculations.

IV.1.1 Single Qubit Dephasing

The dephasing rates due to noise with a flat spectral density can be calculated with Fermi’s golden rule. Within the same approximations the rates calculated with the non-equilibrium Green’s functions must certainly yield the same dephasing rates. For a single qubit with energy splitting $\epsilon$ the free Green’s functions are $G_{0}^{\pm}=(\omega-\epsilon\pm i0)^{-1}$ and $G^{K}_{0}=-2\pi iF(\omega)\delta(\epsilon-\omega)$. The rate function takes then the form

[TABLE]

The Green’s function for a single qubit is strongly peaked for energies close to the qubit energy splitting. As long as the spectral density $S(\omega)$ is not strongly peaked around $\omega=0$ and we have no static component, i.e. $J(0)\to 0$, we can approximate the self-energy with its value at $\omega\approx\epsilon$. This yields the single-qubit dephasing rate

[TABLE]

This result corresponds to the golden-rule single-qubit dephasing rate due to a bath characterized by a flat spectral density.

IV.1.2 Coupled System

Now, we calculate the self-energy for coupled but non-interacting electrons with $H_{f}=\sum_{ij}t_{ij}c_{i}^{\dagger}c_{j}^{\phantom{\dagger}}$. We can diagonalize the Hamiltonian with a transformation

[TABLE]

where $U$ is a unitary transformation matrix and the Hamiltonian is diagonal in $a_{k}$, $H_{f}=\sum_{k}\epsilon_{k}a_{k}^{\dagger}a_{k}{\phantom{\dagger}}$. We express the Green’s functions of the fermions in the original basis with the Green’s functions $\tilde{G}_{kk^{\prime}}$ of the diagonalized system as

[TABLE]

For the on-site ideal Green’s functions this relation simplifies to

[TABLE]

where the Green’s functions in the diagonal basis take the same form as the single-qubit Green’s functions with $\epsilon\to\epsilon_{k}$. We find the rate function

[TABLE]

For a flat spectral density, for which $J(\omega)\coth\frac{\omega}{2T}\approx\mathrm{const.}$, the rate is frequency independent and we expect the Lindblad master equation to yield identical results as the Keldysh calculations. If the spectral density varies on the frequency range determined by the rate function, the rate function will be frequency dependent. Consequently, the constant-rate master equation does not yield the same results anymore. In this situation, we use a Bloch-Redfield master equation which evaluates the spectral density at the systems energies instead of the simple Lindblad equation.

To compare the results of the fermionic theory with master-equation calculations, we choose a simple system: a linear chain with on-site energies $\epsilon$ and nearest neighbor hopping $t_{i,i+1}=g$,

[TABLE]

with eigenenergies $\epsilon_{k}=\epsilon+g\cos k$, $k=2\pi n/N$. For an ohmic bath, the condition of a flat spectral-density is fulfilled for high temperatures $T\gg\epsilon+g$. In Fig. 2, we compare master-equation calculations with the fermionic Green’s functions obtained from Eq. (1). In case of the high temperature system we use a Lindblad master equation with the single qubit dephasing rates $\Gamma_{2}$ from Eq. (47). The fermionic self-energy is given by Eq. (51). We compare results for two different bath temperatures $T\gg\epsilon+g$ and $T\ll\epsilon+g$. As expected, the Lindblad master equation and Keldysh theory are identical for high temperatures while the description with single qubit dephasing rates fails for a low temperature bath, i.e. frequency dependent self-energy. Thus, for the low temperature situation (upper plot) we used QuTip’s Bloch-Redfield implementation to obtain the master equation results. Deviations between Bloch-Redfield and Keldysh results are due to the fact, that the spectral density vanishes very fast for $\omega\approx 0$ and the Bloch-Redfield solutions become unreliable. We conclude, that the mapping of dephasing noise to fermionic operators works as expected and compares well to usual master-equation calculations.

IV.1.3 Effect of dephasing on large systems

Another interesting question is the spectral resolution of a noisy quantum simulator. We expect intuitively that the resolution is determined by the decoherence rates $\Gamma_{2}$ of the single qubits. For larger systems the eigenenergies become more and more dense and we expect an increasing loss of information due to decoherence. To check the validity of these assumptions we calculate the retarded component of the perturbed Green’s function for different dephasing rates $\Gamma_{2}$ and increasing system size. In Fig. 3, we plot the spectral weight of a tight-binding chain for increasing single-qubit dephasing rates in the top graph. As expected, the perturbed quantum simulator can only resolve features further apart than the dephasing rates and the resolution decreases with increasing dephasing rate. In the bottom graph, we show the influence of system size on the loss of information by comparing the perturbed spectral function (dots) with the ideal one (solid) for different system sizes. Here, we also confirm the expectation that while the disturbed quantum simulator can clearly resolve most of the peaks for the smaller system ($N=20$), more peaks vanish for a larger system ($N=40$).

IV.2 Relaxation

Here, we show how the mapping to fermions works when describing decay for a specific type of noise. Namely, we use the Jordan-Wigner transformation to calculate relaxation due to a bath of two level systems at low temperature. Although academic at first glance, this describes a relevant problem. Experimentally and theoretically it is well understood, that two-level systems (TLS) are responsible for a large part of decoherence in superconducting devices.Martinis et al. (2005) For example, the omnipresent $1/f$ noise can be attributed to a bath of TLS with homogeneously distributed properties.Schriefl et al. (2006) TLS are low energy excitations such as tunneling atoms, dangling bonds or trapped electrons inside amorphous and dielectric layers of superconducting devices.Cole et al. (2010); Lisenfeld et al. (2016) The TLS couple via their dipole moment to the electrical field inside the qubit junction and the coupling takes the form

[TABLE]

where $\sigma^{\pm}$ are qubit operators and $\tau^{\pm}$ are pauli matrices in TLS space. Mapping only qubit operators to fermions, we explicitely violate fermion number conservation, see Eq. (8). For this particular model we can circumvent this problem by mapping both sets of spin-$\frac{1}{2}$ operators, i.e. qubits and TLS, to fermions according to Eq. (8). As we will show below, this approach ensures particle conservation and we can use standard non-equilibrium perturbation theory. In order to use a Jordan-Wigner transformation on both types of Pauli operators, we have to order the TLS in some way together with the qubits. We choose the following convention: We label the first qubit as fermion number one, $\sigma_{1}^{-}\to c_{1}$. Then follow all $N_{1}$ TLS that couple to the first qubit according to $\tau^{-}_{1,s}\to a_{1,s}$. The next fermion is the second qubit $c_{2}$ at overall fermion position $1+N_{1}+1$ and so on. With this numbering, the JWT of the different parts of the Hamiltonian reads

[TABLE]

where $N_{i}$ is the number of TLS coupling to qubit $i$. The JWT introduces many-particle interactions between TLS fermions and qubit fermions. Both qubit-qubit hopping as well as interactions between qubits and TLS are affected by this interaction and the entire problem becomes intractable. However, TLS responsible for decay have energies $\epsilon_{\rm TLS}$ close to typical qubit energies and are thus at low temperatures in typical experimental situations. For TLS at low temperatures $T_{B}\ll\epsilon_{\rm TLS}$ we can simplify the transformed Hamiltonian because the expectation values for the TLS occupation number $\hat{n}_{is}$ is negligible small and we can safely let $\hat{n}_{is}\to 0$ in the Hamiltonian. With this we recover the original qubit Hamiltonian with an additional effective interaction of the form

[TABLE]

This Hamiltonian describes tunneling of particles between two fermionic leads, one comprised of the TLS fermions $a_{is}$ and one of qubit fermions $c_{i}$ with tunneling amplitude $g_{is}$. To obtain the last equality we introduced vectors ${\bf c}=(c_{1},\dots,c_{N})^{\rm T}$ and a tunneling matrix $\hat{T}$.

Because the Hamiltonian is quadratic in fermionic operators we can exactly calculate the Green’s functions. Extending the tunneling matrix to Keldysh space we find that the Keldysh component of the tunneling matrix $\hat{T}$ vanishes, i.e. $T^{K}=0$. Expanding the qubits Green’s function in the tunneling matrix we find

[TABLE]

with self-energy $\Sigma=\hat{T}G_{B}\hat{T}^{\dagger}$, where $G_{B}$ is the free TLS Green’s function. Here, $G$, $G_{B}$, $\hat{T}$, and $\Sigma$ are matrices in Keldysh space. This type of self-energy is known as an embedding self-energy in the context of small systems such as quantum dots, coupled to large fermionic systems, e.g., electronic leads.Stefanucci and van Leeuwen (2013) Using that the tunneling matrix is diagonal in qubit space we find the Dyson equation

[TABLE]

where $a=\pm,K$ denotes the components in Keldysh space. Defining the TLS spectral density $J_{i}(\omega)=2\pi\sum_{s}|g_{is}|^{2}\delta(\omega-\omega_{is})$ we find the retarded and advanced Green’s functions

[TABLE]

In the last step, we have neglected the real part of the self-energy. We identify the decay rate $\hat{\Gamma}_{i}(\omega)=J_{i}(\omega)$. To calculate the Keldysh component of the Green’s function we employ the fluctuation-dissipation relation

[TABLE]

with the TLS distribution function $F_{\rm tls}(\omega)=\tanh(\beta_{B}\omega/2)\approx 1$. We calculate the Keldysh component of the Green’s function with the relation

[TABLE]

For uncoupled qubits, e.g. $t_{ij}=\epsilon_{i}\delta_{ij}$ we find

[TABLE]

with $\Gamma_{i}=J_{i}(\epsilon_{n})$. The rates $\Gamma_{i}$ correspond to Fermi’s golden-rule decay-rates for single qubits coupled to a bath of TLS. Looking at Eq. (66), we see that the qubit system equilibrates to the bath temperature $T\approx 0$ and obtains a finite lifetime due to interactions with the bath.

V Relaxation - Transient behavior

In the previous section, we analyzed stationary-state properties of a quantum simulator coupled to the environment. The stationary state analysis yields information such as decoherence rates of the system. However, many experiments on quantum simulators are far away from stationary states. In a typical experiment, we prepare the quantum simulator in a certain initial state that is not necessarily a stationary state of the quantum-simulator Hamiltonian. For a general initial state $\rho_{0}=|\psi\rangle\langle\psi|$ the Green’s functions depend on the initial state configuration and on both time arguments independently. We use the full contour without the vertical branch according to the Kadanoff-Baym equations (KBE), Eqs. (25-30).

For the model discussed in section IV.2, a quantum simulator coupled to a bath of non-interacting TLS at zero temperature, we have already derived an expression for the self-energy, Eq. (61),

[TABLE]

Since the TLS are in thermal equilibrium at, $T_{B}=0$, we can use the fluctuation-dissipation relation to express the Keldysh component of the self-energy as

[TABLE]

Numerical integration of the KBE with this self-energy is straightforward. Assuming that the relevant energies of the quantum simulator lie well within the excitation spectrum of the TLS, we can use the wide band approximation for the TLS spectral density, $J_{i}(\omega)=\sum_{s}|g_{is}|^{2}\delta(\omega-\epsilon_{is})=\Gamma_{i}=\mathrm{const.}$. In this approximation the bath of TLS corresponds to a white noise with time-local self-energies

[TABLE]

where the rate $\Gamma_{i}$ coincides with the decay rate $\Gamma_{1,i}$ used in a master-equation approach. This approximation for the self-energy, known as the Markov approximation, is valid as long as internal TLS decoherence times are fast compared to typical system time scales. Thus, we can easily compare the results obtained via KBE with master equation calculations. For a time-independent Hamiltonian $\hat{H}_{f}$ retarded and advanced components depend only on the time difference and we can define the spectral weight $\mathcal{A}(\omega)=i(G^{+}(\omega)-G^{-}(\omega))$. The retarded and advanced Green’s function are identical to the equilibrium functions (65). For time-independent Hamiltonian we find the equal time kinetic function

[TABLE]

which is directly related to the (equal time) time dependent distribution function $F(t,t)$ of the simulator. Here, $[\hat{\Gamma}]_{ij}=\delta_{ij}\Gamma_{i}$ is the matrix of decay rates. This result is quite intuitive. For $t=0$ the qubit system is described by the initial distribution $F(\hat{H}_{\rm ini})$ specified by the preparation of the simulator. Subsequently, the simulator evolves with $\hat{H}_{f}$ under the influence of the environment. The environment exchanges excitations with the simulator leading to a decay of the initial distribution to a stationary value. For long times $t\gg\text{min}(\hat{\Gamma})^{-1}$ the simulator reaches a steady state with distribution $F_{\rm tls}$. At this point the quantum simulator reaches equilibrium with the TLS environment at the temperature $T$ provided by the environment. We note that for the steady state we can define $G^{K}(\omega)$ for which the fluctuation-dissipation theorem (66) holds. For $t\geq t^{\prime}$ and $[\hat{H}_{f},\hat{\Gamma}]_{-}=0$ we find the analytical solution

[TABLE]

If both time arguments are large compared to the smallest decay rate, $t,t^{\prime}\gg\Gamma_{\rm min}^{-1}$, only the first term in the square brackets survives and the Green’s function does no longer depend on the initial state. We find the stationary state solution $G^{K}=-i\exp\{(-i\hat{h}-\frac{\hat{\Gamma}}{2})(t-t^{\prime})\}$ which depends only on the time difference and corresponds to a stationary solution with distribution function $1-2f=1$. This behavior resembles the expected behavior of a system in an arbitrary initial state subject to decay into an environmental bath. During a time $\sim\Gamma^{-1}$ the information of the initial state decays into the environment and the Keldysh Green’s function depends only on the difference of time arguments. The decay of the initial state/ distribution function into an effective distribution with $\beta\to\infty$ describes the decay of particles into the bath of TLS at zero temperature.

To demonstrate the validity of the method and compare it to master-equation calculations, we solve the KBE numerically for a linear chain of five qubits coupled to a bath of TLS. Initially, the zeroth qubit was excited, i.e. $n_{0}(0)=\langle c_{0}^{\dagger}(0)c_{0}(0)\rangle=1$, while all remaining qubits were in the ground state. This corresponds to an initial density matrix with Hamiltonian $H_{\rm ini}=-c_{0}^{\dagger}c_{0}+\sum_{i=1}^{N-1}c_{i}^{\dagger}c_{i}$. In Fig. 4 we plot the time dependent occupation numbers $n_{i}(t)=\frac{1}{2}[G^{K}_{ii}(t,t)+i]$, and the total number of particles $n_{\rm tot}(t)$ in the system. We clearly see the expected exponential decay of the total number of particles with a rate dominated by the decay rate of the most populated qubit. As expected for the wide-band limit, master equation and KBE yield the same results (bottom), while for a bath consisting of ten individual TLS with corresponding spectral density, KBE and master equation differ (top).

VI Conclusions

We have applied fermionic non-equilibrium Green’s function methods to a noisy analog quantum simulator. We used the Jordan-Wigner transformation to map the noisy qubit system to a coupled fermion-fermion (relaxation) and fermion-boson (dephasing) system. Dephasing maps to a simple on-site interaction, $\sim c_{i}^{\dagger}c_{i}^{\phantom{\dagger}}$, and is equivalent to an electron-phonon interaction. Using this connection we calculated disturbed Green’s functions of noisy systems with diagrammatic many-body methods. For small systems we compared the results with master-equation calculations. In the correct limit, both methods agree with each other. We found that the spectral resolution of large quantum simulators is limited by the spectral width of the qubits.

Contrary to dephasing, the Jordan-Wigner mapping to fermionic degrees of freedom becomes troublesome for transverse coupling, $\propto\sigma_{x}$, due to the non-conserved fermion number. In this work, we discussed the special case of a bath of TLS at low temperatures, which circumvented these problems.

The mapping to fermionic degrees of freedom comes with the advantage of the advanced and powerful many-body toolbox, ranging from diagrammatic expansions to full field-theoretical approaches. These methods help to establish a connection between the ideal Green’s function that should be the result of quantum simulations and the perturbed Green’s function obtained from a noisy quantum simulator.