Relaxation of glassy systems

Relaxation of glassy systems
Submitted by: Ofek Asban
Advisor: Prof. Moshe Schechter
Department of Physics
Faculty of Natural Sciences
Ben-Gurion University of the Negev
April 28, 2016
Relaxation of glassy systems
Submitted by: Ofek Asban
Advisor: Prof. Moshe Schechter
Department of Physics
Faculty of Natural Sciences
Ben-Gurion University of the Negev
April 28, 2016
Author:
Date:
Supervisor:
Date:
Chairman of M.Sc. degree committee:
Date:
Abstract
In this work we study the relaxation process of glassy systems and how it is affected by the
interactions. These systems are represented by two similar minimal models: (1) Electronglass (EG), composed of N localized states with random energies and M < N electrons
interacting via unscreened Coulomb interaction. (2) TLS-glass (TG), composed of N Twolevel systems (TLSs), each representing an atom or a configuration of atoms which can have
two possible states with a random energy difference. Also, a pair of TLSs interact via dipolar
interactions. In the pseudospin form the model can be identified as spin-glass Ising model
with a random field.
We start with a review on a work done on the electron-glass model and apply a similar
method for the interacting TLSs model, with more rigorous derivations, while addressing
the common and different characteristics of the two models. The method is conducted as
follows. By solving the linearised rate equations of the sites occupancies in the mean-field
approximation, we obtain the time dependence of the occupation vector, the total relaxation
to equilibrium of the system is then calculated using the norm of the occupation vector. In
the limit where the decay rates are continuous, the norm takes a form of an exponential
integral which for short enough time scale is approximated as a logarithm.
After retrieving the logarithmic time dependence of the occupations for both models,
we then continue to study how the interactions, disorder and system size affect the relaxation,
showing there are also differences between the models. Under several different schemes of
parameter tuning, we conclude that in the electron-glass, stronger interactions and also larger
system size slow down the relaxation. Finally, we apply the same schemes for the interacting
TLSs showing same effects, except in a specific case where the opposite effect is obtained,
i.e., the interactions actually accelerate the relaxation process.
3
Contents
1 Introduction
6
2 Theoretical framework
2.1
2.2
10
The Electron-glass model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.1.1
Mean-Field approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.1.2
Relaxation process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.1.3
Effect of the disorder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
Two-level systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
2.2.1
Single two-level system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
2.2.2
The standard tunnelling model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
2.2.3
Thermodynamic properties of TLSs . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
2.2.4
The TLS-glass model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
3 Numerical solution of the TLSs self-consistent equations
28
4 Relaxation process of TLSs
31
4.1
Rate equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
4.2
The Linear term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
5 Effect of the interactions disorder and system size
37
5.1
TLS-Glass model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
5.2
Electron-Glass model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
6 Summary and conclusions
47
A Fermi-Dirac occupations
51
4
B Single two-levels Hamiltonian
53
C Transition rate for the STM
56
D TLS-TLS interactions
60
E Analytical derivation of the dipole gap
64
F Expansion of the rate equation
67
G Perturbation treatment of the quadratic term
69
H External perturbation and initial deviations
72
I
75
Additional graphs of TG and EG models
J Numerical routines
80
J.1
DistanceMatrix.m file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
J.2
EnergyIteration.m file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
82
J.3
HcTLS.m file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
5
Chapter 1
Introduction
Amorphous solid is a material that lacks long range order characteristic of crystalline solids.
Interactions, in such disordered media, lead to a rough energy landscape which naturally
give rise to a multiplicity of low-lying metastable states that are close in energy. In the
low temperature regime, T < Tg , where Tg is the glass transition temperature, we get
frustration, a situation in which the competition between all the interactions cannot be
satisfied simultaneously in order to reach the equilibrium state. Systems that consist of
these properties are called glasses. Three of the main features characterising glasses are:
1. Memory: the system’s behaviour depends on the previous perturbation applied on it
([1, 2]).
2. Aging: the relaxation process depends on the time duration of the applied perturbation
([3, 4]).
3. Slow relaxation: non-exponential relaxation with a broad distribution of time scale
([5, 6]).
Therefore, many of the glassy systems in low temperatures sharing these universal features
can be captured by minimal models, two of such known models are the Electron-glass and
6
the Two-level systems.
The TLS-glass (TG). A TLS can be pictured as two local minima with an energy
difference (asymmetry energy) separated by a potential barrier, as depicted in Fig. 1.1. Each
TLS acts as an electric or strain dipole source, and interact with other TLSs via dipole-dipole
interactions. The TG Hamiltonian is given in Eq. 2.23.
Figure 1.1: Two-level systems.
Degrees of freedom represented by
displacements in an amorphous material, at low temperatures. Atoms A,
B and C show transverse, longitudinal and rotational displacements respectively. The figure is taken from
Heidelberg University’s website.
The Electron-glass (EG), depicts localized electrons induced by the disorder, which
interact via unscreened Coulomb interaction as can be seen in Fig. 1.2. The EG Hamiltonian
is given in Eq. 2.1).
Figure 1.2:
glass.
Electron-
Electrons at the region of the
Fermi energy that are subject to
a random potential caused by the
bulk. The figure is taken from Ariel
Amir’s website
In the temperature range of interest, both systems experience population transitions
through a weak interaction with a phonon bath, and the dipolar and Coulomb interactions
can be disregarded as the main transition mechanism (which give rise to an energy diffusion
in the system). The similarity of such models naturally leads to the same qualitative slow
logarithmic relaxation as presented in various experiments [5, 6, 3, 4], see also Fig. 1.3.
7
(a)
(b)
Figure 1.3: Relaxation measurements.
Measurements of the logarithmic relaxation for two example amorphous solids after an external electric DC field was applied
for time tw . (a) Relaxation of the relative dielectric constant, δ/, in Mylar sample, depicted effectively by the TLS model in
low temperatures. Where F is the external bias field, ts is the swipe time and T is the temperature. The figure was taken from
Ref. [5]. (b) A typical relaxation of the conductance, G, in In2 O3−x film, effectively depicted by the Electron-glass model in
low temperatures, for IR field. The figure was taken from Ref .[6]. Note how the relaxation process is dependant on tw (aging).
The purpose of this work is: (1) derive, using the mean-field model, the relaxation
to local equilibrium of the TG model. (2) to study how the interactions and system size
affect the relaxation in both the electron-glass and the TLSs models. This is obtained by
using a similar method conducted on the EG model [7]. The method is as follows. The
local-equilibrium of the system is defined by Fermi-Dirac occupation at each site with an
energy obeying the mean-field equations at finite temperature (see Mean-field equations,
Eq. 2.27 for the TG model and Eq. 2.4 for the EG model). The dynamics is represented
by the time dependence of the site occupations (occupations vector) which is obtained by
solving the Pauli master-equation. For small enough perturbation, the occupations vector is
pushed slightly out of the its equilibrium value, allowing one to expand the rate equation to a
linear order. The solution is obtained by diagonalizing numerically the linear prefactor (rate
8
matrix) resulting in a set of N decoupled differential equations for each decay rate, with a
decaying exponential solution. The total time dependence of the system’s relaxation is then
defined as the ”l = 1” norm of the occupations vector. In the limit where the decay rates are
continuous, the norm takes a form of an exponential integral which for short enough times
is approximated as a logarithm.
The work is divided in to two parts: (1) theoretical background which consists a brief
introduction of the Electron-glass model and a review on its relaxation properties as done
in [7], followed by a full derivation of the TG model. (2) new results obtained in this work
which is the derivation of the relaxation process of the TLS model and the dependence of
the relaxation of both model on interactions, disorder and system’s size.
9
Chapter 2
Theoretical framework
2.1
The Electron-glass model
As mentioned in the introduction, there are two competing energy terms; the disorder energy
W and interactions J. In the case of the Electron-glass model (EG), W is the energy range
of site energy distribution (usually uniform distribution) and the interactions are Coulombic
in nature. In addition the sites locations are distributed randomly. The derivation of EG
Hamiltonian from the general many-body Hamiltonian is achieved under few approximations
relevant to the typical material at hand. The main approximations are: (1) Direct tunnelling
term is neglected since the localization length is much smaller than the typical nearest
neighbouring distance ξ rnn , as a result the dominant transport mechanism is phononassisted tunnelling. (2) Strong disorder causes the electrons to be localized and therefore
there is no screening, causing the electrons to interact with a bare long-range Coulomb
interaction. Additional simplifications are taken by neglecting the Hubbard and the exchange
terms: The Hubbard energy is high compared to the characteristic energy scale of the system,
thus, only single electron occupation is allowed at each site. Also, the exchange interaction
is much smaller than the Coulomb interactions making the electron’s spin irrelevant.
10
The final simplified Hamiltonian of the system, for N sites and M < N electrons, is
given by:
Heg =
N
X
i=1
N X 2
X
e
i (ni − K) +
(ni − K) (nj − K)
r
ij
i=1 j>i
(2.1)
where i are the random site energies of the system in absence of interactions, rij is the
distance between sites i and j,
i and j, K =
M
N
e2
rij
is the Coulomb interaction between the electrons and sites
is the background charge which is determined as the ratio of the number
of electrons to the number of sites and ni , nj ∈ [0, 1] are site occupations. The sites are
distributed uniformly on a square area with periodic boundary conditions.
The system, 2.1, is then coupled to a phonon bath:
Hel−ph =
XX
ij
Mq a−q + a†q c†i cj
(2.2)
q
the coupling to phonons is taken as perturbation and accounts for the transitions of electrons
between the sites (see below 2.6). Mq is the strength of the interaction, c†i , cj are the creation
and annihilation operators of electrons in sites i and j and a†q , a−q are the creation and
annihilation operators of phonon with momentum q and −q.
We review the work on electron-glass model, Ref. [7], in three sections. In Section
2.1.1 the single particle density of states (DOS) at equilibrium is obtained from the numerical calculation of the mean-field self-consistent equations of EG Hamiltonian, 2.1, we then
continue to section 2.1.2, using the MF energies to obtain the distribution of the relaxation
rates relevant for the calculation of the total relaxation of the system. Finally, the chapter
ends with the section 2.1.3, showing how the density of states and the rates depends on the
interactions by changing the disorder.
2.1.1
Mean-Field approximation
Taking into account the long-range nature of the Coulomb interactions (even in 2D where
each site effectively interacts with less sites) the mean-field approximation can give a good
11
qualitative estimate. The single particle density of states (DOS) of Electron-glass model
is then found. Starting from the variational derivative of the EG Hamiltonian, 2.1, after
adding a positive background charge K =
1
2
and setting the chemical potential at zero, for
electro-neutrality (half filling), the obtained site energies are:
N
X e2
δHeg
Ei =
= i +
δni
rij
j6=i
1
nj −
2
(2.3)
where i is distributed uniformly in the range − W2 , W2 for disorder W , the interaction
strength parameter is defined as the average Coulomb interaction of nearest neighbours
J = e2 /rnn , which is controlled by scaling the electron charge, and number of sites N
(system size).
In equilibrium, the site occupations obey the Fermi-Dirac statistics (i.e. n0i =
1
eEi /T +1
,see Appendix A), which lead to the self-consistent equations:
1 X e2
Ei = i −
tanh
2 j6=j rij
Ej
2kB T
(2.4)
where Ei are the mean-field energies and kB is Boltzmann constant (set to unity from now
on).
The DOS obtained from the self-consistent Eqs. 2.3, show a gap around the chemical
potential, known as the Coulomb-gap, first predicted by Pollack [8] and then by Efros and
Shklovskii [9]. In two dimensions, the DOS have a linear dependence on the energy around
the chemical potential.
Starting from random distributed values in each realization, the mean-field (MF) energies Ei are found by iterations according the to algorithm detailed in Chap .3 and a a similar
code to the one presented in Appendix J. The solution is then presented in Fig. 2.1.
As can be seen from Fig. 2.1, a linear gap is obtained which is consistent with Ref. [10].
The finite value of the density of states (DOS) at the chemical potential (µ = 0) is due to
the finite temperature. Furthermore, it has been shown, [11, 12, 13], that the DOS at the
12
Probability density
0.3
Figure 2.1: Single particle Density of states
of the Electron-glass.
0.2
The normalized histogram of average site energies Ei for halffilling and N = 10000 sites, Ref. [7]. The energies i where uniformly distributed in the interval [− W
,
2
is set to unity, W = 1, and temperature
0.1
W
2
]. Disorder energy
e2
rnn T
= 20. All pa-
rameters are in units of average nearest neighbours interaction
J=
e2
,
rnn
where rnn is the average nearest neighbours distance.
The sites are distributed uniformly on a square with periodic
−4
−3
−2
−1
0
1
2
E
3
4
boundary condition and averaged over 300 realizations.
chemical potential in 2D has a linear temperature dependence. For large enough values the
gap disappears completely.
2.1.2
Relaxation process
Following the same steps and notation as in [7]. The solution of rate equation, i.e. the
time dependence of the electron occupations, n(t), is obtained, using the MF energies found
numerically in the last section.
The rate equation, (see also 4.1) is given by:
dni X
γji − γij
=
dt
j6=i
(2.5)
where ni is the electron occupation at site i and γij the the transition rate between site i
and site j.
In the Electron-glass model γij takes the form of Miller-Abrahams transitions, [14],

−rij

Γ0ij ni (1 − nj ) e ξ [N (|∆E|) + 1] , Ei > Ej
γij =
(2.6)
−r

