Supplementary Material for:
Nonharmonic contributions to the high-temperature phonon thermodynamics of Cr
C. M. Bernal-Choban
§
,
1,
∗
H. L. Smith
§
,
2,
†
C. N. Saunders,
1
D. S. Kim,
3
L. Mauger,
4
D. L. Abernathy,
5
and B. Fultz
1,
‡
1
Applied Physics and Materials Science, California Institute of Technology, Pasadena, CA 91125, USA
2
Physics and Astronomy, Swarthmore College, Swarthmore, PA 19081, USA
3
Department of Materials Science and Engineering,
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
4
Jet Propulsion Laboratory, Pasadena, CA 91109, USA
5
Neutron Scattering Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
I.
MULTIPLE SCATTERING
The experiment was designed to minimize the effects
of multiple scattering. This was achieved by choosing a
sufficiently thin sample of Cr. The ratio of double (mul-
tiple) scattering to single scattering is
1
/
2 (
t/τ
)
,
(1)
where
t
is the sample thickness and
τ
is the characteristic
scattering length of
τ
= 1
/
(
ρσ
)
.
(2)
Here
σ
= 3
.
5 barns for the total scattering of a Cr nucleus
and
ρ
= 8
.
3
×
10
22
nuclei/cm
3
. For the Cr pieces in
this experiment, which were
∼
0.15 cm thick, the ratio of
single to double scattering is 0.022, or 2%. Despite some
directions in the sample where the neutron paths were
longer, the sample is still thin enough to neglect multiple
scattering.
II.
AB INITIO MODELING
A. Quasi-harmonic
Phonon calculations with the quasi-harmonic (QH)
approximation used the Phonopy package [1]. The fi-
nite difference method was used to generate forces on
atoms. Force constant matrices were constructed and
transformed into the dynamical matrix,
D
(
q
)
e
q
s
=
ω
2
q
s
e
q
s
,
(3)
where
D
(
q
) is the dynamical matrix,
q
is the phonon
wave vector,
s
is the phonon band index,
ω
q
s
and
e
q
s
are the frequency and polarization vector of the phonon
∗
cmbchoban@gmail.com
†
hsmith@swarthmore.edu
‡
btf@caltech.edu
§
Equally contributing authors.
mode with
q
, s
. With all
{
ω
q
s
}
, the vibrational contri-
bution to the free energy is
F
vib
=
1
2
X
q
s
ℏ
ω
q
s
+
k
B
T
X
q
s
ln
1
−
exp (
−
ℏ
ω
q
s
/k
B
T
)
.
(4)
Temperatures from 0-1500 K were considered us-
ing a grid of 0 K supercells scaled in volume by
±
0
.
5%
,
±
1%
,
±
1
.
5%. The Helmholtz free energy for each
volume was calculated and minimized at different tem-
peratures using the Birch-Murnaghan equation of state.
Minima of the fitted free energy curves gave the QH vol-
umes at the temperatures of interest. Quasi-harmonic
phonon properties were calculated with these volumes.
The calculated phonon density of states was convoluted
with a Gaussian of 1.0 meV to approximate the broaden-
ing due to instrumental resolution.
B. Anharmonic
Phonon calculations with the anharmonic (AH) ap-
proximation used the stochastic Temperature Dependent
Effective Potential Method (sTDEP) [2]. Supercells with
displaced atoms were created by stochastic sampling. For
a cell of
N
a
atoms with mass
m
i
, the atomic positions,
{
u
i
}
, were created using a harmonic normal-mode trans-
formation,
u
i
=
3
N
a
X
s
=1
e
is
⟨
A
is
⟩
p
−
2 ln
ξ
1
sin (2
πξ
2
)
,
(5)
where
ξ
n
are uniformly distributed numbers between 0
and 1 (the Box-Muller transform). The thermal ampli-
tude,
⟨
A
is
⟩
, of normal mode
s
with eigenvector
e
is
is
⟨
A
is
⟩
=
1
ω
s
s
ℏ
ω
s
(
n
s
+
1
2
)
m
i
≈
1
ω
s
r
k
B
T
m
i
,
(6)
where
n
s
= (
e
ℏ
ω
s
/k
B
T
−
1)
−
1
is the thermal occupation
of phonon mode
s
and
ℏ
ω
≪
k
B
T
denotes the classical
limit at high temperatures.
Energies and forces from calculations of an ensemble
of supercells with thermally displaced atoms were ob-
tained with a least squares fit to a model Hamiltonian.
2
A grid of four temperatures,
{
0
,
330
,
1000
,
1500
}
K, and
nine volumes was created. For each volume-temperature
point, ten thermally displaced supercells were stochas-
tically generated using force constants from a model
pair potential. Static calculations were performed on
each supercell to obtain energy-force-displacement data.
Quadratic and cubic force constants were re-generated
for each temperature using these data sets. The free en-
ergy surface,
F
(
V,T
) =
U
0
(
V,T
) +
F
vib
(
V,T
)
,
(7)
for each volume-temperature point was calculated using
sTDEP. The baseline,
U
0
(
V,T
), and the free energy from
lattice vibrations,
F
vib
=
Z
∞
0
g
(
ω
)
n
k
B
T
ln
h
1
−
e
−
ℏ
ω
k
B
T
i
+
ℏ
ω
2
o
dω,
(8)
depend explicitly on volume and temperature. These
free energies were fit to the Birch-Murnaghan equation
of state to find the optimized volume for a given temper-
ature.
The phonon DOS,
g
(
ω
) =
X
s
δ
(
ω
−
ω
s
)
,
(9)
also depends explicitly on volume and temperature.
Shifts, ∆
s
, and linewidths, Γ
s
, of phonons arise from
AH, or phonon-phonon, interactions. To calculate these
renormalized frequencies, results from many-body per-
turbation theory were used. The self-energy correction
of phonons is
X
(Ω) = ∆(Ω) +
i
Γ(Ω)
,
(10)
where
E
=
ℏ
Ω is a probing energy. Linewidths are ob-
tained from the imaginary component,
Γ
s
(Ω) =
ℏ
π
16
X
s
′
s
′′
|
Φ
ss
′
s
′′
|
2
×{
(
n
s
′
+
n
s
′′
+ 1)
δ
(Ω
−
ω
s
′
−
ω
s
′′
) + (
n
s
′
−
n
s
′′
)
×
[
δ
(Ω
−
ω
s
′
+
ω
s
′′
)
−
δ
(Ω +
ω
s
−
ω
s
′′
)]
}
.
(11)
This sum is taken over all possible three-phonon interac-
tions, where Φ
ss
′
s
′′
is the three-phonon matrix element
obtained from the cubic force constants.
A Kramers-Kronig transformation was used to get the
frequency shifts,
∆(Ω) =
1
π
Z
Γ(
ω
)
(
ω
−
Ω)
dω.
(12)
As shown in Eq. 10, these shifts are the real part of the
phonon self-energy.
Anharmonic phonon densities of states curves were cal-
culated using the real and imaginary parts of the phonon-
self energy,
g
anh
(
ω
) =
X
s
2
ω
s
Γ
s
(
ω
)
[
ω
2
−
ω
2
s
−
2
ω
s
∆
s
(
ω
)]
2
+ 4
ω
2
s
Γ
2
s
(
ω
)
.
(13)
In the limit where ∆, Γ
→
0 Eq. 13 reduces to Eq. 9.
III.
BORN-VON K
́
ARM
́
AN MODELING
The Born-von K ́arm ́an formalism relates interatomic
force constants of crystalline solids to phonon frequen-
cies. To solve for phonon frequencies, the interatomic
interactions in the dynamical matrix are truncated to
the first few nearest-neighbor interactions. Below is a
description of how the force constants are extracted from
the DOS of body-centered cubic (bcc) Cr using a genetic
algorithm to obtain temperature-dependent tensorial and
radial force constants [3].
A. Genetic algorithm
A genetic algorithm optimization was performed to
obtain interatomic force constants from experimental
phonon DOS employing the open-source package mystic
[4]. To begin, the algorithm was used to generate many
candidate solutions called a population. Each popula-
tion consisted of a potential set of force constants that
was randomly generated within a set of reasonable force
constant bounds, where
K
1
NN
=
1XX 1XY 1XY
1XY 1XX 1XY
1XY 1XY 1XX
,
K
2
NN
=
2XX
0
0
0
2YY
0
0
0
2YY
(14)
are the first and second nearest-neighbor tensorial force
constant matrices.
The dynamical matrix was con-
structed from each population (set of force constants)
and used to generate a DOS. These DOS are compared
with the experimental DOS, and solutions that best re-
produce the experimental DOS (by minimizing the mean
squared error) are selected as the parents of the next
generation. The parent solutions seeded the next pop-
ulation of potential force constants by selecting random
combinations of parent parameters and also introducing
random changes, or mutations. Least squares fitting and
comparison to experiments were repeated until the popu-
lation converged on a set of force constants that provided
the best fit to the experimental density of states (Fig. 1).
Throughout this process, solutions that generated nega-
tive frequencies were discarded due to their implication of
dynamical instability in the lattice. Fully optimized force
constants were then used to generate phonon dispersions
with
q
-space resolved information.