T HE U NIVERSITY OF T OKYO
M ASTER T HESISTransport Processes and Quantum Critical Phenomena
in Heat Transport via a Two-State System

Author:Tsuyoshi YAMAMOTO

Supervisor:
Takeo K ATO

A thesis submitted in fulfillment of the requirements
for the degree of Master of Physics
in theKato Group
Department of Physics

March, 2019

iii

THE UNIVERSITY OF TOKYO

Abstract
Graduate School of Science
Department of Physics
Master of Physics
Transport Processes and Quantum Critical Phenomena
in Heat Transport via a Two-State System
by Tsuyoshi YAMAMOTO
Many-body effects and non-equilibrium properties in heat transport via a nanoscale object have attracted much interest for a few decades. A number of theoreticalworks have clarified that there appear various kinds of transport processes despite
the simpleness of the target system. Recent development in experimental techniques
for measurement and sample fabrication has now enabled us to access informationof heat current flowing via a nano-scale circuit. In near future, experimental studies
on heat transport via a two-state system (a qubit) will be possible.
In this thesis, we systematically investigate heat transport via a two-state system,
using the spin-boson system near thermal equilibrium. We derive several analyticalformulas for the thermal conductance within the linear response theory. We also perform an accurate numerical simulation with the quantum Monte Carlo method andconfirm these predictions for arbitrary types of thermal baths. The present numerical calculation classifies transport mechanism and shows that the noninteractingblip approximation quantitatively describes thermal conductance in the incoherent
transport regime.
In addition, we study quantum critical behavior in heat transport via a two-statesystem with sub-Ohmic reservoirs theoretically. We calculate temperature dependence of thermal conductance near quantum phase transition and discuss its criticalexponents. We also propose superconducting circuits to realize the sub-Ohmic spinboson model, which can be used for observation of quantum critical phenomena.

This is a version uploaded to the public web. To avoid the problem of copyrights,
the figures taken from previous studies are omitted.

v

Acknowledgements
The work presented in this thesis could not have been done without the support ofa lot of people. I am deeply grateful to Professor Takeo Kato. His broad knowledge
and great intuition were very helpful for proceeding our research.
I also want to thank the other people in the Kato group. Assistant Professor Rui
Sakano taught me how to make the logical presentation. Tomohiro Tamaya gave me
various advices as a academic researcher. Masahiro Hasegawa and Marc Ken Ritterare excellence and I have stimulated from them at seminars. Hiroko Eguchi often
encouraged me. Without their support, I could have not had my enjoyable graduate
life.
Next, I would like to thank Professor Keiji Saito, who came all the way from Keio
University in Yokohama. Discussions with him have inspired me. I am also gratefulto Masanari Kato. Our work could not proceed smoothly without his work.
Finally, I want to thank my parents. The continued encouragement and financial
support of them have helped me.

vii

List of Publications• T. Yamamoto, M. Kato, T. Kato, and K. Saito, New J. Phys. 20, 093014 (2018).
• T. Yamamoto and T. Kato, Phys. Rev. B 98, 245412 (2018).

ix

Contents
Abstractiii

Acknowledgementsv

1

Introduction
1
1.1 Electronic transport vs. Heat transport . . . . . . . . . . . . . . . . . . . 1
1.2 Spin-boson model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Theories on heat transport via a nano-scale object . . . . . . . . . . . . 4
1.3.1 Analytical formulas . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.2 Kondo effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Experiments on heat transport via a nano-scale object . . . . . . . . . . 7
1.5 Quantum critical phenomena . . . . . . . . . . . . . . . . . . . . . . . . 8
1.6 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.7 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2

Formulation
13
2.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Thermal conductance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3

Transport Processes and Quantum Critical Phenomena
3.1 Effective tunneling amplitude . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Sequential tunneling . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Co-tunneling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Incoherent tunneling . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Noninteracting-blip approximation (NIBA) . . . . . . . . . . . . . . . .
3.6 Quantum critical phenomena for the sub-Ohmic case . . . . . . . . . .

21
21
25
28
29
30
31

4

Continuous-time Quantum Monte Carlo Method
4.1 The kink representation of the spin-boson model . . . . . . . . . . . . .
4.2 Calculation algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Analytical continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Binder analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33
33
34
36
37

5

Numerical Results and Comparison with Analytical Formulas
5.1 Thermal conductance . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1 The Ohmic case (s = 1) . . . . . . . . . . . . . . . . . . . . . . .
5.1.2 The sub-Ohmic case (0 < s < 1) . . . . . . . . . . . . . . . . . .
5.1.3 The super-Ohmic case (1 < s < 2) . . . . . . . . . . . . . . . . .
5.1.4 The super-Ohmic case (2 ≤ s) . . . . . . . . . . . . . . . . . . . .
5.2 Quantum critical phenomena . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 Phase diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 Thermal conductance . . . . . . . . . . . . . . . . . . . . . . . .

39
39
39
41
43
43
44
44
44

x
6

Experimental Realization of Sub-Ohmic Reservoirs
6.1 Superconducting circuit model . . . . . . . . . . . . . . . . . . . . . . .
6.1.1 Transmission line . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1.2 Flux qubit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1.3 Charge qubit . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 The spectral densities . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.1 Flux qubit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.2 Charge qubit . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.3 Charge qubit (0 < s ≤ 0.5) . . . . . . . . . . . . . . . . . . . . .

47
47
47
49
51
53
54
55
56

7

Summary57

A Instanton Method
59
A.1 Exponential factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
A.2 Pre-exponential factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
B The Influence function of the Spin-boson Model

67

C The Noninteracting-blip Approximation (NIBA)

71

D Impedances of the Transmission Line
D.1 The flux qubit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.2 The charge qubit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.3 The charge qubit (0 < s ≤ 0.5) . . . . . . . . . . . . . . . . . . . . . . .

77
77
78
79

1

Chapter 1

Introduction
In this chapter, we summarize the background of this thesis. After we discuss similarities between heat transport and electronic transport in Section 1.1, we introduce
the spin-boson model as a minimal model for describing heat transport via a nanoscale object in Section 1.2. We review existing theoretical and experimental studies on heat transport via a nano-scale object in Section 1.3 and Section 1.4, respectively. We also summarize theoretical studies on quantum phase transition in the
sub-Ohmic spin-boson model in Section 1.5. Based on these previous studies, we
describe the purpose of this thesis in Section 1.6. Finally, we give the organization ofthis thesis in Section 1.7.
In this chapter, we give just a brief summary and detailed descriptions of analytical results are given in the subsequent chapters. Throughout this thesis, we employ
the unit of h̄ = 1.

1.1

Electronic transport vs. Heat transport

In macroscopic scale, it is known that electronic transport is well described by the
Ohm’s law, which states the charge current is proportional to the negative gradientin the electronic potential, where the linear coefficient is the conductivity. The corresponding macroscopic law in heat transport is the Fourier’s law, which states thatthe time rate of the heat current density j through a material is proportional to thenegative gradient in the temperature T:
j = −κ ∇ T,

(1.1)

where κ is the thermal conductivity.
For quantum transport via nano-scale objects such as quantum point contact,
quantum wires, and quantum dots, electronic transport is described by the Landauer
formula [1] based on the scattering theory. This formula predicts quantum transportphenomena, e.g., the quantization of the electronic conductance [2]. It is remarkablethat there exists the corresponding formula for heat transport, that is the Landauerlike formula for photons (or phonons) [3]:
J=

1
2π

Z ∞
0

dω ω T (ω ) [n L (ω ) − n R (ω )] ,(1.2)

where T (ω ) is the transmission probability of photons (or phonons) and nν (ω ) is theBose-Einstein distribution function in the reservoir ν. Then, the quantized thermal
conductance is predicted from Eq. (1.2) as in electronic transport and has indeed
been observed in heat transport through one-dimensional wires [4].

2Chapter 1. Introduction
TABLE 1.1: Electronic transport vs. heat transport

Macroscopic law
Microscopic law
Nano-scale object
Quantum phasetransition.
Electronic Transport
Ohm’s law
Landauer’s formula
Quantum dot
Multi-channel Kondo model

Heat Transport
Fourier’s law (1.1)
Photon (or phonon) version of
Landauer formula (1.2)
Two-state system (a qubit)
Sub-Ohmic spin-boson model

As indicated by these two examples, there exist strong similarities between electronic transport and heat transport (see Table 1.1). It is meaningful to consider counterparts of electronic transport phenomena in heat transport. Electronic transport
via a quantum dot has been attracting considerable attention because the system exhibits various intriguing phenomena, e.g., the Kondo effect [5]. It is quite natural
to ask what happens if one considers heat transport via a zero-dimensional system.
Indeed, as a counterpart of electronic transport via quantum dots, we can consider
heat transport via a two-state system (or a qubit). This system is described by the
spin-boson model (see the next section) and is proved to be related to the Kondo
model [6], which exhibits the Kondo effect. Recently, a number of theoretical studies on heat transport via nano-scale objects including two-state systems have shownthat various transport phenomena can occur in this system (see Section 1.3). Heattransport via a zero-dimensional system is also important in designing thermal devices such as thermal rectification [7] and thermal transistors [8].
Realizing quantum critical phenomena (QCP) near quantum phase transition ina controlled way is one of the central topics in condensed matter physics. Recently,
QCP have been studied for the multi-channel Kondo effect realized in artificial nanostructures [9, 10, 11], and quantum critical behavior has been observed via electronictransport properties in good agreement with theories [12, 13]. It is possible to consider even its counterpart in heat transport. The spin-boson model with sub-Ohmic
reservoirs is known to display quantum phase transition at zero temperature whena system-reservoir is tuned at a critical value (for details, see Section 1.5). Heat transport via a two-state system coupled to the sub-Ohmic reservoirs has, however, notbeen discussed in previous works.
While the similarities between electronic and heat transport via a nano-scale object have been discussed for a few decades, it was generally difficult to study the corresponding heat transport experimentally since high-accuracy calorimetry for nanoscale devices was challenging. Recent development of experimental technology has,
however, changed the experimental situation greatly. For example, measurement ofheat current via a harmonic oscillator (a circuit resonator) and its control by an external bias voltage have been realized quite recently [14] (for details, see Section 1.4).
We expect that experimental studies on heat transport via nano-scale objects will be
accelerated in near future.

1.2

Spin-boson model

We consider the spin-boson model as a minimum model describing heat transport
via nano-scale objects. Despite simpleness of the model, it exhibits fruitful fundamental physics such as the Kondo effect and quantum phase transition. In this section, we introduce the spin-boson model and discuss its properties briefly.

1.2. Spin-boson model

3

F IGURE 1.1: Schematics of the spin-boson model composed of a twostate system and two thermal reservoirs.

A two-state system coupled to thermal reservoirs is called as the ‘spin-boson
system’ (see Fig. 1.1). The Hamiltonian is composed of three parts
H = HS +∑ HB,ν + ∑ HI,ν ,ν= L,R(1.3)

ν= L,R

where HS , HB,ν , and HI,ν describe a two-state system, a bosonic thermal reservoir ν,
and a coupling between the two-state system and the reservoir ν, respectively. Eachterm of the Hamiltonian is given by as follows:∆
HS = − σx − ϵσz ,2†
HB,ν = ∑ ων,k bν,kbν,k ,

(1.4)
(1.5)

k

HI,ν = −



σz
†
λ
b
+
b
ν,k
ν,k
ν,k ,
2 ∑
k

(1.6)