Γ0 ni (1 − nj ) e ξij N (|∆E|)
, Ei < Ej
ij
where N =
1
e|∆E|/T −1
is the phonon occupation, ∆E is the energy difference of sites i and
j, rij is the distance between sites i and j, ξ is the localization length of the electron, the
13
prefactor Γ0ij '
2π
|Mq |2 ν,
~
given ν as the phonon density of states and |Mq | as the electron-
phonon coupling strength for a phonon wave number q. Furthermore, Γ0ij depends on ∆E
and rij (see [15]) and approximated as constant, i.e. Γ0ij ∼ 1Hz ([7]).
Small enough perturbation causes the sites occupations, ni , to deviate not too far from
the local equilibrium values n0i . By Expanding the rate equation up to first order in the
small parameter δni = ni − n0i one obtains the linearised rates equation, ([7]):
dδni
dni X
=
=
Aij δnj
dt
dt
j
Aij =
γij0
n0j
X
X e2 γ 0 1
1
ik
;
A
=
−
Aij
−
−
ii
T
r
r
1 − n0j
ij
jk
i
k6=j,i
(2.7)
where the diagonal elements, Aii , are chosen to fulfil the requirement of particle number
P
conservation i Aij = 0 and the superscript X 0 indicates equilibrium values. Notice that
due to the term nj (1 − nj ) the matrix Aij is not symmetric which means that some of its
eigenvalues might be complex. Furthermore, given the Onsanger reciprocity relations, [16],
the rate matrix takes a symmetric form (real eigenvalues). This can be accomplished by
obtaining a linearised rate equation of conjugate variables, i.e., the external perturbation is
represented as a change of the chemical potential (rather than a change of site occupation
as in Eq. 2.7) and a response of the system which is the change of the site occupations,
P
dδni
= j Ãij δµi .
dt
Even though the interaction term (the second term of the off-diagonal element of Aij )
can posses large values at low temperatures, it gives a negligible contribution compared to
the exponential behaviour of the first term. This is confirmed also by a numerical calculation
(see Fig. 6 in [7]). Therefore, after neglecting the interaction term, the rate matrix Aij has
a definite positive off-diagonal elements. The resulted absolute value of the eigenvalues of
the off-diagonal matrix are bounded by unity (see Perron-Frobenius theorem [17]), with the
addition of vanishing sum of each column, the stability of the equilibrium point is guaranteed,
i.e., the real part of all the eigenvalues is negative or zero.
14
After neglecting the direct interactions term (as discussed above) and numerically
diagonalizing the first term of Aij given in Eq. 2.7, one obtains ([7]) the complex eigenvalues
(decay rates) with a negative real part that are distributed as
1
|λ|
(see Fig. 2.2).
0.1
1000
Probability density
Probability density
0.08
800
600
0.06
0.04
0.02
400
1
1.4
1.8
2.2
−25
2.6
−20
(a)
−15
−9 −8
−5
ln(-λ)
−4
x 10
−λ
(b)
Figure 2.2: Decay rates distribution.
Normalized histograms of the real part of the decay rates (originally done in Ref. [7]), obtained by numerical diagonalizing the
rate matrix Aij , 2.7, while neglecting the direct interactions term. This is done for N = 1000,
e2
rnn T
= 10 and
rnn
ξ
= 10.
The disorder energy and density of sites are the same as in Fig. 2.1. The graph is averaged over 1000 realizations. (a) normal
plot with
1
|λ|
fit, (b) log plot of the rates, showing the cut-offs values around the region with the most probable events in the
realizations
Generally, the rate matrix can take a simpler exponential form. For most mean-field
energy values the inequality T < |Ei |, |Ej |, |Ei − Ej | holds. Under this conditions Aij , is
approximated as follows:
Aij ∼ e−rij /ξ e−(|Ei |−|Ej |+|Ei −Ej |)/2T
(2.8)
Roughly speaking, the fact that the rate matrix has an exponential form with a broad
distributed power, Eq. 2.8, give rise to the
1
|λ|
distribution of the rates which leads to the
logarithmic decay. This occurs even without the presence of interactions and can be seen
by further simplifying the rate matrix, as conducted for a in [7]: for low density, rnn ξ,
the rate matrix can be approximated by exponential random euclidean matrix Aij ∼ e−rij /ξ ,
15
saving only the near neighbour terms one can obtain analytically the same distribution of
eigenvalues. However, the next section is devoted to the effect of interaction, therefore the
energy terms in the rate matrix Aij are needed.
The solution to the rate equation (Eq. 2.7) is straightforward. Applying a unitary
transformation that diagonalizes the rate matrix, Aij , results in a set of N decoupled differential equations:
X
d X
Uki δni =
Uki Aij Ujl−1 Ulj δnj
dt i
ijl
X
dδnk
=−
λl δkl δnl = −λk δnk
⇒
dt
l
(2.9)
⇒ δnk = ck e−λk t
where U is the unitary transformation, the deviation from equilibrium of the k’th eigenvector
P
P
is δnk (t) = i Uki δni (t), the diagonalization of the rate matrix matrix is i,j Uki Aij Ujl−1 =
−λk δkl , δkl is Kronecker delta and ck ≡ δnk (0) is the initial deviation of k’th eigenvector.
The total time dependence of the relaxation process can be taken as a summation of
all the occupation deviations (i.e. ”l = 1” norm):
|δn| =
X
|δnk (t)| =
k
X
|ck |e−λk t
(2.10)
k
For the continuous limit, substituting p(λ) = λ1 , the final form of relaxation is obtained:
Z
λmax
|δn| =
c(λ)p(λ)e
−λt
λmin
Z
λmax
dλ ' c
λmin
e−λt
dλ
λ
(2.11)
where the assumption is that the rate eigenvectors are excited roughly in a uniform probability c(λ) ' c, except for the eigenvector associated with the zero eigenvalue, which can not
be excited since the total particle number is conserved, [7]. One way to explain the uniform
distribution is that any two typical rate eigenvectors in the site location basis consists each
a set of linear combination of sites, the two sets of coefficients typically are distributed with
16
a similar broad distributions (given the uniform distribution of sites and relatively broad
distribution of energies).
For
1
λmax
<t<
1
λmin
Z
the integral can be estimated as Exponential-integral:
λmax
c
λmin
e−λt
dλ ' cEi(λmin t) ' c(γE − log(λmin t))
λ
where γE ≈ 0.577 is Euler constant and Ei (a) ≡
R∞
a
e−x
dx
x
(2.12)
is the Exponential integral. For
the evaluation of the Exponential integral see [18], p.252.
2.1.3
Effect of the disorder
As shown in the last section, given the relevant low site density, the interactions generally do
not affect the functional shape of the rates distribution. This can be seen in the exponential
form on the rates matrix, Eq. 2.8, the broad distributed power is caused by the distances part
rather than the energy. Nevertheless, the interactions affect the relaxation of the system,
as shown in [7]. This is carried out by varying directly the disorder W while setting the
interactions J = 1 and the ratio W/T = 10.
As can be seen numerically in Fig. 2.3, as the disorder grow the width of the DOS also
grows. As a result there is a shift of the rates to higher values, see Fig. 2.4. Meaning, as
W/J grows, the interactions have a smaller effect and the rates are faster.
Note that the main goal is to show the relation between the interaction (indirectly
though the variation of disorder) and the shift of the rate distribution. For this purpose the
cut-off range at which the rates are plotted, as seen in Fig. 2.2.b, can be chosen arbitrarily,
as long as the same range of rates is taken for all the cases, in order to see the full shift
of rates. In Fig. 2.4 the cut-off, for each disorder strength, was extracted by plotting the
distribution of the logarithm of the rates and choosing, symmetrically around the peak, the
rate values that are confined in a region of size unity. The same method can be viewed in
Fig. 2.2.
17
Probability density
0.2
Figure 2.3: Density of states of the EG
system for different disorder values.
0.1
The disorder is varied W = 1 (blue stars), W = 5 (green
squares) and W = 10 (red triangles) for W/T = 10, constant interactions J = 1 and N = 1000. The density of
−8
−6
−4
−2
0
2
4
6
8
sites is set as in Fig .2.1.
E
1200
Probability density
Probability density
0.12
800
400
0
0.5
1
1.5
2
2.5
3
3.5
0.04
−25
4
−20
−15
−10
−5
0
ln(−λ)
−3
x 10
−λ
0.08
(b)
(a)
Figure 2.4: Decay rates distribution of the EG system for different disorder values.
The distribution of eigenvalues of the rate matrix, Eq. 2.7, see Ref. [7]. The parameters are the same as Fig. 2.3.
(a) Rates in the regular scale with a fit to ∝
2.2
1
x
function, (b) Rates in the natural log scale.
Two-level systems
This chapter covers two known models: (1) The Standard TLS model (STM), originally
developed to explain thermodynamic anomalies that were raised in the 1970’s, [19], in noncrystalline solids at low temperatures (2) The TLS-glass model, an extension of the (STM)
which includes a TLS-TLS dipolar interactions. The addition of interactions is done naturally
to include some effects that the STM cannot predict such as a gap in the TLS DOS and a
18
correction to the duration of the relaxation process.
The section is divided as follows. Section 2.2.1 starts with a single TLS Hamiltonian
(derived from a general potential in appendix B) in which the energies and states are found
in terms of the asymmetry energy ∆ and tunnelling amplitude ∆0 , as seen in Fig. 2.5. The
following subsection, 2.2.2, addresses an ensemble of TLSs coupled to phonons, also known
as the standard TLS model (STM), which is used to obtain the observed temperature dependence of heat capacity and heat conductance of amorphous materials. Also, the obtained
temperature dependence explain the long standing question concerning the a priory unexpected differences between the amorphous materials and ordered solids. In Section 2.2.4 we
continue to the full Hamiltonian of interacting TLSs (TLS-glass), derived in appendix D,
which ultimately boils down to the addition of the TLS-TLS dipole interactions to the STM.
2.2.1
Single two-level system
Starting from the two-level system depicted in Fig. 2.5 we derive in Appendix B the single
TLS Hamiltonian:
Figure 2.5: A schematic representation of the Two-level
system.
V0 is the height of the barrier, d is the distance between the two minima and ~Ω is the
ground state energy of the harmonic oscillator (see Appendix B for details).

Htls = ∆S z + ∆0 S x =
19

1  ∆ ∆0 
2
∆0 −∆
(2.13)
1
Where ∆ is the asymmetry energy and ∆0 ≈ 2~Ωe− 2
q
2mV0
d
~2
≡ 2~Ωe−λ is the tunnelling
amplitude (see Eq. B.7).
Applying a rotation transformation, U , on Htls :
H̃ =U HU −1 U |ii =

cos θ
1
= 
2
− sin θ

E 0
1
= 
2
0 −E
sin θ


∆
∆0

cos θ − sin θ

∆0 −∆
sin θ

|ψ i
  +  = 1 E S̃z |ψi
2
|ψ− i
cos θ

cos θ


cos θ
sin θ
− sin θ cos θ


|Li
|Ri


(2.14)
where |Li and |Ri are the left and right state of the TLS in the real space,
|ψ+ i = cos θ|Li + sin θ|Ri, |ψ− i = − sin θ|Li + cos θ|Ri are the eigenstates, E± = ± 12 E =
p
± 12 ∆2 + ∆20 are the energies and tan(2θ) = ∆∆0
Detailed predictions of the TLS model are strongly dependent on the approximated
distributions of energies ∆ and ∆0 . For ∆, it is argued that negative and positive values are
equally likely and therefore are distributed symmetrically. Together with the need of a ”soft”
∆2
cut-off value, we assume a Gaussian distribution for the asymmetry energy, p(∆) = P0 e 2W 2 ,
where W 2 is the variance and the density of TLSs and P0 is a proportionality constant that
depends on the material. The tunnelling amplitude, ∆0 , on the other hand, depends on the
material and should be determined by experiment. Nevertheless, a qualitative result can be
obtained. Since ∆0 depends exponentially on λ, a small range of values of λ gives a wide
range of ∆0 . Thus, the distribution of λ can be approximated as uniform. This leads the
distribution of the tunnelling amplitude to take the form of: p(∆0 ) =
1
.
∆0
Assuming that ∆ and ∆0 are independent parameters, the resulted density of states
per unit volume for the non-interacting TLSs is:
p(∆, ∆0 ) =
20
P0 − ∆22
e 2σ
∆0
(2.15)
Usually in the literature [20, 21, 22, 23, 24], under the framework of the standard
tunnelling model, the relevant scale of temperatures is in a small region around 1K. In this
range the Gaussian distribution can be approximated as uniform. As a result the distribution
of the TLS system takes the form:
p(∆, ∆0 ) =
P0
∆0
(2.16)
Also, in the framework of the mean-field (MF) approximation done in chapter 3, calculating the density of states of the interacting TLSs system with the distribution 2.16, rather
than 2.15, gives the same qualitative results.
2.2.2
The standard tunnelling model
The Standard Tunnelling Model (STM) usually refers to a system of TLSs coupled to a
phonons with the given distribution (see 2.16). This simplified picture can describe many
properties of amorphous solids in low temperatures (below 1K). Some of these properties,
such as heat capacity and heat conductivity, are addressed at the end of this section.
At low temperatures only acoustic phonons with long wave lengths are excited. In this
case, where the wave length is much larger than the average inter-atomic spacing, one can
treat the amorphous material as a continuous media. This leads to interaction between TLS
and the local strain, in the deformation potential approximation.
Starting from the expansion of the TLS Hamiltonian (2.13) in powers of local strain
field, one can show how the TLSs interact with phonons:
X X X ∂ n Htls µν
i
n
H=
µν n (uis )
∂(u
)
is
n
i µνs
(2.17)
∂uµ
is
where Hitls is the TLS Hamiltonian, 2.13, at the site i, uµν
is ≡ ∂xν is the local strain field,
uµis = uµs is a component of the local displacement vector at the location of site i, µ, ν =
x=xi
1, 2, 3 are Cartesian indices and s = l, t1, t2 indicates the polarisation, l for longitudinal
orientation and t1, t2 for transverse orientations.
21
Assuming the coupling to the local strain, uµν
is , is small and expanding up to first order
X X ∂Htls i
uµν
µν is
∂u
x=x
i
is
i
i µνs
X
X X ∂∆i ∂∆0i z
x
z
x
=
(∆i Si − ∆0i Si ) +
Si − µν Si uµν
µν is
∂uis x=xi
∂uis x=xi
i
i µνs
H=
X
Hitls +
(2.18)
As can be seen, the TLS Hamiltonian 2.13 is obtained, with an addition of a linear coupling
of TLS to the local strain field, uµν
is .
Considering a small coupling between the tunnelling amplitude (∆0 ) to the strain field
[25, 26, 22], Eq. 2.18 will take the form:
H=
X
X X µν z µν
∆i Siz − ∆i0 Six +
γis Si uis
i
where
µν
γis
=
∂∆i ∂uµν
is
i
(2.19)
µνs
is called the deformation potential.
x=xi
Using the golden rule, the obtained transition rate of TLS at site j due to coupling to
phonons is (see Appendix C):
τi−1
≡ λi =
γil2 2γit2
+ 5
c5l
ct
∆20 Ei
coth
2πρ~4
Ei
2T
=
ai ∆20i Ei
coth
Ei
2T
(2.20)
where nqs is the occupation number of phonon with momentum q and polarization s, Ei =
p
∆2i +∆20i is the
energy difference between the upper and lower states of the TLS at site i,
hai i ≡
2
γil
γ2
+ it5
c5
2ct
l
2πρ~4
≈ 108 K −3 s−1 is the averaged prefactor of the transition rate (see [27]), γis
is the coupling constant averaged over all orientations µ, ν and Boltzmann constant kB = 1.
2.2.3
Thermodynamic properties of TLSs
The thermodynamic properties, such as specific heat and thermal conductivity, of amorphous and disordered solids, at temperatures under 1K, have shown different dependence on
22
temperature than their ordered counterparts. In the framework of the STM, the heat capacity turns out to be proportional to the temperature C ∝ T , and the thermal conductivity
κ ∝ T 2 , see experiment [19], rather than C, κ ∝ T 3 dependence, as Debye’s model predicts.
The STM achieves these temperature dependencies under a simple consideration, which is:
the TLSs are the major excitations in this temperature range.
Following [21], the two most known results of the model are presented. For TLS with
energy difference E and a given partition function Z = 2 cosh 2kEB T we can calculate
the specific heat and the thermal conductivity of a system of TLSs with density of states
P (E) = P0 . The total specific heat per unit volume is:
Z ∞
Z ∞
π2 2
P0 (E/4kB T )2 cosh−2 (E/2kB T ) dE = kB
P (E)C1 dE =
Cv =
P0 T
6
0
0
(2.21)
where C1 is the specific heat at constant volume of a single TLS. The upper energy cut-off is
evaluated as infinity since for low temperatures (than the typical TLS energy) the integrand
is exponentially small. The thermal conductivity:
1X
κ=
3 s
Z
0
∞
3
ρkB
Cs (ω)vs lph dω =
6π~2
X vs
P0 γs
s
!
T2
(2.22)
where ω is the phonon frequency, lph is the phonon mean free path (see Eq. C.11
derived at the end of Appendix C), γs is a typical value of the deformation potential, s is the
polarization of the phonon and vs is the speed of sound which we denote as cs from here on.
In both calculation of Cv and κ the dependence on ∆0 is neglected, leading approximately
to uniform distribution of energies, P (E) = P0 .
After briefly showing the significance of the STM which explains the equilibrium thermodynamic properties, the next section addresses the addition of the TLS-TLS interactions
to the STM and how it might effect the relaxation process.
23
2.2.4
The TLS-glass model
Starting from the Standard Tunnelling Model (STM), shown in Eq. 2.19. After using two
transformations, approximating a weak TLS-phonon coupling and summing over the bath
degrees of freedom, we obtain the model of interest, see Eq. D.6 in Appendix D for derivation.
The final form of the Hamiltonian is:
H = Htg =
X
(∆i Siz + ∆0i Six ) −
i
X uij
ij
3
rij
Siz Sjz +
X X ∆0i gik i
k
~ωk
a†−k − ak iSiy
(2.23)
≡ Htg + Htls−ph
where the tilde sign denoting representation in the diagonal TLS base (see Eq. 2.14) and the
pseudospin operator is proportional to the Pauli matrices, S̃iz = 12 σiz .
The Hamiltonian contains three parts, each plays a similar role as its counterpart in
the Electron-glass model, Eq. 2.1. Using spin terminology, the first term is a random field
that accounts for disorder W , the second term is the interaction term (Ising spin-glass) and
the third term is the coupling between spin system and a phonon bath, which is responsible
for the transitions.
As explained in Appendix D, by treating H̃tls−ph as perturbation we obtained the
transition rate which takes the same form as the rate in the STM, 2.20. Therefore, for all
practical purposes we treat the TLS-glass model (TG), 2.23, as the STM with the addition
of Ising spin-glass interactions.
Usually, given the random nature of the TLS orientations, the average over angles
q
between a pair of interacting TLSs is zero huij i = 0 with standard deviation U0 ≡ hu2ij i ≈
γ2
,
ρc2
where the average deformation potential γ = hγs is is defined in 2.19 and c = hcs is is the
average speed of sound. Since uij is symmetric and falls-off rapidly above some typical value
(U0 ) its distribution is modelled as Gaussian.
24
The interactions are also related to the known unitless universal constant,
χ = P0 U0 =
P0 γ 2
ρc2
(2.24)
which ranges roughly between 10−4 to 10−3 in all known glasses, and has been shown to
be almost material independent, even when the parameters P0 , ρ, c and γ can vary greatly.
Furthermore, χ ≈ λ/lph where lph and λ are the mean free path (see Eq. 2.20) and wave
length of thermal phonons, which are responsible for heat conductivity. This equality was
one of the first clues which showed that the interactions must be playing an important role
in the low temperature regime of amorphous and disordered solids (see [28, 29]).
Just as in the Electron-glass case, in the next section we obtain a gap in the TLS
density of states (DOS) caused by the TLS-TLS interactions, and also address the interplay
between the relaxation process and the interactions, disorder and system size.
The validity and relevance of the model is determined mainly by two conditions:
1. TLS transitions due to phonons is the major mechanism of relaxation.
2. The only excitations in the system are TLSs and phonons.
Conditions (1) sets a lower bound for the temperature with a cross-over temperature T 0 ≈
q
3 2k~aB χ3/2 , in which, TLS-TLS induced transitions take a major part (see [30]), where a is
a typical constant related to the TLS-phonon coupling. The range of T 0 is roughly between
the values ∼ 10mK − 0.3K (given the range of χ).
Condition (2) sets the upper bound. In the case of higher temperature range (around
1K − 3K) excitations come into play such as Fractons, which could interact with TLSs [31],
and also strongly interacting TLSs [32].
As in the Electron-glass, we evaluate the energy spectrum of the unperturbed Hamiltonian Htg = Htls + Htls−tls and use a weak coupling to phonons Htls−ph , as the mechanism
of transitions. Also, under the same consideration as before, given the long range nature of
25
Htg , the MF model could give a rather good qualitative estimate of the energy spectrum, in
a relatively low computational cost.
TLSs Mean-Field model
Starting from the TLS-glass Hamiltonian in the location basis:
Htg =
X
(∆i Siz + ∆0i Six ) −
i
X uij
ij
3
rij
Siz Sjz
(2.25)
We obtain the MF equations by applying a variational derivative with respect to Siz ,
∆0i =
X uij
δHtg
= ∆i −
Sz
z
3 j
δSi
r
j6=i ij
(2.26)
where ∆i is the asymmetry energy of TLS i and is normal distributed with variance W which
is the disorder, ∆0i is the new asymmetry energy, the interaction strength parameter is defined
3
which is controlled by
as the average dipolar interaction of nearest neighbours J = uij /rnn
scaling U0 which is the standard deviation of uij (setting rnn = 1 in the numerical calculations
we obtain J = U0 ), and number of sites N (system size).
After thermal averaging (see A.4), the obtained self-consistent equations are:
1 X uij
Ej
0
∆i = ∆i +
tanh
3
2 j6=i rij
2T
(2.27)
where we used hSiz i = 21 hσiz i. Accordingly, the single TLS shifted Hamiltonian is
HT0 LS =
X
(∆0i Siz + ∆0i Six )
(2.28)
i
and the energy difference of the i’th TLS is
q
Ei = sgn(∆ ) ∆0 2i + ∆20i
0
(2.29)
since the interactions can give negative contributions, we use the sign function in order to
maintain half of the solutions of the self-consistent equations (negative branch of the DOS).
26
Note that the TLS-Phonon interaction term in Eq. 2.23 is invariant under the rotation
transformation U 0 , which diagonalizes the mean-filed single TLS Hamiltonian, Eq. 2.28.
Meaning, U 0 S y U 0−1 = S y where U 0 is the same rotation transformation an in Eq. 2.14 only
with the new mean-field asymmetry energies, tan(2θ) =
27
∆0
.
∆0
Chapter 3
Numerical solution of the TLSs
self-consistent equations
Previous analytical work (see Appendix. E and [10]) have shown that the DOS of the TLS
system in 3D has a logarithmic gap which results from the dipole interactions. In this section
we present the numerical solution of the self-consistent equations 2.26, which gives the same
logarithmic dependence, and in addition the behaviour for larger energy values far from the
gap. By doing so, two goals are accomplished: (1) we verify the validity of the mean-field
approximation by comparing the numerical solution to the DOS obtained in the Monte-Carlo
method, [33, 34], and (2) we gain an ability to control the parameters of the model such as
disorder, interaction strength and system size.
By following the algorithm presented in [35], we calculate numerically the solution of
the Mean-field (MF) equations for finite temperature which is given by Eq. 2.26. We start
by setting the initial values for the TLS asymmetry energy vector ∆0i (see 2.27) to uniform
distribution, and then perform a three steps procedure 1. Choosing a TLS randomly and update its asymmetry energy according the MF equations, 2.27. We then repeat this procedure to update the entire energy vector (∆0i ).
28
2. Repeating step 1 until the maximum difference of TLS’s energy between successive iterations, ∆Emax , is smaller than the predetermined minimum value ∆E0 . This iterative
procedure converges relatively fast (in a few iterations). ∆Emax usually has characteristic values much smaller than the chosen minimum value, therefore the resulted
energies showed no dependence on ∆E0 for a wide range of values (up to ∆E0 = 10−8 ,
in units of near neighbour average interaction value J).
3. Repeating steps 1 and 2 while heating and then cooling back the system in 10 equally
spaced temperature values. By doing so the system relaxes to lower energetic configuration, see [35].
A Normalized histogram (Density of states) is then made from the final energy difference
p
vector of each TLS, Ei = sgn(∆0 ) ∆0 2i + ∆20i , as can be seen in Fig. 3.1. The code is
presented in Appendix J.
Probability density
0.04
Figure 3.1: Density of states of the TLS
system.
0.03
The normalized histogram of average TLS energies Ei =
q
sgn(∆0 ) ∆0 2i + ∆20i obtained by solving the self consistent equations 2.26, for N=10000, uij is normal dis0.02
tributed with standard deviation, U0 = 1, the disorder
energy W = 1, the average nearest neighbours interaction J =
0.01
U0
3
rnn
where rnn is the average distance of
nearest-neighbours,
J
T
= 20, and
W
J
= 1. The sites
where distributed uniformly in a cube. The graph is the
−20
−15
−10
−5
0
5
10
15
20
E
average of 300 realizations.
As can be seen from Fig. 3.1, a logarithmic gap is obtained which is consistent with
Eq. E.7 given for small interaction and in Ref. [10] for a broader range of interaction values.
As in the electron glass DOS, the logarithmic gap disappears completely for large enough
temperatures. Note that ∆0 do not evolve with the iterations since their coupling to phonons
is neglected (2.18, 2.19).
29
In all the numerical calculations, the parameters of Eq. 2.26 are measured in units
of average nearest-neighbours interaction J, χ = 10−3 and homogeneous spatial density of
sites. In addition, we initially set ∆ and ∆0 by the distribution written in Eq. 2.15, where
the disorder energy W is defined as the standard deviation of ∆, and ∆0 is in the range
[10−7 , 10−1 ], (see [36, 22] for cut-off values of ∆0 ). Also, no dependence on the initial values
of ∆0 , (the mean-field asymmetry energy, 2.27), have been observed for a wide range of
values, i.e., for ∆0 ∈ [−1, 1] up to [−100, 100].
30
Chapter 4
Relaxation process of TLSs
As mentioned in the introduction, like most systems with disorder and interactions, there
is a glassy behaviour leading to slow relaxations. In particular, we address the relaxation
to local minima of the free energy (local-equilibrium) of the excited interacting TLSs given
the TLS-glass model. Meaning, a specific state of the system, one of many states which
lay close in energy, that minimizes the Two-TLS excitation (given in Eq. E.1). Also, as in
the EG model, we treat the phonon degrees of freedom as a bath. This picture explains
the slow relaxation in short time scales, roughly smaller than the inverse of the minimal
relaxation rate, i.e., t < 1/λmin , where usually in most systems 1/λmin is of the order of
hours. For the longer relaxation time scales there are two typical dominating mechanisms
that are not addressed in this work: (1) TLSs with non-adiabatic initial excitation (See [36])
(2) transitions between many-particle metastable-states of the system.
In sub section 4.1 we start with a short discussion on the general master-equation and
derive its explicit form for the TLS system. Then, we expand around the equilibrium point
up to second order (quadratic term). In Sec. 4.2 we neglect the quadratic term and show
the results for the linear master equation, obtaining the logarithmic relaxation, just as in
the electron-glass case.
31
4.1
Rate equation
Starting from the general Pauli Master-equation
dpi (t) X
=
ωij pj (t) − ωji pi (t)
dt
j6=i
(4.1)
where we denote pi (t) as the occupation of state i at time t which can take the values
0 ≤ pi (t) ≤ 1, and ωij is transition rate from state j to state i.
P
Eq. 4.1 already contains probability conservation, i.e.
i pi (t) = const, which for the
P i
TLS system m pm (t) = pi1 + pi2 = 1 where i indicates the site of the TLS and p1 , p2 are
the occupations of the lower and upper levels of the TLS, respectively. Particularly, for the
P
electron glass the probability conservation takes the form of i pi (t) = M where M is the
number of electrons.
In the TLS model, only transitions between the two states of a given TLS are possible,
therefore, for the i’th TLS occupations, Eq. 4.1 takes the form of two coupled rate equations:
dpi1 (t)
i i
i i
= ω−
p2 (t) − ω+
p1 (t)
dt
dpi2 (t)
i i
i i
= ω+
p1 (t) − ω−
p2 (t)
dt
(4.2)
i
i
are TLS-phonon upward and downward transition rates respectively, see
and ω−
where ω+
Eqs. C.8 and C.9.
Finally, we can reduce Eq. 4.2 to one parameter. In the pseudospin representation, i.e.
σi ≡ hσiz iT = pi2 − pi1 and using pi1 + pi2 = 1, the rate equation will take the form:
dσi
1
1
= −2Ai σi ni +
+
= −λi σi − Ai
dt
2
2
where ni =
1
eβEi −1
(4.3)
i
i
is the equilibrium phonon occupation at energy Ei , β = T1 , λi ≡ w−
+w+
=
Ai coth (Ei /2T ) TLS-phonon relaxation rate and Ai = ai ∆20i Ei (see 2.20 for derivation).
σ ) and
Eq. 4.3 has a simple form but has a hidden complexity, the relaxation rate λi (σ
σ ) at site i depends on the interactions and also on the out-of-equilibrium
the constant Ai (σ
32
occupation values of all the TLSs in the system. Therefore, one of the approaches is to
expand the rate equation near the equilibrium point and obtain the r.h.s of Eq. 4.3 in terms
of the known equilibrium values, A.4.
For sufficiently small deviations from equilibrium values of the TLS occupations, δσ 1, caused by the initial perturbation, we can expand the r.h.s. of 4.3 near the equilibrium
point up to 2nd order, see Appendix F. The resulted rate equation for TLS i is:
dσi
∆0i X uij
' −Ai (2ni + 1) δσi + 2Ai βni (ni + 1)
δσj δσi
3
dt
Ei j6=i rij
(4.4)
and the general solution of Eq. 4.4 is:
δσi (t) = ci e−λi t e
where fij ≡ 12 Ai β sinh−2
Ei
2T
i uij
3
Ei rij
∆
P
j6=i fij
Rt
0
δσj (t0 )dt0
(4.5)
and ci ≡ δσi (0) is the initial deviation of TLS i at the
moment the external strain driving force has stopped.
As can be seen in Eq. 4.4, only the second order (quadratic term) in the rate equation
depends explicitly on the interactions, where in the electron-glass a similar dependence is
introduced already in the first order. This is a result of probability conservation that takes
place in each TLS separately (p1 + p2 = 1), while for the electron-glass the total probability
P
of all the sites is conserved i pi = const. The direct interaction terms presented in the rate
equation, Eq. 4.4, are discussed in Appendix G.
4.2
The Linear term
As mentioned in the previous section, the deviations from equilibrium of spin occupations are
small, this means that the linear term gives the main contribution to the relaxation dynamics.
By neglecting the second term in the exponential power of Eq. 4.5, which corresponds to
neglecting the quadratic term in the rate equation, Eq. 4.4, we obtain the solution of the
33
linear rate equation:
δσi (t) = ci e−λi t
(4.6)
where λi is the transitions rate which depends on the MF energies 2.29.
Comparing the linear rate equation in the EG model, Eq. 2.7, to the linear rate equation
of TLS model (Eq. 4.5 without the second power) we see that in the TLS model the rate matrix in the pseudospin representation is diagonal. This results from probability conservation
for each TLS separately.
As before, having an exponential dependence with a broad power results in
1
|λ|
dis-
tribution of decay rates, as can be seen in Fig. 4.1. Repeating the calculation done in
Eqs. 2.10-2.12, we obtain the total relaxation of the TLS system:
Z λmax
X
e−λt
c(λ)
|δσ| =
|δσi | '
' c (γE − log (λmin t))
λ
λ
min
i
(4.7)
As in the EG model, the linear part of the rate equation is responsible for the slow
logarithmic relaxation, however with several distinctions. First, the minimum rate λmin has
a different dependence on the MF energies which in turn depends on the dipolar rather
than Coulomb interactions. Secondly, in the pseudospin representation the rate-matrix is
already diagonal, which means that the sites basis is also the eigen-basis of the rate matrix.
Also, The uniform probability of the initial TLS values (c(λ) = c) is achieved under more
strict considerations (see Appendix H). Furthermore, the quadratic term gives only small
corrections (for small deviations), and is addressed as perturbation in appendix G.
Comparing the rate distribution of the two models, the TG (Fig. 4.1) and the EG
(Fig. 2.2), we can see that the former has a wider cut-off range. This can be explained as
follows: In the TG model, following the conventional parametrization, we set the tunnelling
amplitude as a parameter which is distributed in a given wide region ∆0 ∈ [10−7 , 10−1 ]
(see [23]), this results in a wide rates distribution. In the EG model, on the other hand,
the tunnelling amplitude is given by the distances exponential, e−rij /ξ (2.6) and leads to a
34
4
x 10
0.04
Probability density
Probability density
3
2
0.03
0.02
0.01
1
2
3
4
5
6
−15
−10
−5
0
5
10
15
ln(−λ)
−6
x 10
−λ
(a)
(b)
Figure 4.1: Decay rates distribution.
The decay rate distribution of rates λi presented in Eq. 4.6 for N = 10000 and
J
T
= 10. The disorder energy, tunnelling
splinting and interaction strength are the same as Fig. 3.1. The graph is averaged over 1000 realizations.
(a) Rates in the regular scale with a fit to ∝
1
x
function. Only the lower cut-off is shown. The distribution extends up to the
upper cut-off value 106 Hz (b) Rates in the natural log scale, showing the full region of values. Just as in Fig .2.2, the cut-off
values are defining the region of the rates bulk which distributes as
1
λ
and is the flat region in the logarithmic plot.
smaller cut-off region. This is perhaps because the largest matrix elements (which gives the
main contribution to the eigenvalues) corresponds to all the pairs of sites that are close in
distance. Thus, a reasonable measure for the rates cut-off region of the matrix e−rij /ξ , in
the low density limit, ξ rnn , is the width of the nearest-neighbours distance distribution
V ar(rnn ) and possibly it would be relevant to add also the next-nearest neighbours variance
to improve the accuracy of the measure. In order to obtain a broader rates distribution in
the EG model we could change the system’s dimensionality or site density in such a way
that the variance of rnn will increase.
Having a simple expression for the rates in the TG model, under a reasonable approximation, we can derive the 1/λ distribution analytically. Given the joint distribution
35
p(∆0 , ∆0 ), (see Eq. 2.15), we can find the distribution of rates by change of variables:
0
0 p(λ, ∆ ) = p(∆0 , ∆ ) where we substitute p(∆0 ) =
dependant on ∆0 ) and
dp(∆0 )
d∆0
dp(∆0 )
,
d∆0
∂∆0
∂∆0
1
,
∆0
∂∆0
∂λ
∂∆0
∂∆0
∂∆0
∂λ
∂∆0
∂∆0
p(∆0 ) ∂∆0 =
∆0 ∂λ (4.8)
approximating p(∆0 , ∆0 ) ' p(∆0 )p(∆0 ) (for ∆0 weakly
= 0. Also, assuming p(∆0 ) varies much slower than p(∆0 ), i.e.
which is valid in most regions in the DOS, we can approximate p(∆0 ) as
constant. After substituting the rate Eq. 2.20, we obtain:
p(λ, ∆0 ) '
where C =
p(∆0 ) √ 2 ∆0 =P
C
λ
(4.9)
and P is some point on p(∆0 ).
As in the Electron glass model, the 1/λ shape of the rates distribution does not depends
on the interactions. Therefore, the next chapter addresses question of how the interactions
affect the rates.
36
Chapter 5
Effect of the interactions disorder and
system size
5.1
TLS-Glass model
In this chapter we apply schemes of parameter variation for the TLS model. One intriguing
outcome is obtained; the interactions speed up the relaxation process rather than slow it
down as obtained for the EG model. A possible explanation can be given through the rates
Ei
term, λi = ai ∆20i Ei coth 2T
. The rates depend monotonically on the values of MF energies.
Meaning, as long as the inequality Ei /T & 1 holds, increase of the energies leads to an
increase of the rates. However, this argument does not hold for Ei /T . 1 since rates do not
depend on the energies but rather on the temperature, i.e. λi ≈ 2ai ∆20i T .
The different schemes are:
• Varying the disorder, W , for constant interactions J and ratio W/T , where T is the
temperature. For higher disorder the DOS broadens and the rates are shifted to higher
values (see Fig. 5.1 and Fig. 5.2.)
• Varying the size of the ratio J/W while holding the sum of the variances constant
37
,W 2 + J 2 = 2. This is done in order to change the strength of the interactions while
preserving the energy variance, (see Eq. 5.1). For decreasing interaction strength the
rates are shifted to lower values (see Fig. 5.3 and Fig. 5.4).
• Varying the number of sites, N , while keeping the density of sites constant. For
increasing N the DOS is narrower (including the gap) and the resulting rates are
shifted to lower values, see Fig. 5.5 and Fig. 5.6
Note that unlike in the EG model, in the TG model the shift of the rates is expressed
only through the lower cut-off while the upper cut-off is constant. This is due to the fact
that the rates distribution is much wider in the TG model, as can be seen by comparing
the plateau regions in the rates log scale distribution of the two models, see for example
Fig. 2.2 for the EG model and Fig. 4.1 for the TG model. Nevertheless, what determines
the relaxation time scale are the low rates and therefore the lower cut-off is the relevant one.
Probability density
0.035
0.025
Figure 5.1: Density of states for different
disorder values.
0.015
The disorder is varied W = 1 (blue stars), W = 10 (green
0.005
squares) and W = 100 (red triangles) for W/T = 10, constant interactions J = 1 and N = 1000. The density of
−20
−15
−10
−5
0
5
10
15
20
sites are the same as in Fig. 3.1.
E
By examining the DOS, particularly at the region of the gap and not too far from it,
we explain, at least qualitatively, the change of the typical small energies, and as a result
(as explained above) also the shift of the low rates (i.e. lower cut-off which governs the
relaxation). Also, through the variance of the DOS, we can explain the typical change in
the energies. The variance does not contain information about the height of and slope of the
38
0.04
Probability density
Probability density
7000
5000
3000
0.03
0.02
0.01
1000
2
4
6
8
−15
−10
−5
0
5
10
15
ln(−λ)
−5
x 10
−λ
(a)
(b)
Figure 5.2: Decay rates distribution for different disorder values.
Decay rates for W = 1 (blue stars), W = 10 (green squares) and W = 100 (red triangles) for W/T = 10, constant interactions
J = 1 and N = 1000. The density of sites are the same as in Fig. 3.1.
(a) Rates distribution in the regular scale (notice that the intermediate case of W = 10 is close to the value of W = 1 and
therefore was discarded), (b) Rate distribution in natural log scale.
Prbability density
0.14
0.1
Figure 5.3: Density of states for different
interaction strengths.
The disorder interactions ratio number is varied. W/J = 1
0.06
(blue asterisks), W/J = 3.1 (green squares), and W/J =
9.8 (red triangles) for W 2 + J 2 = 2 and N = 1000. The
0.02
temperature and density of sites are the same as in Fig. 3.1.
Notice how the gap gets narrower and higher for greater
−20
−15
−10
−5
0
5
10
15
20
values of W/J.
E
dipole-gap, which is also responsible for the shift of the low rates. Nevertheless, the energy
variance can give an additional perspective about how the parameters such as disorder W ,
interaction J and system size N affect the width of DOS. Thus, the approximated variance
39
4
x 10
0.04
Probability density
Probability density
8
6
4
2
0.4
0.8
1.2
1.6
0.03
0.02
0.01
−15
2
−10
−5
0
5
10
15
ln(−λ)
−5
x 10
−λ
(a)
(b)
Figure 5.4: Decay rates distribution for different interaction strengths.
Decay rates for J/W = 1 (blue asterisks), J/W = 1/3.1 (green squares), and J/W = 9.8 (red squares). The rest of the
parameters are the same as in Fig. 5.3.
(a) Rates in the regular scale, (b) Rate in the natural log scale.
Probability density
0.04
0.03
Figure 5.5: Density of states for different number of sites.
The number of N is varied, 100 (red triangles), 1000
0.02
(green squares) and 10000 (blue asterisks). The temperature, asymmetry energy, disorder energy and den0.01
sity of sites are the same as in are the same as in Fig. 3.1.
Notice how the DOS getting wider and the gap getting
−20
−15
−10
−5
0
5
E
10
15
20
shorter as the number of sites getting smaller.
of TLS i is:
X U02 2 Ei2 ≈ W 2 + ∆20i +
Sj
6
r
ij
j
(5.1)
where rij is the distance between sites i and j and U0 is the interaction (see, 2.23) and is equal
to near neighbour average interaction J (for rnn = 1). For changing the disorder scheme,
the energies are shifted to higher values for higher disorder, this can be seen in the linear
40
4
x 10
0.04
Probability density
Probability density
3
2
0.03
0.02
0.01
1
2
3
4
−λ
5
−15
6
(a)
−10
−5
0
5
10
15
ln(−λ)
−6
x 10
(b)
Figure 5.6: Decay rates distribution for different site number.
Decay rates for N = 100 (red triangles), N = 1000 (green squares) and N = 10000 (blue asterisks). Disorder, interactions,
temperature and TLS density are the same as Fig. 4.1.
(a) Rates in the regular scale, (b) Rate in the natural log scale.
dependence between the variance of the site energy and the disorder, Eq. 5.1. Meaning,
higher disorder gives wider DOS and higher rates. Also, the same effect is obtained for
constant temperature, see Fig. I.1 and Fig. I.2. The same logic holds for changing the
interaction scheme. Each interaction term in Eq. 5.1 gives a smaller contribution which
overcomes the broadening of the DOS, caused by increasing disorder, also, naturally the
gap height increases which suggest that there are additional low energy TLSs. Therefore,
a typical low energy values decreases and the rates shift to lower values. Also, typically a
greater change occurs in the DOS when only the interactions are varied for constant disorder
(see Fig. I.3 and Fig. I.4). Finally, for increasing the number of sites, each TLS gains an extra
interaction terms, which have smaller energy values since the density od TLSs is constant.
Despite this, the resulted DOS has a narrower gap and width (and not wider as expected
from Eq.5.1). Therefore, the overall effect is a shift of the rates to lower values.
Note that we discuss the indirect influence of the interactions through the MF energies,
while neglecting the quadratic term in the rate equation 4.4. By treating the quadratic
41
term as perturbation we obtain corrections to the rates given in the linear approximation,
see appendix. G, thus, introducing a direct dependence on the interactions and a possible
additional effects to the relaxation process.
5.2
Electron-Glass model
In order to compare between the models, the same parameter-varying schemes are conducted
for the Electron-glass model. As in the previous section, we study how the change of the
DOS affects the rates. This is done by studying the typical change of rate matrix element,
Aij = e−
rij
ξ
(n0i /n0j )N (∆E) as defined in Eq. 2.8, and use it as a measure for the shift of the
rates, i.e., a typical increase of the matrix elements will cause an increase of the eigenvalues
(rates) and vice versa.
The different schemes are:
• Varying the disorder, W , for constant interactions J and ratio W/T , where T is the
temperature. Originally done in Ref. [7] and repeated in Subsec. 2.1.3. For higher
disorder the DOS broadens and the rates are shifted to higher values (see Fig. 2.3 and
Fig. 2.4).
• Varying the size of the ratio J/W while holding the sum of the variances constant
,W 2 + J 2 = 2. This is done in order to change the strength of the interactions while
preserving the energy variance, (see Eq. 5.2). For decreasing interactions the rates are
shifted to higher values, which is the opposite from the TLS case (see Fig. 5.7 and
Fig. 5.8.
• Varying the number of sites, N , while keeping the density of sites constant. For
increasing N the DOS is narrower (including the gap) and the rates are shifted to
lower values (see Fig. 5.9 and Fig. 5.10).
42
Probability density
0.7
0.5
Figure 5.7: Density of states for different
interaction values.
0.3
The ratio of disorder to interaction is varied W/J = 1
√
√
(blue stars), W/J = 17 (green squares) and W/J = 97
0.1
(red triangles) for constant sum of variances W 2 +J 2 = 2,
−4
−3
−2
−1
0
1
2
3
temperature T = 0.1, and N = 1000. The density of sites
4
E
is set as in Fig .2.1.
1200
0.12
Probability density
Probability density
1000
800
600
400
0.08
0.04
200
0
0
0.5
1
1.5
2
−25
2.5
−20
−15
(a)
−10
−5
0
ln(−λ)
−3
x 10
−λ
(b)
Figure 5.8: Decay rates distribution for different interaction values.
The distribution of eigenvalues of the rate matrix Eq. 2.7 the rest of the parameters are the same as Fig. 5.7.
(a) Rates in the regular scale with a fit to ∝
1
x
function, (b) Rates in the natural log scale.
By examining the DOS, we explain the change of the site energies leading to the change
of the typical rate matrix element which ultimately is observed by the shift of the rates. In
turn, the change of the site energies is explained using the qualitative shape of the DOS and
the energies variance (see the discussion at the end of section 5.1 about the contribution of
the variance). Note though, that in the EG model, given the transformation from N 2 − N
43
Probability density
0.3
Figure 5.9: Density of states for Different site number.
0.2
The number of sites N is varied, 10 (red triangles), 100
(green squares), 1000 (blue asterisks) and
e2
rnn T
= 10.
The temperature, disorder energy, density of sites and
0.1
localization length are the same as in Fig. 2.2.
Notice how the DOS is wider and the gap is shorter
for smaller number of sites.
The same quantitative
behaviour can be seen for changing the interactions
−4
−3
−2
−1
0
1
2
3
4
strength to smaller values.
E
0.14
800
Probability D ensity
Probability density
1000
600
400
0.1
0.06
200
0.02
1
2
−25
3
−20
−15
−10
−5
0
ln(−λ)
−3
x 10
−λ
(b)
(a)
Figure 5.10: Decay rates distribution for different site number.
Rate distribution of the decay matrix 2.7, for N = 10 (red triangles), N = 100 (green squares)and N = 1000 blue asterisks. i
energies, disorder energy, density of sites and localization length are the same as in Fig. 2.2. The solid line is a fit of 1/x curve.
(a) Rates in the regular scale with a fit to ∝
1
x
function, (b) Rates in the natural log scale.
transition rates to N relaxation rates (i.e. diagonalization procedure), the typical matrix
element can not generally be evaluated by a simple limit (as conducted for a specific cases
below, Eq. 5.3).
The approximated variance of site i is:
X e4 n2j
Ei2 ≈ W 2 +
2
r
ij
j
44
(5.2)
where W is the disorder and rij is the distance between sites i and j. For changing the
parameter J/W scheme, Each interaction term in Eq. 5.2 gives a smaller contribution which
overcomes the broadening of the DOS, caused by increasing disorder. Also the gap height
increases which is a by-product of the periodic boundary conditions (for open boundary
conditions the height increases for smaller number of sites). Unlike in the TLS case, in the
EG model the opposite effect is obtained, i.e., the rates shift to higher values. A simple
approximation of the matrix element Aij can be taken in this case. For any two sites
close in energy (which is the majority case, as can be seen from the DOS in Fig. 5.7), the
approximated energy part of the matrix element will be:


e−|∆E|/2T , |∆E| > 2T
aij ∼

 T
, |∆E| < 2T
|∆E|
(5.3)
As can be seen from Eq. 5.3, the largest matrix elements, which are the ones that are
generally close in energy, show an increase for decreasing energy difference. Also a greater
change occurs in the DOS when only the interactions are varied (see Fig. I.7 and Fig. I.8. For
increasing the number of sites scheme. The resulted DOS has a narrower gap and the overall
effect is a shift of the rates to lower values. This effect can be explained by the fact that
changing the number of sites affects two different terms in the rate matrix. The first term
depends only on the distances between the sites, e−rij /ξ , and the second term depends only on
the site energies, which can be also approximated as in Eq. 5.3. The two terms compete and
the observed overall effect is dominated by the first term. This is confirmed numerically, see
Fig. I.9. Note that the same qualitative effect (resulted from changing the system’s size) is
obtained also in the TLS model, but from different origins: (1) the tunnelling amplitude, ∆0 ,
is distributed independently from the site-distribution of TLSs (which is used to calculate the
TLS DOS), (2) the rates have an opposite dependence on the energies. Finally, for increasing
the disorder, the rates shifts to higher values which may infer that the typical matrix element
increases also. As stated above, given the diagonalization procedure, a simple qualitative
45
argument can not be given, since the typical matrix element can not be obtained by a simple
limit as done for the other two cases. Furthermore, the same effect is obtained for constant
temperature, which shows constant gap height (see Fig. I.5 and Fig. I.6).
46
Chapter 6
Summary and conclusions
In this work we examined the relaxation to local-equilibrium and its dependence on interactions and system’s size, for both the TLS-glass and the Electron-glass models. This is
achieved by solving the rate equation and obtaining the time evolution of the localized states
occupancies. The occupancies eventually are proportional, in the linear response regime, to
a physical measurable quantities such as electrical conductance in the EG model, [7], and
sound or heat conductance in the TLS model, Appendix H.
Following a similar work done on the Electron-glass model, [7]. We apply the meanfield approximation at finite temperature to both of the models and use it to calculate the
density of states. We then expand the rate equation to linear order and obtain the
1
λ
rate
distribution which leads to the observed logarithmic relaxation. Finally, we conclude with
the effect of interactions and system size. We conduct a comparison between the models by
varying the disorder, W , interaction, J and system size, N , which causes the rates to shift
on the
1
λ
curve. Specifically, for both models, changing the disorder to higher values causes
a shift of the rates distribution to higher values and increasing the number of sites causes a
shift of the rates distribution to lower values. On the other hand, changing the interactions
strength directly gives an opposite effect, i.e., in the EG model, increasing the interactions
47
shifts the rates to lower values while in the TG model rates are shifted to higher values.
Generally, the logarithmic shape of the relaxation process and the dependence on the
interaction and system size can be traced back to the basic similarities and differences between models.
The similarities:
1. The disorder and long range interactions, found in both models, give rise to similar
MF equations. As a result the interactions lead to a gap in the DOS. Also, in both
models, the DOS have the same qualitative dependence on disorder, interactions and
system size.
2. The tunnelling amplitude ∆0 in the TG model plays an identical role as the distancesexponential given by the Miller-Abrahams rates in the EG case. Therefore, in both
models the rates are distributed as λ1 . Also, the interactions can only shift the rates
on the
1
λ
curve.
The differences:
1. The dipole interactions in the TG model are naturally weaker than the Coulomb interactions present in the EG model. This leads, among other ingredients, to a narrower
gap in the TG DOS. In addition, the weaker interactions in the TG model generally
causes a smaller shift of the rates distribution.
2. Transitions between states of different TLSs are not allowed. Therefore, the probability
is conserved for each TLS separately which is the cause for the block-diagonal shape of
the TLS rate matrix (or diagonal in the pseudospin representation). This is in contrast
to the non-diagonal rate matrix in the EG model. As a result the interactions appear
only in the quadratic term of rate equation, in the TG model, while in the EG the
interactions enter in the linear term. Also, the direct interactions term in the EG model
48
found to be negligible while in the TG it might give a small correction (see Appendix
G).
3. The total occupation in each TLS is unity, therefore the Pauli-exclusion term is not
added to the rates as in the EG model. Together with difference (2) the resulting rates
of each model have a different dependence on the MF energies. Meaning, in the TLS
model the relaxation rates have a monotonic dependence on the DOS, whereas in the
EG model the dependence is not as direct given the transformation of the EG rate
matrix from N 2 − N transition rates to N relaxation rates, i.e. the diagonalization
procedure.
49
Appendix A
Fermi-Dirac occupations
Generally, given the interactions, it is not trivial to obtain Fermi-Dirac distribution at equilibrium. Following [37], one obtains the Fermi-Dirac distribution of the average electronic
occupation, assuming only that site i has a given energy i which in turn can be obtained
by any approximation method. The average occupation of site i with energy i and chemical
potential µ is defined as:
P
T r [ni exp (−β (H − µ i ni ))]
P
hni i =
T r [exp (−β (H − µ i ni ))]
where H = Heg is the electron glass Hamiltonian 2.1, β =
1
T
(A.1)
and the trace is taken over all
the states of the system n1 = ±1, ..., nN = ±1.
Decomposing the sum over the sites into two parts Heg = i ni + H0 and substituting
in A.1:
P
T r [ni exp (−β (i ni − µ)) exp (−β (H0 − µ i ni ))]
P
hni i =
T r [exp (−β (i ni − µ)) exp (−β (H0 − µ i ni ))]
P
T rni [ni exp (−β (i ni − µ))] T r0 [exp (−β (H0 − µ i ni ))]
P
=
T rni [exp (−β (i ni − µ))] T r0 [exp (−β (H0 − µ i ni ))]
T rni [ni exp (−β (i ni − µ))]
1
=
=
T rni [exp (−β (i ni − µ))]
exp (β (i − µ)) + 1
51
(A.2)
For the specific case of the mean-field model, one obtains the single-particle energies which
also results in the Fermi-Dirac occupation.
Under the same assumptions we can obtain the TLSs average occupations. For TLS i
the upper state occupation pi1 with energy 21 Ei is:
pi1 =
1
exp
1
βEi
2
(A.3)
+1
where we set the chemical potential µ = 0.
Using probability conservation pi2 + pi2 = 1 we can represent the thermal average of
TLS at site i as a pseudospin:
σi0
≡
hσiz iT
=
pi1
−
pi2
=
2pi1
− 1 = − tanh
1
βEi
2
(A.4)
where σ z is the z component Pauli matrix which can take values in the range [−1, 1]. The
superscript X 0 represents equilibrium value.
52
Appendix B
Single two-levels Hamiltonian
Following a similar method as used in [21] and [24]. Given a general one dimensional Hamiltonian with a local potential one can derive the known shape of the effective two-level
Hamiltonian. Two main assumptions are used: (1) the local potential is a smooth function
that changes slowly in the relevant length scale, (2) there is non zero overlap of left and right
states which is small compared to the on-site energies. It is worth mentioning that another
way to derive the TLS effective Hamiltonian is the WKB method which does not requires
some of the assumptions and approximations done here.
Holding in mind the same potential as in Fig. 2.5, the Hamiltonian of the atom or
molecule center of mass is:
H=
p2
+ U (x)
2m
(B.1)
where U (x) is the local double-well potential and x is the location measured from the left
well.
Using the fact that the potential U (x) is a smooth function and varies slowly in x, one
can expand the potential locally near the left and right minimum points, x = 0 and x = d
respectively. Separating the total Hamiltonian in to two parts, an harmonic oscillator, Hho ,
53
and a residual term:
HL =
p2
d2 U (x) + U (0) +
x2 + V1 (x) = Hho + U (0) + V1 (x)
2m
dx2 x=0
(B.2)
d2 U (x) p2
+ U (d) +
(x − d)2 + V2 (x) = Hho + U (d) + V2 (x)
2m
dx2 x=d
P
P
d2 U (x) d2 U (x) 2
where V1 (x) = n=2 dx2 x and V2 (x) = n=2 dx2 (x − d)2
x=0
x=d
d2 U (x) d2 U (x) For simplicity, the same frequency is set for both wells dx2 = dx2 HR =
x=0
(B.3)
x=d
= ~Ω,
and in addition, treating V1 and V2 as perturbations. The resulted approximated energies of
both wells are:
1
1
1
hL|H|Li ≈ hL|U (0) + Hho |Li ≈ ∆ + ~Ω ≈ ∆
2
2
2
(B.4)
Same for the right well1
hR|H|Ri ≈ − ∆
2
2
(x−d)
x2
√
√
1 −
1 −
2
2
where H = HL + HR , hx|Li = ( πx0 )− 2 e 2x0 , hx|Ri = ( πx0 )− 2 e 2x0 are the
1
~ 2
harmonic oscillator ground states for left and right wells respectively, x0 = ( mΩ
) , U (0) = 21 ∆
and U (d) = − 12 ∆. Furthermore, it was assumed that 12 ~Ω |∆| and the next excited states
of harmonic oscillator on both wells are frozen out at low enough temperatures.
In order to evaluate the tunnelling splitting,
1
∆
2 0
= hL|H|Ri, the Hamiltonian is
expanded around the center of the barrier x = d2 .
p2
H=
+U
2m
where
d2 U (x) dx2 x= d2
2
d
d
− ~Ω x −
+ V3 (x) = Hho + V0 + V3 (x)
2
2
= −~Ω and denoting U
d
2
(B.5)
= V0 for the barrier height.
After some algebra the obtained tunnelling amplitude is:
"
#
2
mΩ 2
1
d
+ V0 e− 4~ d
∆0 = 2~Ω − mΩ2
2
2
54
(B.6)
Finally, setting V0 = 21 mΩ2
d 2
2
one gets the known exponential form of the tunnelling
splitting:
1
∆0 ≈ 2~Ωe− 2
where λ =
1
2
q
2mV0
d
~2
q
2mV0
d
~2
≡ 2~Ωe−λ
(B.7)
denotes the strength of the overlap between the right and left wave
functions.
Using B.4 and B.6 the TLS Hamiltonian takes the final form:




1  hL|H|Li hL|H|Ri  1  ∆ ∆0 
≈
2
2
hR|H|Li hR|H|Ri
∆0 −∆
55
(B.8)
Appendix C
Transition rate for the STM
In this appendix the transition rate of a given TLS in the standard tunnelling model is
calculated. This is done by using the TLS-phonon interaction as a perturbation as done in
[38, 23].
Using the rotation transformation, U (see Eq. 2.14) on the Hamiltonian 2.19,
0
H =
HT0 LS
+
0
Htls−ph
=
X
i
Ei Siz
+
XXX
s
i
µν
µν
γis
∆0i x ∆i z
S +
S uµν
is
Ei i
Ei i
(C.1)
where H0 = U HU −1 , i = 1, ..., N is the TLS site index, µ, ν = 1, 2, 3 are Cartesian orientations, s = l, t1, t2 is the three phonon polarizations (one longitudinal and two transverse).
Using the Fourier series expansion of the displacement vector to calculate the strain
tensor at some site i of TLS Siz ,
s
µ
X
∂u
∂
~
†
is
iqxi
uµν
≡
=
a
+
a
=
qs
−qs eqsµ e
is
∂xν
∂xν q
2V ρcs |q|
s
X
~
†
=
iqν aqs + a−qs eqsµ eiqxj
2V
ρc
|q|
s
q
(C.2)
where V is the volume of the sample, ρ is the mass density, q ∈ [−∞, ∞], eqsµ is the µ
component of the direction of mode with wave number q and polarization s.
56
Substituting Eq. C.2 in the second term of Eq. C.1 one obtains:
s
X ∆0i
~
µν
H0 tls−ph =
γis
iqν aqs + a†−qs eqsµ eiqxi Six
Ei 2V ρcs |q|
iqsµν
(C.3)
After averaging over orientations, µ, ν, (or equivalently, approximating an isotropic
µν
media, γis
= γis δ µν , which is a good approximation for an amorphous solid) and using the
identities:
• for the longitudinal mode, s = l,
eqlµ =
qµ
|q|
(C.4)
• for the two transverse modes, s = t1, t2.
eq,t1 · q = eq,t2 · q = eq,t1 · eq,t2 = 0
X
qµ qν
eq,s,µ eq,s,ν = δµν − 2
q
s=t1,t2
the resulted interaction term will take the form of:
s
∆
X ∆0i
X ~
0i x
H0 tls−ph =
i
qγis aqs + a†−qs eiqxi Six =
gik ak + a†−k
Si
E
2V
ρω
E
i
q
i
iqs
jk
(C.5)
(C.6)
Where ωq = cs |q| is the angular frequency of phonon with wave number size q. For future
reference we define k ≡ (q, s) as general state of momentum vector q and polarization s,
q
and gik ≡ i 2V ~ρωq qeiqxi γis .
P
Calculating the matrix elements of Eq. C.6, i,j,n0 ,n00 hψ 00 , n00 |Htls−ph |ψ 0 , n0 i, where |ψ 0 , n0 i
is a state with n0 phonons and a TLS in state ψ 0 . Also, using the Rotating wave approximation (see [39], and for arbitrary coupling [40]) the two resonant contributions are left, for
upward and downward transitions of TLS at site j, respectively:
r
q
∆0i √
0(i)
iqxi
hψ+ ; nqs |Htls−ph |ψ− ; nqs + 1i = ie
γis
nqs
2V ρcs
Ei
r
q
∆0i p
0(i)
iqxi
hψ− ; nqs + 1|Htls−ph |ψ+ ; nq i = ie
γis
nqs + 1
2V ρcs
Ei
57
(C.7)
where nqs = e~ωqs /T + 1
−1
is the occupation of a phonon with energy ~ωqs = ~cs |q|.
Substituting the upward transition matrix element, given in Eq. C.7, in Fermi’s golden
rule, one obtains the upward rate by summing over all phonon modes:
Z
2π X
0(i)
i
ω+ =
dEq |hψ+ ; nqs |Htls−ph |ψ− ; nqs + 1i|2 δ(Eqs − Ei ) =
~ s
X γ 2 ∆2 Ei
X 2π γ 2 ∆20j Z
d3 q
is
is
0i
|q|n
δ(~c
|q|
−
E
)
=
=
nE
qs
s
i
2
3
5
~ 2V ρcs Ei
(2π) /V
cs 2πρ~4 i
s
s
(C.8)
R d3 q
P
where the continuous limit is done, for small phonon level spacing, i.e.
→
=
q
(2π)3 /V
R
V
dqq 2 , also, isotropic orientation is taken since the integrand depends only on the size
2π 2
of the phonon wave vector q. The same calculation for the downward rate gives:
i
ω−
=
X γ 2 ∆2 Ei
is
0i
(nEi + 1)
5 2πρ~4
c
s
s
(C.9)
Finally, the relaxation rate can be defined by the TLS rate equation at equilibrium
(see 4.3) and take the form of the sum of the upward and downward rates:
2
X γ 2 ∆20j Ei
γil 2γit2 ∆20 Ei
Ei
is
−1
i
i
λi ≡ τi = ω+ + ω− =
(2nEi + 1) =
+ 5
coth
≈
c5s 2πρ~4
c5l
ct
2πρ~4
2kB T
s
Ei
γi2 ∆20 Ei
coth
≈ 5
c 2πρ~4
2kB T
(C.10)
where the summation over the polarizations s is done. Also, the averaged coupling constant
and speed of sound, γi = hγis is , c = hcs is respectively, are taken.
Furthermore, following the derivation in [21], one can obtain the phonon transition
rate in a similar manner as conducted above. Given a phonon with energy E = ~ω that
scatters on a TLS at site i, one obtains the phonon rate equations by using the relation
g(E)ṗph = 21 σ̇i , where g(E) is the phonon DOS. After summing over the TLS modes and
using the DOS of TLSs P =
τω−1
P0
∆0
the resulted phonon transition rate is:
X πω
~ω
−1
2
=
γ P tanh
= clph
2 s 0
ρc
2k
T
B
s
s
58
(C.11)
where lph is the phonon mean free path and c is the speed of sound averaged over the three
polarization orientations.
59
Appendix D
TLS-TLS interactions
There are several techniques to derive the TLS-TLS interactions, such as: dipole strain source
in a continuous media [41], second-order perturbation theory with a canonical transformation
to a new equilibrium positions of the normal oscillator coordinates [42, 43] and polaron
transformation [44, 45]. For semantic reasons we derive the TLS-TLS interaction with the
polaron transformation approach.
Starting from the Hamiltonian in the TLS location basis as seen in 2.19, averaged over
orientations µ, ν, with the addition of phonon bath term, we obtain:
Hsb =
X
∆j Sjz
+
∆0j Sjx
+
j
X
gjk ak +
a†−k
Sjz +
X
~ωk a†k ak
(D.1)
jk
jk
where gjk and the index k are defined in Eq. C.6. In more general context this Hamiltonian
is called the spin-boson (SB) Hamiltonian.
Applying the polaron transformation:
"
#
X gik0 Up = exp −
ak0 − a†−k0 Siz
0
~ω
k
ik0
which is simply a displacement operator Up = e−
P
ik
Siz pk ak
, where ak ≡
(D.2)
N γis
,
V ρ|q|c2s
N is the
number of units cells, V is the total volume, ρ is the mass density and cs is the speed
60
of sound for phonons with polarization s. Meaning, the Polaron transformation displaces
distance of ak the phonon wave function at the position xj , for each phonon state (q, s)
differently depending on the occupation state of the TLS Sjz .
Using also Baker-Hausdorff identity:
eB̂ Âe−B̂ = Â + [B̂, Â] +
1
[B̂, [B̂, Â]] + ...
2!
(D.3)
we obtain:
Up Sjz Up† = Sjz
P gjk
P gjk
1
1
−
a−a†
a−a†
Up Sjx Up† = Sj+ e k ~ωk ( −k ) + Sj− e k ~ωk ( −k )
2
2
X g∗
ik
Up ak Up† = ak −
Siz
~ω
k
i
∗
X
gik
Up a†−k Up† = a†−k −
Siz
~ω
k
i
∗
X
X X gik gjk
1
∗ †
Up a†k ak Up† = a†k ak −
gik ak + gik
SzSz
ak Siz +
2 i j
~ωk i
(~ω
)
k
i
j
(D.4)
h
i
where we used the spin and bosonic commutation relations Sil , Sjm = ilmn Sin δij , ak , a†k0 =
h
i
†
†
∗
0
0
δk,k , [ak , ak ] = ak , ak0 = 0 respectively, gi,−k = gi,k
and a symmetric phonon spectrum
ω−k = ωk .
The resulted form of the Hamiltonian is:
H0 = Up Hsb Up† =
∗
X
X
X ∆0j
X X gik gjk
Sj+ Bj− + Sj− Bj+ −
Siz Sjz +
~ωk a†k ak
=
∆j Sjz +
2
~ωk
ij
j
j
k
k
where
Sj±
=
Sjx
±
Sjy
and we define
Bj±
±
=e
P
gjk
k ~ωk
(D.5)
(a−a†−k ) .
For weak TLS-Phonon coupling we can expand B ± keeping zero and first orders, Bj± =
P g
1 ± k ~ωjkk (ak − a†−k ), and retain the original Hamiltonian but with the addition of TLS-TLS
61
interactions:
H0 =
X
j
∗
X X gik gjk
X
X X i∆0j gjk
∆j Sjz + ∆0j Sjx −
(ak − a†−k )Sjy −
Siz Sjz +
~ωk a†k ak
~ω
~ω
k
k
j
ij
k
k
k
(D.6)
By using transformation in Eq. 2.14 we rotate Eq. D.6 to the diagonal basis of the
non-interacting TLSs:
H̃ = U H0 U −1 =
X
X X i∆0j gjk ak − a†−k S̃jy −
=
Ej S̃jz −
~ωk
j
j
k
X
∗ X X gik gjk
∆i z ∆0i x
∆j z ∆0j x
−
S̃i −
S̃i
S̃j −
S̃j +
~ωk a†k ak
~ωk
Ei
Ei
Ej
Ej
ij
k
k
where we used the following identities,
U S x U −1 =
∆0 z
S̃
E
+
∆ x
S̃
E
∆
E
= cos(2θ),
∆0
E
= sin(2θ), U S z U −1 =
∆ z
S̃
E
(D.7)
−
∆0 x
S̃ ,
E
and U S y U −1 = S̃ y .
For temperatures T & 0.03K and g E2 E2
,
∆ ∆0
we can approximate the direct TLS-PH
interaction as the dominating term for the TLS transitions (see Sec. 2.2.4 and [30]), thus,
neglecting TLS-TLS interactions of the form Six Sjx , Siz Sjx and Six Sjz and keeping the diagonal
interaction term Siz Sjz as a correction to the non-interacting TLSs energies. Consequently
the approximated Hamiltonian is:
H̃ = U H0 U −1 =
∗
X
X X i∆0j gjk X X gik gjk
∆i ∆j z z X
=
Ej S̃jz −
ak − a†−k S̃jy −
S̃i S̃j +
~ωk a†k ak
~ω
~ω
E
E
k
k
i j
j
j
ij
k
k
k
(D.8)
Since we are treating the second term in Eq. D.8 as perturbation the summation over
the phonon degrees of freedom is performed only on the unperturbed Hamiltonian, giving
62
rise to the final shape of the model Hamiltonian:
H̃ =
X
Ej S̃jz −
j
X uij ∆i ∆j
ij
3
rij
Ei Ej
S̃iz S̃jz +
X X i∆0j gjk j
k
~ωk
a†−k − ak S̃jy
(D.9)
where rij = |ri − rj | is the distance between two TLSs i and j. When performing the
summation one uses the phonon spectral function J(ωk ) ∝ ωk3 (A super ohmic bath) with
the Debye’s approximation ωk ≡ ωqs = cs |q| as in [46] and also we dropped a constant.
To get the full shape of the TLS-TLS interactions as seen in [41, 42] one must sum over
µν ∗µν
P P gik
gjk
the exact coupling constants which are not averaged over orientations, i.e.
k
µν
~ωk
µν
where gik
∝ γiµν and γiµν is the deformation potential at site i, as defined in Eq. 2.19, (see
[44] for more details).
By using the third term in Eq. D.9 as perturbation and repeating the calculation of the
transition rates performed in Appendix C an important realization follows, we obtain the
same transition rate as in the standard tunnelling model given in Eq. C.10, the reason for
this is that under the Rotating wave approximation the matrix elements of the perturbation
π
have an extra phase e−i 2 which Fermi’s golden rule is not sensitive to.
Note that the diagonal basis of a single TLS is generally the natural choice when
dealing with dissipation. However, for greater accuracy and numerical compatibility, we
start in section 2.2.4 with the Hamiltonian presented in the location basis, given in Eq. D.6,
and transform the Hamiltonian after making the mean-field approximation.
63
Appendix E
Analytical derivation of the dipole gap
Following Burin’s derivation ,[36]. A gap appear in the TLS DOS for small energy values.
The gap takes a logarithm shape due to the characteristics of the dipole interactions and the
dimensionality of the system.
Given the MF equations, 2.29, and the temperature range, (see the discussion in
Sec. 2.2.4), the elementary excitations in the system can be approximated as two-TLS excitation (see [47] for the TLS-glass and[48] for the Electron-glass):
Eij ≈ Ei + Ej − 2Jij
where Jij =
uij
3 .
rij
(E.1)
For simplicity one approximates Jij = ± rU30 with probability
ij
1
2
to be with
+ or − sign. Eq. E.1 can be explained as follows:
In the pair approximation, the system is defined to be in the ground state when all
two-particle excitations are positive (for weak interactions other multi-particle excitations
can be neglected). Accordingly, the density of DOS of the non interacting TLSs, P (2.16)
takes an additional factor:
Pi → P
Y
Θ(Ei + Ej − 2Jij )
j
where Θ is a step function.
64
(E.2)
In terms of the TLS energies, P takes the following form:
P0
E
p
Θ(E − ∆0 )Θ(∆ − ∆0min )Θ(W − E)
(E.3)
2
∆0 E − ∆20
p
where the addition of step functions imposes a cut-off for E = ∆2 + ∆20 and ∆0 , ∆0min ∼
P ≈
10−7 K is the lower cut-off value for ∆0 and W ≥ 1K is upper cut-off value for both ∆ and
∆0 .
For finite temperature, only pairs with interactions that exceed the temperature, 2Jij >
T , are effecting the density of states,
Pi → P
Y
Θ(Ei + Ej − 2Jij )
Y
Θ(2Jij 0 − T )
(E.4)
j0
j
Expanding up to the first correction in the small parameter χ = P0 U0 1 (see 2.24):
Pi ≈ P
Y
Θ(Ei + Ej − 2Jij )Θ(2Jij 0 − T ) = P
Y
[1 − Θ(−Ei − Ej + 2Jij )]Θ(2Jij − T ) (E.5)
j
j
After averaging over the sites j and expanding again in the parameter χ, a decrease in
the density of states is obtained:
*
+
X
P 0 (E, ∆0 ) ≈ 1 −
Θ(−Ei − Ej + 2Jij )Θ(2Jij − T )
j
j
P0
E
p
(E.6)
Θ(E − ∆0 )Θ(∆0 − ∆0min )
2
∆0 E − ∆20
"
#
Z
Z W
Z
P0
d∆00 W
2U
2U
0
0
× 1−
dE 0 Θ −E − E 0 + 3 Θ
dr
−T
0
3
0
2
∆
r
r
∆0min
∆0
0
=
where the factor
1
2
accounts for half of the neighbours with the positive coupling (for Jij < 0
the step function in Eq. E.6 gives zero). Also the factor √
E0
E 2 −∆20
in the integral over E 0 is
approximated as unity since for the relevant energies E 0 ≥ T the limit ∆00 E 0 holds.
After performing the integrals, the final shape of the gap is obtained:
P0
E
4π
W
W
p
P (E, ∆0 ) ≈
Θ(E − ∆0 ) 1 −
P0 U0 log
log
∆0 E 2 − ∆20
3
E+T
∆0min
0
65
(E.7)
where the additional term log(W/T )log(W/∆0min ) ∼ 10−1 was neglected.
The weak interaction assumption taken above, is valid if the correction to the density
of states (E.3) is small,
P − P0
≈ P0 U0 log
P
W
T
log
W
∆0min
1
(E.8)
which applies in most amorphous solids.
Note that doing the same calculation for the DOS where the asymmetry energy is
Gaussian distributed (2.15) does not change the logarithmic shape of the gap, it only yields
a different constant prefactor for the logarithm. Furthermore, the dipole gap and coulomb
gap can be derived without the assuming weak interactions, see [10].
66
Appendix F
Expansion of the rate equation
Starting from the rate equation 4.3. For small deviation from the equilibrium values of the
TLSs (δσ 1) we can take only the first two orders of the expansion:
dσi X ∂F 1 X ∂ 2 F '
δσj +
δσj δσk
dt
∂σj σ =σσ0
2 j,k ∂σj ∂σk σ =σσ0
j
(F.1)
σ )σi − Ai (σ
σ ), the zero’th-order is by definition zero, i.e. F(σ
σ 0 ) ≡ 0, δσ
σ≡
where F ≡ −λi (σ
Ei
σ − σ 0 is the TLSs vector of deviations from equilibrium values, and σi0 = − tanh 2T
(see
A.4). Also, we expand up to second order to include also a direct term of the interactions.
Using the equality:
∂ni
∂ni ∂Ei ∂∆0i
∆0i uij
=
·
=
−βn
(n
+
1)
·
i
i
3
∂σj
∂Ei ∂∆0i ∂σj
Ei rij
(F.2)
where ni is the phonon occupation, Ei and ∆0i are the MF energy and asymmetry energy
defined in Eq. 2.27 and Eq. 2.26 respectively, and uij is the interaction constant defined in
appendix D.
67
The first and second derivatives are:
0 P
P ∂F 0 ∆i
δσ
=
−A
(2n
+
1)
δσ
+
2A
βn
(n
+
1)
σ
j
i
i
i
i
i
i
i
j ∂σj j6=i
E
i
0
σ =σ
σ
And
P
1
∆0i
∂2F δσ
δσ
=
2A
βn
(n
+
1)
j
k
i
i
i
jk ∂σj ∂σk
2
Ei
0
σ
σ =σ
2 P
β∆0i
uij uik
+Ai βni (ni + 1) Ei
3 r 3 δσj δσk
j,k6=i rij
ik
uij
3 δσj
rij
uij
3 δσj δσi +
j6=i rij
P
(F.3)
where Ai ≡ ai ∆20 Ej and ai is a constant defined in Eq. 2.20.
After keeping only the terms that are proportional to δσi we end with the N-species
(or N-TLSs) Lotka-Volterra equations:
dσi
∆0 X uij
' −Ai (2ni + 1) δσi + 2Ai βni (ni + 1) i
δσj δσi
3
dt
Ei j6=i rij
(F.4)
As can be seen in Eq. F.4, the prefactor of the linear term is the transition rate
Ei
derived at Appendix C. The linear order coincide with the
λi = Ai (2ni + 1) = Ai coth 2T
relaxation time approximation given in literature [23], although in this case it is calculated
directly from the rate equation instead of being assumed, and the energies are the mean-field
energies.
Finally, the general solution to Eq. F.4 takes an exponential form:
P
δσi (t) = bi e−λi t e
j6=i fij
Rt
0
δσj (t0 )dt0
(F.5)
where bi ≡ δσi (0) is the initial value of TLS i right after the external strain driving force has
0
Ei ∆i uij
.
stopped and fij ≡ A2i sinh−2 2T
Ei T r3
ij
68
Appendix G
Perturbation treatment of the
quadratic term
The purpose of this appendix is to evaluate the contribution of the direct interaction term
(the quadratic prefactor in Eq. F.4) to the exponential decays that are given from the linear
term. This is done as an analogy to the electron-glass model in which the interactions term
is introduced already in the linear expansion (second term of Aij in Eq. 2.7).
Since the expansion of the rate equation holds for small deviations from equilibrium,
the quadratic term in Eq. F.4 must be small compared to the linear term. This, with
the addition of studying the low rates λi t 1 allow us to treat the quadratic term as
perturbation. Accordingly, the given inequality is:
X
t
Z
δσj (t0 )dt0 λi t 1
fij
(G.1)
0
j6=i
By approximating the left hand side we can express inequality G.1 in terms of the small
parameter bj ≡ δσj (0):
1X
fij
t j6=i
Z
t
0
0
δσj (t )dt ≤
0
X
bj fij '
sX X
j6=i
j6=i
69
j 0 6=i
bj bj 0 fij fij 0 ≈
sX
j6=i
b2j fij2
(G.2)
where the equality G.2 is done under the assumption of stability of the linear solution,
which means that the deviation from equilibrium value is maximal at time t = 0 when the
perturbation is applied (δσj (t0 ) ≤ bj ), also in the last approximation we used the large N
limit and the fact the fij is Gaussian distributed with zero mean value.
Expanding the exponent in Eq. F.5 results in:
!
Z tZ t
Z t
X
X
δσi (t) ' bi e−λi t 1 +
σj (t0 )dt0 +
bj bj 0 fij fij 0
σj (t0 )σj 0 (t00 )dt0 dt00 + ...
bj fij
j6=i
0
0
j,j 0 6=i
0
(G.3)
where time ordering scheme is not applied since [σj (t), σj 0 (t0 )] = 0, ∀j, j 0 .
0
As a first approximation we evaluate σj (t0 ) ≈ bj e−λj t , after substitution in Eq.G.3 we
obtain:
"
−λi t
δσi (t) ' bi e
1+
X
j6=i
#
X
fij
fij fij 0
−λj t
−λj t
−λj 0 t
bj
+ ...
1−e
+
bj bj 0
1−e
1−e
λj
λj λj 0
j,j 0 6=i
(G.4)
The corrections of the quadratic term, to the exponential decay of the TLSs given in
the linear term, does not effect the stability of the solution, since the they are smaller than
unity, which ultimately results in a decay of the spin deviation, δσi (t), as expected.
To evaluate roughly the direct interaction terms in Eq. G.4 we calculate leading order of
fij for several limits which are relevant to the low rates (i.e. we disregard the limit ∆0i ∆0 i
or ∆0i & T ):

λi
−

Ei − ETi
λi
ai ∆2 T

0i
e
≈
e

2
2T
2ai ∆0i T


q



 1 1 − ∆0i 2
u
2
T
ij
fij ' ai ∆20i 3
rij 


1





 ∆0i
2T
, Ei T, ∆0i ∆0 i
, Ei ∼ T
(G.5)
0
, Ei T, ∆0i ∆ i
, Ei T, ∆0i ∼ ∆0 i
As can be seen from Eq. G.5, the direct interaction terms may give, in some limits
and parameter values, a small but relevant correction to the exponential decay given in the
70
linear term in Eq. G.4. This is in contrary to the Electron-Glass model where the direct
interaction term is found to give a negligible corrections, [7].
71
Appendix H
External perturbation and initial
deviations
An applied external stress perturbation shifts the asymmetry energy ∆0 i of TLS i by the
amount φi , the resulted new equilibrium energy is:
q
Eiφ = (∆0 i + φi )2 + ∆20i
(H.1)
In the linear response regime the energy shift can be represented by the coupling
between the external stress-field and the local strain-dipole moment caused by the TLS:
φi = −2
X
µν
v0 Gµν
i σi
(H.2)
µν
where σiµν is the external stress, Gµν
i is the local strain dipole, v0 is a typical unit cell volume
and µ, ν indicate the Cartesian orientations.
The energy shift can be expressed also in terms of the local strain and deformation
potential (defined in Eqs. 2.17 and 2.19 respectively):
φi = −2
X
µν
72
v0 γiµν uµν
i
(H.3)
where we used the relations Gµν
i =
P
αβ
S αβµν γiαβ and σiµν =
P
αβ (S
−1 αβµν µν
)
ui ,
Siαβµν is the
elastic susceptibility (also known as the compliance tensor) which also relates the applied
stress to local strain in the linear regime (see [49], [50] for more information). Furthermore,
the compliance tensor is generally complex, its real part is related to a deviation in the sound
velocity ∆c/c and the complex part is the dissipation of the system which takes the form
of internal friction coefficient Q−1 . Both quantities are measured in experiments (see for
example [51]).
In most cases, the largest strain that could be applied without heating the sample
is umax ≈ 1.2 × 10−5 (see [23]), taking into account also the typical TLS-phonon coupling
γ ∼ 1eV we get an upper bound for the energy scale of the perturbation, φ . 0.14K.
After an external perturbation φi have been applied we show that the initial deviations
depend only weakly on the rates and thus can be approximated as uniformly distributed
over the rate space. To show that, we approximate the deviations as the difference between
the equilibrium occupations of TLS i before and after the applied field (i.e. the adiabatic
approximation):
0
E
E
i
iφ
0
ci ≡ δσi (t = 0) ≈ σi − σiφ ≈ tanh
− tanh
2T
2T (H.4)
Given the upper bound of the external stress perturbation, φmax ≈ 0.14K and temperature range of T ≈ 0.03K − 0.1K, we evaluate leading order of Eq. H.4 for two relevant
limits, φi T and φi ∼ T :


φi




T
ci ∼ 0.1




E

e− Ti
, φi T, ∀Ei
, φi ∼ T Ei
(H.5)
, φi ∼ T Ei
where two limits for the mean-field asymmetry energy are taken: ∆0 i ∆0i and ∆0 i ∼ ∆0i .
The third scale, ∆0 i ∆0i is neglected, since we are interested in TLSs that have low
73
relaxation rates. In other words, a TLS with the maximal rate for a given energy Ei and
temperature T has ∆0i ' Ei , which is another way to express the limit ∆0 i ∆0i .
As can be seen from Eq. H.5, for the first two limits ci is independent of rates (see
the rates dependence on MF energies, Eq. 2.20), up to leading order. For the third limit we
have:
ci ∼

∆0 i


e− T


e−
2λi
ai T 3
, ∆0 i ∆0i
1
3
'e
1/3
λ
i
− 10
(H.6)
, ∆0 i ∼ ∆0i
where in the first case, ci depends weakly on the rates as argued before (see the derivation of
Eq. 4.9). For the second case, in order to be consistent with the condition T Ei (as seen
in Eq. H.5) the rates must be λi > 103 1/s which are fast relaxing TLSs and therefore this
limit is not relevant. Meaning, for a given external adiabatic perturbation, φi ∼ T , there are
no slow relaxing TLSs with ∆0 i ∼ ∆0i .
74
Appendix I
Additional graphs of TG and EG
models
For all the schemes presented in this appendix we set the number of sites, N = 1000,
−d
the temperature, T = 0.1, density of sites, rnn
= 1, where d = 2, 3 is the dimensions of the
Electron-glass and the TLS-glass systems respectively. Also, all the parameters are presented
3
for the TLS-glass and
in units of average nearest-neighbours interactions J which is U0 /rnn
e2 /rnn for the Electron-glass.
Probability density
0.03
0.02
Figure I.1: TG model: Density of states
for different disorder values.
0.01
The disorder is varied W = 1 (blue stars), W = 10 (green
squares) and W = 100 (red triangles) for interactions J = 1
−20
−15
−10
−5
0
E
5
10
15
20
and constant temperature T = 0.1.
75
0.04
Probability density
Probability density
15000
10000
5000
0.03
0.02
0.01
0
0
1
2
3
4
−15
5
(a)
−10
−5
0
5
10
15
ln(−λ)
−5
x 10
−λ
(b)
Figure I.2: TG model: Decay rates distribution for different disorder values.
The distribution of eigenvalues Eq. 2.20, the rest of the parameters are the same as Fig. I.3.
(a) Rates in the regular scale with a fit to ∝
1
x
function. The W = 10 curve (green triangles) is omitted since the differences
from W = 1 are small, (b) Rates in the natural log scale.
Probability density
0.18
0.14
0.1
0.06
Figure I.3: TG model: Density of states
for different interaction values.
0.02
The interaction is varied J = 1 (blue stars), J = 1/3 (green
−20
−15
−10
−5
0
E
5
10
15
20
squares) and W/J = 1/10 (red triangles) for disorder W = 1.
76
4
x 10
Probability density
Probability density
7
5
3
0.03
0.02
0.01
1
0.2
0.4
0.6
0.8
−15
1
(a)
−10
−5
0
5
10
15
ln(−λ)
−5
x 10
−λ
(b)
Figure I.4: TG model: Decay rates distribution for different interaction values.
The distribution of eigenvalues Eq. 2.20, the rest of the parameters are the same as Fig. I.3.
(a) Rates in the regular scale with a fit to ∝
1
x
function, (b) Rates in the natural log scale.
Probability density
0.3
0.2
Figure I.5: EG model: Density of states
for different disorder values.
0.1
The disorder is varied, W = 1 (blue stars), W = 5 (green
squares) and W = 10 (red triangles) for constant interaction
−10
−5
0
E
5
10
strength J = 1 and constant temperature T = 0.1.
77
0.1
Probability density
Probability density
1600
1200
800
400
0.08
0.06
0.04
0.02
0
1
2
3
4
5
6
−35
7
(a)
−30
−25
−20
−15
−10
−5
0
ln(−λ)
−4
x 10
−λ
(b)
Figure I.6: EG model: Decay rates distribution for different disorder values.
The distribution of eigenvalues of the rate matrix Eq. 2.7 the rest of the parameters are the same as Fig. I.5.
(a) Rates in the regular scale with a fit to ∝
1
x
function, (b) Rates in the natural log scale.
Probability density
1
0.8
0.6
0.4
Figure I.7: EG model: Density of states
for different interaction values.
0.2
The interaction is varied, J = 1 (blue stars), J = 1/3 (green
squares) and J = 1/9 (red triangles) for constant disorder
−4
−3
−2
−1
0
E
1
2
3
4
W = 1.
78
1600
Probability density
Probability density
0.12
1200
800
400
0.1
0.08
0.06
0.04
0.02
0
1
2
3
0
−25
4
−20
−15
x 10
(a)
−10
−5
0
ln(−λ)
−3
−λ
(b)
Figure I.8: EG model: Decay rates distribution for different interaction values.
The distribution of eigenvalues of the rate matrix Eq. 2.7 the rest of the parameters are the same as Fig. I.7.
(a) Rates in the regular scale with a fit to ∝
1
x
function, (b) Rates in the natural log scale.
0.16
0.12
Probability desity
Probability densiy
30
20
10
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.04
−20
0.08
−15
−10
−5
0
5
ln(-λ)
−λ
(b)
(a)
Figure I.9: EG model: Decay rates distribution for different site number for no energy
dependence.
The eigenvalues distribution of the rate matrix, with no energy dependence, Aij = erij /ξ . N = 1000 (blue stars), N = 100
(green squares) and N = 10 (red triangles) The average nearest-neighbours distance rnn = 1 (constant density of sites) and
localization length ξ = 0.1rnn .
(a) Rates in the regular scale with a fit to ∝
1
x
function, (b) Rates in the natural log scale.
79
Appendix J
Numerical routines
Three Matlab routines are displayed:
1. DistanceMatrix.m file: Builds a three dimensional sample, or two dimensional for
the EG model, with uniform distributed sites and returns the sites distance matrix
Rij = |ri − rj |.
2. EnergyIteration.m file: Finds the solution to the self-consistent equation of the asymmetry energies Eq. 2.27 (or site energies, Eq. 2.4, for the EG).
3. HcTLS.m file: Heats and cools the system. In each temperature step this file uses
EnergyIteration.m file to find the solution of the asymmetry energies in the current
q
temperature value. Finally, the energies are calculated Ei = ∆20i + ∆0i 2 .
J.1
DistanceMatrix.m file
function [R_ij, r_nn] = DistanceMatrix3D(N)
80
L=1*((N^(1/3))-1); %N is the number of sites and L^3 is the volume of the
system. Making the density of sites constant by setting average nearest
neighbours distance r_nn=1. This expression is accurate enough for uniform
distributed site in the large N limit.
c = L*(rand(N,3)-0.5); %Creating a uniform distributed array for the locations
of the sites.
R_ij = zeros(N); %Calibrating distances matrix
for i = 1:N
for j = i+1:N
%Difference of distance between sites i and j.
dx = c(i,1)-c(j,1);
dy = c(i,2)-c(j,2);
dz = c(i,3)-c(j,3);
%%%%%%%%%Periodic boundary conditions%%%%%%%%%
if abs(dx) > L/2
dx = dx - L*sign(dx);
end
if abs(dy) > L/2
dy = dy - L*sign(dy);
end
if abs(dz) > L/2
dz = dz - L*sign(dz);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81
R_ij(i,j) = sqrt(dx^2+dy^2+dz^2); %Calculating the shortest euclidean
distance between sites i and j
R_ij(j,i) = R_ij(i,j) %Symmetrizing the distances matrix
end
end
J.2
EnergyIteration.m file
function [Delta] = EnergyIterations(N, delta, Delta, delta_0, R_ij, T, r_nn,
u_ij)
rng(’shuffle’) %For cluster usage, changing the seed at each run of the File to
produce new random numbers on different CPUs. Needed for randperm function.
maxdDelta = 1;
PrevDelta = zeros(1,N);
diffLimit = 10^-4; % Convergence limit for the iterations
count = 1;
%%%%%%%%%%%%Iterative evolution of the initial asymmetry energies%%%%%%%%%%%%
while (maxdDelta > diffLimit)
v = randperm(N); %Random vector with no repetition of values. Needed for
random permutation in the Iterative evolution process.
82
for i=(1:N)
%Calculating the new TLS energy at site v(i) using the Mean-Field self
consistent equations
Delta(v(i)) = delta(v(i)) +
2*sum((r_nn^3)*((u_ij(v(i),1:v(i)-1)).*((1./((1+exp((sign(Delta(1:v(i)-1)).*sqrt(Delta(1:
+ delta_0(1:v(i)-1).^2))/T)))-0.5)./((R_ij(v(i),1:v(i)-1)).^3)))) +
2*sum((r_nn.^3)*((u_ij(v(i),v(i)+1:N)).*((1./((1+exp((sign(Delta(v(i)+1:N)).*sqrt(Delta(v
+ delta_0(v(i)+1:N).^2))/T)))-0.5)./((R_ij(v(i),v(i)+1:N)).^3))));
end
if(count) %Saving the maximum change between following iterations (for the
exit condition).
DDelta = Delta - PrevDelta;
maxdDelta = max(DDelta);
end
PrevDelta = Delta;
count = count + 1;
end
J.3
HcTLS.m file
function [] = HcTLS(RealNum, N)
rng(’shuffle’) %Needed for the use of random numbers of DistanceMatrix3D file.
83
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Calibrating variables%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_low = 0.1; %Initial and final temperature in the heating-cooling process.
T_high = 1; %Highest temperature in the heating-cooling process
%The temperature is defined though nearest neighbours average interaction ,J_nn,
for example: J_nn/T=1,10,100.
delta = normrnd(0,1,1,N); %Normal distributed TLS asymmetry energy.
Delta = 20*(rand(1,N)-0.5); %Initial mean-field TLS assymetry energy vertor,
Delta = delta for no interactions.
cd /fastspace/users/asban/parallelTLS/WithDiffDistanceMatrices %Changing
directory.
[R_ij, r_nn] = DistanceMatrix3D(N); %Creating new distance matrix in each
realization.
u_ij = normrnd(0,1,N,N); %u_ij is normal distributed prefactor of the
interaction that comes from the random relative orientation between TLSs i
and j.
u_ij = u_ij - tril(u_ij,-1) + triu(u_ij,1)’; %Symmetrizing u_ij.
%%%%%%%%%%%%%%%%%%%%%%%%%%%Random tunnelling Splitting%%%%%%%%%%%%%%%%%%%%%%%%%%
84
%Creating a random vector size N that represents delta_0 (tunnelling splitting)
with p/delta_0 distribution by using an exponent with a power of uniform
distributed vector v. delta_0 is in range [10^-7, 10^-1] in units of J_nn.
Keeping the universal relation p*U_0=10^-3, where (U_0)^2 = <(u_ij)^2>.
p=10^(-3);
n=10^(p*6); %The base of the log which together with parameter ’a’ gives the
range of the distribution of 1/delta_0.
a=10;
y1=(log(1/(a*(n^(1/p)))))/(log(n)); % the range of the random vector of uniform
distribution that with it I create the vector with 1/delta_0
y2=(log(1/a))/(log(n));
v = (abs(y2-y1))*rand(1,N)+min(y1,y2); %Uniform distributed vector in the range
[1/a*n^{1/p}, 1/a]
delta_0 = n.^(v); %Tunnelling splitting vector with p/delta_0 distribution.
%%%%%%%%%%%%%%%%%%%%%%%%Heating and then Cooling procedure%%%%%%%%%%%%%%%%%%%%%%%
T = T_low; %Initializing the temperature
Num_Of_T_Steps = 10; %Number of heating and cooling steps.
size_Of_T_Step = 2*(T_high-T_low)/(Num_Of_T_Steps); %Size of each
Heating/Cooling step .
for k = 1:(Num_Of_T_Steps)
[Delta] = EnergyIterationsSlave(N, delta, Delta, delta_0, R_ij, T, r_nn,
u_ij); %Finding the Mean-Field energies for the current temperature, T.
85
if(k <= Num_Of_T_Steps/2) %Heating the system
T = T + size_Of_T_Step;
end
if(k > Num_Of_T_Steps/2) %Cooling back the system
T = T - size_Of_T_Step;
end
end
E = sign(Delta).*sqrt(Delta.^2 + delta_0.^2); %Energy difference vector of all
the TLSs.
cd /fastspace/users/asban/output_mat %Changing directory
fname = sprintf(’Eanddelta0_%d.mat’, RealNum); %Adding a realization number to
the file name.
save(fname, ’E’, ’delta_0’); %Saving E and delta_0 vectors.
86
Bibliography
[1] D. J. Salvino, S. Rogge, B. Tigner, and D. D. Osheroff, “Low temperature ac dielectric
response of glasses to high dc electric fields,” Phys. Rev. Lett., vol. 73, pp. 268–271, Jul
1994.
[2] Z. Ovadyahu and M. Pollak, “Disorder and magnetic field dependence of slow electronic
relaxation,” Phys. Rev. Lett., vol. 79, pp. 459–462, Jul 1997.
[3] K. Fritsch, J. Friedrich, and B. M. Kharlamov, “Nonequilibrium phenomena in spectral
diffusion physics of organic glasses,” The Journal of Chemical Physics, vol. 105, no. 5,
pp. 1798–1806, 1996.
[4] A. Vaknin, Z. Ovadyahu, and M. Pollak, “Aging effects in an anderson insulator,” Phys.
Rev. Lett., vol. 84, pp. 3402–3405, Apr 2000.
[5] S. Ludwig, P. Nalbach, D. Rosenberg, and D. Osheroff, “Dynamics of the destruction
and rebuilding of a dipole gap in glasses,” Phys. Rev. Lett., vol. 90, p. 105501, Mar
2003.
[6] Z. Ovadyahu, “Optical excitation of electron glasses,” Phys. Rev. B, vol. 83, p. 235126,
Jun 2011.
[7] A. Amir, Y. Oreg, and Y. Imry, “Mean-field model for electron-glass dynamics,” Phys.
Rev. B, vol. 77, p. 165207, Apr 2008.
88
[8] M. Pollak, “Effect of carrier-carrier interactions on some transport properties in disordered semiconductors,” Discuss. Faraday Soc., vol. 50, pp. 13–19, 1970.
[9] A. L. Efros and B. I. Shklovskii, “Coulomb gap and low temperature conductivity of
disordered systems,” Journal of Physics C: Solid State Physics, vol. 8, no. 4, p. L49,
1975.
[10] S. Baranovskii, B. Shklovskii, and A. Efros, “Elementary excitations in disordered systems with localized electrons,” Soviet Journal of Experimental and Theoretical Physics,
vol. 51, p. 199, 1980.
[11] J. H. Davies, P. A. Lee, and T. M. Rice, “Properties of the electron glass,” Phys. Rev.
B, vol. 29, pp. 4260–4271, Apr 1984.
[12] E. Levin, V. Nguen, B. Shklovsii, and A. Efros, “Coulomb gap and hopping electric
conductioncomputer simulation,” Soviet Physics JETP, vol. 65, no. 4, pp. 842–848,
1987.
[13] F. G. Pikus and A. L. Efros, “Critical behavior of density of states in disordered system
with localized electrons,” Phys. Rev. Lett., vol. 73, pp. 3014–3017, Nov 1994.
[14] A. Miller and E. Abrahams, “Impurity conduction at low concentrations,” Phys. Rev.,
vol. 120, pp. 745–755, Nov 1960.
[15] B. Shklovskii and A. Efros, Electronic Properties of Doped Semiconductors. Springer,
Berlin, 1984.
[16] L. Onsager, “Reciprocal relations in irreversible processes. i.,” Phys. Rev., vol. 37,
pp. 405–426, Feb 1931.
[17] D. Cox and H. Miller, The Theory of Stochastic Processes. Science paperbacks, Taylor
& Francis, 1977.
89
[18] C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and
Engineers: Asymptotic methods and perturbation theory. Springer, 1999.
[19] R. C. Zeller and R. O. Pohl, “Thermal conductivity and specific heat of noncrystalline
solids,” Phys. Rev. B, vol. 4, pp. 2029–2041, Sep 1971.
[20] W. Phillips, “Tunneling states in amorphous solids,” Journal of Low Temperature
Physics, vol. 7, no. 3-4, pp. 351–360, 1972.
[21] W. A. Phillips, “Two-level states in glasses,” Reports on Progress in Physics, vol. 50,
no. 12, p. 1657, 1987.
[22] P. w. Anderson, B. I. Halperin, and c. M. Varma, “Anomalous low-temperature thermal
properties of glasses and spin glasses,” Philosophical Magazine, vol. 25, no. 1, pp. 1–9,
1972.
[23] P. Esquinazi, Tunneling Systems in Amorphous and Crystalline Solids. Berlin: Springer,
1998.
[24] K. Binder and W. Kob, Glassy Materials and Disordered Solids: An Introduction to
Their Statistical Mechanics. Singapore: World Scientific Publishing Company, 2005.
[25] A. Burin Soviet J. Low Temp Phys, vol. 17, p. 456, 1991.
[26] S. Hunklinger and A. Raychaudhuri, Thermal and elastic anomalies in glasses at low
temperatures, vol. 9. Amsterdam: Elsevier, 1986.
[27] P. Esquinazi, R. Knig, and F. Pobell, “Acoustic properties of amorphous sio2 and
pdsicu, and of crystalline ag, nbti and ta at very low temperatures,” Zeitschrift fr
Physik B Condensed Matter, vol. 87, no. 3, pp. 305–321, 1992.
[28] J. J. Freeman and A. C. Anderson, “Thermal conductivity of amorphous solids,” Phys.
Rev. B, vol. 34, pp. 5684–5690, Oct 1986.
90
[29] C. Yu and A. Leggett, “Low temperature properties of amorphous materials: Through
a glass darkly,” Comments on Condensed Matter Physics, vol. 14, no. 4, pp. 231–251,
1988.
[30] A. Burin and Y. Kagan, “Low-energy collective excitations in glasses. new relaxation mechanism for ultralow temperatures,” Journal of Experimental and Theoretical
Physics, 1994.
[31] S. Alexander, C. Laermans, R. Orbach, and H. M. Rosenberg, “Fracton interpretation
of vibrational properties of cross-linked polymers, glasses, and irradiated quartz,” Phys.
Rev. B, vol. 28, pp. 4615–4619, Oct 1983.
[32] M. Schechter and P. C. E. Stamp, “Inversion symmetric two-level systems and the lowtemperature universality in disordered solids,” Phys. Rev. B, vol. 88, p. 174202, Nov
2013.
[33] E. Cuevas, R. Chicón, and M. Ortuño, “Density of states for a disordered system of
interacting dipoles,” Physica B: Condensed Matter, vol. 160, no. 3, pp. 293–296, 1989.
[34] A. Churkin, D. Barash, and M. Schechter, “Nonhomogeneity of the density of states
of tunneling two-level systems at low energies,” Phys. Rev. B, vol. 89, p. 104202, Mar
2014.
[35] M. Grunewald, B. Pohlmann, L. Schweitzer, and D. Wurtz, “Mean field approach to
the electron glass,” Journal of Physics C: Solid State Physics, vol. 15, no. 32, p. L1153,
1982.
[36] A. Burin, “Dipole gap effects in low energy excitation spectrum of amorphous solids.
theory for dielectric relaxation,” Journal of Low Temperature Physics, vol. 100, no. 3-4,
pp. 309–337, 1995.
91
[37] A. Efros, “Thermodynamics and transport properties of interacting systems with localized electrons,” in Phase Transitions and Self-Organization in Electronic and Molecular
Networks (M. Thorpe and J. Phillips, eds.), Fundamental Materials Research, pp. 247–
262, Springer US, 2001.
[38] J. Jäckle, “On the ultrasonic attenuation in glasses at low temperatures,” Zeitschrift
für Physik, vol. 257, no. 3, pp. 212–223, 1972.
[39] C. Gerry and P. Knight, Introductory Quantum Optics. Cambridge University Press,
2004. Cambridge Books Online.
[40] E. K. Irish, “Generalized rotating-wave approximation for arbitrarily large coupling,”
Phys. Rev. Lett., vol. 99, p. 173601, Oct 2007.
[41] J. L. Black and B. I. Halperin, “Spectral diffusion, phonon echoes, and saturation
recovery in glasses at low temperatures,” Phys. Rev. B, vol. 16, pp. 2879–2895, Sep
1977.
[42] M. Schechter and P. C. E. Stamp, “What are the interactions in quantum glasses?,”
Journal of Physics: Condensed Matter, vol. 20, no. 24, p. 244136, 2008.
[43] J. Joffrin and A. Levelut, “Virtual phonon exchange in glasses,” Journal de Physique,
vol. 36, no. 9, pp. 811–822, 1975.
[44] K. Kassner and R. Silbey, “Interactions of two-level systems in glasses,” Journal of
Physics: Condensed Matter, vol. 1, no. 28, p. 4599, 1989.
[45] K. Agarwal, I. Martin, M. D. Lukin, and E. Demler, “Polaronic model of two-level
systems in amorphous solids,” Phys. Rev. B, vol. 87, p. 144201, Apr 2013.
92
[46] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger,
“Dynamics of the dissipative two-state system,” Rev. Mod. Phys., vol. 59, pp. 1–85, Jan
1987.
[47] A. Burin, “Universal properties of glasses in a pseudospin model,” JETP Lett, vol. 54,
no. 6, 1991.
[48] S. Baranovskii, B. Shklovskii, and A. Efros, “Elementary excitations in disordered systems with localized electrons,” Soviet Journal of Experimental and Theoretical Physics,
vol. 51, p. 199, 1980.
[49] A. S. Nowick, Anelastic relaxation in crystalline solids, vol. 1. Elsevier, 1972.
[50] G. Leibfried and N. Breuer, Point defects in metals I: introduction to the theory. Springer
tracts in modern physics, Springer, 1978.
[51] P. Esquinazi, R. Knig, and F. Pobell, “Acoustic properties of amorphous sio2 and
pdsicu, and of crystalline ag, nbti and ta at very low temperatures,” Zeitschrift fr
Physik B Condensed Matter, vol. 87, no. 3, pp. 305–321, 1992.
[52] M. O. Michael Pollak and A. Frydman, The Electron Glass. New York: Cambridge
University Press, 2013.
93