† ) is the annihilation (creation)
where σα (α = x, y, z) is the Pauli matrix and bν,k (bν,koperator of bosonic excitation with the wavenumber k in the reservoir ν. The parameters, ∆, ϵ, ων,k , and λν,k , denote the tunneling amplitude, the detuning energy, theenergy dispersion of bosonic excitation, and the system-reservoir coupling strength,
respectively. In the spin-boson model, the properties of the reservoirs are characterized by the spectral density:Iν (ω ) ≡ ∑ λ2ν,k δ(ω − ων,k ).

(1.7)

k

Assuming that the number of modes is so large to form a continuous spectrum, thespectral density becomes a smooth function of ω, and is usually taken as [15, 6]
Iν (ω ) = αν Ĩ (ω ),  s −1ωĨ (ω ) = 2ωe−ω/ωc ,
ωc

(1.8)
(1.9)

where αν is the dimensionless coupling strength between the two-state system andthe reservoir ν. We introduced the exponential cutoff function e−ω/ωc , where thehigh-frequency cutoff frequency ωc is assumed to be considerably larger than othercharacteristic frequencies, e.g., ∆, ϵ, and kB T. The exponent s of the spectral density is crucial for determining the properties of the reservoirs. The case s = 1 iscalled ‘Ohmic,’ whereas the cases s > 1 and s < 1 are called ‘super-Ohmic’ and‘sub-Ohmic,’ respectively. Note that effect of high-frequency modes (ω ≳ ωc ) isrenormalized into the tunneling amplitude ∆.

4

Chapter 1. Introduction

F IGURE 1.2: Schematics of (a) a superconducting transmission-line
resonator with a flux qubit and (b) a molecular junction.

The spin-boson model has been used in theoretical studies on decoherence phenomena due to the environments [15, 6]. Recently, this model has been used also
for the study of heat transport due to phonons or photons via a nano-scale object.
By coupling the two-state system to the two reservoirs with different temperatures(TL > TR ), the heat current flows from one reservoir to the other one (see Fig. 1.1).
Since the spin-boson model is a simple model for describing heat transport, it has
been used for a variety of systems, e.g., superconducting circuits [16, 17, 18], photonic waveguides with local two-state systems [19], and molecular junctions [20, 21](see Fig. 1.2).
In the following sections, we briefly review existing theoretical and experimental
studies on heat transport via a two-state system using the spin-boson system. The
explicit expression of the heat current in this model will be given in Chapter 2.

1.3

Theories on heat transport via a nano-scale object

1.3.1

Analytical formulas

In a number of theoretical studies, heat transport via a two-state system has been investigated by analytical methods [7, 22, 23, 24, 25, 26, 27]. In the pioneering work by
Ruokola and Ojanen [23], heat current via a two-state system was calculated using
perturbation calculation with respect to the coupling Hamiltonian HI . This calculation gives transition rates between the eigenstates of the two-state Hamiltonian, |±⟩.
Solving the stochastic equation, we obtain the heat current as
J = J+ + J− ,
4

J± =

×

π sin θ8

(1.10)Z ∞0dω ωIL (ω ) IR (ω ) [n L (ω ) − n R (ω )]1ω − ω0 ± 2i Γ∓−

1ω + ω0 ∓ 2i Γ∓,(1.11)√
where tan θ = ∆/ϵ, ω0 = ϵ2 + ∆2 , nν (ω ) is the Bose-Einstein distribution function
in reservoir ν, and Γ± is a transition rate out of state |±⟩ calculated by the Fermi’s
golden rule.
At the high-temperatures limit (kB T ∼ ∆), since resonant energy transfer at ω =ω0 is dominant, the heat current is obtained asJ=

π sin2 θ ω0 IL (ω0 ) IR (ω0 ) [n L (ω0 ) − n R (ω0 )].
2 IL (ω0 ) [2n L (ω0 ) + 1] + IL (ω0 ) [2n L (ω0 ) + 1]

(1.12)

1.3. Theories on heat transport via a nano-scale object

5

−4

10

−6

10

−8

10

−10

10

−2

10

−1

10

0

10

F IGURE 1.3: Schematics of (a) the sequential tunneling and (b) cotunneling processes. (c) The left-reservoir temperature dependenceof the total heat current for the symmetric Ohmic case (s = 1, ϵ =0) with α L = α R = 0.1/π [23]. Right-reservoir temperature is
TR = 0.95TL . The solid line represents the numerical solution of
Eq. (1.11), the green dash-dotted line represents the high-temperatureapproximation (1.12), and the blue dashed line represents the lowtemperature approximation (1.13).

This formula indicates a heat transport process induced by boson absorption andemission (see Fig. 1.3 (a)). This transport process is called the ‘sequential tunneling’.
In contrast, at the low-temperature limit (kB T ≪ ∆), the frequency dependence ofthe energy denominators in Eq. (1.11) can be ignored and the excitation of the twostate system is exponentially suppressed. Then, the expression of the heat current
has been obtained as
π sin4 θ
J=
2ω02

Z ∞
0

dω ωIL (ω ) IR (ω ) [n L (ω ) − n R (ω )] .

(1.13)

This expression indicates heat transport via a virtual excitation in the two-state system (see Fig. 1.3 (b)). This transport process is called the ‘co-tunneling’.
Roukola and Ojanen have also shown that the crossover between the sequential
tunneling and co-tunneling is quite sharp as shown in Fig. 1.3 (c). This means that
one of the two analytical formulas holds well except for a narrow transition region.
We should note that the heat current (1.13) in the co-tunneling regime obtained by
the perturbative method does not include higher-order processes, which is important in the Kondo effect.
There are many other theoretical approaches to treat heat transport via a twostate system. For weak system-reservoir coupling, one can consider perturbative
calculation for general quantum systems using the quantum master equation [7]and the nonequilibrium Green’s function method [25]. For strong system-reservoir
coupling, the noninteracting-blip approximation (NIBA) [22, 24] and the polarontransformed Redfield equation (PTRE) [26, 27] have been applied. Although theseanalytical methods are powerful in analyzing heat transport, the condition for which

6

Chapter 1. Introduction

these approximations hold well is non-trivial.

1.3.2

Kondo effect

In the spin-boson model, the many-body effect emerges due to the strong interaction
between the local system and the environments. In particular, it is known that thesymmetric Ohmic spin-boson model (s = 1, ϵ = 0) can be related to the anisotropic
Kondo model [5] using the bosonization technique [6, 15, 28]. The Kondo model describes a single magnetic impurity with spin 1/2 interacting with conduction electrons. The ground state forms a singlet state between the spin of the impurity andspins of conduction electrons, which is called the ‘Kondo singlet’. The temperature
at which the Kondo singlet is formed is characterized by the Kondo temperature TK .
The anisotropic Kondo Hamiltonian takes the form
HK = ∑ ϵk c†k,s ck,s +
k,s



J∥
J
σz ∑ sc†k,s ck′ ,s + ⊥ ∑ σ+ c†k,↓ ck′ ,↑ + σ− c†k,↑ ck′ ,↓ ,4L k,k′ ,s
2L k,k′(1.14)

where ck,s (c†k,s ) is the annihilation (a creation) operator of the conduction electronwith wavenumber k and spin polarization s =↑ or ↓, L is the system size, andσ± = (σx ± iσy )/2. The J∥ -term describes the scattering of the impurity in which
the spin polarization is conserved while the J⊥ -term describes spin-flip scattering.
The correspondence of the parameters between the Ohmic spin-boson model andthe anisotropic Kondo model is [6, 15, 28]



πρJ∥
2
∆/ωc = ρJ⊥ cos arctan,(1.15)4πρJ∥ 22,(1.16)
α = 1 − arctan
π4where ρ is the density of states of conduction electrons. In the Ohmic spin-bosonmodel, the characteristic energy scale (the Kondo temperature) is given by [6, 15, 29]  α/(1−α)∆ ∆, (0 < α < 1),(1.17)kB TK =ωc
0,(1 ≤ α ).
We show the “Kondo temperature” as a function of the system-reservoir coupling α
in Fig. 1.4. As shown Fig. 1.4, at low temperature T < TK , the Kondo-like many-bodystate is realized in the spin-boson model.
Saito and Kato have studied the Kondo-like effect in heat transport in the Ohmic
spin-boson model using the quantum Monte Carlo simulation [29]. They have clarified the Kondo signature in the universal temperature dependence of the linear
thermal conductance proportional to T 3 at low temperatures (T ≪ TK ) (see Fig. 1.5).
They have also shown that at high temperatures the thermal conductance is proportional to T 2α−1 . While the detailed features of the thermal conductance havebeen clarified for the Ohmic spin-boson model, the non-Ohmic case (the sub-Ohmic
and super-Ohmic cases) has not been studied so much. It is a fundamental question whether the Kondo-like effect emerges in heat transport in the non-Ohmic spinboson model or not.

1.4. Experiments on heat transport via a nano-scale object

7

F IGURE 1.4: The “Kondo temperature” TK as a function of α. This line
indicates a crossover from the Kondo regime (blue) to the incoherent
tunneling regime (orange). At zero temperature, a quantum phasetransition occurs at α = 1.

F IGURE 1.5: The temperature dependence of the thermal conductance for various α [29]. The thermal conductance is proportional toT 2α−1 for high temperatures and T 3 for low temperatures.

1.4

Experiments on heat transport via a nano-scale object

The recent, rapid development in techniques of heat measurement that has enabled
us to experimentally access the heat current in nano-scale objects is remarkable.
In particular, Ronzani et al. have observed heat transport via resonator-systemresonator quite recently (see Fig. 1.6 (a)) [14]. They have considered two models:
the quasi-Hamiltonian (QH) model for γ ≃ g and the non-Hamiltonian (NH) modelfor γ ≪ g, where γ and g are the reservoir-resonator and resonator-system couplings, respectively. For each model, they have measured the total heating powers
for the different source temperatures with high accuracy as shown in Fig. 1.6 (b) and(c).
In addition, it has been theoretically proposed that a superconducting qubit coupled to transmission lines realizes the Ohmic spin-boson model [17, 18], and experimentally demonstrated [31, 30]. The experimental setup composed of a superconducting flux qubit (a two-state system) coupled to transmission lines (reservoirs) isshown in Fig. 1.7 (a). The flux qubit consists of a superconducting loop interrupted
by three (or more) Josephson junctions [32] (for details, see Section 6.1). Magazzú et
F IGURE 1.6: (a) A circuit representing the experimental setup in
Ref. [14]. A transmon qubit, whose tunneling amplitude can be tunable by the external flux Φ, is coupled to two Lr Cr -resonators andeach resonator is terminated at one end by a mesoscopic normalmetal resistor which plays the role as the reservoir. The total heating
powers as a function of the external flux Φ in (b) the QH model and(c) the NH model regimes for different source temperatures TS [14].

8

Chapter 1. Introduction
F IGURE 1.7: (a) Experimental setup in Ref. [30]. A coplanar waveguide running across the chip plays the role as the Ohmic reservoir.
The inset is a scanning electron micrograph showing the flux qubit
and the scale bar is 2 µm. (b)-(d) Experimental and theoretical transmissions for α = 0.007, 0.21, and 0.8, respectively [30].

al. have experimentally investigated an impact of quantum coherence by the systemreservoir coupling strength [30]. When the microwave is injected into the transmission line, the detuning energy is described by
ε(t) = ε 0 + ε p cos(ω p t),(1.18)

where ε p and ω p are the amplitude and frequency of microwave, respectively, andthe static detuning energy ε 0 is given as ε 0 ∝ Φε − Φ0 /2 (Φε : the external flux,
Φ0 = h/2e: flux quantum). When the microwave amplitude is weak, the linear
response theory gives the transmission at the microwave frequency ω p as

T ( ω p ) = 1 − i N ω p χ ( ω p ),(1.19)

where N is a coupling constant determined by fitting and χ(ω ) is the dynamic susceptibility of the qubit. Fig. 1.7 (b)-(d) show the experimental (the left side) andtheoretical (the right side) transmissions as a function of the static detuning energy
ε 0 and the probe frequency ω p for different coupling strengths. The theoretical results have been calculated in the NIBA [15, 6]. For the weakest coupling strength,
α = 0.007, one can see the high quantum coherence as shown in Fig. 1.7 (b). As
the coupling strength is increased, for α = 0.21, the quantum coherence is reduced,
but the hyperbolic curve characteristic of it can be recognized in Fig. 1.7 (c). However, for the strongest coupling strength, α = 0.8, the quantum coherence is broken
completely and the transmission is independent of the prove frequency as shown in
Fig. 1.7 (d).
We expect that heat transport via a two-state system will be measured in near future by combining the accurate heat measurement technique [14] with the controlled
system-reservoir coupling [30].

1.5

Quantum critical phenomena

One of the intriguing phenomena in the spin-boson model is that, at zero temperature, a quantum phase transition occurs when the coupling strength between the
two-state system and thermal reservoirs is tuned to a critical point αc . For α < αc ,the ground state is described by a coherent superposition of two wavefunctions localized at each well, and is called a ‘delocalized phase’. For α > αc , the ground-state
becomes two-fold degenerate because the coherent superposition is completely broken due to the disappearance of quantum tunneling between the two wells. This
phase is called a ‘localized phase’. For the Ohmic case, a Kosterlitz-Thouless-typephase transition occurs at αc = 1 [6, 33, 34]. The sub-Ohmic case shows a secondorder phase transition [35, 36, 37, 38, 39, 40]. The super-Ohmic case does not have a
distinct phase transition.
In this thesis, we focus on quantum critical phenomena in the sub-Ohmic case.
Fig. 1.8 shows schematic of the population ⟨σz ⟩ε as a function of the detuning energy
ϵ near the critical point α = αc . In the delocalized phase (α < αc ), the slope at ϵ = 0

1.5. Quantum critical phenomena

9

F IGURE 1.8: Population, ⟨σz ⟩ε , as a function of the detuning energy,
ϵ near the critical point. For α < αc (blue line), ⟨σz ⟩ε is a continuousfunction of ϵ, and the static susceptibility, χ0 , can be defined by theslope at ϵ = 0. At the critical point (green line; α = αc ), ⟨σz ⟩ε iscontinuous, but the static susceptibility diverges at ϵ = 0. For α > αc
(red line), ⟨σz ⟩ε is discontinuous at ϵ = 0.
TABLE 1.2: Summary of the critical exponents.

Exponentγβ′ηxDefinition
χ0 ∝ ( αc − α ) − γ′
m z ∝ ( α − αc ) βC ( τ ) ∝ τ −η
χ0 ∝ T − x

Condition
α < αc , T = 0
α > αc , T = 0
α = αc , |τ | ≪ 1/kB T
α = αc , T > 0

corresponds to the static susceptibility defined as⟨σz ⟩ε,ϵ →0 ϵχ0 = lim(1.20)

where ⟨· · ·⟩ implies an equilibrium average. The static susceptibility diverges as
the system-reservoir coupling α approaches the critical point αc from below. In thelocalized phase (α > αc ), the population ⟨σz ⟩ε jumps from −mz to +mz at ϵ = 0,
where mz ≡ ⟨σz ⟩ϵ→+0 is the spontaneous magnetization. Information of the system dynamics in thermal equilibrium is included in the imaginary-time correlation
function C (τ ) = ⟨σz (τ )σz (0)⟩. In Table 1.2, we summarize the critical exponents. In
principle, all the exponents can be determined experimentally by measurement of
the population ⟨σz ⟩ε . By using the two exponents related to the fixed point, y∗t andy∗h , the critical exponents in Table 1.2 are expressed as [37]γ=2y∗h − 1
,
y∗t

β′ =

1 − y∗h,y∗t

η = 1 − x = 2 − 2y∗h .

(1.21)For the case of 0 < s ≤ 0.5, the exponents, y∗t and y∗h , are given by the mean-field
theory as [41, 42, 37]
y∗t =

1
,
2

y∗h =

3
,
4

(1.22)

because the transition occurs above the upper critical dimension. Hence, we obtain
γ = 1,

β′ =

1
,
2

η=

1
,
2

x=

1
.
2

(1.23)

10

Chapter 1. Introduction

F IGURE 1.9: The static susceptibility χ0 and magnetization mz close
to the critical point for s = 0.2 and ∆/ωc = 0.1 [37]. (a) χ0 as afunction of the distance from the critical point −δ ≡ −(α − αc )/αc >0 for different temperatures. As the temperature is lowered, the static
susceptibility approaches a line proportional to (−δ)−1 . (b) mz as afunction of δ > 0 for different s. The straight lines are proportional toδ1/2 . (c)-(d) T 1/2 χ and T −1/4 mz collapse onto one point at δ = 0.

Winter et al. have demonstrated using the quantum Monte Carlo simulation and confirmed that these critical exponents for 0 < s ≤ 0.5 agree with the numerical resultsin Fig. 1.9 [37]. For 0.5 < s < 1, the exponents, y∗t and y∗h , are nontrivial functionsof the exponent s of the spectral density. Luijten has calculated the exponents using
the ε-expansion as follows [41]:14A(s) 2y∗t = s + ε −ε + O(ε3 ),
6
9s

y∗h =

1+s 1
A(s) 2+ ε−
ε + O(ε3 ),2
46s(1.24)

where ε = 2s − 1, A(s) = s [ψ(1) − 2ψ(s/2) + ψ(s)], and ψ( x ) is the digammafunction. Vojta et al. have studied quantum critical phenomena using the numerical renormalization group (NRG) for all values of s [35, 36, 38, 40]. The NRG produces the phase diagram consistent with one from the quantum Monte Carlo methods [35, 37]. However, the NRG simulation gives incorrect critical exponents which
do not agree with those from the quantum Monte Carlo method [37] or Eq. (1.22)
for the case of 0 < s ≤ 0.5. In Ref. [40], the origin of this disagreement has been
discussed and attributed to a failure in describing quantum critical phenomena ina naive NRG simulation. Finally, we show a schematic phase diagram of the subOhmic spin-boson model in Fig. 1.10.
So far, experimental realization of the sub-Ohmic spin-boson model has been discussed only in few theoretical works. In Refs. [43, 44], experimental realization of the
sub-Ohmic reservoirs of s = 0.5 was theoretically proposed using transmission lines
with resistance elements. The sub-Ohmic reservoirs with arbitrary s have, however,1.6. Purpose

11

F IGURE 1.10: The phase diagram of the sub-Ohmic case. At zero temperature, a quantum phase transition occurs between the delocalized
and localized phases at the critical point α = αc . In contrast, at finite
temperature, there is a crossover between the co-tunneling and incoherent tunneling regimes.

not been discussed.

1.6

PurposeIn this thesis, we study heat transport using the spin-boson model, which is thesimplest model to describe heat transport via a nano-scale object. Despite its simpleness, the spin-boson model exhibits remarkable many-body phenomena such as
the Kondo effect and the quantum phase transition. Theoretical analysis of the spinboson model not only contributes to the understanding of transport processes andmany-body effect in heat transport but also provides a concrete foundation for designing novel thermal devices.
In this thesis, we focus on heat transport near equilibrium, for which the linear
response theory is applicable. We calculate the thermal conductance, which is thelinear-response coefficient in heat current under a temperature gradient, using bothanalytical and numerical methods. We first derive analytical formulas for three different transport processes, i.e., sequential tunneling, co-tunneling, and incoherent
tunneling. In a previous study [45], Kato has derived a new asymptotically exact
expression including higher-order perturbations, which are not taken into accountin the previous work by Ruokola and Ojanen [23]. This new formula is applicable
also to the non-Ohmic spin-boson model and is useful to understand the Kondolike effect at low-temperatures. We perform the continuous-time quantum Monte
Carlo simulation and compare the numerical results with these analytical formulas. The present Monte Carlo study is a substantial extension of the previous workby Saito and Kato [29] toward arbitrary types of reservoirs; in Ref. [29], only the
Ohmic reservoirs have been considered. We confirm the temperature dependencies predicted by our expressions for the co-tunneling and the sequential-tunnelingtransport regime, and in addition, discuss the accuracy of noninteracting-blip approximation (NIBA) [15] in incoherent tunneling regime. We consider a general
picture to understand transport properties at extremely low temperatures for thewhole regime of spectral densities and also characterize the transport mechanism
for all temperature regimes quantitatively.

12

Chapter 1. Introduction

In this thesis, we also investigate quantum critical phenomena in heat transport
via a two-state coupled to the sub-Ohmic reservoirs. We calculate the temperaturedependence of thermal conductance in the critical regime and discuss how its criticalexponent is determined. We also propose superconducting circuits to realize thesub-Ohmic spin-boson model with an arbitrary value of s. The present study will
provide a useful starting point for observing quantum critical phenomena via heattransport.

1.7

OutlineThis thesis is organized as follows. In Chapter 2, we explain the spin-boson modelin detail and derive the linear thermal conductance using the Keldysh formalism. In
Chapter 3, after introducing an effective tunneling amplitude, we classify the transport mechanism and derive analytical expressions for the thermal conductance ineach process. We also consider quantum critical phenomena in the sub-Ohmic case.
In Chapter 4, we explain the numerical method to calculate the thermal conductancebased on the continuous-time quantum Monte Carlo method. In Chapter 5, we show
numerical results of the thermal conductance and compare the results with analytical approximations derived in Chapter 3. We also determine the critical exponentof the thermal conductance in the sub-Ohmic case. In Chapter 6, we propose superconducting circuits with a qubit to realize the sub-Ohmic reservoirs. Finally, wesummarize our study in Chapter 7.

13

Chapter 2

Formulation
In this chapter, we first introduce a general Hamiltonian coupled to reservoirs andconsider the reduction of it to the spin-boson Hamiltonian. Second, we derive theMeir-Wingreen-Landauer-type formula for the thermal conductance using the standard Keldysh formalism.

2.1

Model

We consider heat transport via a local quantum system coupled to two reservoirs
denoted by L and R. The model Hamiltonian is given by
H = HS +∑ HB,ν + ∑ HI,ν ,

ν=L,R

(2.1)

ν=L,R

p2
+ V ( q ),
2M

pν,k
1
2
2
= ∑
+ mν,k ων,k xν,k ,
2mν,k
2
k

HS =

(2.2)

HB,ν

(2.3)

HI,ν =∑−Cν,k xν,k q +

k

2
Cν,k
2
2mν,k ωk,ν

!
q

2

,

(2.4)

where HS , HB,ν , and HI,ν describe the local system, the reservoir ν (=L, R), and the interaction between them, respectively. The operators p and q are the momentum andposition of the local system, respectively, and V (q) is the potential energy. The reservoirs comprise multiple phonon (or photon) modes, which are described in generalby harmonic oscillators with frequency ων,k and mass mν,k , where the subscript denotes the phonon (photon) wavenumber k in the reservoir ν. The momentum andposition of an individual oscillator are denoted by pν,k and xν,k , respectively. For
simplicity, the system-reservoir coupling HI,ν is consider as a bilinear form of q andxν,k , and the interaction strength is denoted by Cν,k . The second term of HI,ν is acounter term to cancel the potential renormalization due to the reservoirs [46, 15].
In this thesis, the potential energy V (q) of the local system is considered as a
double-well potential as shown in Fig. 2.1:
Mω02
V (q) =
2

q2 − 14 q20
q0

!2
,

(2.5)

where ω0 is the frequency of small oscillation in each well. The double-well potential, V (q), has minimal potential energy V (±q0 /2) = 0 at q = ±q0 /2 and thebarrier height Vb ≡ V (q = 0) = Mω02 q20 /32. Quantum tunneling between the two

14

Chapter 2. Formulation

F IGURE 2.1: Symmetric double-well potential of the local system. An
energy spacing of quantum levels in each well is ω0 (indicated by theblue solid lines), and an energy splitting due to quantum tunneling(indicated by the red dashed lines) is ∆ = Ee − Eg , where Eg and Ee
are the ground-state energy and the first-excited-state energy, respectively.

wells induces small energy splitting ∆0 between the ground-state energy Eg and thefirst-excited-state energy Ee .
We consider the parameter regime∆0 , |ϵ|, kB T ≪ ω0 ≪ Vb ,(2.6)

where ω0 is the energy spacing of quantum levels in each well, ϵ is a detuning energy. The tunneling amplitude ∆0 can be calculated using instanton methods asfollows (see Appendix A) [15, 47]:
s


16Vb
2Vb
ω0 exp −,(2.7)∆0 = 8πω0
3ω0
the local system with the double-well potential is effectively restricted to a twodimensional Hilbert space spanned by the states |↑⟩ and |↓⟩ localized in the rightand left wells, respectively, and described by the two-state Hamiltonian in terms of
the Pauli matrices {σi } in the pseudo-spin form (see Fig. 2.2):∆
HS = − σx − ϵσz .
2

(2.8)

Here, for later convenience, we employ the renormalized tunneling amplitude ∆,which includes the effect of the high-frequency oscillators using adiabatic renormalization (for the detail discussion, see Section 3.1), instead of the bare tunneling amplitude ∆0 . In this thesis, we examine the symmetric double-well potential (ϵ = 0).
For the symmetric case (ϵ = 0), the system Hamiltonian HS describes the tunneling
splitting ∆ between the ground state (σx = +1) and the first excited state (σx = −1).
Introducing annihilation and creation operators of the harmonic oscillators with thewavenumber k in the reservoir ν:rmν,k ων,kpν,kbν,k =
xν,k + i,(2.9)
2mν,k ων,krmν,k ων,kpν,k†bν,k =
xν,k − i,(2.10)
2mν,k ων,k2.1. Model

15

F IGURE 2.2: Schematic of the model comprises a two-state system
coupled to two bosonic reservoirs (L and R) with temperatures TL and
TL , respectively. Symmetric double-well potential of the local system.
(a) σx -basis and (b) σz -basis for the local two-state system.

the Hamiltonians of the reservoir and the system-reservoir interaction can be expressed in terms of these operators:
HB,ν =†bν,k ,∑ ων,k bν,k(2.11)kHI,ν = −



σz
†
,
λν,k bν,k + bν,k∑2 k(2.12)

prespectively, where λν,k = q0 Cν,k / 2mν,k ων,k . Here we have dropped the irrelevant
constant zero-point energies. Consequently, we obtained the spin-boson Hamiltonian:H = HS + HB + HI ,∆HS = − σx − ϵσz ,2†
HB = ∑ HB,ν = ∑ ων,k bν,kbν,k ,
ν=L,R(2.13)(2.14)(2.15)

ν,k

σz





†
.
∑ HI,ν = − 2 ∑ λν,k bν,k + bν,k

HI =

ν=L,R

(2.16)

ν,kThe reservoir effects are characterized by the spectral density:
I (ω ) =∑ Iν (ω ) = ∑ ∑ λ2ν,k δ(ω − ων,k ) = α Ĩ (ω ),

ν=L,R


Ĩ (ω ) = 2ω(2.17)

ν=L,R kωωc

 s −1e−ω/ωc ,(2.18)

where α = αL + αR is the dimensionless coupling strength between the two-statesystem and the reservoirs, and ωc is the cutoff frequency. In this thesis, we focus on

16

Chapter 2. Formulation

the symmetric case (ε = 0) and only use the detuning energy to defined the static
susceptibility
χ0 = lim
ε →0

⟨σz ⟩
,
ε

(2.19)

where ⟨· · ·⟩ denotes an equilibrium average.

2.2

Thermal conductanceIn this section, we derive the linear thermal conductance from the spin-boson model.
First, based on previous studies [23, 29, 3], we derive the Meir-Wingreen-Landauer
formula for the heat current in the Keldysh formalism. The heat current flowing
from the reservoir ν into the local two-state system is defined as follows:Jν ≡ −



dHB,ν
σz
†
= i [ HB,ν , H ] = −i ∑ λν,k ων,k −bν,k + bν,k
.
dt
2 k

(2.20)

The average of the heat current is written as⟨ Jν (t)⟩ = Tr [ρJν (t)] = Tr [ρ(t) Jν ] ,(2.21)

Here ρ is the initial-state density matrixρ = ρS ⊗ ∏ ρ ν ,(2.22)

ν=L,R

ρν =

1 − β ν HB,ν,eZν(2.23)

where ρ(t) = e−iH (t−t0 ) ρeiH (t−t0 ) and Zν is the partition function of the isolated reservoir ν. Hereafter, we consider the average of the observables after realizing thesteady state by taking the limit t0 → −∞. Substituting the definition of the heatcurrent (2.20) into Eq. (2.21), the heat current is rewritten as∗ ων,k λν,k<<⟨ Jν (t)⟩ = ∑Gσ ,b† (t1 , t2 ) + Gσ ,b† (t1 , t2 )(2.24)
z ν,k
z ν,k
2
t1 = t2 = tk"<= Re ∑ ων,k λν,k Gσ ,b† (t1 , t2 ),(2.25)kz

ν,kt1 = t2 = t

where Gσ< ,b† (t1 , t2 ) is the lesser Green function defined as
z

ν,k

′
′
G<
A,B ( t, t ) = −i ⟨ B ( t ) A ( t )⟩ ,(2.26)

where A(t) = eiHt Ae−iHt and B(t) = eiHt Be−iHt . In addition, the retarded and advanced Green functions are defined as
GrA,B (t, t′ ) = −iθ (t − t′ ) ⟨[ A(t), B(t′ )]⟩ ,(2.27)

GaA,B (t, t′ )(2.28)′′= iθ (t − t) ⟨[ A(t), B(t )]⟩ ,respectively, where θ (t) is the Heaviside step function.

2.2. Thermal conductance17

F IGURE 2.3: Keldysh contour C. C+ is the forward branch and C− isthe backward branch.

Next, let us calculate the lesser Green function Gσ< ,b† (t1 , t2 ) in the Keldysh forz

ν,kmalism [48]. We introduce the Green function on the Keldysh contour as
G A,B (τ, τ ′ ) = −i ⟨ TC A(τ ) B(τ ′ )⟩ ,(2.29)

where TC is a time-ordered product on the Keldysh contour. The Keldysh contour
consists of two oriented branches C = C+ ∪ C− , the forward branch C+ extending
from −∞ to ∞ and the backward branch C− extending from ∞ to −∞ (see Fig. 2.3).
Using the formal expansion with respect to λν,k , the Green function can be writtenas†Gσz ,b† (τ1 , τ2 ) = −i ⟨ TC SC σ̃z (τ1 )b̃ν,k(τ2 )⟩ .

(2.30)

ν,kHere Õ(τ ) = ei( HS + HB )τ Oe−i( HS + HB )τ and SC is the S-matrix defined as∞(−i )nSC = ∑
n=0 n!

ZC

du1 · · ·ZC

dun H̃I (u1 ) · · · H̃I (un ),(2.31)

R
where ui is a time variable on the Keldysh contour and C dui implies a line integral
along the Keldysh contour. Thus, we obtain the relation
Gσz ,b† (τ1 , τ2 ) =
ν,k

λν,k
2

Z
C

du Gσz ,σz (τ1 , u) gbν,k ,b† (u, τ2 ),(2.32)

ν,kwhere gbν,k ,b† (u, τ2 ) is the Green function for the isolated reservoir ν. By projection
ν,konto the real-time axis using the Langreth rule [48], the lesser Green function can be
written asZλν,k ∞<r<<aGσ b† (t1 , t2 ) =ds Gσz ,σz (t1 , s) gb ,b† (s, t2 ) + Gσz ,σz (t1 , s) gb ,b† (s, t2 ) .
z ν,kν,k ν,kν,k ν,k2 −∞(2.33)

18

Chapter 2. Formulation

Substituting this expression of the lesser Green function (2.33) into Eq. (2.25), theheat current is then rewritten as⟨ Jν (t)⟩ = Re

Z ∞
−∞



×

ds ∑

λ2ν,k
2

k

ων,k

Gσr z ,σz (t1 , s) gb< ,b† (s, t2 ) + Gσ<z ,σz (t1 , s) gba ,b† (s, t2 )
ν,k
ν,k



ν,k

ν,k

t1 = t2 = t

Z ∞


∂  r
<
a
= Re
ds (−i )Gσz ,σz (t1 , s)Σ<ν ( s, t2 ) + Gσz ,σz ( t1 , s ) Σν ( s, t2 ) t1 =t2 =t ,∂t2−∞(2.34)′a′where Σ<ν ( t, t ) and Σν ( t, t ) are the lesser and advanced self-energies′Σ<ν ( t, t ) ≡

λ2ν,k

i

∑ 2 gb< ,b (t, t′ ) = − 2 ∑ λ2ν,k nν (ων,k )e−iω (t−t )
ν,k

k

= −

i
2

Z ∞
0

†
ν,k

′

ν,k

k

dω ∑ λ2ν,k δ(ω − ων,k )nν (ω )e−iω (t−t )
′

k

Z ∞

′
i
dω Iν (ω )nν (ω )e−iω (t−t ) ,2 0
λ2ν,k a′iΣaν (t, t′ ) ≡ ∑gb ,b† (t, t′ ) = θ (t′ − t) ∑ λ2ν,k e−iων,k (t −t)ν,k ν,k
2
2
k
k

= −

=
=

i
θ (t′ − t)
2

Z ∞

i
θ (t′ − t)
2

Z ∞

0

0

(2.35)dω ∑ λ2ν,k δ(ω − ων,k )e−iω (t −t)′k′dω Iν (ω )e−iω (t −t) ,(2.36)

 −1respectively. Here nν (ω ) = e β ν ω − 1is the Bose-Einstein distribution function
in reservoir ν. We have used the definition of the spectral density (2.17). Performinga Fourier transform on the heat current (2.34), we obtain


Z
 r

1 ∞
i <
⟨ Jν ⟩ =(2.37)dω ωIν (ω ) Im Gσz ,σz (ω ) nν (ω ) − Gσz ,σz (ω ) .
2 0
2The conservation law of energy ⟨ JL ⟩ = − ⟨ JR ⟩ ≡ ⟨ J ⟩ gives⟨ J⟩ ==αR
αL⟨ JL ⟩ −⟨ JR ⟩
αL + αR
αL + αR
Z ∞


αL αR
dω ωIm Gσr z ,σz (ω ) Ĩ (ω ) [nL (ω ) − nR (ω )] ,2 ( αL + αR ) 0(2.38)(2.39)

1r
where Ĩ (ω ) = α−ν Iν ( ω ) (see Eq. (2.17)). Rewriting Gσz ,σz ( ω ) with χ ( ω ), we finallyobtain⟨ J⟩ =

αγ
8

Z ∞
0

dω ωIm [χ(ω )] Ĩ (ω ) [nL (ω ) − nR (ω )] ,(2.40)

where α = αL + αR , γ = 4αL αR /α2 is an asymmetric factor, and χ(ω ) is the dynamical susceptibility of the two-state system defined byχ ( ω ) = −i

Z ∞0

dt eiωt ⟨[σz (t), σz (0)]⟩ .

(2.41)

2.2. Thermal conductance

19

This expression of the heat current is called the Meir-Wingreen-Landauer formula.
The thermal conductance is defined as⟨ J⟩.
∆T →0 ∆Tκ ≡ lim(2.42)

Using the formula (2.40), the linear thermal conductance is given as
αγkB
κ=
8

Z ∞
0

2
βω/2
dω Im [χ(ω )] Ĩ (ω ),sinh( βω/2)(2.43)

where χ(ω ) is evaluated for the thermal equilibrium and β = 1/kB T. Thus, we needto calculate the dynamical susceptibility, χ(ω ), for evaluating the linear thermal conductance.
For convenience of discussion, we also introduce a symmetrized correlation function and its Fourier transformation:
S(t) ≡
S(ω ) =1⟨σz (t)σz (0) + σz (0)σz (t)⟩ ,
Z2
∞

−∞

dt eiωt S(t).

(2.44)(2.45)

From the fluctuation-dissipation theorem [15], the imaginary part of the dynamicalsusceptibility is related to S(ω ) asβωS(ω ) = cothIm [χ(ω )] .
(2.46)
2
Then, the thermal conductance is rewritten in terms of the symmetrized correlation
function, S(ω ), as
κ=

αγkB
8

Z ∞
0


dω tanhβω2S(ω ) Ĩ (ω )

2
βω/2
.
sinh( βω/2)(2.47)

21

Chapter 3

Transport Processes and Quantum
Critical Phenomena
The dynamics of dissipative two-state systems have long been studied using various
approximations [15, 6]. In this chapter, we re-examine such analytical approximations from the viewpoint of heat transport. In Section 3.1, we first consider the effective tunneling amplitude and discuss a quantum phase transition driven by strong
system-reservoir coupling. Next, we consider the three mechanisms, which we call‘sequential tunneling’ (Section 3.2), ‘co-tunneling’ (Section 3.3), and ‘incoherent tunneling’ (Section 3.4). We derive analytical expressions for the thermal conductancein each transport process. We also introduce the noninteracting-blip approximation(NIBA) in Section 3.5.
In Section 3.6, we consider in detail the quantum phase transition for the subOhmic case. In the critical regime near the quantum phase transition, the temperature dependence of thermal conductance becomes power law which is reflected
the nature of the quantum phase transition. We derive the critical exponent of thethermal conductance.

3.1

Effective tunneling amplitude

One important effect of the system-reservoir coupling is renormalization of the tunneling amplitude ∆. In this section, we derive the effective tunneling amplitude via
adiabatic renormalization [15, 6] (see Fig. 3.1).
We assume that high-frequency oscillators (∆0 ≪ pωc < ω, p ∈ (0, 1)) instantaneously adjust their quantum states to the current value of σz . Then, the two lowest
energy eigenstates for the global system are described by| Eg ⟩ =| Ee ⟩ =1√ (|ΨL ⟩ + |ΨR ⟩) ,2
1
√ (|ΨL ⟩ − |ΨR ⟩) ,
2

(3.1)
(3.2)

where |ΨL ⟩ and |ΨR ⟩ are give by−|ΨL ⟩ = |σz = −1⟩ ⊗ ∏ ′ |ψν,k
⟩,

(3.3)

ν,k

+
|ΨR ⟩ = |σz = +1⟩ ⊗ ∏ ′ |ψν,k
⟩,

(3.4)

ν,k

respectively. Here the prime symbol indicates that the product is in the range pωc <±ων,k . |ψν,k⟩ is the ground state of the oscillator with wavenumber k in reservoir ν22

Chapter 3. Transport Processes and Quantum Critical Phenomena

F IGURE 3.1: A conceptual diagram of the adiabatic renormalization.
The spectral density is separated into two parts, Iν0 (ω ) = Iνlf (ω ) +Iνhf (ω ). In the first step, the tunneling amplitude ∆0 is renormalized
into the dressed tunneling amplitude ∆ by high-frequency oscillators(orange arrow). In the second step, the dressed tunneling amplitude∆ is renormalized into ∆( p) by oscillators with pωc < ω < ωc (bluearrow). In general, the spectral density with high-frequency oscillators is complex shapes.

when the state of the local two-state system is located at right (left) well; it is obtained
0 ⟩ for the isolated oscillators as
by translation of the ground-state wavefunction |ψν,k
±
0
|ψν,k
⟩ = e±iδν,k pν,k |ψν,k⟩,Cν,k q0
λν,k1=pδν,k =.
22
ωmν,k ων,k2mν,k ων,k ν,k

(3.5)
(3.6)

Adiabatic renormalization suggests that the tunneling amplitude is renormalized by
the overlap between the ground states of the oscillators for different localized states(σz = ±1):00+−∆( p) = ∆0 ∏ ′ ⟨ψν,k|ψν,k⟩ = ∆0 ∏ ′ ⟨ψν,k| e−2iδν,k pν,k |ψν,k
⟩
ν,k

ν,k

"

=
=
=
=

#

λ
0
†
0
∆0 ⟨ψν,k| exp ∑ ′ ν,k bν,k− bν,k |ψν,k
⟩
ω
ν,k
ν,k
"
#

 n
∞
1
0
′ λν,k†0∆0 ∑⟨ψν,k | ∑bν,k − bν,k
|ψν,k
⟩
n!
ω
ν,k
n =0
ν,k
!n
!
2
2
∞
λ
λ
1
1
1
− ∑ ′ ν,k∆0 ∑= ∆0 exp − ∑ ′ ν,k22n!
2
2 ν,k ων,k
ω
n =0
ν,k
ν,k
"
#
Z ∞
1
Iν0 (ω )∆0 exp − ∑dω.
2 ν pωc
ω2


(3.7)

In general, the spectral density can be separated into two parts (see Fig. 3.1):Iν0 (ω ) = Iνlf (ω ) + Iνhf (ω ),(3.8)

3.1. Effective tunneling amplitude

23

F IGURE 3.2: Schematics of the ground-state wavefunction (a) below
the transition (0 ≤ α < αc ) and (b) above the transition (α > αc ). The
former state is delocalized, whereas the latter is localized at one of thetwo wells. For the localized state, quantum tunneling between thetwo wells is forbidden since the overlap integral between the statesin the two wells vanishes.

where Iνlf (ω ) is the spectral density for the low-frequency oscillators and given by
the power law behavior (2.18) (Iνlf (ω ) = Iν (ω )), and Iνhf (ω ) represents the highfrequency oscillators and has complex behavior. In the first stage of the adiabatic
renormalization procedure, we consider high-frequency oscillators, ω > ωc , withthe spectral density Iνhf (ω ). Then, the bare tunneling amplitude, ∆0 , is changed into
the dressed tunneling amplitude, ∆, by renormalization of the high-frequency oscillators (the orange arrow in Fig. 3.1):
"
#
Z ∞
1
Iνlf (ω ) + Iνhf (ω )∆( p) = ∆0 exp − ∑dω
2 ν pωc
ω2
"
#
"
#
Z ∞
Z ∞
1
1
Iνlf (ω )Iν (ω )
= ∆ exp − ∑= ∆ exp − ∑dωdω
2 ν pωcω22 ν pωc
ω2#"
Z ∞(ω/ωc )s−1 −ω/ωc
e
.
(3.9)
= ∆ exp −αdωωpωcIn the second stage, we consider the low-frequency oscillators with the power law
spectral density Iνlf (ω ) (2.18). If the renormalized tunneling amplitude ∆( p) is lessthan pωc , the adiabatic renormalization can continue by reducing the factor p (theblue arrow in Fig. 3.1). If ∆( p∗ ) = p∗ ωc holds at p = p∗ , adiabatic renormalizationmust be stopped there and the finite effective tunneling amplitude ∆eff = ∆( p∗ ),which is a characteristic energy scale of the global system, is obtained. In contrast,
if ∆( p) < pωc holds for an arbitrary value of p, adiabatic renormalization can becompleted even at p = 0, yielding an effective tunneling amplitude of zero (∆eff =0).
For the Ohmic case (s = 1), the effective tunneling amplitude is obtained asfollows:  α/(1−α)∆ ∆, (0 < α < 1),′∆eff=(3.10)ωc
0,(1 ≥ α ).

24

Chapter 3. Transport Processes and Quantum Critical Phenomena

F IGURE 3.3: Schematic of the sequential tunneling process. Heattransport occurs by a combination of (a) phonon (photon) absorptionand (b) phonon (photon) emission.

In this thesis, we employ a modified effective tunneling amplitude multiplied by adimensionless function of α [15]:′∆eff = [Γ(1 − 2α) cos(πα)]1/2(1−α) ∆eff.

(3.11)

This result indicates a phase transition at zero temperature, for which the criticalvalue of the system-reservoir coupling is α = 1 [49, 50]. For system-reservoir couplings below the transition (0 < α < 1), the ground state is non-degenerate, as
shown in Fig. 3.2 (a), indicating the coherent superposition of the two localized states|↑⟩ and |↓⟩. We call this ground state ‘delocalized’. For strong system-reservoir couplings above the transition (1 < α), the coherent superposition of the two localizedstates is completely broken, leading to the doubly-degenerate ground states shownin Fig. 3.2 (b). We call this ground state ‘localized’. In this localized regime, quantum
tunneling between the wells is forbidden at zero temperature since there is no mixing (∆eff = 0) between the two localized states. Thus, the present quantum phasetransition can be recognized as a ‘localization’ transition that separates the delocalized and localized regimes at zero temperature.
For the sub-Ohmic case (0 < s < 1), the adiabatic renormalization always leads
to an effective tunneling amplitude of zero (∆eff = 0). This is correct in the limit∆/ωc → 0, as discussed in a previous study [6]. However, for a finite value of∆/ωc , the naive adiabatic renormalization procedure yields incorrect results andshould be improved. In subsequent theoretical studies [51, 52] , it was found thatthe localization transition actually occurred at a critical system-reservoir coupling(α = αc ), where the critical value αc depended on both s and ∆/ωc . The existence ofthe localization transition was also confirmed via numerical calculations [36, 37]. Insummary, for the sub-Ohmic case, the ground state is delocalized for 0 < α < αc , asshown in Fig. 3.2 (a), and localized for αc < α as shown in Fig. 3.2 (b).
For the super-Ohmic case (s > 1), the effective tunneling amplitude is always
finite:∆eff = ∆e−αΓ(s−1) ,(3.12)

where Γ(z) is the Gamma function. Therefore, there is no localization transition andthe ground state is always delocalized, as shown in Fig. 3.2 (a).

3.2. Sequential tunneling

3.2

25

Sequential tunneling

For weak system-reservoir couplings (α ≪ 1), the system and the reservoirs arealmost decoupled and the interaction Hamiltonian HI,ν can be regarded as a perturbation. For the second-order perturbation, the system dynamics are describedby a stochastic transition between the ground state (σx = +1) and the excited state(σx = −1), as shown in Fig. 3.3. The transition from the ground state to the excited state involves phonon (photon) absorption, and the inverse transition involves
phonon (photon) emission. A combination of these two processes induces heattransport. We refer to this type of transport process as ‘sequential tunneling’. The
transition rates for processes of absorption and emission can be calculated using
Fermi’s golden rule:Γ j = 2π |⟨ f | HI |i ⟩|2 δ(ϵ f − ϵi ),(3.13)

where |i ⟩ (| f ⟩) and ϵi (ϵ f ) are eigenstate of HS + HB and its eigenvalue, respectively.
Introducing the lowering and raising operatorsσz − iσy
1 1 −1+σx ≡=,(3.14)
2
2 1 −1σz + iσy
1
1
1
−
=
,
(3.15)
σx ≡
−1 −1
2
2
respectively, the interaction Hamiltonian HI,ν is rewritten in terms of σx+ and σx− as
HI,ν = −



1
+
+ †
−
− †
λ
σ
b
+
σ
b
+
σ
b
+
σ
b
ν,k
x ν,k
x ν,k
x ν,k
x ν,k .
2∑ν,k(3.16)

Here the second and third terms on the right side of Eq. (3.16) correspond to theprocess of emission and absorption, respectively. Substituting the expression of theinteraction Hamiltonian (3.16) corresponding to the process of emission into Fermi’s
golden rule (3.13), we obtain the transition rates for process of emission as follows:Γe = 2π ⟨σx = +1| ⊗ ∑ ⟨nν2 ,k2 |ν2 ,k21†− ∑ λν,k σx+ bν,k
2 ν,k

!

2

|σx = −1⟩ ⊗ ∑ |nν1 ,k1 ⟩ν1 ,k1×δ(ων,k − ∆)2λ2†
= 2π ∑ ∑ ∑ ν,k ⟨nν2 ,k2 | bν,k|nν1 ,k1 ⟩ δ(ων,k − ∆)4ν,k ν ,k ν2 ,k2
1= 2π ∑
ν,k

=

1

λ2ν,k
4

[nB (ων,k ) + 1] δ(ων,k − ∆)ππ[nB (∆) + 1] ∑ λ2ν,k δ(ων,k − ∆) = I (∆) [nB (∆) + 1] ,2
2
ν,k(3.17)

26

Chapter 3. Transport Processes and Quantum Critical Phenomena

 −1where I (ω ) = IL (ω ) + IR (ω ) and nB (ω ) = e βω − 1is the Bose-Einstein distribution function. Similarly, the rate for the process of absorption is obtained byΓa = 2π ⟨σx = −1| ⊗ ∑ ⟨nν2 ,k2 |ν2 ,k2=1− ∑ λν,k σx− bν,k
2 ν,k

!

2

|σx = +1⟩ ⊗ ∑ |nν1 ,k1 ⟩×δ(∆ − ων,k )π
I ( ∆ ) nB ( ∆ ).
2

ν1 ,k1(3.18)

For weak system-reservoir coupling, the density matrix of the two-state system
can be described by the quantum master equation in the Lindblad form [53]:


dρ(t)
1 †
†
†
= −i [ HS , ρ(t)] + ∑ Γ j L j ρ(t) L j −(3.19)
L j L j ρ + ρL j L j ,
dt
2
j=a,e
where La = σx− and Le = σx+ are Lindblad operators. Since Hilbert space of the twostate system is spanned by just two states, an excited state |σx = −1⟩ and a groundstate |σx = +1⟩, it is convenient to represent the density matrix as
11
1 + ⟨σz (t)⟩⟨σx (t)⟩ − i ⟨σy (t)⟩ρ(t) = (1 + ⟨⃗σ(t)⟩ ·⃗σ) =. (3.20)⟨σx (t)⟩ + i ⟨σy (t)⟩1 − ⟨σz (t)⟩22With the help of the algebra of the Pauli matrices, the Lindblad equation (3.19) gives
the differential equations:d ⟨σx (t)⟩dtd ⟨σy (t)⟩dtd ⟨σz (t)⟩dt= Γe (1 − ⟨σx (t)⟩) − Γa (1 + ⟨σx (t)⟩) ,(3.21)= ∆ ⟨σz (t)⟩ − Γ ⟨σy (t)⟩ ,(3.22)= −∆ ⟨σy (t)⟩ − Γ ⟨σz (t)⟩ ,(3.23)

where Γ = (Γa + Γe )/2. We observe that the x-component, ⟨σx (t)⟩, given by Eq. (3.21)
decays exponentially with rate 2Γ. The differential equations for y- and z-components(3.22), (3.23) are written in matrix form




d
⟨σz (t)⟩⟨σz (t)⟩= G,(3.24)⟨σy (t)⟩⟨σy (t)⟩dt
where G is a coefficient matrix

G =

−Γ −∆
+∆ −Γ


.

(3.25)

Using the quantum regression theorem [53], which states that the correlation function satisfies the same system of differential equations:d⟨σi (t)σj (0)⟩ = ∑ Gik ⟨σk (t)σj (0)⟩ ,dt
k(3.26)

3.2. Sequential tunneling27

we obtain the differential equations for the correlation functionsd⟨σz (t)σz (0)⟩ = −Γ ⟨σz (t)σz (0)⟩ − ∆ ⟨σy (t)σz (0)⟩ ,dtd⟨σy (t)σz (0)⟩ = +∆ ⟨σz (t)σz (0)⟩ − Γ ⟨σy (t)σz (0)⟩ .
dt(3.27)(3.28)

From Eq. (3.27) and (3.28), under boundary conditions:⟨σz (0)σz (0)⟩ = 1,

(3.29)d⟨σz (t)σz (0)⟩|t=0 = 0,dt

(3.30)

the z-z correlation function is obtained as⟨σz (t)σz (0)⟩ =1 −Γte[∆ cos(∆t) + Γ sin(∆t)] .
∆(3.31)

Since the symmetrized correlation function (2.44) has time-reversal symmetry, it is
written as
S(t) = ⟨σz (|t|)σz (0)⟩ =1 −Γ|t|e[∆ cos(∆|t|) + Γ sin(∆|t|)] .
∆(3.32)

Performing a Fourier transform on it, we finally derive the symmetrized correlation
function as follows:
S(ω ) =

4Γ(∆2 + Γ2 ).
[(ω − ∆)2 + Γ2 ] [(ω + ∆)2 + Γ2 ](3.33)

The correlation function S(ω ) has two peaks at ω = ±∆, reflecting the coherent
system dynamics. Because Γ ≪ ∆ always holds in the weak coupling regime, thecorrelation function is approximated asS(ω ) ≃ π [δ(ω − ∆) + δ(ω + ∆)] ,(3.34)

where is δ( x ) is the delta function. The thermal conductance for the weak coupling
regime is obtained by substituting Eq. (3.34) into Eq. (2.47) as follows
παγkB
κ≃
tanh
8



β∆
2



2
β∆/2
Ĩ (∆).
sinh( β∆/2)(3.35)

This result is identical to the formula derived in previous research [29, 7] using themaster equation approach and is consistent with the perturbation theory [23]. For
actual comparison with the numerical simulation in Section 5, we improve the approximation by replacing ∆ with ∆eff using adiabatic renormalization. Then, thesymmetrized correlation function is improved toS(ω ) =

4Γ(∆2eff + Γ2 ).
[(ω − ∆eff )2 + Γ2 ] [(ω + ∆eff )2 + Γ2 ](3.36)For this function to have two peaks at ω = ∆eff , it is necessary to hold the condition
as followsΓ=π[1 + 2nB (∆eff )] I (∆eff ) ≪ ∆eff .
4(3.37)

28

Chapter 3. Transport Processes and Quantum Critical Phenomena

F IGURE 3.4: Schematic of the co-tunneling process. At kB T ≪ ∆eff ,heat transport via a virtual excitation in the local system is dominant.

For the sub-Ohmic case (s < 1), this condition is never satisfied, indicating the absence of a sequential tunneling regime. For the Ohmic case (s = 1), the conditionis equivalent to α ≪ 1, whereas for the super-Ohmic case (s > 1), the conditionis always satisfied for a moderate temperature (kB T ∼ ∆eff ). At high temperatures(kB T ≫ ∆eff ), the condition is always satisfied for s ≥ 2, whereas for 1 < s < 2, it
becomesωc ∆eff 2−s
∗
T<T =
.
(3.38)
αkB ωcwhere T ∗ is the crossover temperature.
The formula for sequential tunneling (3.35) predicts the exponential decrease in
the thermal conductance as the temperature is lowered. At low temperatures, thethermal conductance behaves as κ ∝ e−∆eff /kB T /T 2 ; this is because the transition fromthe ground state to the excited state is strongly suppressed if the thermal fluctuationis smaller than the effective energy splitting, i.e., when kB T ≪ ∆eff . When the sequential tunneling process is strongly suppressed at low temperatures, Eq. (3.35)
becomes invalid since another process becomes dominant, as discussed in the next
subsection.

3.3

Co-tunneling

At low temperatures, heat transport via the virtual excitation of the local two-statesystem becomes dominant (see Fig. 3.4); this transport process is known as ‘cotunneling’. In a previous study [23], an analytical expression for the thermal conductance was derived using the fourth-order perturbation theory with respect to theinteraction HI,ν . However, in this calculation, the renormalization of the tunneling
amplitude at a low temperature has not been considered.
Here, we derive an asymptotically exact formula for the thermal conductancewithout any approximations [54]. For this purpose, we focus on an asymptoticallyexact relation called the generalized Shiba relation [55, 56]: χ 2
S(ω )
0
,
= πα
ω →0+ Ĩ ( ω )2lim(3.39)

where χ0 is the static susceptibility defined in Eq. (2.19). This exact relation holds
at low temperatures (kB T ≪ ∆eff ) for arbitrary environments and arbitrary systemreservoir couplings. At low temperatures (kB T ≪ ∆eff ), the dominant contribution
to the integral of Eq. (2.47) comes from the low-frequency part (0 ≤ ω ≃ kB T ≪ ∆eff )3.4. Incoherent tunneling29

F IGURE 3.5: Schematic of the incoherent tunneling process. The
wavefunction is localized in the two wells, and a stochastic transition occurs between them.

due to the factor of the Bose-Einstein distribution function. By substituting the lowfrequency asymptotic form S(ω ) ≃ πα (χ0 /2)2 Ĩ (ω ) into Eq. (2.47), we obtain [54]πkB χ20
κ≃
8

Z ∞
0

2
βω/2
dω IL (ω ) IR (ω ).
sinh( βω/2)(3.40)

This expression is similar to the co-tunneling formula in previous studies [23, 57] but
significantly differs in terms of the static susceptibility, χ0 , which considers higherorder processes. Eq. (3.40) can be rewritten as1
kB T 2s+1
πkB αL αR ωc3 χ20
F ( s ),
2
ωc
2

Z ∞
x/2
2s
,
F (s) =dx xsinh( x/2)0
κ ≃(3.41)(3.42)

where F (s) is a dimensionless function of s. Thus, we find that the thermal conductance is proportional to T 2s+1 at low temperatures. The same temperature dependence has been derived by the perturbation theory [23, 57]. However, the perturbation theory cannot treat renormalization effect due to higher-order processes on
the static susceptibility and fails in predicting a correct prefactor including χ0 . In
contrast, the present result given in Eq. (3.40) is asymptotically exact, incorporating
the renormalization effect appropriately.
The co-tunneling formula (3.40) holds universally at low temperatures for an
arbitrary exponent, s, as long as the ground state of the system is delocalized (∆eff >0). In a previous study [29], the thermal conductance in the Ohmic case (s = 1)
was shown to be proportional to T 3 , which is consistent with Eq. (3.40), and this T 3 dependence was discussed in terms of the emergence of the Kondo effect. However,
it is worth noting that the power-law temperature dependences are derived in a
unified way even in the non-Ohmic cases. These temperature dependences result
from non-trivial many-body effects due to strong mixing between the system andthe reservoirs.

3.4

Incoherent tunneling

For a strong reservoir-system coupling, the coherent superposition of the two localized states is completely broken. In such a situation, heat transport is induced
by stochastic dynamics between the two localized states |↑⟩ and |↓⟩, as shown in
Fig. 3.5. We call this transport process ‘incoherent tunneling’.

30

Chapter 3. Transport Processes and Quantum Critical Phenomena

Within the Markov approximation [58], the stochastic dynamics of the system is
described by the master equation:dPL (t)= −ΓPL (t) + ΓPR (t),dtdPR (t)= ΓPL (t) − ΓPR (t),dt

(3.43)

where PL (t) and PR (t)(= 1 − PL (t)) are the probabilities that the wavefunction of
the system is localized at the well on the left-hand side (σz = −1) and that on theright-hand side (σz = +1), respectively, at time t. The transition rate Γ is calculated
via second-order perturbation with respect to the Hamiltonian HS as follows [6]:
∆2
2

Z ∞



′
dτ e−Q (τ ) cos Q′′ (τ ) ,
−∞


Z ∞
βω
I (ω )′coth
Q (τ ) =dω[1 − cos(ωτ )] ,
ω2
2
0


Z ∞
βω
I (ω )′′Q (τ ) =dω
cothsin(ωτ ).
ω2
2
0
Γ =

(3.44)(3.45)(3.46)

Note that this expression for the transition rate of incoherent tunneling is valid whenΓ ≪ kB T [58]. By solving the master equation (3.43), the symmetrized correlation
function is calculated asS(ω ) =4Γ
ω 2 + 4Γ2.

(3.47)

In contrast to sequential tunneling, S(ω ) has only one peak at ω = 0 with a width of
2Γ, indicating the destruction of the superposition of the two localized states.
The long-term dynamics are well described by the Markov approximation [6].
Therefore, one may expect that the thermal conductance in the incoherent tunnelingregime would be well approximated by substituting Eq. (3.44) into Eq. (2.47). However, the results of the Markov approximation show clear deviation from the numerical results, as discussed in Section 5. The reason for this is summarized as follows.
Note that incoherent tunneling occurs when Γ ≪ kB T. Under this condition, theintegrand of Eq. (2.47) is proportional to ω s−2 for Γ ≪ ω ≪ kB T since S(ω ) ∝ ω −2(see Eq. (3.47)). Then, the integral in Eq. (2.47) diverges if the high-frequency cutoff occurring due to the Bose-Einstein distribution function is absent. This indicatesthat the high-frequency part of the integral in Eq. (2.47) makes the dominant contribution to the thermal conductance. Although the Markov approximation yields
reasonable results for the low-frequency behavior of S(ω ), it fails to reproduce theaccurate high-frequency behavior of S(ω ) in general, leading to incorrect results forthe thermal conductance.

3.5

Noninteracting-blip approximation (NIBA)

To study the short-term (high-frequency) dynamics in the incoherent tunneling regime,
we introduce the noninteracting-blip approximation (NIBA), which is a natural extension of the Markov approximation in the previous section [6, 15]. In the NIBA, theFourier transformation of the symmetrized correlation function is expressed as [15](the detail derivation is given in Appendix C)
1S(ω ) = 2Re,(3.48)−iω + Σ(−iω )3.6. Quantum critical phenomena for the sub-Ohmic case31

where Σ(λ = −iω ) is the frequency-dependent self-energy defined as
Σ(λ) = ∆

2

Z ∞
0



′
dτ e−λτ e−Q (τ ) cos Q′′ (τ ) .

(3.49)

Here, Q′ (τ ) and Q′′ (τ ) are given by Eqs. (3.45) and (3.46), respectively. The thermal
conductance is then calculated by substituting Eq. (3.48) and (3.49) into Eq. (2.47).
From the definition, it is easy to check that the NIBA reproduces the Markov approximation if one neglect the frequency dependence of the self-energy and replace
it with the zero-frequency value Σ(0) = 2Γ. Since the NIBA appropriately considers
the non-Markovian properties, it is suitable to describe the thermal conductance in
the incoherent tunneling regime.
The condition for the NIBA is well known [15, 6]. As expected from the fact that
the NIBA is an extension of the Markov approximation, it works well for the incoherent tunneling regime. Roughly, the incoherent tunneling mechanism becomes crucial in a regime wherein both the sequential tunneling formula and the co-tunnelingformula fail. (a) The NIBA holds at moderate-to-high temperatures in the sub-Ohmic(0 < s < 1) and Ohmic cases (s = 1). (b) It holds for T > T ∗ in the super-Ohmic caseof 1 < s < 2, where T ∗ is the crossover temperature defined in Eq. (3.38). Note that
the NIBA never holds for s ≥ 2 since the crossover temperature T ∗ diverges.
Here, the NIBA has been introduced to improve the Markov approximation in
the incoherent regime. This introduction of the NIBA may give an impression to thereaders that the NIBA is a good approximation only in the incoherent regime. However, the NIBA is known to be applicable for a wider parameter region not restricted
to the incoherent regime [15, 6]. The NIBA holds also in the weak coupling regime(α ≪ 1) at arbitrary temperature for the unbiased case (ϵ = 0), where the interblip
interaction is shown to be much weaker than the intrablip interaction. For this reason, the NIBA yields almost the same result as the sequential tunneling formula orthe co-tunneling formula if the system-reservoir coupling is sufficiently weak.
In Section 5, we show that the NIBA is an excellent approximation for reproducing the numerical results for a wide region of the parameter space at moderate-tohigh temperatures. Thus, the short-term (high-frequency) non-Markovian behavior
in the system dynamics is important for calculating the thermal conductance in theincoherent tunneling regime.

3.6

Quantum critical phenomena for the sub-Ohmic case

We explained three tunneling processes (sequential tunneling, co-tunneling, and incoherent tunneling) and derived the analytical expressions for the thermal conductance in each transport process. However, the temperature dependence of the thermal conductance near the quantum phase transition for the sub-Ohmic case (as described in Section 3.1) has not been discussed. In this section, we discuss quantum
critical phenomena induced by the second-order phase transition for the sub-Ohmiccase at zero temperature and determine the critical exponent of the thermal conductance using other critical exponents discussed in Section 1.5. In the quantum criticalregime near the quantum phase transition, the thermal conductance exhibits a distinctive power-law temperature dependence reflecting the nature of the quantum
phase transition:κ ∝ T c , ( α = αc )(3.50)

32

Chapter 3. Transport Processes and Quantum Critical Phenomena

where c is the critical exponent dependent on s.
Let us discuss the critical exponent, c, defined in Eq. (3.50). The static susceptibility defined in Eq. (2.19) is rewritten byχ0 = β ⟨σ̄z2 ⟩eq ,(3.51)

1β(3.52)

σ̄z =

Z β
0

dτ σz (τ ).

where ⟨· · ·⟩eq indicates the thermal average and σz (τ ) is the imaginary time path.
Combining Eq. (3.51) with (3.52), the static susceptibility is expressed as
χ0 =

Z ∞
0

dτ C (τ ),(3.53)

where C (τ ) is the correlation function C (τ ) ≡ ⟨σz (τ )σz (0)⟩eq . At the critical point,
the correlation function exhibits the power-law decay:
C (τ ) = C ( β − τ ) ∼ τ −η , (ωc−1 ≪ τ ≪ β/2),(3.54)

where η is the critical exponent (1.21). Then, the temperature dependence of thestatic susceptibility at the critical point is obtained:χ 0 ∼ β 1− η .

(3.55)

By using the Fourier transformation of the correlation function:
C (iωn ) =Z β0dτ eiωn τ C (τ ),(3.56)

the critical behavior of the imaginary part of the dynamical susceptibility χ(ω ) =C (iωn → ω + iδ) is obtained:
Im [χ(ω )] ∼ ω η −1 .

(3.57)

Substituting this into Eq. (2.43), the thermal conductance at the critical point behavesas κ ∼ T c , where the exponent is given by:
c = s + η.

(3.58)

The critical exponent η is a function of s and has been analyzed in previous theoretical studies [37, 41]. The phase transition for 0 < s ≤ 0.5 belongs to the mean-field
universality class and leads to η = 1/2 (see Eq. (1.22)). For 0.5 < s < 1, η is a
nontrivial function of s and is evaluated by the ε-expansion (see Eq. (1.24)) [41]. Insummary, the exponent of the thermal conductance is given as follows:1s + ,2c=
1 − 1 ε − A(s) ε2 + O(ε3 ),23s(0 < s ≤ 0.5),(3.59)(s > 0.5),

where ε = 2s − 1, A(s) = s [ψ(1) − 2ψ(s/2) + ψ(s)], and ψ(z) is the digamma function. The detail has been discussed using the numerical calculations in Ref. [41].

33

Chapter 4

Continuous-time Quantum Monte
Carlo Method
While the analytical approaches discussed in the previous section are sufficiently
powerful for clarifying the mechanism of heat transport in a two-state system, thedetailed conditions justifying each approximation are not trivial. To understandall features of heat transport, unbiased numerical simulation without any approximation would be helpful. In this thesis, we employ the continuous-time Monte
Carlo (CTQMC) algorithm [59]. According to this algorithm, the partition function
is rewritten in path integral form with respect to an imaginary time path, σz (τ ),(kink representation) and the weight of this path is defined. We describe the way to
calculate the imaginary part of the dynamical susceptibility, Im [χ(ω )], by applying
the Monte Carlo method to the kink representation using the cluster update algorithm [37]. In addition, we explain the way to determine the critical point of thequantum phase transition for the sub-Ohmic case.

4.1

The kink representation of the spin-boson model

The partition function of the spin-boson model (2.13) is written in the imaginary
time path integral form as [15, 37] (for the detained derivation, see Appendix B)hi Z
Z = Tr e− βH = D σz (τ ) e−S[σz (τ )] ,S [σz (τ )] = −

1
4

Z β

dτ

Z τ

0

0

dτ ′ σz (τ )K (τ − τ ′ )σz (τ ′ ),

(4.1)
(4.2)

where σz (τ ) is a spin path defined on the imaginary time axis, D σz (τ ) indicates theintegral measure for all possible paths σz (τ ), S [σz (τ )] is a influence function, andK (τ ) is a kernel defined as
K (τ ) =Z ∞0dω I (ω )cosh [ω ( β/2 − τ )].
sinh( βω/2)(4.3)

The kernel has the symmetry K (τ ) = K ( β − τ ) and asymptotic form K (τ ) ∼ τ −(1+s)
for ωc−1 ≪ τ ≪ β/2. Note that the partition function (4.1) shows the equivalence
between the spin-boson model and the ferromagnetic Ising model with the longrange spin-spin interaction K (τ ) [60].
To apply the CTQMC method, we change over from the spin representation (4.1)
to the kink representation. As shown in Fig. 2.1 (a), the path σz (τ ) is assigned by analternative configuration of kinks (jumps from σz = −1 to σz = +1) and anti-kinks(jumps from σz = +1 to σz = −1) and described by the positions τi (i = 1, 2, . . . , 2n)

34Chapter 4. Continuous-time Quantum Monte Carlo Methodof the kinks (qi = +1) and anti-kinks (qi = −1). Thus, the derivative of the spin pathreads
2ndσz (τ )= ∑ 2qi δ(τ − τi ),
dτ
i =1

(4.4)

where n is the number of the pairs of kinks and anti-kinks. By substituting Eq. (4.4)
into (4.1), we obtain the kink representation of the partition function:


Z τ2n
Z τ2
2n
∞  2n Z β∆dτ2ndτ2n−1 · · ·dτ1 exp  ∑ qi q j W (τi − τj ) , (4.5)
Z= ∑
2
0
0
0
n =0
⟨i,j⟩
where W (τ ) is obtained from the relation W ′′ (τ ) = −K (τ ) asW (τ ) =

4.2

Z ∞
0

dω

I (ω ) cosh ( βω/2) − cosh [ω ( β/2 − τ )].
ω2sinh ( βω/2)(4.6)

Calculation algorithm

We apply the CTQMC method to the partition function for the kink representation (4.5). The present CTQMC algorithm [59] employs a cluster-flip update similar
to that in the Swendsen–Wang cluster algorithm [61]. Employing this algorithm, we
can overcome the problem of the critical slowing down (the correlation time diverges
near the critical point of second-order phase transition). The cluster-flip update isconstructed as follows [37] (see Fig. 4.1):
1. Inserting new vertices with Poisson statics given by P(∆τ ) = Γe−Γ∆τ withmean value Γ−1 = 2/∆.
2. Connecting two segments (the line segments between neighboring vertices), si
and s j , with the probabilityhip[si , s j ] = 1 − δσz (si ),σz (s j ) 1 − e−2A ,(4.7)

A = W (τi−1 − τj−1 ) − W (τi−1 − τj )−W (τi − τj−1 ) + W (τi − τj ),(4.8)

where σz (si ) is the value of σz in the segment si , and positions of the vertices at
the two edges of the segment si are denoted by τi−1 and τi , respectively.
3. Flipping each segment cluster with probability 1/2.
4. Removing the redundant vertices within segments.
The Monte Carlo data presented in this paper typically represent averages over 103 104 updates at low temperatures and 107 –108 updates at high temperatures.
Using the CTQMC sampling method, we evaluate the spin correlation function
C (iωn ) defined in Eq. (3.56) as follows:C (iωn ) =1⟨|ρ(iωn )|2 ⟩ ,
βωn2

(4.9)

4.2. Calculation algorithm

35

(a)
+1
0
−1

(b) (i)

0

(ii)
(iii)
(iv)
(v)
(c)
+1
0
−1
F IGURE 4.1: The spin path on the imaginary time axis for the spinboson model and the CTQMC update process: (a) An initial spin
path, (b-i) the vertex representation of (a), (b-ii) insertion of new cuts
using the Poisson distribution, (b-iii) connection of segments with theprobability (4.7), (b-iv) flipping each cluster with probability 1/2, (bv) removal of redundant cuts, and (c) the final spin path after the
update.

where ⟨· · ·⟩ denotes the average obtained via CTQMC sampling and ρ(iωn ) is theFourier transformation of ρ(τ ) = dσz (τ )/dτ. From Eq. (4.4), ρ(iωn ) is expressed as2n

ρ(iωn ) = ∑ 2(−1) j eiωn τj .

(4.10)

j =1

Then, to evaluate the thermal conductance (2.43) numerically, we need to calculate
as follows:
1. Calculating the Fourier transformation of the correlation function, C (iωn ), from
Eq. (4.9).
2. Performing the analytical continuation
χ(ω ) = C (iωn → ω + iδ).
3. Evaluating the thermal conductance from Eq. (2.43).

(4.11)

36

Chapter 4. Continuous-time Quantum Monte Carlo Method

For using the co-tunneling formula (3.40), we need to calculate the static susceptibility χ0 . Typically, a simple estimate χ0 ≃ 2/∆eff yields quantitatively correct
results. However, for the sub-Ohmic case, χ0 has non-trivial temperature dependence, even at low temperatures. For this case, we numerically calculate the static
susceptibility, χ0 , using the CTQMC method as follows:χ0 = β ⟨σ̄z2 ⟩ ,
σ̄z =

4.3

1
β

Z β
0

dτ σz (τ ) =(4.12)
2 2n−1(−1) j τj + 1.
β j∑
=0

(4.13)

Analytical continuation

To obtain the dynamical susceptibility, χ(ω ), we need to perform the analytical continuation (4.11). If the analytical expression for the Fourier transformation of thecorrelation function, C (iωn ) is clarified, we have only to replace iωn with ω + iδ.
However, since the Fourier transformation of the correlation function includes numerical errors, the analytical continuation is difficult. In this thesis, we usually employ Padé approximation [62].
In Padé approximation, given values of the analytical function for upper halfplane F (z) at z = iωn (n = 0, 1, . . . , N) numerically, F (ω ) is approximated by acontinued fractional function:
GN (ω ) =a ( ω − z N −1 )a1 a2 ( ω − z1 )
··· N
,
1+
1+
1

(4.14)

where ai are coefficients determined by GN (ω ). Using Thiele’s reciprocal difference
algorithm, we can determine coefficients ai as follows. We introduce a function
g1 ( z n ) = F ( z n ) ,(4.15)and define sequentially functions
gn ( z n ) =g n −1 ( z n −1 ) − g n −1 ( ω ).
( ω − z n −1 ) g n −1 ( ω )(4.16)

Then, coefficients ai are obtained by
a i = gi ( z i ) .

(4.17)In this thesis, F (z) and GN (ω ) correspond to the Fourier transformation of the correlation function C (iωn ) and dynamical susceptibility χ(ω ), respectively.
For the weak coupling regime, however, Padé approximation yields poor results
since the imaginary part of the pole nearest to the real frequency axis is small. In thiscase, we employ another approximation based on the fitting [63]. We assume that
the Fourier transformation if the correlation function as
C (iωn ) ≃aω03+ const,(ωn + λ)2 + ω03(4.18)

where a, ω0 , and λ are the fitting parameters determined using the least-squares
method. It is easy to obtain the imaginary part of the dynamical susceptibility
Im [χ(ω )] using the fitting function (4.18) with optimized parameters. Note that

4.4. Binder analysis

37

B

(a)

(c)

B

(b)

F IGURE 4.2: An example of the Binder parameter analysis. The results for s = 0.6 and ∆/ωc = 0.01 are shown. (a) The Binder parameter as a function of the coupling constant α for different temperatures. (b) Enlarged view of (a) around α = 0.08. Crosses represent the intersection points α β,2β between the two neighboring inverse temperature, βωc and 2βωc , respectively. (c) The intersection
points of the Binder parameters. The red solid curve shows the fitted
quadratic function. The dashed horizontal line indicates the criticalvalue αc = 0.0615 obtained via the present analysis.

this fitting method works well for weak couplings since it is compatible with the
dynamic susceptibility for the sequential tunneling process.

4.4

Binder analysis

In this section, we describe how to determine the critical point of the quantum phasetransition for the sub-Ohmic case (0 < s < 1) [37]. To study this quantum phasetransition, we utilize the finite size scaling form for observable O close to the criticalpoint∗⟨O⟩ ( T, α) = β xO fO ( βyt δ),(4.19)

where δ = (α − αc )/αc denotes the distance from the critical point, xO and fO are
the scaling exponent and scaling function of the observable O, respectively. The
exponent y∗t is related to the correlation length exponent ν, ξ ∼ δ−ν at δ ≃ 0, as

38

Chapter 4. Continuous-time Quantum Monte Carlo Method

follows [37, 41, 42]:

1

 ,
y∗t = ν
 1 + 1 − s,
ν 2(0 < s ≤ 0.5),(4.20)(s > 0.5).

We introduce the Binder parameter, which is defined as follows:!
1⟨σ̄z4 ⟩
B=,3−2
⟨σ̄z2 ⟩2

(4.21)

R
where σ̄z = β−1 dτ σz (τ ) indicates the average obtained via the Monte Carlo sampling. Since the Binder parameter has x B = 0, the critical point δ = 0 (α = αc )is determined as the point for which the Binder parameter is independent of thetemperature at sufficiently low temperatures. In Fig. 4.2, we show an example of
the Binder analysis for s = 0.6 and ∆/ωc = 0.01. The curve of the Binder parameter for different temperatures has intersection points around α = 0.08, as shownin Fig. 4.2 (a). To accurately determine the critical point, we consider the intersection points α β,2β between the two neighboring inverse temperatures, β and 2β (seeFig. 4.2 (b)), and plot the intersection points as a function of ( βωc )−1 , as shown in
Fig. 4.2 (c). By extrapolating α β,2β in the limit ( βωc )−1 → 0 using fitting to thequadratic function of ( βωc )−1 , the critical value αc = 0.0615 is obtained for thisparameter set.

39

Chapter 5

Numerical Results and Comparison
with Analytical Formulas
In this chapter, we perform numerical simulations based on the CTMQC method
and compare the simulation results with the analytical formulas introduced in Section 3. We separately consider the Ohmic (Section 5.1.1), sub-Ohmic (Section 5.1.2),
and super-Ohmic cases (Sections 5.1.3 and 5.1.4). The dynamics of the spin-bosonmodel has been studied by using various numerical methods. However, no systematic comparisons between analytical approximations and numerical simulations
have been performed in the context of heat transport near thermal equilibrium. This
comparison allows us to discuss the validity of various approximations critically.
In Section 5.2, we show numerical results near the quantum phase transition
point for the sub-Ohmic case and compare them with the analytical temperaturedependence of the thermal conductance discussed in Section 3.6.

5.1

Thermal conductance

5.1.1

The Ohmic case (s = 1)

In Fig. 5.1, we show the thermal conductances for α = 0.05, 0.1, 0.5, and 0.7 as functions of temperature. We plot the graph using the normalized temperature kB T/∆eff
and the normalized thermal conductance κ/(kB γ∆eff ), where ∆eff is the effective
tunneling amplitude defined in Eq. (3.11). As shown in Fig. 5.1, the numerical results fall on a universal scaling curve at each value of α regardless of the ratio ∆/ωc(∆/ωc ≪ 1) obtained via this normalization. This universal behavior is characteristic of the Kondo-like effect [29]. In Fig. 5.1 (c), we also show the exact solution (theToulouse point) for α = 0.5 (indicated by the brown dotted–dashed line) [15, 29, 63].
The agreement between the numerical results and the exact solution indicates thecorrectness of the CTQMC simulation.
At low temperatures (kB T ≪ ∆eff ), the numerical results agree well with thoseof the approximate formula for the co-tunneling process (Eq. (3.41); indicated byblue dashed lines in Fig. 5.1). In this regime, the thermal conductance is always
proportional to T 3 (= T 2s+1 ), which is consistent with the previous study [29].
At moderate (kB T ∼ ∆eff ) and high temperatures (kB T ≫ ∆eff ), the numerical results deviate from the co-tunneling formula and agree well with the NIBA (indicatedby black solid lines in Fig. 5.1). Note that the thermal conductance obtained by the
NIBA is proportional to T 3−2α at low temperatures, as shown in Fig. 5.1. The NIBA
agrees well even with the low-temperature numerical results for the weak systemreservoir coupling (α ≪ 1), whereas it deviates from these results as this coupling
becomes large. It is remarkable that the NIBA agrees well with the numerical results
at arbitrary temperatures for α ≪ 1, as shown in Fig. 5.1 (a).

40

Chapter 5. Numerical Results and Comparison with Analytical Formulas
(a)

(b)

(c)

(d)

F IGURE 5.1: The temperature dependence of the thermal conductance for (a) α = 0.05, (b) 0.1, (c) 0.5, and (d) 0.7. The symbols indicate
the numerical results obtained using the CTQMC method. The blacksolid, green dotted-dashed, blue dashed, and orange dotted lines represent the NIBA, sequential tunneling formula, co-tunneling formula,and Markov approximation for incoherent tunneling, respectively. In(c), the exact solution for the Toulouse point (α = 0.5) is indicated by
the brown dotted line.

In Fig. 5.1 (a) and (b), we also show the approximate formula for sequential tunneling (indicated by green dotted-dashed lines). As shown in this figure, the sequential tunneling formula at moderate temperatures (kB T ∼ ∆eff ) agrees with thenumerical results of the weak system-reservoir coupling (α ≪ 1). This agreement is
consistent with the previous study [29]. However, note that the NIBA agrees with
the numerical results for a wider temperature region than the sequential tunnelingformula.
The Markov approximation for incoherent tunneling, indicated by orange dotted lines in Fig. 5.1, clearly deviates from the numerical results for α = 0.05, 0.1,
and 0.7, indicating the importance of the non-Markovian properties of the system.
The Toulouse point (α = 0.5) is an exception, as shown in Figure 5.1 (c); the NIBA
coincides with the Markov approximation since at this point the self-energy in the
NIBA becomes independent of the frequency for the unbiased case [15]. A detailed
discussion on the failure of the Markov approximation is given in Section 5.1.2.
As described in Section 3.1, quantum phase transition occurs at αc = 1 for the
Ohmic case. For αc ≥ 1, the effective tunneling amplitude ∆eff becomes zero, indicating complete destruction of the superposition of the two localized states. Therefore, heat transport is induced by incoherent tunneling at arbitrary temperatures. In
Fig. 5.2, we show the thermal conductance for α = 1.0, 1.5, and 2.0 as a function oftemperature. As indicated by the black solid lines in the figure, the numerical results

5.1. Thermal conductance41

F IGURE 5.2: Temperature dependence of the thermal conductance.
The symbols with error bars indicate the numerical results obtainedusing the CTQMC method for α = 1.0, 1.5, and 2.0 from top to bottom. The black solid and orange dashed lines represent the NIBA andthe Markov approximation for incoherent tunneling.

(a)

(b)

F IGURE 5.3: (a) The thermal conductance calculated using the
CTQMC method for s = 0.9, ∆/ωc = 0.01 and α = 0.1. (b) The symmetrized correlation function calculated using the CTQMC method
at kB T = ωc /64 for parameters same as those considered in (a). The
black solid, blue dashed, and orange dotted lines represent the NIBA,
co-tunneling formula, and Markov approximation for incoherent tunneling, respectively.

agree well with the NIBA formula for arbitrary temperatures. Note that for α ≥ 1,
the condition for the co-tunneling regime kB T ≪ ∆eff is never satisfied. In Fig. 5.2,
we also show the Markov approximation for incoherent tunneling (indicated by theorange dashed line). For α ≥ 1, the difference between the NIBA and the Markov
approximation is not considerably large.

5.1.2

The sub-Ohmic case (0 < s < 1)

We first discuss the thermal conductance for the sub-Ohmic case wherein the systemreservoir coupling is below the critical value for the quantum phase transition. In
Fig. 5.3 (a), we show the thermal conductance as a function of the temperature fors = 0.9, ∆/ωc = 0.01, and α = 0.1, for which the ground state is delocalized(α < αc (s, ∆)). At moderate and high temperatures, the numerical results agree well

42

Chapter 5. Numerical Results and Comparison with Analytical Formulas
(a)

(b)

F IGURE 5.4: Temperature behavior of the thermal conductance calculated by a Monte Carlo simulation. The data represent results for
s = 0.6, ∆/ωc = 0.01, (a) α = 0.02, and (b) α = 0.1. The blacksolid lines and the blue dashed line represents the NIBA and the cotunneling formula, respectively.

with the NIBA, which is shown by the black solid line. We note that the sequentialtunneling formula cannot be applied to the sub-Ohmic case. At low temperatures(kB T ≪ ∆eff ), the numerical results agree well with the co-tunneling formula, showing T 2s+1 -dependence.
We also show the results of the Markov approximation for incoherent tunnelingby the orange dotted line in Fig. 5.3 (a). The Markov approximation clearly deviatesfrom the numerical results. To understand the failure of the Markov approximation,
we show the numerical and analytical result of the symmetrized correlation function
S(ω ) as a function of ω/ωc for kB T = ωc /64 in Fig. 5.3 (b). While the Markov approximation for the incoherent tunneling process agrees with the numerical results
at a low frequency, clear deviation is observed at higher frequencies; the numericalresult indicates that the high-frequency decay of S(ω ) is much faster than that of theMarkov approximation, which is proportional to ω −2 (see Eq. (3.47)). We note thatthe numerical result of S(ω ) is well reproduced by the NIBA at arbitrary frequencies. These observations indicate that the non-Markovian properties of the system
dynamics are important for obtaining correct thermal conductance results for the
sub-Ohmic case.
The quantum phase transition remarkably affects the temperature dependenceof the thermal conductance. In Fig. 5.4, we show the thermal conductance as a function of the temperature for s = 0.6 and ∆/ωc = 0.01, for which a quantum phasetransition occurs at α = αc = 0.0615. The detailed procedure for the determinationof the critical point is given in Section 4.4. Fig. 5.4 (a) shows the temperature dependence in the delocalized regime (α = 0.02 < αc ), for which ∆eff remains finite.
The numerical results agree well with the co-tunneling formula at low temperatures
and with the NIBA at moderate-to-high temperatures. This feature is the same asthat shown in Fig. 5.3. Fig. 5.4 (b) shows the temperature dependence in the localized regime (α = 0.1 > αc ), for which ∆eff = 0. Reflecting the quantum phasetransition, the numerical results agree with the NIBA at arbitrary temperatures, as
shown in Fig. 5.4 (b). Since the condition for the co-tunneling regime, kB T ≪ ∆eff ,is never satisfied for ∆eff = 0, the thermal conductance does not show a universal
T 2s+1 -dependence due to the co-tunneling process at low temperatures.

5.1. Thermal conductance
(a)

43
(b)

F IGURE 5.5: Temperature behavior of the thermal conductance calculated by Monte Carlo simulation. The data represent results for
s = 1.5, ∆/ωc = 0.05, (a) α = 0.1, and (b) α = 0.5. In both figures,
the black solid, blue dashed, green dotted-dashed, and orange dotted lines represent the NIBA, co-tunneling formula, sequential tunneling formula, and Markov approximation for incoherent tunneling,respectively.

5.1.3

The super-Ohmic case (1 < s < 2)

In Fig. 5.5, we show the numerical thermal conductance results obtained using the
CTQMC method as a function of temperature for s = 1.5. Here, the horizontal andvertical axes are the normalized temperature kB T/∆eff and the normalized thermal
conductance κ/(kB γ∆eff (∆eff /ωc )2s−2 ), respectively, where ∆eff is the effective tunneling amplitude defined in Eq. (3.12). Note that there is no quantum transition forthe super-Ohmic case (s > 1); ∆eff is finite for arbitrary system-reservoir couplings.
At low temperatures (kB T ≪ ∆eff ), the numerical results agree with the co-tunnelingformula (indicated by blue dashed lines) and show T 2s+1 -dependence, regardlessof the strength of the system-reservoir coupling. As shown in Fig. 5.5 (a), the numerical results for α = 0.1 agree with the sequential tunneling formula at moderatetemperatures (kB T ∼ ∆eff ) and with the NIBA at high temperatures. However, from
Fig. 5.5 (b), it is evident that the numerical results for α = 0.5 agree better with the
NIBA than with the sequential tunneling formula at moderate-to-high temperatures(kB T ≳ ∆eff ). This change can be explained by the crossover temperature T ∗ , which
separates the sequential (T < T ∗ ) and incoherent (T > T ∗ ) tunneling regimes (seeEq. (3.38)). As the system-reservoir coupling α increases, the temperature region
for which the numerical results agree with the NIBA is widened since the crossover
temperature T ∗ is lowered.
The Markov approximation for incoherent tunneling is indicated by orange dotted lines in Fig. 5.5. The incoherent tunneling formula clearly deviates from numerical results, indicating the importance of the non-Markovian properties of the system
dynamics. The origin of this disagreement is the same as that for the sub-Ohmic case(see Section 5.1.2).

5.1.4

The super-Ohmic case (2 ≤ s)

In Fig. 5.6, we show the numerical results of the thermal conductance obtained using
the CTQMC method as a function of the temperature for s = 2.0. The normalizationof the horizontal and vertical axes as well as the line types of the analytical formula are same as those in Fig. 5.5. At low temperatures, the numerical results agree

44

Chapter 5. Numerical Results and Comparison with Analytical Formulas

10
5

1

0.5

1

F IGURE 5.6: Temperature dependence of the thermal conductancecalculated using CTQMC simulation. The data represent the result fors = 2.0, ∆/ωc = 0.05, and α = 0.5. The linetypes are same as those in
Fig. 5.5. The inset shows an enlarged graph in the high-temperature
region.

well with the co-tunneling formula and show T 2s+1 -dependence, regardless of thestrength of the system-reservoir coupling. In contrast to the case of 1 < s < 2, thenumerical results agree with the sequential tunneling formula at moderate-to-hightemperatures. This is reasonable since the crossover T ∗ becomes of the order of ωcfor s = 2.0.

5.2

Quantum critical phenomena

For the sub-Ohmic case (0 < s < 1), a quantum phase transition occurs at zerotemperature when the reservoir-system coupling reaches the critical value αc , where
αc is a function of s and ∆/ωc as described in Section 3.1. In this section, we discuss
quantum critical phenomena near the transition point.

5.2.1

Phase diagram

For α < αc , the ground state is described by a coherent superposition of two wave
functions localized at each well (delocalized phase). For α > αc , the ground state
becomes two-fold degenerate because the coherent superposition is completely broken owing to the disappearance of quantum tunneling between the two wells (localized phase). The phase diagram of the spin-boson model determined by the
CTQMC simulations for ∆/ωc = 0.1 is shown in Fig. 5.7 (for details on determining the critical value, αc , see Section 4.4). The transition separating the two phases
is of second-order for the sub-Ohmic case (the empty squares) or of the KosterlitzThouless-type [49, 50] for the Ohmic case (the filled circle). This phase diagram isconsistent with previous numerical studies [37, 35].

5.2.2

Thermal conductanceIn Fig. 5.8, we show the temperature dependence of the thermal conductance fors = 0.5 and ∆/ωc = 0.1, where the critical system-reservoir coupling is αc = 0.1074.
Figs. 5.8 (a) and (b) show the delocalized-phase side (α ≤ αc ) and the localized-phaseside (α ≥ αc ), respectively. At the critical point, the thermal conductance exhibits
distinctive power-law behavior κ ∼ T c , where c is the critical exponent, determined5.2. Quantum critical phenomena

45

F IGURE 5.7: The phase diagram of the sub-Ohmic case for ∆/ωc =0.1. The solid line indicates the second-order transition line separating the delocalized and localized phases. The empty squares indicate
the critical system-reservoir coupling that is numerically determined
for the sub-Ohmic case (0 < s < 1), whereas the filled circle represents the known transition point αc = 1 for the Ohmic case (s = 1).

F IGURE 5.8: Temperature dependence of the thermal conductance for(a) α ≤ αc and (b) α ≥ αc . The plots represent the CTQMC simulation
results for s = 0.5 and ∆/ωc = 0.1, for which the critical systemreservoir strength is αc = 0.1074.

by the nature of the quantum phase transition in Section 3.6. From Eqs. (3.50) and(3.59), the critical exponent c = 1 is obtained for s = 0.5. This conclusion is consistent
with the critical exponent c = 1 obtained by the CTQMC simulation for s = 0.5 (seeFig. 5.8).
As the system-reservoir coupling is reduced below the critical value (α < αc ),the temperature dependence of the thermal conductance deviates from one at thecritical point. For a sufficiently small system-reservoir coupling (e.g., α = 0.07 in
Fig. 5.8 (a)), the thermal conductance becomes proportional to T 2s+1 at low temperature, presumably for heat transport due to co-tunneling (3.41). The temperaturedependence of the thermal conductance also deviates as the system-reservoir coupling is increased above the critical value (α > αc ). Its temperature dependencecannot be explained by a simple formula such as the NIBA up to α = 0.13.

47

Chapter 6

Experimental Realization of
Sub-Ohmic Reservoirs
In this chapter, we discuss an experimental realization of sub-Ohmic reservoirs using a superconducting circuit. A previous theoretical study [43] has shown that aspatially-uniform transmission line can realize a sub-Ohmic reservoir with s = 0.5.
For a controlled experiment of the quantum phase transition, however, it is favorable
to realize a sub-Ohmic reservoir with an arbitrary value of s. We propose a superconducting circuit to realize sub-Ohmic reservoirs for arbitrary s by introducing spatial
dependence to the circuit elements.

6.1

Superconducting circuit model

In this section, we show that a qubit coupled to two transmission lines reduces tothe spin-boson model (see Fig. 6.1 (a)). We consider two kinds of qubits, a flux qubit(Section 6.1.2) and a charge qubit (Section 6.1.3). For simplicity, we first consider
uniform LC transmission lines and confirm that the qubits coupled to the transmission lines reduce to the Ohmic spin-boson model. To treat arbitrary transmission
lines, we derive general linear-response relations between the spectral density andan impedance of the transmission lines for both types of qubits.

6.1.1

Transmission line

First, we explain that the transmission line, which consists of resistances, inductances, and capacitances (see Fig. 6.1 (b)), is modeled as a set of harmonic oscillators.
For simplicity, we consider a uniform transmission line made from constant capacitances C and inductances L without resistances. It is described by the Hamiltonian
N

H=∑

j =1

Q2j
2C

( ϕ j − ϕ j −1 )2
.
2L
j =2
N

+∑

(6.1)

Here we introduced a flux operator ϕj (t) by dϕj (t)/dt ≡ Vj (t), where Vj is a voltage drop at the jth capacitance. The charge operator Q j is the canonical conjugate
variable to the flux operator and these operators satisfy the canonical commutation
relation: [ϕj , Qk ] = iδj,k . In the continuous limit, N → ∞, the Hamiltonian (6.1) isapproximated as follows:
"

 #
Z d
Q ( x )2
1 ∂ϕ( x ) 2
H=
dx
+
,
(6.2)
2c
2l
∂x
0

Chapter 6. Experimental Realization of Sub-Ohmic Reservoirs=

48

F IGURE 6.1: (a) A superconducting circuit composed of (a-i) flux or(a-ii) charge qubits and two transmission lines. (b) The circuit of
the transmission lines proposed to realize the sub-Ohmic spin-bosonmodel, consisting of resistances R j , inductances L j , and capacitances
Cj .

where d, c, and l are the length of the transmission line, the capacitance and inductance per length, respectively. Expanding ϕ( x ) in the modes {ϕk }ϕ( x ) = ∑ ϕk f k ( x ),

(6.3)

k

where f k ( x ) satisfies
Z d0dx f k ( x ) f k′ ( x ) = dδk,k′ ,

Z d
0

dx

∂ f k ( x ) ∂ f k′ ( x )= dk2 δk,k′ ,
∂x
∂x

(6.4)

the Hamiltonian is rewritten follows in terms of ϕk and its canonical conjugate variable Qk as!
1 2 dcωk2 2
H=∑
Q +
ϕk ,
(6.5)
2dc k
2
k

6.1. Superconducting circuit model

49√
where ωk = k/ lc. Since ϕk and Qk are expressed by ϕk = αk ak + α∗k a†k and Qk =β k ak + β∗k a†k , where ak (a†k ) is an annihilation (a creation) operator of boson, respectively, the effective Hamiltonian is written as
!
"
!
#

2
2
β
dcω
1
| β k |2 dcωk2 |αk |2kk 2a†k ak ++∑H̃ = ∑++
αk a2k + h.c. . (6.6)
dc
2
2
2dc
2
k
k
Choosing appropriate αk and β k under the condition αk β∗k − α∗k β k = i from the canonical commutation relation, the effective Hamiltonian (6.6) reduces to a set of the harmonic oscillators with the frequency ωk .

6.1.2

Flux qubit

We consider a flux qubit coupled to two transmission lines (or two junction arrays),as shown in Fig. 6.1 (a-i). We first consider a uniform transmission line with constantcapacitance and inductance (Cj = C, L j = L) while neglecting resistance (R j = 0).
The Hamiltonian of the present circuit is given by
H = HS + HB + HI ,
" 2
#
3
Q J,k
HS = ∑− E J,k cos(ϕ J,k /ϕ0 ) ,2C J,k
k =1
" 2
#
N
Qν,j
(ϕν,j+1 − ϕν,j )2
HB = ∑ ∑+,2C
2L
ν j =1HI =(ϕa − ϕL,1 )2 (ϕR,1 − ϕb )2
+
,
2L1
2L1

(6.7)
(6.8)
(6.9)
(6.10)

where HS , HB (= ∑ν HB,ν ), and HI (= ∑ν HI,ν ) describe the flux qubit, the transmission lines, and the system-reservoir coupling, respectively, and ϕ0 = h̄/2e is theflux quantum. The flux qubit comprises three Josephson junctions with Josephsonenergies E J,k (k = 1, 2, 3), and the charge and flux operator of the k-th Josephson
junction are denoted by Q J,k and ϕ J,k , respectively. Similarly, the charge and flux
operators of the transmission line (see Fig. 6.1 (b)) are denoted by Qν,j and ϕν,k , respectively. These operators satisfy the exchange relations [ϕ J,k , Q J,k′ ] = iδk,k′ and[ϕν,j , Qν′ ,j′ ] = iδj,j′ δν,ν′ , respectively. The flux operators at the two sides of the flux
qubit are expressed by ϕa and ϕb (see Fig. 6.1 (a-i)).
To make the flux qubit, the area of one junction is reduced by a factor of α (E J,1 =E J,3 = E J , C J,1 = C J,3 = C J , E J,2 = αE J , and C J,2 = α−1 C J ). Then, the Hamiltonian of
the flux qubit Hamiltonian (6.8) can be rewritten as [17, 32]
Hqb =

Q2J,+

+

Q2J,−

+ V (ϕ J,+ , ϕ J,− ),2C J,+ 2C J,−
V (ϕ J,+ , ϕ J,− ) = − E J [2 cos(ϕ J,+ /2ϕ0 ) cos(ϕ J,− /2ϕ0 )(6.11)+ α cos ((Φext − ϕ J,− )/2ϕ0 )] ,(6.12)

where ϕ J,± = (ϕ J,1 ± ϕ J,3 )/2, its conjugate operator is denoted by Q J,± , and V (ϕ J,+ , ϕ J,− )is the Josephson energy that plays the role of the potential energy. When the magnetic flux through the loop is tuned to be half of the flux quantum (Φext = ϕ0 /2),
the Josephson energy, V (ϕ J,+ , ϕ J,− ), has two energy minima on the line ϕ J,+ = 0.
Due to quantum tunneling effects, there is an energy splitting ∆ between the ground50

Chapter 6. Experimental Realization of Sub-Ohmic Reservoirs

state and the first-excited state. Since these lowest two eigenstates are well separated from the other eigenstates, we can truncate the system into the lowest two
eigenstates, thus leading to the two-state system Hamiltonian (6.8). The √wavefunctions of the lowest two states are described as |σx = ±1⟩ = (|↑⟩ ± |↓⟩)/ 2, where|↑⟩ and |↓⟩ are the two-dimensional wavefunctions localized at the two potential
energy minima, respectively. Introducing the new variables ϕ± = ϕR,1 ± ϕL,1 andΦ± = ϕb ± ϕa and using ϕ J,+ ∝ Φ+ ≃ 0, the system-reservoir coupling (6.10) is
rewritten as HI = −ϕ− Φ− /2L1 . After truncation into the two-state system, we obtain:
HI = −

ϕ−
ϕ0 ⟨ φ− ⟩ σz ,2L1

(6.13)

where ⟨↑| Φ− |↑⟩ ≡ ϕ0 ⟨ φ− ⟩, ⟨↓| Φ− |↓⟩ ≡ −ϕ0 ⟨ φ− ⟩, and ⟨↑| Φ− |↓⟩ = ⟨↓| Φ− |↑⟩ =0.
For simplicity, we consider the continuous limit ∆x → 0 while keeping the lengthof the transmission line, d = N∆x, constant, where ∆x is the size of each elementaryisland. Then, the system-reservoir coupling can be rewritten by [17]HI = −1 ∂ϕ( x )ϕ0 ⟨ φ− ⟩ σz .
l ∂x x=0(6.14)The transmission line Hamiltonian (6.9) can be The flux, ϕ( x ), can be expressed by


1
eikx
ϕ( x ) = ∑ √ √bk + bk† ,Lt 2cωkk(6.15)

where bk and bk† are bosonic annihilation and creation operators, respectively. Then,the Hamiltonians for the transmission lines and the system-reservoir coupling canbe rewritten as follows:∑ ωk bk† bk ,(6.16)

HI = −(6.17)

λk =



σz
λk bk† + bk ,∑2 kr
2ϕ0 ⟨ φ− ⟩ ωk
√
,
c
vl Lt

(6.18)HB =k√where v = 1/ lc is the speed of light in the transmission line. This model corresponds to the spin-boson model with an Ohmic reservoir.
Now, we discuss the general linear response relation. The electric current operator at the position x is defined by I( x ) = l −1 ∂ϕ( x )/∂x and is calculated at x = 0:iλk  †
bk + bk .
2ϕ0 ⟨ φk ⟩kI0 ≡ I( x = 0) = ∑(6.19)

From Eqs. (2.17), (6.17), (6.18), and (6.19), the spectral density is written as
I (ω ) = −hi
4ϕ02 ⟨ φ− ⟩2
Im GIR0 (ω ) ,π(6.20)

where GIR0 (ω ) is the Fourier transformation of the current-current correlation function defined by GIR0 (t) = −iθ (t) ⟨[I0 (t), I0 (0)]⟩.

6.1. Superconducting circuit model51

A current I0 (ω ) is induced by applying an electric field E(r, ω ) which can betreated as a perturbation
Hext (ω ) = −i
ωZdr i (r ) · E(r, ω ),(6.21)

where i (r, ω ) is the current density operator. Using the linear response theory, thecurrent can be expressed as a linear response to the perturbation as follows [64]:I0 ( ω ) =Z

dσr er · i (r, ω ) =Zdσr dσr′ dr ′ er ·i R′G′ ( ω ) er ′ E (r , ω ),
ω i (r )i (r )(6.22)

where GiR(r )i(r ′ ) (t − t′ ) = −iθ (t − t′ ) ⟨[i (r, t), i (r ′ , t′ )]⟩. Note that we have employed acoordinate system (r, σr ), where r is a coordinate parallel to the electric field, E(r, ω ) =er E(r, ω ), and σr are coordinates on the plane perpendicular to the r-direction. Since
the current is not dependent on r due to the current conservation, the linear response
can be rewritten by

I0 ( ω ) =

Z

dr ′

i R
G ( ω ) E (r ′ , ω ),
ω I0

(6.23)In addition, the current through a circuit with an impedance Z (ω ) can be also written as

I0 ( ω ) = G ( ω )V ( ω ) =Z

dr G (ω ) E(r, ω ).

(6.24)

Comparing Eq. (6.23) and Eq. (6.24), we obtain the linear response relation between
the current-current correlation function and the total impedance of the transmission
lines:
i1= GIR0 (ω ).
Z (ω )
ω(6.25)

Substituting Eq. (6.25) into (6.20), the spectral density can be expressed by the joint
impedance of the two transmission lines (Zs (ω ) = ∑ν Zν (ω )) as follows:
4ϕ02 ⟨ φ− ⟩2I (ω ) = ∑ Iν (ω ) =I0 (ω ),πνhi
I0 (ω ) = ωRe Zs (ω )−1 .

(6.26)(6.27)

Although we have derived them for a special case, i.e., the case of uniform transmission lines without damping, Eqs. (6.26) and (6.27) hold for arbitrary circuits of thetransmission lines because we expect that the coupling form I0 Φ− does not change
and the details of the transmission lines are included only in its impedance Zν (ω ).

6.1.3

Charge qubit

We consider a charge qubit coupled to two transmission lines as shown in Fig. 6.1 (aii). As in the case of the flux qubit (Section 6.1.2), we first consider a uniform transmission line with constant capacitance and inductance.

52

Chapter 6. Experimental Realization of Sub-Ohmic Reservoirs

F IGURE 6.2: Eigeneneries of the charge qubit Hamiltonian as a function of the number of the offset charge n g . At the charge sweet spots(the vertical dashed lines), the two lowest energy states are separated
from others and the energy splitting between them is E J .

The Hamiltonian of the circuit composed of the charge qubit and transmission
lines is given by
H = HS + HB + HI ,2

HS = EC (n − n g ) − E J cos(ϕ/ϕ0 ),
CJ
QV0 ,
HI =
C J + Cg

(6.28)(6.29)(6.30)

where HS and HI (= ∑ν HI,ν ) describe the charge qubit and the system-reservoir
coupling, respectively, and HB (= ∑ν HB,ν ) is the Hamiltonian of the transmission
lines given in Eq. (6.9) (or Eq. (6.16)). The charge qubit comprises a Josephson junction with Josephson energies E J and the flux drop ϕ, and has the Coulomb energy
EC = 2e2 /(C J + Cg ) and the offset charge 2en g . The operator n is the number of
Cooper pairs transferred through the Josephson junction (Q = 2en). The number
operator and flux operator satisfy the exchange relation [ϕ, n] = iϕ0 . Considering asmall Josephson junction with E J ≪ EC and setting the number of the offset charge
as n g = 0.5 (see Fig. 6.2), two lowest charge states n = 0 and n = 1 are well separated from others, and this system is reduced to the two-state system (2.8) withenergy splitting ∆ ∼ E J . The wavefunctionsof the lowest two states are described√
as |σx = ±1⟩ = (|n = 0⟩ ± |n = 1⟩)/ 2. After truncation into the two-state system,
the system-reservoir coupling Hamiltonian is rewritten as
HI = −σz C J
2|e|V0 .
2 C J + Cg

(6.31)

Considering the continuous limit, the voltage operator V0 is given by
V0 ≡∂ϕ( x = 0, t)
.
∂t
t =0

(6.32)

6.2. The spectral densities

53

Using the expression of the flux in the transmission line:


i
1
†
√
√
b
(
t
)
−
b
(
t
)
∑ L 2cω k
k
t
k
k


1
i
= ∑√ √
bk e−iωk t − bk† eiωk t ,Lt 2cωkk

ϕ( x = 0, t) =the voltage operator can be rewritten asrωk 
V0 = ∑
bk + bk† .
2Lt ck(6.33)(6.34)

Then, the Hamiltonian for the system-reservoir coupling can be rewritten asσz(6.35)
HI = − ∑ λk bk + bk† ,
2 k
r
CJ
ωk
2| e |.
(6.36)λk =C J + Cg
2Lt cThis model corresponds to the Ohmic spin-boson model.
Now, we consider the voltage-voltage correlation function defined by GVR0 (t) =−iθ (t) ⟨[V0 (t), V0 (0)]⟩ and obtain the relation between the spectral density and theimpedance of the transmission lines in the same way as the flux qubit case. From
Eq. (6.34), the Fourier transformation of the voltage-voltage correlation function isrelated with the spectral density as follows:4| e |2
I (ω ) = −πCJC J + Cg2hi
Im GVR0 (ω ) .

(6.37)

Using the linear response relation (6.25) and GVR0 (ω ) = Z (ω )2 GIR0 (ω ), we obtain the
spectral density in terms of the impedance
4| e |2∑ Iν (ω ) = πνI0 (ω ) = ωRe Z p (ω ) ,I (ω ) =CJC J + Cg2
I0 (ω ),

(6.38)(6.39)

where Z p (ω )−1 = ∑ν Zν (ω )−1 .

6.2

The spectral densities

To realize a sub-Ohmic reservoir with an arbitrary exponent, s, we propose a superconducting circuit, as shown in Fig. 6.1 (b). The circuit comprises resistances R j ,inductances L j , and capacitances Cj (j = 1, 2, · · · , N). For simplicity, we assume thatthe two transmission lines are constructed by the same circuit. The joint impedanceof the two transmission lines is then calculated as Zs (ω ) = 2Z1 (ω ) for the case of theflux qubit and Z p (ω )−1 = 2Z1 (ω )−1 for the case of the charge qubit, where Zj (ω )(j = 1, 2, · · · , N) includes from jth to Nth islands and is given by a recurrence
relation:
Zj (ω ) = R j + iωL j +1
Z j +1(ω )−1 + iωC,j(6.40)

54

Chapter 6. Experimental Realization of Sub-Ohmic Reservoirs

F IGURE 6.3: The spectral density of the superconducting circuit fors = 0.5 and 0.25, corresponding to (n, p, m) = (2, 0, 2) and (6, 0, 6). The
circuit parameters are set to N = 104 , R0 = 1 kΩ, L0 = 13 nH, andC0 = 1 pF. Each circuit element has 1% relative randomness.

with ZN +1 (ω )−1 = 0. In the following, we assume that circuit elements have spatial
dependence:
 n
 p
 m
j
j
j
R j = R0
, L j = L0
, Cj = C0.
(6.41)N
N
N
where n, p, and m are non-negative real numbers.

6.2.1

Flux qubit

We assume that the exponents of the circuit elements (see Eq. (6.41)) are set asn ≥ 0,

p = 0,m ≥ 0.

(6.42)We show the spectral density, I0 (ω ) = ωRe[ Zs (ω )−1 ] (see Eq. (6.27)), of the superconducting circuit in Fig. 6.3 for (n, p, m) = (2, 0, 2) and (6, 0, 6). The parameters are
set to R0 = 1 kΩ, L0 = 13 nH, C0 = 1 pF, and N = 104 and referred to experimental
studies on Josephson junction arrays [65]. In Fig. 6.3, we added 1% relative randomness for each circuit element to introduce tolerance to circuit parameter fluctuations.
We determined that the spectral density is approximately proportional to ω s ina certain range of the frequency with the exponent 0 < s < 1. This indicates thatthe present circuit can realize a sub-Ohmic reservoir with an arbitrary value of s.
Certainly, the analytical calculation concludes:
I (ω ) ∝ ω 2/(m+2) ,( ω ∗ ≪ ω ≪ ωc ) .

(6.43)The detailed calculation is given in Appendix D.1. This result is in good agreementwith Fig. 6.3; m = 2 and 6 correspond to s = 0.5 and 0.25, respectively. The lower
frequency limit for the sub-Ohmic spectral density, ω ∗ , is calculated as follows:
"
∗

ω =

 m 2nR0m+2

2N

C0n L0m+n+2#1/(m+2n+2).

(6.44)

6.2. The spectral densities

55

F IGURE 6.4: The spectral density of the superconducting circuit fors = 0.5 and s = 0.25, corresponding to (n, p, m) = (4, 2, 0) and (12, 6, 0).
The circuit parameters are set to N = 104 , R0 = 0.1 kΩ, L0 = 13 nH,
and C0 = 0.01 pF. Each circuit element has 1% relative randomness.

Therefore, the exponent n for the resistance (see Eq. (6.41)) controls the lower limitof the sub-Ohmic spectral density. In contrast, the higher frequency limit, ωc , is acomplex function of the circuit parameters.

6.2.2

Charge qubit

We set the exponents of the circuit elements asn ≥ p ≥ 0,

m = 0.

(6.45)In Fig. 6.4, we show the spectral density, I0 (ω ) = ωRe[ Z p (ω )] (see Eq. (6.39)), of
the superconducting circuit for (n, p, m) = (4, 2, 0) and (12, 6, 0) with R0 = 0.1 kΩ,
L0 = 13 nH, C0 = 0.01 pF, and N = 104 , and added 1% relative randomness for eachcircuit element to introduce tolerance to circuit parameter fluctuations.
From Fig. 6.4, the transmission line realizes a sub-Ohmic reservoir with arbitrary
exponent 0 < s < 1 in a certain range of the frequency. Calculating the transmission line with spatial dependent elements (6.45) analytically (for detail, see Appendix D.2), we obtain the frequency dependence asI (ω ) ∝ ω 2/( p+2) ,( ω ∗ ≪ ω ≪ ωc ) .

(6.46)

This analytical result is consistent with Fig. 6.4; p = 2 and 6 correspond to s = 0.5
and 0.25, respectively. The lower frequency limit, ω ∗ , is calculated as follows:
"
∗

ω =

p +2

 p  2( n − p )2NR0n − p n +2
L0C0

#1/(2n− p+2).

(6.47)The lower frequency limit is controlled by the exponent n for the resistance. As in
the case of the flux qubit, the higher frequency limit, ωc , is a complex function of thecircuit parameters.

56

Chapter 6. Experimental Realization of Sub-Ohmic Reservoirs

F IGURE 6.5: The spectral density function of the superconducting circuit for s = 0.25 corresponding to (n, p, m) = (2, 0, 0). The circuit parameters are set as N = 102 , R0 = 0.1 mΩ, and C0 = 0.1 mF. The redand blue lines represent cases that each circuit element has 0% and
10% relative randomnesses, respectively.

6.2.3

Charge qubit (0 < s ≤ 0.5)

There is another way to realize a sub-Ohmic reservoir with the restricted exponent
0 < s ≤ 0.5 using the charge qubit. To realize the transmission line, we set theexponent of circuit elements asn ≥ 0,

p = 0,

m = 0.

(6.48)In Fig. 6.5, we show the spectral density, I0 (ω ), of the transmission line for (n, p,m) = (2, 0, 0) with 0% and 10% relative randomness for each circuit element. The
parameters are set to N = 102 , R0 = 0.1 mΩ, and C0 = 0.1 mF.
As shown in Fig. 6.5, the present circuit realizes the sub-Ohmic reservoir in acertain range of the frequency. Compared with the above two cases (Sections 6.2.1
and 6.2.2), the present circuit (6.48) can realize a sub-Ohmic reservoir with a smaller
number of islands, N, and has a higher tolerance for randomness. By the analytical
calculation in Appendix D.3, we obtain the frequency dependence of the spectral
density:
I (ω ) ∝ ω 1/(n+2) ,( ω ∗ ≪ ω ≪ ωc ) .

(6.49)

This result is in good agreement with Fig. 6.5; n = 2 corresponds to 0.25. However,from analytical calculation (6.49), the exponent of the spectral density is restricted to
0 < s ≤ 0.5. In this case, the lower frequency limit, ω ∗ , and the higher frequency
limit, ωc , can be calculated as follows:
ω∗ =
ωc =

√ ! n +2
1  n 2
2 2
1+,
R0 C0 2N
n


1
2N n
.
R0 C0
n

(6.50)(6.51)

57

Chapter 7

Summary
We systematically studied heat transport via a local two-state system for all types ofreservoirs, i.e., for the Ohmic case (s = 1), sub-Ohmic case (s < 1), and super-Ohmiccase (s > 1) near equilibrium. We used the exact expression for the thermal conductance obtained from the Keldysh formalism and studied it using both analytical andnumerical methods within the linear response.
We showed that heat transport via a two-state system can be described by three
types of transport processes: sequential tunneling, co-tunneling, and incoherent tunneling. Sequential tunneling is a transport process occurred by a combination ofphoton (or phonon) absorption and emission, and the temperature dependence ofthe thermal conductance is the exponential decrease as the temperature is lowered.
Co-tunneling is a transport process via a virtual excitation in the local two-state system, and the thermal conductance behaves as κ ∝ T 2s+1 at low temperatures. Incoherent tunneling is a transport process induced by a stochastic transition between
two localized wavefunctions. In this thesis, we gave the condition for the sequential
tunneling using the quantum master equation and showed that the sub-Ohmic casedoes not satisfy this condition. We also pointed out that the Markov approximation
yielded incorrect results for the thermal conductance in the incoherent tunnelingregime since the non-Markovian properties are important. However, for the incoherent tunneling regime, the NIBA yielded correct results.
We used a continuous-time Monte Carlo algorithm and systematically compared
the numerical results with those of the analytical approximation formulas. We foundthat all numerical results were well reproduced by one of three formulas, i.e., thesequential tunneling formula, co-tunneling formula, or NIBA. We summarize therelevant transport mechanisms for each regime in Table 7.1. We also showed thatfor 0 < s ≤ 1, the quantum phase transition between the delocalized and localized
phases strongly affected the temperature dependence of the thermal conductance.
For the delocalized phase (α < αc ), the thermal conductance is well described by theco-tunneling formula at low temperatures and by NIBA at moderate-to-high temperatures. On the contrary, for the localized phase (α ≫ αc ), the NIBA holds arbitrary
temperatures.
Furthermore, we studied quantum critical phenomena (QCP) in the sub-Ohmiccase from the viewpoint of heat transport. By the continuous-time quantum Monte
Carlo simulations, we showed that the thermal conductance at the critical point(α = αc ) has a characteristic power-law temperature dependence determined by
the nature of the quantum phase transition. We also clarified the means by which
the critical exponent of the thermal conductance is related to other critical exponents
discussed in previous theoretical studies.
Finally, we proposed a superconducting circuit with the flux or charge qubits thatrealize the sub-Ohmic spin-boson model. Giving the spatial dependence to circuit elements in transmission lines (or Josephson junction arrays), the sub-Ohmic spectral

58

Chapter 7. Summary
TABLE 7.1: Summary of the relevant transport process. The lastcolumn shows the temperature dependences of the thermal conductance. The temperature dependence of the NIBA is complex in general, and the symbol (∗ ) indicates the high-temperature limit. Near
the quantum critical point for the sub-Ohmic case, the thermal conductance shows a quantum critical phenomena (QCP). The criticalexponent c is a function of s (see Eq. (3.59)).

Exponent

Condition

Transport Process

Dependence

0<s<1(sub-Ohmic)

α < αc , kB T ≪ ∆eff
α < αc , kB T ≳ ∆eff
α > αc , arbitrary temp.
α = αc , kB T ≪ ∆eff
α < 1, kB T ≪ ∆eff
α < 1, kB T ≳ ∆eff
α > 1, arbitrary temp.
kB T ≪ ∆eff
∆eff ≲ kB T ≲ kB T ∗kB T > kB T ∗kB T ≪ ∆eff
kB T ≳ ∆eff

Co-tunneling
Incoherent tun.
Incoherent tun.
QCP
Co-tunneling
Incoherent tun.
Incoherent tun.
Co-tunneling
Sequential tun.
Incoherent tun.
Co-tunneling
Sequential tun.

T 2s+1s=1(Ohmic)1<s<2(super-Ohmic)
2≤s(super-Ohmic)

Tc
T3
T 2α−1 (∗ )T 2α−1 (∗ )T 2s+1e−∆eff /kB T /T 2
T 2s+1e−∆eff /kB T /T 2

density with an arbitrary value of the exponent s is realized in a certain frequency
range.
The present study provides a theoretical basis for describing heat transport via
nano-scale objects. The present study also provides a new platform for experiments
attempting to access quantum phase transitions directly upon measuring heat transport in mesoscopic devices.

59

Appendix A

Instanton Method
In this appendix, we derive the tunneling amplitude for the double-well potential
systems, Eq. (2.7), using an instanton method. In the imaginary-time path-integral
formalism, the transition amplitude from q(−τ0 /2) = qi to q(τ0 /2) = q f is writtenas⟨q f | e− HS τ0 |qi ⟩ = N

Z

where S(E) [q(τ )] is the Euclidean action:
(
Z
S

(E)

[q(τ )] =+τ0 /2−τ0 /2dτ

D q ( τ ) e−S(E) [ q ( τ )]

,

(A.1)

)


M q(τ ) 2
+ V [q(τ )] ,
2
dτ

(A.2)

R
where N is a normalization factor, and D q(τ ) denotes integration over all pathq(τ ). Now suppose the semi-classical approximation, in which the action takes large
values than the Plank constant (note that h̄ is set as unity in this thesis). In thissituation, the integral (A.1) is dominantly determined by the stationary solution andfluctuation of the paths around it. The stationary solution X (τ ) is obtained from thestationary condition δS(E) [ X (τ )] = 0 as
−M

d2 X (τ ) dV [ X (τ )]
+
= 0.
dτ 2
dX (τ )

(A.3)

This equation corresponds to the classical equation of motion of particle under theupside-down potential, −V (q). If there is one stationary point, we can estimate theexponential factor of the transition amplitude as follows:⟨q f | e− HS τ0| qi ⟩ = NZ

D q ( τ ) e−S(E) [ q ( τ )](E)

∼ e − S0 ,

(A.4)

(E)

where S0 = S(E) [ X (τ )] is the action of the stationary solution.
Next, we consider the pre-exponential factor of the transition amplitude. The
path q(τ ) of the integral in Eq. (A.1) is expanded around the stationary solution asq ( τ ) = X ( τ ) + ∑ c n x n ( τ ),

(A.5)

n

where xn (τ ) is a complete set of orthonormal functions that vanish at the boundary:Z +τ0 /2
−τ0 /2dτ xn (τ ) xm (τ ) = δm,n , xn (±τ0 /2) = 0.

(A.6)

60

Appendix A. Instanton Method

Thus, the integral measure can be written in the formZ

D q(τ ) = ∏

Z

n

dc
√ n .
2π

(A.7)

Expanding the Euclidean action up to the second order around the stationary point"#
S(E) X (τ ) + ∑ cn xn (τ ) − S0 [ X (τ )]
(E)

n

=

Z +τ0 /2
−τ0 /2dτ∑ cn xn (τ )!n

1 d2 V [ X (τ )]M d2−+
2 dτ 2 2 dX (τ )2∑ cn xn (τ )

!
.

n

(A.8)
We note that the first-order terms with respect to xn (τ ) vanish as X (τ ) is the stationary solution. Here, we assume that { xn (τ )} are eigenfunctions of the differential
operator − Md2 /dτ 2 + V ′′ [ X (τ )]:
d2 V [ X (τ )]d2
x n ( τ ) = ϵn x n ( τ ),−M 2 +dτdX (τ )2

(A.9)

where {ϵn } is the corresponding eigenvalues. Then, the Euclidean action is transformed into the diagonal form
"
#
1
(E)
(E)
S
X (τ ) + ∑ cn xn (τ ) = S0 + ∑ ϵn c2n .
(A.10)
2 nn
Performing Gaussian integration, we obtain the transition amplitude(E)⟨q f | e− HS τ0 |qi ⟩ = e−S0 N ∏ ϵn−1/2(A.11)

n

= e

(E)

− S0

d2
d2 V [ X (τ )]N det − M 2 +dτdX (τ )2

−1/2.

(A.12)In the following, we calculate the exponential and pre-exponential factors concretely
in the double-well potential (2.5).

A.1

Exponential factor

Let us first consider the stationary solution for the double-well potential system. The
classical equation of motion (A.3) has not only trivial solutions
X (τ ) = ±q0,
2(A.13)

but also nontrivial solutions which connect the points ±q0 /2. These nontrivial solutions are crucial to calculate the quantum tunneling amplitude. In the limit τ0 → ∞,these nontrivial solutions of Eq. (A.3) areq0
ω0 (τ − τc )X± (τ ) = ± tanh,
(A.14)
22

A.2. Pre-exponential factor

61

F IGURE A.1: Schematic of the instanton and anti-instanton. The instanton is the imaginary time path from the point −q0 /2 to +q0 /2 inthe upside-down potential V [q(τ )] = −V [q(t)]. The anti-instanton is
the imaginary time path opposite to the instanton.

where τc is the arbitrary parameter. These solutions, X+ (τ ) and X− (τ ) , are called‘instanton’ and ‘anti-instanton’, respectively (see Fig. A.1). Using the instanton or(E)
anti-instanton solutions, we obtain the action S0(E)

S0= S(E) [ X± (τ )] =

16Vb
.
3ω0

(A.15)

This exponential factor is consistent with one derived from the WKB approximation [66]. Note that τc may take an arbitrary time, and the action of the instanton,(E)
S0 , does not depend on the τc , reflecting the invariance with respect to time translation. Each concrete solution (A.14) has a definite position with respect to the origin,and thus there exits an infinite family of solutions distributed along the imaginarytime axis.

A.2

Pre-exponential factor

First, we consider the single-instanton contribution to the transition amplitude, Eq. (A.1),
and calculate the pre-exponential factor of it. The transition amplitude (A.12) isrewritten as−1/2(E)
d2− HS τ0− S0
2⟨±q0 /2| e|∓q0 /2⟩single = Ne
det − M 2 + ω0dτ( )−1/2
det − Md2 /dτ 2 + V ′′ [ X (τ )]×(A.16).
det − Md2 /dτ 2 + ω02
We multiplied and divided by the determinant for harmonic oscillator, which is calculate easily as [66]s−1/2d2
Mω0
2
N det − M 2 + ω0=
(A.17)dτ
2π sinh(ω0 τ0 )r
Mω0 −ω0 τ0 /2∼e, τ0 → ∞.
(A.18)π62

Appendix A. Instanton Method

Substituting the explicit expression X+ (τ ) (A.14) into V ′′ [ X (τ )], we obtain the eigenvalue equation for determinant in the numerator:
"
(
#)
3
1
d2
2
− 2 + ω0 1 −
xn (τ ) = ϵxn (τ )(A.19)dτ
2 cosh2 (ω0 (τ − τc )/2)
It can be regarded as a certain Schorödinger equation. Under the boundary condition xn (±τ0 /2) = 0, τ0 → ∞, Eq. (A.19) has discrete spectrum for ϵ < ω0 andcontinuous spectrum for ϵ > ω0 [67]. There are two discrete levels; one has the
eigenvalue ϵ1 = (3/4)ω02 , and the other ϵ0 = 0. The eigenfunction of the latteris [67]r
3ω0
1x0 ( τ ) =.
(A.20)
8 cosh2 (ω0 (τ − τc )/2)
The naive integration with respect to the coefficient c0 (A.7) (the zero mode) gives
divergence since the action is independent of c0 . However, we can overcome thisdivergence by changing the integration over c0 into τc . Note that the eigenfunction
x0 (τ ) is proportional to dX (τ )/dτ. This means that the change of c0 can be absorbed
into the translational shift of the instanton solution along the imaginary-time axis.
When c0 is changed into c0 + ∆c0 , the change of the path, ∆q(τ ), is given as ∆q(τ ) =x0 (τ )∆c0 (see Eq. (A.5)). On the other hand, the change ∆q(τ ) under a shift ∆τc isqdq(E)∆q(τ ) = ∆X (τ ) =(A.21)∆τc = − S0 x0 (τ )∆τc .
dτcComparing the two increments, we obtainq(E)
dc0 = S0 dτc .

(A.22)

Then, considering separately the zero mode in the ratio of the determinant, )−1/2
det − Md2 /dτ 2 + V ′′ [ X (τ )]det − Md2 /dτ 2 + ω02s( )−1/2(E)
det′ − Md2 /dτ 2 + V ′′ [ X (τ )]S0=
ω0 dτc2πdet′ − Md2 /dτ 2 + ω02s( )−1/2(E)
det′ − Md2 /dτ 2 + V ′′ [ X (τ )]S0=,
ω0 dτc2π
ω0−2 det − Md2 /dτ 2 + ω02((A.23)

where det′ [· · · ] denoted the reduced determinant with the zero mode removed.
We now consider the non-zero modes. It is easiest to deal with the second discrete level, whose eigenvalue is (3/4)ω02 . If we denote by Φ the ratiodet′ − Md2 /dτ 2 + V ′′ [ X (τ )] ,Φ=(A.24)
ω0−2 det − Md2 /dτ 2 + ω02
the contribution of this level to Φ as τ0 → ∞ is 3/4. We turn to other modes, withϵ > ω02 . In this region, Eq (A.19) has a continuous spectrum without the boundary condition. Therefore, let us forget the boundary condition for a moment. The

A.2. Pre-exponential factor

63

solution in the continuous spectrum has the following properties: First, the solutions are labeled by a continuous index p, which is related to the eigenvalue ϵ by√p = ϵ p − ω02 and ranges over the entire interval (0, ∞). Second, the continuum
wavefunctions determined by the Schrödinger equation, Eq. (A.19), include no reflection. Namely, choosing one of the linearly independent solutions in such a waythatx p (τ ) = eipτ ,τ → +∞,(A.25)

we have in the other asymptotic region the same exponential:x p (τ ) = eipτ +iδp ,τ → −∞.

(A.26)The entire dynamical effect is reduced to the phase [67]eiδp =

1 + ip/ω0 1 + 2ip/ω0.
1 − ip/ω0 1 − 2ip/ω0(A.27)

Since the second linearly independent solution can be chosen in the form x p (−τ ),the general solution is Ax p (τ ) + Bx p (−τ ), where A and B are arbitrary constants.
Considering the boundary condition x (±τ0 /2) = 0, the arbitrary constants satisfy
Ax p (±τ0 /2) + Bx p (∓τ0 /2) = 0.

(A.28)In order to have nontrivial solutions for Eq. (A.28), it is necessary to satisfy
x p (τ0 /2)= ±1,
x p (−τ0 /2)(A.29)

which gives the equation for p:pτ0 − δp = πn,n = 0, 1, . . . .

(A.30)

If we denote the solutions of Eq. (A.30) by p̃n = pn + δp /τ0 where pn = πn/τ0 , theratio Φ of the determinant in ϵ > ω0 is written as∞

ω02 + p̃2n.
22n =0 ω0 + p n

Φcon. ≃ ∏(A.31)

Since the difference p̃n − pn is small in the limit τ0 → ∞, Φcon. can be written as
"
#
"
#
∞
∞
ω02 + p2n2pn ( p̃n − pn )Φcon. = exp ∑ ln 2
≈ exp ∑.
(A.32)
ω0 + p2n
ω02 + p2nn =0n =0Changing from the summation over n into integration over pn and using Eq. (A.30),
we obtain

 Z ∞


 
Z
2δp p
1
p2
1 ∞ dδp
dp 2
ln 1 + 2 dp
Φcon. = exp= exp −π 0
π 0 dpp + ω02
ω0




Z ∞

2
1
2
1
= exp −dy+ln 1 + y2 = .
(A.33)
22
π 0
1+y
1 + 4y9

64

Appendix A. Instanton Method

F IGURE A.2: Dilute instanton gas configuration.

Here, we introduced the dimensionless variable y = p/ω0 . Finally, combining thecontribution of the second discrete level and the continuous spectrum, we obtain
Φ=

3 1
1
× = ,
4 9
12

(A.34)

and from Eqs. (A.18) and (A.23), the single-instanton contribution to the transition
amplitude (A.16):⟨±q0 /2| e− HS τ0 |∓q0 /2⟩singless −1/2
(E)
(E)
S0
Mω0
1
− S0
ω0 dτc=e
2π sinh(ω0 τ0 ) 2π
12
! r q
!
r
Mω0 −ω0 τ0 /2
6(E) −S0(E)e
S0 e
ω0 dτc .
=ππ(A.35)

Note that the factor in the first bracket corresponds to a simple harmonic operator, the exponential factor e−S0 is classical contribution accounting for the action of
the instanton, and the remaining represents for the fluctuation around the bouncing
segments of the path. This result can be trusted as long as
q
(E)
(E)
S0 e−S0 ω0 τ0 ≪ 1.
(A.36)The energy of the lowest state is determined by the transition to the limit τ0 →∞. At large τ0 , paths constructed of many instantons and anti-instantons are important. Considering n instantons or anti-instantons with the arbitrary parametersτ1 , τ2 , . . . , τn (see Fig. A.2). If the characteristic intervals satisfy |τi − τj | ≫ ω0−1 , thenindividual instantons and anti-instantons do not interact each other. Therefore, theA.2. Pre-exponential factor

65

transition amplitudes ⟨−q0 /2| e− HS τ0 |+q0 /2⟩ is obtained by summation over n:⟨−q0 /2| e− HS τ0 |+q0 /2⟩
! Z
r
Z τn
Z τ2
+τ0 /2Mω0 −ω0 τ0 /2edn
ω0 dτnω0 dτn−1 · · ·ω0 dτ1= ∑π−τ0 /2−τ0 /2−τ0 /2n∈odd!
rMω0 −ω0 τ0 /2(ω0 τ0 )n
= ∑
e
dn
π
n!
n∈odd
r
Mω0 −ω0 τ0 /2esinh (ω0 τ0 d) ,(A.37)=
π
where d is the instanton density
r q
(E)
6
(E)
d=
S 0 e − S0 .
π
Similarly, ⟨+q0 /2| e− HS τ0 |+q0 /2⟩ is obtained asrMω0 −ω0 τ0 /2− HS τ0⟨+q0 /2| eecosh (ω0 τ0 d) .
|+q0 /2⟩ =π(A.38)(A.39)

Finally, we consider how the transition amplitude (A.37), (A.39) can be understood in physical terms. Considering the two-state system (2.8), the transition amplitude can be expressed as⟨±q0 /2| e− HS τ0 |+q0 /2⟩= ⟨±q0 /2| e(∆0 /2)σx τ0 |+q0 /2⟩hi= ⟨±q0 /2| | Eg ⟩ e(∆0 /2)τ0 ⟨ Eg | + | Ee ⟩ e−(∆0 /2)τ0 ⟨ Ee | |+q0 /2⟩ .

(A.40)

Using ⟨+q0 /2| Eg ⟩ = ⟨−q0 /2| Eg ⟩ = ⟨+q0 /2| Ee ⟩ = −⟨−q0 /2| Ee ⟩,⟨+q0 /2| e− HS τ0 |+q0 /2⟩ ∝ cosh(∆0 τ0 /2),⟨−q0 /2| e− HS τ0 |+q0 /2⟩ ∝ sinh(∆0 τ0 /2).

(A.41)(A.42)

Comparing this expression with Eq. (A.37), the tunneling amplitude ∆0 is expressedass(E)
6S0 −S(E)∆0 = 2ω0 d = 2ω0e 0 .
(A.43)π
Rewriting it in terms of the potential barrier height Vb , we obtain Eq. (2.7).

67

Appendix B

The Influence function of theSpin-boson ModelIn this appendix, following Ref. [15], we derive the influence function of the spinboson model (4.2). Let us start with the Caldeira-Leggett model (2.1). For simplicity, we consider the local system coupled with one reservoir. Using a standard
imaginary-time path integral, the density matrix at β = 1/kB T can be written as
′′

′′

′

′

Wβ (q , x ; q , x ) = Zt−1

Z q( β)=q′′q(0)=q′D q(τ )

Z x( β)= x′′x(0)= x′D x ( τ ) e−S(E) [ q ( τ ),x ( τ )]

,

(B.1)

where Zt is the partition function of the global system. As explicitly indicated, thefunctional integrations run over all the paths taken by the coordinates of the local
system q(τ ) and reservoir x(τ ) with endpoints q(0) = q′ , x(0) = x′ , q( β) = q′′ , andx( β) = x′′ . The Euclidean action S(E) is given by the sum of the contributions fromthe system, the reservoir, and the interaction as follows:
S

(E)

(E)
(E)
(E)
= SS + SB + SI =

Z β
0

h
i
(E)
(E)
(E)
dτ LS (τ ) + LB (τ ) + LI (τ ) ,(B.2)

where the Euclidean Lagrangians of the Caldeira-Leggett model (2.1) are given byMq̇(τ )2 + V [q(τ )],
2
N

mα 
(E)
LB ( τ ) = ∑ẋα (τ )2 + ωα2 xα (τ )2 ,
2
α =1


N
1 Cα2 q(τ )2(E).
LI (τ ) = ∑ −Cα xα (τ )q(τ ) +2 mα ωα2
α =1(E)

LS ( τ ) =

(B.3)
(B.4)
(B.5)

Now, we focus on a reduced description in which the reservoir coordinates, x(τ ),are eliminated. To this end, we introduce the reduced density matrix
′′

′



′′

′′

′

′



ρ β (q , q ) = TrB Wβ (q , x ; q , x ) =Z ∞−∞dx′ Wβ (q′′ , x′ ; q′ , x′ ),(B.6)

and the reduced partition function
Z ( β) =Zt ( β),ZB ( β)

(B.7)

68

Appendix B. The Influence function of the Spin-boson Model

where ZB is the partition function of the isolated reservoir described by the collectionof N harmonic oscillators
N

N

1
2
sinh
(
βωα /2)α =1

ZB = ∏ ZB = ∏
(α)

α =1

(B.8)

Thus, the reduced density matrix can be written in the form′′′

ρ β (q , q ) = Z−1Z q( β)=q′′q(0)=q′(E)

D q(τ ) e−SS [q(τ )] F (E) [q(τ )],(B.9)

where F (E) is the influence functional(E)

F (E) [q(τ )] ≡ e−Sinfl [q(τ )]
I
(E)
(E)
1
=
D x(τ ) e−SB [q(τ ),x(τ )]−SI [q(τ ),x(τ )] .
(B.10)ZB
H
The integration D x(τ ) indicates the sum over all periodic paths with period β
taken by the coordinates of the reservoir x. In the following, we calculate the Eu(E)clidean influence action Sinfl [q(τ )].
We choose periodic continuation of the paths q(τ ) and x(τ ). Therefore, we can
expand the coordinates as Fourier seriesxα (τ ) =1 ∞ iνn τ∑ e x̃α,n ,β n=−∞(B.11)q(τ ) =

1 ∞ iνn τ∑ e q̃n ,β n=−∞(B.12)∗∗
where x̃α,n = x̃α,−n , q̃n = q̃−n , and νn = 2πn/β, ( n = 0, ±1, ±2, . . . ) is a Matsubara
frequency. Substituting the Fourier expansion (B.11) and (B.12) into Eq. (B.2), we
obtain
(E)

(E)

(E)

SB,I [q(τ ), x(τ )] ≡ SB [q(τ ), x(τ )] + SI [q(τ ), x(τ )]1 ∞ mα
= ∑
∑ 2
β n=−
∞
α =1
N

Cα q̃n 2
νn2 | x̃α,n |2 + ωα2 x̃α,n −mα ωα2!
. (B.13)

Next, we split x̃α,n into the classical part x̄α,n and the quantum part yα,n asx̃α,n = x̄α,n + yα,n =Cα q̃n+ yα,n ,mα (νn2 + ωα2 )(B.14)

where x̄α,n is the stationary solution determined from the equation of motionmαd2 x̄α (τ )− mα ωα2 x̄α (τ ) + Cα q(τ ) = 0.
dτ 2

(B.15)

Appendix B. The Influence function of the Spin-boson Model

69

Substituting Eq. (B.14) into (B.13), we obtain
(E)

(E)

(E)

SB,I [q(τ ), x(τ )] = SB [y(τ )] + Sinfl [q(τ )](B.16)

1 ∞ mα 2∑ β ∑ 2 νn + ωα2 |yα,n |2
n=−∞
α =1
N

(E)

SB [y(τ )] =

N

mα 

Z β

∑ 0 dτ 2

=

α =1
N

ẏα (τ )2 + ωα2 yα (τ )2

C2 1 ∞= ∑ α∑
2mα β n=−
∞
α =1

(E)
Sinfl [q(τ )]



1
1
− 2
ωα2
νn + ωα2

(B.17)| q n |2 .

(B.18)The path sum of all β-periodic quantum fluctuations y(τ ) directly yields the partition function of the reservoir,
ZB =

I

(E)

D y(τ ) e−SB [y(τ )] .

(B.19)

Hence this term cancels the factor ZB−1 in Eq. (B.10). The influence action (B.18) may
be concisely written as(E)

Sinfl [q(τ )] =
ξn =

M1 ∞
ξ n | q n |2 ,2 β n∑=∞(B.20)

1 N Cα2νn22=∑
M α=1 mα ωα2 (νn2 + ωα2 )
Mq20

Z ∞
0

dω

I (ω )
νn2, (B.21)ω (νn2 + ω 2 )where I (ω ) is the spectral density defined in Eq. (2.17). Performing the inverse
Fourier transformation, the influence action (B.20) becomes(E)
Sinfl [q(τ )]

Z β

=

dτ

0

Z τ
0

dτ ′ k (τ − τ ′ )q(τ )q(τ ′ ),M ∞ iνn τ∑ e ξn.
β n=−∞k(τ ) =

Note that the kernel k(τ ) satisfies k (τ ) = k ( β − τ ) = k (τ + β) andUsing these properties, the influence function (B.22) is rewritten as(E)

Sinfl [q(τ )] = −1
2(B.22)

Z β

dτ

Z τ

0

0


2
dτ ′ k (τ − τ ′ ) q(τ ) − q(τ ′ ) .

(B.23)

Rβ0dτ k (τ ) = 0.

(B.24)(E)

Evidently, the influence action Sinfl [q(τ )] is fully nonlocal. Alternatively, the influence action is expressed in terms of the kernelk̃ (τ ) = µ : δ(τ ) : −k (τ ),(B.25)

where : δ(τ ) : is the periodically continued δ-function,(E)

Sinfl [q(τ )] =

1
2

Z β
0

dτ

Z τ
0

2

dτ ′ k̃ (τ − τ ′ ) q(τ ) − q(τ ′ ) .

(B.26)

70

Appendix B. The Influence function of the Spin-boson Model

With choiceN

2
Ca α 2= 2
2m ωq0
α =1 α α

µ= ∑

Z ∞

I (ω ),
ω(B.27)

ωdω I (ω ) 2
eiνn τ
νn + ω 2(B.28)

cosh[ω ( β/2 − τ )].
sinh(ωβ/2)(B.29)dω0

the kernel can be expressed ask̃ (τ ) ==∞2∑
q20 β n=−∞

1
q20

Z ∞
0

Z ∞
0

dω I (ω )Considering the spin-boson model, the path q(τ ) is piecewise constant and it
occasionally jumps between the positions +q0 /2 and −q0 /2, which are positions
give minimal potential energy of the local system. It is convenient to putq(τ ) =q0
σz (τ ),
2

(B.30)

where σz (τ ) is a spin path. Thus, the spin path, σz (τ ), jumps forth and back between
the variables +1 and −1, and is expressed by2mσz (τ ) = 1 + 2 ∑ (−1) j θ (τ − s j ),(B.31)

j =1

which implies a path of 2m alternating flips with centers {s j }. Introducing the kernel
K (τ ) ≡ q20 k̃ (τ ), the influence function (B.10) takes form


Z
Z τ


1 β
(E)
′
′
′ 2
F [σz (τ )] = exp −dτdτ K (τ − τ ) σz (τ ) − σz (τ )(B.32)8 00Z ∞cosh[ω ( β/2 − τ )].
(B.33)K (τ ) =dω I (ω )sinh(ωβ/2)
0From these results, it is straightforward to derive the partition function (4.1) with thekernel (4.3).

71

Appendix C

The Noninteracting-blip
Approximation (NIBA)
In this appendix, based on Ref. [15], we derive the symmetrized correlation function (3.48) in the NIBA and discuss the conditions for the NIBA. The density matrixof the global system (2.8), Ŵ (t), is written in the real-time path integral form as⟨q f , x f | Ŵ (t) |q′f , x′f ⟩=Z

dqi dqi′ dxi dxi′ K (q f , x f , t; qi , xi , 0) ⟨q f , x f | Ŵ (0) |q′f , x′f ⟩×K ∗ (q′f , x′f , t; qi′ , xi′ , 0),(C.1)

where qi/ f and xi/ f are the coordinates of the local system and reservoir, respectively,
and K is the coordinate representation of the time evolution operator
K (q f , x f , t; qi , xi , 0) = ⟨q f , x f | e−iHt |qi , xi ⟩=Here, the functional integration,
pointsq (0) = q i ,RZ

D q(t)D x(t) eiS[q(t),x(t)] .

(C.2)
(C.3)

D q(t)D x(t), extends over all paths with end-

q(t) = q f ,x (0) = x i ,x(t) = x f ,(C.4)

and S[q(t), x(t)] is the action of the global system. The density matrix of the global
system describes the dynamics of the system-plus-environment complex as a whole.
However, in most cases of interest, the only information one wishes to have is thesystem’s dynamics under the reservoir’s influence. Then, the quantity one is reallyinterested in is the reduced density matrix defined as
ρ(q f , q′f ; t) =Z

dx f ⟨q f , x f | Ŵ (t) |q′f , x f ⟩ .

(C.5)

Assuming that the density matrix of the global system at t = 0 is in a productinitial stateŴ (0) = ρS (0) ⊗1 − βHB
e
,
ZB

(C.6)

72

Appendix C. The Noninteracting-blip Approximation (NIBA)F IGURE C.1: Graphical representation of sojourns and blips. Λ j,k isthe interaction between blip j and blip k, and X j,k is the interaction
between blip j and sojourn k.

where ρS is the density matrix of the local system, e.g., ρS (0) = |↑⟩ ⟨↑|, the reduced
density matrix is expressed by
ρ(q f , q′f ; t)=

JFV (q f , q′f , t : qi , qi′ , 0) =Z

dqi dqi′ JFV (q f , q′f , t : qi , qi′ , 0)ρ(qi , qi′ ; 0),Z

D qD q′ ei(SS [q]−SS [q ]) FFV [q, q′ ],

(C.7)

′

(C.8)

where FFV [q, q′ ] is the Feynman-Vernon influence functional. In the two-state system, it can be expressed in terms of the spin path [68]:

FFV [σz , σz′ ]ZZ ′
1 t ′ t ′′ dt
dt σz (t′ ) − σz′ (t′ )= exp −
4 0
0


′
′′
× Q̈(t − t )σz (t′′ ) − Q̈∗ (t′ − t′′ )σz′ (t′′ ) ,



(C.9)

where σz (τ ) = q(τ )/(q0 /2), σz′ (τ ) = q′ (τ )/(q0 /2), and Q(t) is the kernel given as




Z ∞
βω
I (ω )
coth
1
−
cos
(
ωt
)]
+
i
sin
(
ωt
)
.
(C.10)
Q(t) =dω[ω20
We note that Q(t) is related with W (τ ) in Eq. (4.6) via the analytical continuation as
Q(t) = W (τ = it). Introducing antisymmetric and symmetric spin paths
σz (t) − σz′ (t),2σz (t) + σz′ (t)η (t) =,
2

ξ (t) =

(C.11)(C.12)

respectively, the influence functional (C.9) is rewritten as

FFV [σz , σz′ ] = exp

Z t
0

′

′

dt ξ̇ (t )

Z t′
0

dt

′′



′

′

′′

′′

′′

′

′′

′′

Q (t − t )ξ̇ (t ) + iQ (t − t )η̇ (t ),
(C.13)

where Q′ (t) ≡ Re[ Q(t)] and Q′′ (t) ≡ Im[ Q(t)].
In the two-state system, the reduced density matrix is 2 × 2 matrix. Intervalsof the spin path dwelling in the diagonal and off-diagonal states have been termed‘sojourn’ and ‘blip’, respectively [6]. During a sojourn, the function ξ (t) is zero,
whereas during a blip interval the function η (t) is zero. There are two sojourn states

Appendix C. The Noninteracting-blip Approximation (NIBA)73(ηi = ±1) and two blip states (ξ i = ±1). A general sojourn-to-sojourn path making
2n transitions at intermediate times t j , ( j = 1, 2, . . . , 2n) is parametrized byn∑ ηj θ (t − t2j ) − θ (t − t2j+1 ) ,η (n) ( t ) =(C.14)∑ ξ j θ (t − t2j−1 ) − θ (t − t2j ) .

(C.15)

j =0n

ξ (n) ( t ) =





j =0

Denoting the intervals spent in blip and sojourn states by τj and s j , respectively, we
have (see Fig. C.1)
τj = t2j − t2j−1 ,s j = t2j+1 − t2j .

(C.16)

Thus, the influence functional (C.13) can be written as

F (n) = Gn Hn ,"(C.17)#nn j −1"#′
Gn = exp − ∑ Q2j,2j−1 exp − ∑ ∑ ξ j ξ k Λ j,k ,j =1"(C.18)

j =2 k =1

n j −1#"n −1

Hn = exp i ∑ ∑ ξ j ηk X j,k = exp i ∑j =1 k =0n#∑ ξ j ηk Xj,k .

(C.19)

k =0 j = k +1

Here, Λ j,k is an interaction between blip j and blip k, and X j,k is an interaction between blip j and sojourn k (see Fig. C.1):′′′′
Λ j,k = Q2j,2k−1 + Q2j−1,2k − Q2j,2k − Q2j−1,2k−1 ,(C.20)′′′′′′′′
Q2j,2k+1 + Q2j−1,2k − Q2j,2k − Q2j−1,2k+1 ,(C.21)

X j,k =

where Q j,k = Q(t j − tk ). The blip interactions are bundled up in the real-valued
function Gn . The first exponential factor in Eq. (C.18) contains the intrablip interactions, and the second exponential factor represents the interblip correlations. Theinteractions between the sojourns and the blips are in the phase factor Hn . However,the sojourns do not interact with each other. The function Gn is a filtering function
which suppresses long blips. Hence, it is favored that the system dwells in sojournstates. Physically, this is because continuous measurement of σz due to the environment suppresses quantum interference between the eigenstates of σz , and results inthe decrease of occupation lengths of blip states.
Using the influence functional (C.17), we obtain the dynamics of the population⟨σz ⟩t as follows:∞⟨σz ⟩t = 1 + ∑ (−1)

m

Z t

m =1

1
D2m,0 {t j } m ∑ Fm ,2 { ξ ±1}0(C.22)

=

where
Z t
0

D2m,0 {t j } =Z t
0dt2m · · ·Z t20

dt1 ∆2m ,(C.23)

74

Appendix C. The Noninteracting-blip Approximation (NIBA)

and the effect of the reservoir is included in the influence function
!
m −1

m

Fm = Gm ∏ cos∑ ξ j Xj,k

k =0

.

(C.24)

j = k +1

In addition, one can derive the generalized master equation for the population:d ⟨σz ⟩t
=−
dt

Z t
0

dt′ Kz (t, t′ ) ⟨σz ⟩t ,(C.25)

Compared Eq. (C.22) with (C.25), the kernel Kz (t, t′ ) is expressed as∞

Kz (t, t′ ) = ∆2 F1 (t, t′ ) + ∑ (−1)n−1
n =2

∆2n
2n

Z t
t′

dt2n−1 · · ·

Z t3
t′

dt2

∑ F̃n ,

(C.26){ξ j =±}

where F̃n is irreducible influence function defined by subtraction of the reducible
parts in Fn . Here, ‘irreducible’ means that it cannot be separated into two or more
uncorrelated parts at sojourns.
Now, we turn to the discussion of the dynamics within the noninteracting-blip
approximation (NIBA). The simple assumption underlying the NIBA is that the average time ⟨s⟩ spent by the system in sojourn states is very large compared to theaverage time ⟨τ ⟩ spent in blip states. In detail, the NIBA assumption leads to two
simple prescriptions regarding the sojourn-blip correlation X j,k and interblip interactions Λ j,k .
1. Set the sojourn-blip correlation(
0,
X j,k =
Q′′ (τk+1 ),( j ̸ = k + 1),( j = k + 1).

(C.27)

2. Set all interblip interactions Λ j,k equal to zero.
With these specifications, the influence functional (C.17) reduces to a factorized formof intrablip correlations in which the sign of the individual blip phase depends on
the label of respective preceding sojourn,
n


(n)
FNIBA = ∏ exp − Q′ (τj ) + iξ j η j−1 Q′′ (τj ) .

(C.28)

j =1

Generally, the NIBA can be justified in three cases discussed in Section 3.5 of themain text.
In the NIBA, since the irreducible influence functions F̃n become zero for n > 1,
the kernel Kz (t, t′ ) is approximated as′′
Kz,NIBA (t, t′ ) = ∆2 F1 (t, t′ ) = ∆2 e−Q (t−t ) cos Q′′ (t − t′ ) .

(C.29)

Since this kernel depends only on the relative time t − t′ , one can perform the Laplacetransformation for it:
Kz,NIBA (λ) = ∆2

Z ∞
0



′
′
dτ e−λτ e−Q (t−t ) cos Q′′ (t − t′ ) .

(C.30)

Appendix C. The Noninteracting-blip Approximation (NIBA)

75

From Eq. (C.25), the Laplace transformation of the kernel is related to the population
under the condition ⟨σz ⟩t=0 = +1 as follows:⟨σz (λ)⟩ =1.
λ + Kz ( λ )(C.31)The symmetrized correlation function S(t) defined in Eq. (2.44) is related to the
population as follows [56]:S(t) = ⟨σz ⟩|t| ,(C.32)

under the condition ⟨σz ⟩t=0 = +1. Performing the Fourier transformation for Eq. (C.32)
and using Laplace transformation of the population (C.31),
S(ω ) = 2ReZ ∞0

dt eiωt ⟨σz ⟩t(C.33)
= 2Re [⟨σz (−iω )⟩] = 2Re1.
−iω + Kz (−iω )(C.34)

Rewriting the kernel Kz (λ) with the self-energy Σ(λ) (3.49), we finally derive Eq. (3.48).

77

Appendix D

Impedances of the Transmission
Line
We analyze the frequency dependence of the spectral density for the circuit model
discussed in Chapter 6. Assuming |ωCj Zj+1 (ω )| ≪ 1, the recurrence relation (6.40)
is given approximately:
Zj (ω ) ≃ R j + iωL j + Zj+1 (ω ) − iωCj Zj+1 (ω )2 .

(D.1)

In the continuous limit N → ∞, where N is the number of islands in the transmissionline, this recurrence relation reduces to the differential equation:∂Z (ω, x )≃ −r ( x ) − iωl ( x ) + iωc( x ) Z (ω, x )2 ,
∂x

(D.2)

where r ( x ), l ( x ), and c( x ) (0 ≤ x = j/N ≤ 1) are the resistance, inductance, andcapacitance per unit length, respectively. From Eq. (6.41), they are given asr ( x ) = r0 x n ,

(D.3)

p

l ( x ) = l0 x ,

(D.4)

m

(D.5)

c ( x ) = c0 x ,where r0 = R0 N, l0 = L0 N, and c0 = C0 N. We note that Z (ω ) = Z (ω, x → 0)/2.

D.1

The flux qubit

For the choice of the parameters given in Eq. (6.42), the differential equation for the
impedance is written as∂Z (ω, x )≃ −r0 x n − iωl0 + iωc0 x m Z (ω, x )2 .
∂x

(D.6)

Under the condition
x∗ ≡



n
√

2ω l0 c02/(m+2)



≪x≪

ωl0
r0

1/n
,

(D.7)

Ż (ω, x ) = ∂Z (ω, x )/∂x and r0 x n are sufficiently small compared with other terms in
the right-hand side of Eq. (D.6). Then, we obtains
l0 −m/2
Z A (ω, x ) =
x
,
(D.8)
c0

78

Appendix D. Impedances of the Transmission Line

In contrast, for x ≃ 0, we can neglect r ( x ) and c( x ), and obtain the following:
ZB (ω, x ) = −iωl0 x + A(ω ).

(D.9)

The constant of integration, A(ω ), can be determined by the equation Z A (ω, x ∗ ) =ZB (ω, x ∗ ). Thus, we arrive at the approximate expression of Z (ω ) as follows:
ZB (ω, x → 0)2s

Z (ω ) ∼12=iωl0 x ∗ +l0 ∗ −m/2
(x )
c0

!
.

(D.10)

From Eq. (6.27), we obtain the following the spectral density:hi
I (ω ) ∝ I0 (ω ) = ωRe Z (ω )−1 ∝ ω 2/(m+2)(D.11)

Note that since m is a non-negative real number, this circuit realizes sub-Ohmic reservoirs with arbitrary sub-Ohmic exponent s (s < 1). This frequency dependence appears for ω ∗ ≪ ω ≪ ωc , where the lower bound, ω ∗ , is obtained by considering thecondition (D.7):
"
∗

ω =

D.2

 m 2nr0m+2

2

c0n l0m+n+2#1/(m+2n+2).

(D.12)

The charge qubit

For the choice of the parameters given in Eq. (6.45), the differential equation for the
impedance is written asŻ (ω, x ) ≃ −r0 x n − iωl0 x p + iωc0 Z (ω, x )2 .

(D.13)The spatial dependence of the elements is different from the case of the flux qubit.
However, we can analyze the behavior of the impedance in the same way. In thiscase, neglecting Ż (ω, x ) and r0 x n give the impedance ass
l0 p/2Z A (ω, x ) =x ,(D.14)
c0
for
∗

x ≡



p
√

2/( p+2)2ω l0 c0



≪x≪

ωl0
r0

1/(n− p).

(D.15)

In contrast, for x ≃ 0, we can neglect terms that elements have the spatial dependence and obtain
ZB (ω, x ) = −1.
iωc0 x + A(ω )(D.16)

D.3. The charge qubit (0 < s ≤ 0.5)

79

From the condition Z A (ω, x ∗ ) = ZB (ω, x ∗ ) and Eq. (6.39), the frequency dependenceof the spectral density is obtained as
I (ω ) ∝ I0 (ω ) = ωRe[ Z (ω )] ∝ ω 2/( p+2) .

(D.17)

Note that since p is a non-negative real number, this circuit realizes sub-Ohmic reservoirs with arbitrary sub-Ohmic exponent s (s < 1). This frequency dependence forI (ω ) appears for ω ∗ ≪ ω ≪ ωc . The lower bound, ω ∗ , is obtained from the condition (D.15):
"
ω∗ =

D.3

#1/(2n− p+2)

p +2

 p  2( n − p )r0.

n − p n +2
l0

2

c0

(D.18)

The charge qubit (0 < s ≤ 0.5)

We can analytically investigate the power of the impedance in the same way as theabove two cases. For the choice of the parameters given in Eq. (6.48), the differential
equation for the impedance is written asŻ (ω, x ) ≃ −r0 x n + iωc0 Z (ω, x )2 .

(D.19)

For
∗

x ≡



n2
4ωr0 c01/(n+2)≪ x,(D.20)

r0 n/2x .
ωc0(D.21)

we can neglect Ż (ω, x ) and obtain:
Z A (ω, x ) = e

−iπ/4

r

On the other hand, for x ≃ 0, we can neglect the first terms of the right hand side in
Eq. (D.19). It gives the impedance for x ≃ 0 as
ZB (ω, x ) = −1,iωc0 x + A(ω )(D.22)

where A(ω ) is the constant of integration and determined by the equation Z A (ω, x ∗ ) =ZB (ω, x ∗ ). Therefore, from Eq. (6.39), we obtain the spectral density as follows
I (ω ) ∝ I0 (ω ) = ωRe [ Z (ω )] ∝ ω 1/(n+2) .

(D.23)

Note that since n is a non-negative real number, this circuit realizes sub-Ohmic reservoirs with 0 < s ≤ 0.5. This frequency dependence appears for ω ∗ ≪ ω ≪ ωc ,where the lower bound, ω ∗ , and cutoff frequency, ωc , are given as
ω∗ =
ωc =
respectively.

√ ! n +2
2 2
1  n 2
1+,
r0 c0 2n
N 2 2N n,r0 c0n(D.24)(D.25)

81

Bibliography[1] R. Landauer, IBM J. Res. Dev. 1, 223 (1957).
[2] S. Datta, Electronic Transport in Mesoscopic Systems (Cambridge University Press,
Cambridge, 1995).
[3] K. Saito, Europhys. Lett. 83, 50006 (2008).
[4] K. Schwab, E. A. Henriksen, J. M. Worlock, and M. L. Roukes, Nature 404, 974(2000).
[5] A. C. Hewson, The Kondo Problem to Heavy Fermions (Cambridge University
Press, Cambridge, 1997).
[6] A. J. Leggett et al., Rev. Mod. Phys. 59, 1 (1987).
[7] D. Segal and A. Nitzan, Phys. Rev. Lett. 94, 034301 (2005).
[8] T. Ojanen and A. P. Jauho, Phys. Rev. Lett. 100, 155902 (2008).
[9] R. M. Potok, I. G. Rau, H. Shtrikman, Y. Oreg, and D. Goldhaber-Gordon, Nature 446, 167 (2007).
[10] Z. Iftikhar et al., Nature 526, 233 (2015).
[11] Z. Iftikhar et al., Science 360, 1315 (2018).
[12] D. L. Cox and A. Zawadowski, Adv. Phys. 47, 599 (1998).
[13] R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
[14] A. Ronzani et al., Nat. Phys. 14, 991 (2018).
[15] U. Weiss, Quantum Dissipative Systems, 4th ed. (World Scientific, Singapore,1999).
[16] Y. Makhlin, G. Schön, and A. Shnirman, Rev. Mod. Phys. 73, 357 (2001).
[17] B. Peropadre, D. Zueco, D. Porras, and J. J. García-Ripoll, Phys. Rev. Lett. 111,
243602 (2013).
[18] J. Leppäkangas et al., Phys. Rev. A 97, 052321 (2018).
[19] J. Bourassa et al., Phys. Rev. A 80, 032109 (2009).
[20] A. Nitzan, Chemical Dynamics in Condensed Phases: Relaxation, Transfer, and Reactions in Condensed Molecular Systems (Oxford University Press, Oxford, 2006).
[21] J. Ren, P. Hanggi, and B. Li, Phys. Rev. Lett. 104, 170601 (2010).
[22] D. Segal, A. J. Millis, and D. R. Reichman, Phys. Rev. B 82, 205323 (2010).

82

BIBLIOGRAPHY[23] T. Ruokola and T. Ojanen, Phys. Rev. B 83, 045417 (2011).
[24] T. Chen, X. B. Wang, and J. Ren, Phys. Rev. B 87, 144303 (2013).
[25] Y. Yang and C. Q. Wu, Europhys. Lett. 107, 30003 (2014).
[26] C. Wang, J. Ren, and J. Cao, Sci. Rep. 5, 11787 (2015).
[27] J. Liu, H. Xu, B. Li, and C. Wu, Phys. Rev. E 96, 012135 (2017).
[28] F. Guinea, V. Hakim, and A. Muramatsu, Phys. Rev. B 32, 4410 (1985).
[29] K. Saito and T. Kato, Phys. Rev. Lett. 111, 214301 (2013).
[30] L. Magazzù et al., Nat. Commun. 9, 1403 (2018).
[31] P. Forn-Díaz et al., Nat. Phys. 13, 39 (2017).
[32] J. E. Mooij et al., Science 285, 1036 (1999).
[33] P. W. Anderson and G. Yuval, J. Phys. C: Solid State Phys. 4, 607 (1971).
[34] J. M. Kosterlitz, Phys. Rev. Lett. 37, 1577 (1976).
[35] R. Bulla, N. H. Tong, and M. Vojta, Phys. Rev. Lett. 91, 170601 (2003).
[36] M. Vojta, N. H. Tong, and R. Bulla, Phys. Rev. Lett. 94, 070604 (2005).
[37] A. Winter, H. Rieger, M. Vojta, and R. Bulla, Phys. Rev. Lett. 102, 030601 (2009).
[38] M. Vojta, N.-H. Tong, and R. Bulla, Phys. Rev. Lett. 102, 249904 (2009).
[39] A. W. Chin, J. Prior, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 107, 160601(2011).
[40] M. Vojta, Phys. Rev. B 85, 115113 (2012).
[41] E. Luijten, Interaction range, universality and the upper critical dimension, PhD
thesis, Delft University of Technology, 1997.
[42] E. Luijten and H. W. Blöte, Phys. Rev. B 56, 8945 (1997).
[43] N. H. Tong and M. Vojta, Phys. Rev. Lett. 97, 016802 (2006).
[44] L. B. Yu, N. H. Tong, Z. Y. Xue, Z. D. Wang, and S. L. Zhu, Sci. China Phys. 55,
1557 (2012).
[45] M. Kato, Master Thesis (Graduate School of Science, The Universe of Tokyo,
2017).
[46] A. O. Caldeira and A. J. Leggett, Physica A 121, 587 (1983).
[47] A. I. Vaı̆nshteı̆n, V. I. Zakharov, V. A. Novikov, and M. A. Shifman, Sov. Phys.
Usp. 25, 195 (1982).
[48] J. Maciejko, An Introduction to Nonequilibrium Many-Body Theory .
[49] A. J. Bray and M. A. Moore, Phys. Rev. Lett. 49, 1545 (1982).
[50] S. Chakravarty, Phys. Rev. Lett. 49, 681 (1982).

BIBLIOGRAPHY

83[51] S. K. Kehrein, A. Mielke, and P. Neu, Z. Phys. B 99, 269 (1996).
[52] S. K. Kehrein and A. Mielke, Phys. Lett. A 219, 313 (1996).
[53] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford
University Press, Oxford, 2002).
[54] M. Kato, Theory of heat transport via a local two-state system, Master’s thesis,
The University of Tokyo, 2016.
[55] H. Shiba, Prog. Theor. Phys. 54, 967 (1975).
[56] M. Sassetti and U. Weiss, Phys. Rev. Lett. 65, 2262 (1990).
[57] B. K. Agarwalla and D. Segal, New J. Phys. 19, 043030 (2017).
[58] U. Weiss and H. Grabert, Phys. Lett. A 108, 63 (1985).
[59] H. Rieger and N. Kawashima, Eur. Phys. J. B 9, 233 (1999).
[60] H. Spohn and R. Dümcke, J. Stat. Phys. 41, 389 (1985).
[61] R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 58, 86 (1987).
[62] G. A. Baker Jr., Essentials of Padé Approximants (Academic Press, New York,1975).
[63] K. Völker, Phys. Rev. B 58, 1862 (1998).
[64] H. Bruus and K. Flensberg, Many-Body Quantum Theory in Condensed Matter
Physics (Oxford University Press, Oxford, 2004).
[65] H. Miyazaki, Y. Takahide, A. Kanda, and Y. Ootuka, Phys. Rev. Lett. 89, 197001(2002).
[66] A. Altland and B. Simons, Condensed Matter Field Theory, 2nd ed. (Cambridge
University Press, Cambridge, 2010).
[67] L. D. Landau and E. M. Lifshitz, Quantum Mechanics: Non-Relativistic Theory,
3rd ed. (Butterworth-Heinemann, Oxford, 1981).
[68] R. P. Feynman and F. L. Vernon, Ann. Phys. (N. Y). 24, 118 (1963).