# Operational experience and evolution of the ATLAS Tile Hadronic Calorimeter
Read-Out Drivers

A. Valero, on behalf of the ATLAS Tile Calorimeter Collaboration A. Valero is
with Instituto de Física Corpuscular (Universidad de Valencia-CSIC), 46980
Valencia, Spain (e-mail: alberto.valero@cern.ch).This work was supported by
the Ministerio de Ciencia e Innovación under projects RTI2018-094270-B-I00 and
RTI2018-100863-B-I00.Copyright 2020 CERN for the benefit of the ATLAS
Collaboration. CC-BY-4.0 license

###### Abstract

TileCal is the central hadronic calorimeter of the ATLAS experiment at the
Large Hadron Collider (LHC). It is a sampling detector where scintillating
tiles are embedded in steel absorber plates. The tiles are grouped forming
cells, which are read-out on both sides by photomultiplier tubes (PMTs). The
PMT digital samples are transmitted to the Read-Out Drivers (ROD) located in
the back-end system for the events accepted by the Level 1 trigger system. The
ROD is the core element of the back-end electronics and it represents the
interface between the front-end electronics and the ATLAS overall Data
AcQuisition (DAQ) system. It is responsible for energy and time
reconstruction, trigger and data synchronization, busy handling, data
integrity checking and lossless data compression. The TileCal ROD is a
standard 9U VME board equipped with DSP based Processing Units mezzanine
cards. A total of 32 ROD modules are required to read-out the entire TileCal
detector. Each ROD module has to process the data from up to 360 PMTs in real
time in less than 10 µs. The commissioning of the RODs was completed in 2008
before the first LHC collisions. Since then, several hardware and firmware
updates have been implemented to accommodate the RODs to the evolving ATLAS
Trigger and DAQ conditions adjusted to follow the LHC parameters. The initial
ROD system, the different updates implemented and the operational experience
during the LHC Run 1 and Run 2 are presented.

###### Index Terms:

Calorimetry, Digital filters, Digital signal processors, Field programmable
gate arrays

##  I Introduction

TileCal [1] is the central hadronic calorimeter of the ATLAS experiment [2] at
the Large Hadron Collider (LHC). It is a sampling detector where scintillating
tiles are embedded in steel absorber plates. The tiles are grouped forming
cells, which are read-out on both sides by photomultiplier tubes (PMTs). The
electric signals produced by the approximately 10000 PMTs are digitized at 40
MHz synchronously with the LHC bunch crossing. The digital samples are stored
in pipelined memories located in the front-end electronics [3]. Seven samples
are transmitted to the Read-Out Drivers (ROD) [4] located in the back-end
system if the Level 1 (L1) trigger accepts the event (Figure 1). The
commissioning of the RODs was completed in 2008 before the first LHC
collisions. Since then, several hardware and firmware updates have been
implemented to accommodate the RODs to the evolving ATLAS Trigger and DAQ
conditions [5] adjusted to follow the evolution of the LHC parameters.

##  II The TileCal read-out driver modules

The ROD is the core element of the TileCal back-end electronics. It represents
the interface between the front-end electronics and the ATLAS overall Data
AcQuisition (DAQ) system. The TileCal ROD module is a standard 9U VME board
equipped with DSP based Processing Units mezzanine cards (Figure 2) and
connected to a Rear Transition Module (RTM) through the VME backplane. The
RODs are responsible for the energy and time reconstruction, trigger and data
synchronization, busy handling, data integrity checking and lossless data
compression. A total of 32 ROD modules are required to read-out the entire
TileCal detector. Each ROD module has to process the data from up to 360 PMTs
in real time in less than 10 µs established by the maximum L1 trigger rate
[6]. A digital Finite Impulse Response (FIR) filter called Optimal Filtering
(OF) is used in the DSPs to reconstruct the amplitude and phase of the pulses
(Figure 3). It exploits the knowledge of the pulse shape and noise of the
electronics and the amount of expected pileup to reduce the contribution of
noise and determine the time of deposition.

Figure 1: Sketch of the TileCal readout electronics.

The input data events and the output-reconstructed fragments are stored in the
DSPs in elastic buffers with configurable depth. A veto signal (busy) is
generated in the DSP and propagated to the ATLAS Central Trigger Processor to
stop the L1 trigger generation when the input buffer is almost full. This
mechanism prevents the overwriting of events in the input buffer with the
consequent data loss but it introduces undesired deadtime in the detector. The
processing time, the output data bandwidth, the depth of the buffers and the
almost full flag are the key parameters to reduce the deadtime introduced by
the RODs.

Figure 2: Block diagram of a ROD system equipped with two Processing Units and
two HOLA cards.

###  II-A The Optimal Filtering method

Optimal Filtering [7] is a relatively simple algorithm used to reconstruct in
real time the energy and phase of the PMT pulses in the DSPs of the RODs. The
OF method relies on the usage of a digitization clock synchronized with the
trigger, thus the signal pulses and the samples have always a fixed phase with
small variations [7]. This feature permits the discrimination of out-of-time
pileup energy depositions.

The amplitude and phase of the real pulse is obtained as a linear combination
of the digital samples ($S_{i}$) and the weights ($a_{i}$, $b_{i}$). There are
two versions of the OF algorithm; OF1 subtracts the pedestal ($p$) from the
samples before applying the filter (1), whereas OF2 includes an additional
constraint in the weights ($\sum_{i=1}^{n}a_{i}=0$) which implies that any
common variation in all the samples is cancelled (2):

$A=\sum_{i=1}^{n}a_{i}(S_{i}-p)\;\;,\;\;\;A\tau=\sum_{i=1}^{n}b_{i}(S_{i}-p),$
(1)

$A=\sum_{i=1}^{n}a_{i}S_{i}\;\;,\;\;\;A\tau=\sum_{i=1}^{n}b_{i}S_{i}.$ (2)

The weights are computed from the known pulse shape, the expected phase and
electronic and pileup noise. The expected phase and electronics noise is
measured with calibration data and they are independent of the LHC conditions.
On the other hand, the pileup noise depends on the LHC beam conditions and the
cell position. Thus, any pulse shape distortion strongly affects the
performance of the Optimal Filtering method. In particular, signal pileup
different than expected deforms the signal of interest and biases the results
of the signal reconstruction. In addition, there is an iterative version of OF
where the weights are selected according to the phase reconstructed in the
previous iteration. This iterative method optimizes the reconstruction when
the trigger is not synchronous with the digitizing clock, like in cosmic runs,
and when the bunch spacing is larger than the 7 samples window ($\pm$75 ns)
which ensures the absence of out-of-time pileup.

Figure 3: Representation of the TileCal pulse shape with the Optimal Filtering
reconstructed magnitudes.

##  III Evolution of data processing algorithms during Run 1

The first LHC beams collided in ATLAS in 2009. Since then, the LHC parameters
like bunch spacing and number of interactions per bunch crossing evolved,
gradually reaching 75% of the nominal instantaneous luminosity in 2012. The
TileCal data acquisition system followed this evolution to optimize its
performance and efficiency. The read-out format and the reconstruction
algorithms were optimized with various firmware upgrades.

###  III-A Commissioning and early Run 1 operation

During the first year the LHC was operated with larger bunch spacing and lower
instantaneous luminosity than the nominal design parameters [8]. Under these
conditions the ATLAS L1 trigger rate was reduced, which permitted the usage of
large time consuming reconstruction algorithms and big output data fragments
in the RODs with a data taking efficiency close to 100%. The iterative version
of the OF algorithm was used to minimize the effect of large pulse shifts.
This method is capable of detecting the peak of the pulse and apply different
filter weights according to its position. The output data fragment included
both the online reconstructed magnitudes and the front-end digital samples for
all the channels. Then, the digital samples were reconstructed offline to
validate and certify the online reconstructed magnitudes.

###  III-B ROD firmware optimizations and performance during Run 1

During the second half of Run 1 there was an increase of the LHC instantaneous
luminosity and a reduction of the bunch spacing and the consequent increase in
the L1 trigger rate forced a change in the read-out strategy. In order to
minimize the effect of out-of-time pileup and to reduce the computing time,
the OF iterative method was replaced by a non-iterative algorithm using pre-
calibrated filter weights for each channel. The increase in the L1 trigger
rate reduced the output bandwidth available for each event fragment. A
lossless data compression algorithm was implemented to reduce the output
fragment size to avoid link saturation and the resulting undesirable deadtime.
Figure 4 shows the ROD output fragment size for each of the 64 ROD output
Read-Out Links (ROL) (32 RODs, 2 ROL each) for a 2012 run with a peak
luminosity of $3.5\times 10^{33}$ cm-2s-1 for the legacy (a) and the
compressed (b) output dataformats. The plots show the limit for a free of
deadtime operation for a L1 trigger rate of 100 kHz and 75 kHz. The legacy
dataformat (Figure 4a) exceeds the limit for 100kHz and thus the data taking
efficiency was affected. With the compressed dataformat (Figure 4b) it is
possible to run at 100 kHz without deadtime from the RODs. The data
compression algorithm relies on the fact that the majority of the channels
have a pedestal-like structure and only few bits are needed to pack the
samples.

ATLAS recorded in Run 1 21.3 fb-1 of data with a data taking efficiency of
93%. After these optimizations, TileCal data taking efficiency at the end of
Run 1 was around 95%, mainly due to failures in the front-end low voltage
power supply system. The TileCal data quality efficiency was 98.7% [9].

Figure 4: ROD output fragment size per Read-Out Link for the old (a) and new
(b) dataformat [10].

##  IV ROD system upgrade for Run 2

In view of the expected LHC instantaneous luminosity increase for Run 2, the
different ATLAS sub-systems used the Long Shutdown 1 (2013–2014) to
consolidate and upgrade their read-out systems. In particular, TileCal
repaired and consolidated the front-end electronics components and replaced
the low voltage power supplies with a more radiation tolerant version. As
presented in the next sections, the ROD system was also upgraded to operate
with the new LHC conditions. The LHC Run 2 started in 2015 reaching a peak
luminosity of two times the nominal value ($2\times 10^{34}$ cm-2s-1 ) in
2016. The number of proton–proton interactions per bunch crossing exceeded the
nominal values thus producing very high detector occupancy, signals pileup and
high L1 trigger rates.

###  IV-A Motivation for a ROD upgrade

During the Long Shutdown 1 TileCal undertook a major hardware upgrade in the
ROD system in order to cope with the expected evolution of the LHC parameters
and in particular the increase of the ATLAS L1 trigger rate while keeping a
high data taking efficiency. The ROD processing power and the output data
bandwidth were doubled. Two Processing Units and two output High-speed Optical
Link for Atlas (HOLA) cards per ROD were installed in the available empty
slots (Figure 5). A total of 64 HOLA cards dismantled from the ATLAS Muons
system were installed in the ROD RTMs. The Processing Units cards were
produced and certified in the laboratory test-bench during the first year of
the Long Shutdown 1. The second year of the Long Shutdown 1 was used to
install the 64 Processing Units in the empty ROD slots and make the
commissioning of the new upgraded ROD system. The ROD firmware and data
acquisition control software was adapted to run with the new hardware
conditions.

Figure 5: Picture of the upgraded ROD system (ROD motherboard and RTM)
equipped with 4 Processing Units and 4 HOLA cards.

###  IV-B ROD upgrade and performance during Run 2

The signal reconstruction algorithm was updated to avoid a feature of the OF2
method which reconstructs negative energies for large out-of-time pulses. The
increase of the pileup above the nominal value enlarged this effect which was
specially affecting high level trigger algorithms based on the total
transverse energy in the calorimeter. Thus, the OF2 method was replaced with
OF1 which requires a periodical pedestal calibration but does not provide
negative energies when it is properly measured. In addition, the pileup noise
was introduced in the calculation of the OF weights. Since the pileup noise
depends on the cell position, different sets of weights per cell are now used
and the database and procedure for loading the values into the DSPs were
updated accordingly.

These ROD hardware upgrades and firmware optimizations permitted a successful
operation during Run 2 and it should provide stable and smooth operation
during the coming Run 3 that should start in 2021 after the Long Shutdown 2
(2019–2020). The ATLAS overall data taking efficiency was 94% during Run 2 and
the data quality losses from TileCal were below 0.4%.

##  V Conclusions and prospects for Run 3 and beyond

The TileCal ROD system has evolved following the requirements imposed by the
evolution of the LHC during Run 1 and Run 2. Overall the TileCal data
acquisition system has operated successfully during this period with data
taking and quality efficiencies compatible with the rest of the ATLAS sub-
systems. Moreover, the ROD system is ready to operate smoothly during Run 3
without major modifications.

The ATLAS read-out strategy will be radically changed for the High Luminosity
LHC where an increase of up to 7.5 times the LHC nominal instantaneous
luminosity is expected. TileCal will replace the read-out electronics and most
of the front-end functionalities will be moved to the off-detector
PreProcessor (PPr) modules which are the natural evolution of the current RODs
[11] . The PMT signals will be digitized and transmitted to the PPr modules
before any event selection is applied (Figure 6). The PPr will be the
interface with the trigger and the ATLAS overall data taking system [11].

Figure 6: Sketch of the TileCal readout electronics for the HL-LHC.

## References

  * [1] ATLAS Collaboration, “ATLAS Tile Calorimeter Detector Technical Design Report,” CERN, Geneva, Technical Design Report ATLAS, Switzerland, Rep. CERN-LHCC 96-42, 1996.
  * [2] ATLAS Collaboration, “The ATLAS Experiment at the CERN Large Hadron Collider,” Journal for Instrumentation, vol. 3, no. S08003, 2008, DOI. 10.1088/1748-0221/3/08/S08003.
  * [3] K. Anderson et al., “Design of the front-end analog electronics for the ATLAS Tile Calorimeter,” Nucl. Instrum. Meth., vol. A, no. 551, 2005, DOI. 10.1016/j.nima.2005.06.048.
  * [4] A. Valero et al., “ATLAS TileCal Read Out Driver production,” Journal for Instrumentation, vol. 2, no. P05003, 2007, DOI. 10.1088/1748-0221/2/05/P05003.
  * [5] ATLAS Collaboration, “ATLAS high-level trigger, data-acquisition and controls: Technical Design Report,” CERN, Geneva, Technical Design Report ATLAS, Switzerland, Rep. CERN-LHCC-2003-022, 1998.
  * [6] ATLAS Collaboration, “ATLAS level-1 trigger: Technical Design Report,” CERN, Geneva, Technical Design Report ATLAS, Switzerland, Rep. CERN-LHCC-98-014, 2003.
  * [7] A. Valero et al., “DSP Online algorithms for the ATLAS TileCal Read-Out Drivers,” Nuclear Science, IEEE Transactions on Nuclear Science., vol. 55, no. 1, p.158-164, Feb. 2008, DOI. 10.1109/RTC.2007.4382840.
  * [8] ATLAS Collaboration et al., ”Readiness of the ATLAS Tile Calorimeter for LHC collisions,” Eur. Phys. J. C., no. 2010-024, vol. 70, p. 1193-1236, Jul. 2010, DOI. 10.1140/epjc/s10052-010-1508-y.
  * [9] ATLAS Collaboration, ”Operation and performance of the ATLAS Tile Calorimeter in Run 1,”Eur. Phys. J. C., vol. C, no. 78, Dec. 2018, DOI. 10.1140/epjc/s10052-018-6374-z
  * [10] ATLAS Collaboration, ”ATLAS Experiment - Public Results,” [Online]. Available: https://twiki.cern.ch/twiki/bin/view/AtlasPublic/TileCaloPublicResultsSignalReadOut. Accessed on: October 21, 2020.
  * [11] ATLAS Collaboration, “Technical Design Report for the Phase-II Upgrade of the ATLAS Tile Calorimeter,” CERN, Geneva, Technical Design Report ATLAS, Switzerland, Rep. CERN-LHCC-2017-019. ATLAS-TDR-0282, Sept. 2017.


# Understanding mechanisms of thermal expansion in PbTiO3 thin-films from
first principles: role of high-order phonon-strain anharmonicity

Ethan T. Ritz  Sibley School of Mechanical and Aerospace Engineering, Cornell
University, Ithaca, New York 14853, USA  Nicole A. Benedek
<EMAIL_ADDRESS>Department of Materials Science and Engineering, Cornell
University, Ithaca, New York 14853, USA

###### Abstract

The thermal properties of materials are critically important to various
technologies and are increasingly the target of materials design efforts.
However, it is only relatively recent advances in first-principles
computational techniques that have enabled researchers to explore the
microscopic mechanisms of thermal properties, such as thermal expansion. We
use the Grüneisen theory of thermal expansion in combination with density
functional calculations and the quasiharmonic approximation to uncover
mechanisms of thermal expansion in PbTiO3 thin-films in terms of elastic and
vibrational contributions to the free energy. Surprisingly, we find that
although the structural parameters of PbTiO3 thin-films evolve with
temperature as if they are dominated by linear elasticity, PbTiO3 thin-films
are strongly anharmonic, with large changes in the elastic constants and
Grüneisen parameters with both misfit strain and temperature. We show that a
fortuitous near-cancellation between different types of anharmonicity gives
rise to this behavior. Our results illustrate the importance of high-order
phonon-strain anharmonicity in determining the temperature-dependent
structural parameters of PbTiO3 thin-films, and highlight the complex manner
in which thermal expansion, misfit strain and elastic and vibrational
properties are intertwined.

Suggested keywords

###### pacs:

Valid PACS appear here

††preprint: AIP/123-QED

##  I Introduction

Recent advances in theoretical and experimental techniques have promoted
increasing interest in the phonon and thermal properties of materials by
making it possible to more easily access the information required to elucidate
atomic-scale mechanisms.Tucker  _et al._ (2005); Mounet and Marzari (2005);
Kim  _et al._ (2018); Ritz and Benedek (2018) Particularly exciting progress
has been made in the family of metal-organic frameworks, where it now not only
appears possible to design materials with unconventional thermal responses
(such as negative thermal expansion), but also to tune those responses in a
highly controllable way Goodwin  _et al._ (2008a, b, 2009). In contrast,
negative thermal expansion is both rarer and apparently more difficult to tune
in inorganic framework materials, such as complex oxides Lind (2012); Miller
_et al._ (2009); Barrera  _et al._ (2005). The most common approach is
chemical substitution, although it is generally challenging to change either
the sign or the magnitude of the thermal expansion significantly in this way
(recent work provides an impressive exception Yanase  _et al._ (2011); Ablitt
_et al._ (2019)). A significant barrier to designing effective strategies for
highly tunable thermal expansion is that, with few exceptions, the microscopic
factors controlling the temperature evolution of structural parameters are
poorly understood.

The ferroelectric $P4mm$ tetragonal phase of PbTiO3 exhibits negative
volumetric thermal expansion from approximately room temperature up to 760 K,
where it undergoes a structural phase transition to a cubic $Pm\bar{3}m$ phase
Chen  _et al._ (2005); Shirane and Hoshino (1951). Thermal expansion in the
ferroelectric phase is characterized by a shrinking $c$-axis and an expanding
$a$-axis – the $c$-axis contracts faster than the $a$-axis expands, so the net
effect is a decrease in volume with increasing temperature. Interestingly,
recent theoretical work has shown that negative thermal expansion in PbTiO3 is
actually _driven_ by the expanding $a$-axis Ritz and Benedek (2018). However,
strong elastic (Poisson-like) coupling between the $a$\- and $c$-axes pulls
the $c$-axis down as the $a$-axis expands.

It is possible to grow high-quality thin-films of PbTiO3 on various
substrates, and a wealth of temperature-dependent structural and ferroelectric
polarization data are available in the literature. Thin-film systems are
excellent platforms for probing mechanisms of thermal expansion because
different substrates will impart different misfit strains (which will change
with temperature, depending on the thermal expansion coefficient of the
substrateJanolin  _et al._ (2007); Janolin (2009)) and hence different elastic
boundary conditions to the growing film. Hence, it becomes possible to
understand how changes in the unit cell dimensions affect thermal expansion
behavior, all while keeping the chemical composition constant.

In previous work, we used first-principles density functional theory (DFT) and
the quasiharmonic approximation (QHA) to investigate the thermal expansion of
PbTiO3 thin-films on SrTiO3, DyScO3 and
(La0.29(5)Sr0.71(5))A-site(Al0.65(1)Ta0.35(1))B-siteO3 (LSAT) substrates. Our
key findings were that: 1) the misfit strain does not just change the
properties of PbTiO3 thin-films, it changes how they change with temperature,
and 2) the mismatch in thermal expansion coefficients between PbTiO3 and the
substrates had a striking effect on thin-film properties. For a given misfit
strain, we found qualitative changes in thermal expansion behavior and
ferroelectric transition temperatures depending on the thermal expansion
coefficient of the substrate.

In this work, we use the Grüneisen theory of thermal expansion in combination
with DFT calculations and the quasiharmonic approximation to elucidate the
microscopic mechanisms underlying the structural properties of PbTiO3 thin-
films in terms of elastic and vibrational contributions to the lattice
parameters at finite temperatures. Surprisingly, we show that the $c$-axis
lattice parameter (perpendicular to the substrate) evolves with temperature as
if it is almost entirely dominated by linear elastic effects, despite the fact
that PbTiO3 thin-films exhibit the hallmarks of strong anharmonicity. We find
that this behavior is due to a fortuitous near-cancellation between the
contributions to thermal expansion from anharmonic elasticity and from
vibrational thermal stress, the details of which depend both on the mismatch
between the lattice parameters of the film and substrate, as well as the
mismatch between their rates of thermal expansion. We use our findings to
explain the temperature evolution of the structural, elastic, vibrational and
thermal expansion properties of PbTiO3 thin-films on experimentally available
substrates.

##  II First-Principles Calculations

All calculations were performed using density functional theory, as
implemented in Quantum Espresso 6.5.0. Giannozzi  _et al._ (2009) We used the
Wu-Cohen exchange-correlation functional Wu and Cohen (2006) with Garrity-
Bennett-Rabe-Vanderbilt ultrasoft pseudopotentials.Garrity  _et al._ (2014)
The following states were included in the valence for each element: 5d106s26p2
for Pb, 3s23p64s23d1 for Ti, and 2s22p4 for O. Zero temperature unit cell
lattice parameters and atomic positions of $P4mm$ PbTiO3 were converged with
respect to the plane wave cutoff energy and $\mathbf{k}$-point mesh density to
within 0.001 Å. Unless otherwise mentioned, structural parameters were found
to be converged at a force cutoff threshold of $3.0\times 10^{-5}$ Ry/bohr
using a 6 $\times$ 6 $\times$ 6 Monkhorst-Pack (MP) mesh and a 100 Ry plane
wave cutoff energy, compared with MP meshes up to 12$\times$12$\times$12 and
plane wave cutoffs up to 120 Ry. Phonon dispersion calculations were performed
using density functional perturbation theory on an 8 $\times$ 8 $\times$ 8
$\mathbf{q}$-point grid.

Finite temperature structural parameters and elastic constants were predicted
using a quasiharmonic approximation (QHA) to the Helmholtz free energy – in
this study, we closely follow the framework outlined in Refs. 20 and Ritz and
Benedek, 2020 regarding the application of the QHA. Our grid of strained
systems consisted of a 6 $\times$ 12 grid spanning -1.5 to +2.5% strain along
the $a$-axis of the $P4mm$ tetragonal phase, and -4.0% to +5.1% strain along
$c$, augmented with 20 additional points with strains spanning -1.1% to -0.3%
strain in $a$ and -4.0% to +5.1% strain in $c$. These points were chosen to
ensure good convergence of finite-temperature elastic constants and to make
sure that for bulk PbTiO3 and all strained films studied, the lattice
parameters as a function of temperature would be bounded by these values ins .
Some structures at very large values of strain exhibited phonons with
imaginary frequencies; these points were discarded from the data set. Note
that the strain values above are defined with respect to the 0 K lattice
parameters of the $P4mm$ phase, including zero-point energy corrections from
vibrational degrees of freedom (we find $a$ = 3.891 Å and $c$ = 4.174 Å).

We use the same definitions of misfit strain and thermal expansion mismatch as
in Ref. Ritz and Benedek, 2020; we reproduce the definitions here for
convenience. For the calculations of epitaxially strained PbTiO3, we clamp the
in-plane ($a$ and $b$) lattice parameters of PbTiO3 to the temperature-
dependent lattice parameters of the substrate. That is, our misfit strain is
defined as,

$\varepsilon_{a}(T)=\frac{a_{\mathrm{substrate}}(T)-a_{\mathrm{PTO}}(T)}{a_{\mathrm{PTO}}(T)},$
(1)

where $a_{\mathrm{substrate}}(T)$ is the lattice parameter of the substrate at
some temperature $T$, and $a_{\mathrm{PTO}}(T)$ is the in-plane lattice
parameter of bulk PbTiO3 at the same temperature from our QHA calculations. We
can re-write this (as in Ref. 21) as

$\varepsilon_{a}(T)=\varepsilon^{300}_{a}+\varepsilon^{thermal}_{a}(T),$ (2)

where $\varepsilon^{300}_{a}$ is the strain at a reference temperature of 300
K (the misfit strain is usually defined at room temperature in the thin-film
literature by convention), and the additional strain that arises from changes
in temperature $\varepsilon_{a}^{\mathrm{thermal}}$ is defined as

$\begin{split}\varepsilon_{a}^{\mathrm{thermal}}(T)=\int_{\tau=300}^{T}\alpha_{a}(\tau)-\alpha_{a}^{\mathrm{bulk}}(\tau)\mathrm{d}\tau\\\
=\int_{\tau=300}^{T}\Delta\alpha_{a}(\tau)\mathrm{d}\tau,\end{split}$ (3)

The misfit strain rate is defined as
$\Delta\alpha_{a}(T)=\alpha_{a}(T)-\alpha_{a}^{\mathrm{bulk}}(T)$, where
$\alpha_{a}(T)$ is the thermal expansion coefficient of the substrate, and
$\alpha^{\mathrm{bulk}}_{a}(T)$ is the rate of thermal expansion along $a$ for
unconstrained bulk PbTiO3 (since we are considering PbTiO3 films coherently
clamped to the substrate, the thermal expansion coefficient of PbTiO3 films
along the $a$-axis is just that of the substrate).

##  III  Grüneisen Theory of Thermal Expansion

In our previous workRitz and Benedek (2020) on thermal expansion in
ferroelectric PbTiO3 thin-films, we showed that the temperature of the
transition ($T_{c}$) to the cubic $Pm\bar{3}m$ phase depends _qualitatively_
not only on the misfit strain but also on the thermal expansion coefficient of
the substrate ($\alpha_{a}$). For example, Figure 1 (reproduced from Ref. Ritz
and Benedek, 2020) shows that for a tensile misfit strain of 1% at 300 K,
corresponding approximately to growth of PbTiO3 on the substrate DyScO3,
$T_{c}$ can be either higher or lower than bulk PbTiO3, depending on the value
of $\alpha_{a}$.

Figure 1: Variation in ferroelectric transition temperature and $c$-axis
behavior of ferroelectric PbTiO3 thin-films as a function of misfit strain and
substrate coefficient of thermal expansion ($\alpha_{a}$) from our first-
principles QHA calculations reported in Ref. Ritz and Benedek, 2020. The color
chart for $T_{c}$ has white set to the bulk PbTiO3 transition temperature of
760 K. Blue (red) squares indicate a strain and $\alpha_{a}$ combination that
produces films with a lower (higher) $T_{c}$ than bulk. Triangles indicate
combinations of strain and $\alpha_{a}$ that produce films in which the
$c$-axis continues to grow with temperature and $T_{c}$ is suppressed. All
other combinations of strain and $\alpha_{a}$ produce films in which the
$c$-axis shrinks with temperature and $T_{c}$ is finite. Note that the films
with the highest transition temperatures may exceed 1600 K or even be
completely suppressed. The letters ‘L’, ‘S’ and ‘D’ denote the strain and
$\alpha_{a}$ conditions corresponding to growth on LSAT, SrTiO3 and DyScO3
substrates. See Ref. Ritz and Benedek, 2020 for further details.

There is an intimate and complex connection between the transition
temperature, the misfit strain, the thermal expansion coefficient of the
substrate, and the thermal expansion coefficient of the out-of-plane $c$-axis
($\alpha_{c}$) of the PbTiO3 thin-film. Rather than continuing to focus on
transition temperatures, in this work we shift our attention to the
fundamental relationship between the misfit strain, $\alpha_{c}$ and the
thermal expansion coefficient of the substrate. Our goal is to elucidate the
microscopic origins of this relationship in terms of the elastic and
vibrational properties of PbTiO3 thin-films. That is, we aim to answer the
question, how can we understand differences in the thermal expanion properties
of PbTiO3 thin-films compared to bulk in terms of fundamental elasticity and
vibrational properties?

###  III.1 The thermal expansion and Grüneisen tensors

We begin by describing the thermodynamic framework through which we will
explore and answer the question posed above. In the Grüneisen theory of
thermal expansion Wallace (1972); Horton and Maradudin (1974), the thermal
expansion can be related to a thermal stress proportional to the mean (or
bulk) Grüneisen parameters. The generalized mode Grüneisen parameter Wallace
(1972) for phonon mode $s$ at wave vector $\mathbf{q}$ is defined as the
change in phonon frequency $\omega$ with respect to an infinitesimal strain
$\epsilon$ along Cartesian directions $i$ and $j$:

$\gamma^{ij}_{s,\mathbf{q}}\equiv-\frac{1}{\omega_{s,\mathbf{q}}}\frac{\partial\omega_{s,\mathbf{q}}}{\partial\varepsilon_{ij}}.$
(4)

The mean Grüneisen parameter Choy  _et al._ (1984); Ashcroft and Mermin (1976)
is then essentially the sum of the individual mode Grüneisen parameters
weighted by their specific heats, $c_{s,\mathbf{q}}$:

$\gamma^{ij}\equiv\frac{\sum_{s,\mathbf{q}}\gamma^{ij}_{s,\mathbf{q}}c_{s,\mathbf{q}}}{\sum_{s,\mathbf{q}}c_{s,\mathbf{q}}}.$
(5)

At a given temperature, the mean Grüneisen tensor can be related to the
thermal strain tensor, $\alpha_{j}$, through the elastic stiffness tensor
$C_{ij}$,Wallace (1972); Munn (1972)

$\sum_{j}C_{ij}\,\alpha_{j}=\frac{C_{\eta}}{V}\gamma^{i},$ (6)

where $C_{\eta}$ is the bulk heat capacity at constant configuration, $V$ is
the equilibrium volume and $\gamma^{i}$ is the mean Grüneisen parameter of
Equation 5 in Voigt notation. Both $\alpha_{j}$ and $\gamma^{i}$ are
6$\times$1 vectors and $C_{ij}$ is a 6$\times$6 tensor, again employing Voigt
notation. We emphasize that all of the materials properties in Equation 6 are
temperature dependent, including the elastic constants.

Since thermal strains are forbidden by symmetry to induce shear strains in
tetragonal systems, for this work we only need to consider the upper left
3$\times$3 block of the thermal expansion tensor. Hence, the thermal expansion
tensor can be written as a 3$\times$1 vector where the first two indices
correspond to the strain per degree Kelvin of the in-plane $a$ and $b$ lattice
parameters (necessarily equal in a tetragonal system), and the third index
corresponds to strain per degree Kelvin of the out-of-plane $c$ lattice
parameter:

$\alpha_{j}=\begin{bmatrix}\alpha_{a}\\\ \alpha_{a}\\\
\alpha_{c}\end{bmatrix}.$ (7)

The mean Grüneisen tensor can similarly be represented as a 3$\times$1 vector
where the first two indices correspond to the mean Grüneisen parameter with
respect to $a$ and $b$ (again, necessarily equal by symmetry) and the third
index corresponds to the mean Grüneisen parameter with respect to $c$:

$\gamma^{i}=\begin{bmatrix}\gamma^{a}\\\ \gamma^{a}\\\
\gamma^{c}\end{bmatrix}.$ (8)

These relations help us understand how changes in the mechanical boundary
conditions applied to PbTiO3 affect its vibrational properties, and in this
study allow us to quantify changes in the thermal expansion properties of the
PbTiO3 film as a function of misfit strain and substrate thermal expansion
coefficient.

We compute the materials properties discussed above using the results from the
QHA. For $\alpha_{j}$, we compute a numerical derivative of strain at each
value of $T$ using the predicted lattice parameters and a central finite
difference method. The $C_{ij}$ relevant to the Grüneisen framework of thermal
expansion are calculated using derivatives of the Helmholtz free energy
surface fit to the points in the quasiharmonic grid with respect to strain
along the $a$ and $c$ axes (see Supplementary Information for details). Note
that since the films are under a nonzero state of in-plane strain, the elastic
stiffness tensor is no longer symmetric, i.e. $C_{13}\neq C_{31}$, and
requires both second derivatives of energy with respect to strain as well as
first derivatives (stress) to calculate Barron and Klein (1965); Wallace
(1965); Wang and Li (2012).

After the $\alpha_{j}$ and $C_{ij}$ are calculated, the Grüneisen parameters
can be found by solving the system of linear equations described by Equation
6. In this study, we are often only concerned with analyzing thermal expansion
along $c$, as $\alpha_{a}$ is defined by the substrate. Thus, the only
Grüneisen parameter necessary is $\gamma_{c}$, which can be solved for using
the following relation derived from 6:

$2C_{31}\alpha_{a}+C_{33}\alpha_{c}=\frac{C_{\eta}}{V}\gamma^{c}.$ (9)

###  III.2 Relationship between linear elastic and anharmonic effects in the
Grüneisen theory

Having established some basic definitions, what we would now like is an
expression that relates changes (compared to bulk) in the thermal expansion
properties of PbTiO3 thin-films to corresponding changes in the elastic
constants and vibrational properties (Grüneisen parameters). Let
$\Delta\alpha_{i}(T)\equiv\alpha_{i}(T)-\alpha_{i}^{\mathrm{bulk}}(T)$ such
that $\alpha_{i}^{\mathrm{bulk}}(T)+\Delta\alpha_{i}(T)=\alpha_{i}(T)$.
Similarly, let $\Delta C_{ij}(T)\equiv C_{ij}(T)-C_{ij}^{\mathrm{bulk}}(T)$,
where C${}_{ij}(T)$ is the temperature-dependent elastic stiffness tensor of
the PbTiO3 thin-film and $C_{ij}^{\mathrm{bulk}}(T)$ is the elastic stiffness
tensor of bulk PbTiO3. Finally, let
$\Delta\gamma_{i}(T)\equiv\gamma_{i}(T)-\gamma_{i}^{\mathrm{bulk}}(T)$, where
$\gamma_{i}(T)$ is the $i^{th}$ element of the Grüneisen tensor for a PbTiO3
thin-film and $\gamma_{i}^{\mathrm{bulk}}(T)$ is the Grüneisen tensor of bulk
PbTiO3. Thus, $\Delta\alpha_{i}(T)$, $\Delta C_{ij}(T)$, and
$\Delta\gamma_{i}(T)$ respectively denote the difference in thermal expansion,
elastic properties, and Grüneisen parameters in the strained film from that of
bulk PbTiO3 at the same temperature.

In the Grüneisen theory of thermal expansion, the difference in thermal
expansion between a PbTiO3 film and bulk PbTiO3, $\Delta\alpha(T)$, is
entirely accounted for by the departure of the Grünesien parameters and
elastic constants from the values they would have in the bulk system (with
higher-order phonon-phonon coupling terms ignored Wallace (1972); Horton and
Maradudin (1974)). Thus, using Equation 6, for the strained film

$\sum_{j}C_{ij}(T)\,\alpha_{j}(T)=\frac{C_{\eta}(T)}{V(T)}\gamma_{i}(T),$ (10)

and therefore,

$\sum_{j}(C_{ij}^{\mathrm{bulk}}(T)+\Delta
C_{ij}(T))\,(\alpha^{\mathrm{bulk}}_{j}(T)+\Delta\alpha_{j}(T))=\\\
\frac{C_{\eta}(T)}{V(T)}(\gamma^{\mathrm{bulk}}_{i}(T)+\Delta\gamma_{i}(T)).$
(11)

Expanding Equation 11 for $i=3$ for a tetragonal system and subtracting from
it the relation in Equation 6 for the bulk system,

$\displaystyle\begin{split}2C^{\mathrm{bulk}}_{31}\Delta\alpha_{a}+C^{\mathrm{bulk}}_{33}\Delta\alpha_{c}+2\Delta
C_{31}\alpha^{\mathrm{bulk}}_{a}+\Delta C_{33}\alpha^{\mathrm{bulk}}_{c}+\\\
2\Delta C_{31}\Delta\alpha_{a}+\Delta
C_{33}\Delta\alpha_{c}=\frac{C_{\eta}}{V}\Delta\gamma^{c}.\end{split}$ (12)

We remind the reader again that since the film is under non-zero stress along
the $a$-axis, the elasticity tensor is no longer necessarily symmetric, and
thus $C_{31}\neq C_{13}$. For notational simplicity, we have also dropped the
$(T)$ from each of these quantities, though they all remain functions of
temperature. Solving Equation 12 for $\Delta\alpha_{c}$ and rearranging (see
Supplementary Information for a more detailed derivation of Equations 11–13),
we obtain,

$\begin{split}\Delta\alpha_{c}=-2\frac{C^{\mathrm{bulk}}_{31}}{C^{\mathrm{bulk}}_{33}}\Delta\alpha_{a}+\\\
\left[\left(2\frac{C^{\mathrm{bulk}}_{31}}{C^{\mathrm{bulk}}_{33}}\Delta\alpha_{a}-\alpha^{\mathrm{bulk}}_{c}\right)\frac{\Delta
C_{33}}{(C^{\mathrm{bulk}}_{33}+\Delta C_{33})}\right.\\\
\left.-2(\alpha^{\mathrm{bulk}}_{a}+\Delta\alpha_{a})\frac{\Delta
C_{31}}{(C^{\mathrm{bulk}}_{33}+\Delta C_{33})}\right]+\\\
\frac{C_{\eta}\Delta\gamma^{c}}{V(C^{\mathrm{bulk}}_{33}+\Delta
C_{33})}.\end{split}$ (13)

Though it appears complex, Equation 13 establishes the sought after connection
between $\Delta\alpha_{c}$, the _deviation_ of $\alpha_{c}$ of the strained
film from that of bulk PbTiO3, and the elastic and vibrational properties of
strained PbTiO3 films and how _they_ depart from that of bulk PbTiO3. By
understanding the physical meaning of each term in Equation 13, we can explain
the microscopic mechanism underlying $\Delta\alpha_{c}$ in strained film
systems in terms of changes in elasticity and vibrational properties:

  1. 1.

We refer to the first term,
$-2\frac{C^{\mathrm{bulk}}_{31}}{C^{\mathrm{bulk}}_{33}}\Delta\alpha_{a}$, as
_linear elastic_ coupling. This is the elastic coupling between in-plane
($a$-axis) and out-of-plane ($c$-axis) thermal strain from linear elasticity.
It accounts for the change in the rate of thermal expansion along $c$ of the
film that would be predicted if the materials properties $C_{ij}$ and
$\gamma^{c}$ did not change at all with strain.

  2. 2.

We refer to the second and third terms (grouped by square brackets in Equation
13) as _elastic anharmonicity_. These terms capture how changes in the elastic
properties in the strained film affect $\Delta\alpha_{c}$.

  3. 3.

We refer to the last term,
$-\frac{C_{\eta}\Delta\gamma^{c}}{V(C^{\mathrm{bulk}}_{33}+\Delta C_{33})}$,
as _anharmonic thermal stress_. This term captures how changes in phonon
anharmonicity (due to phonon-strain coupling), and the thermal stress it
induces in the strained film, affect $\Delta\alpha_{c}$.

It is clear that the difference in the $c$-axis thermal expansion of strained
PbTiO3 thin-films compared to bulk arises from three sources – the mismatch
between the thermal expansion coefficients of the substrate and bulk PbTiO3
($\Delta\alpha_{a}$), the difference between the elastic constants of the
strained film compared to those of bulk PbTiO3 ($\Delta C_{31}$ and $\Delta
C_{33}$), and the difference between the Grüneisen parameters along the
$c$-axis of the strained film compared to those of bulk PbTiO3
($\Delta\gamma_{c}$). Note that the misfit strain at a given temperature gives
rise to changes in the elastic constants and the Grüneisen parameters
($\Delta\gamma^{c}$, $\Delta C_{31}$ and $\Delta C_{33}$) at that temperature,
whereas $\Delta\alpha_{a}$ determines the misfit strain rate. In the following
section, we consider how each of $\Delta C_{31}$, $\Delta C_{33}$
$\Delta\gamma^{c}$, and $\Delta\alpha_{a}$ affect the thermal expansion of
strained PbTiO3 thin-films along the $c$-axis, and how they depend on misfit
strain and the thermal expansion coefficient of a given substrate.

##  IV Results

###  IV.1 Thermal expansion in strained thin-films – Importance of
anharmonicity

Figure 2 shows the rate of thermal expansion along the $c$-axis for PbTiO3
thin-films as a function of strain and two different substrate thermal
expansions compared with bulk PbTiO3. When $\alpha_{a}$ = 1$\times$10-5, the
negative thermal expansion along $c$ is suppressed compared to bulk for each
of the strained films. When the rate of thermal expansion of the substrate is
increased ($\alpha_{a}$ = 2.5$\times$10-5), the strained films initially show
an enhanced rate of negative thermal expansion along $c$ compared to bulk, up
to about 400 K. Above this temperature, bulk PbTiO3 once again shows the
largest rate of negative thermal expansion along $c$. The behavior of the
lattice-matched films (zero misfit strain at 300 K) is always ‘between’ that
of the tensile and compressively strained films. In order to prevent the
discussion below from becoming overly complex, we do not explicitly discuss
the lattice-matched films.str

Figure 2: Comparison between $c$-axis thermal expansion coefficients
($\alpha_{c}$) as a function of temperature for bulk PbTiO3 and a series of
strained films from our QHA calculations for growth on two different
substrates with $\alpha_{a}$ = 1.0$\times$10-5 (left) and 2.5$\times$10-5
(right) – compressive ($\varepsilon_{a}^{300}$ = -0.75%, tensile
($\varepsilon_{a}^{300}$ = +0.75%) and lattice-matched
($\varepsilon_{a}^{300}$ = 0.00%). The $\alpha_{c}$ data are denoted by closed
squares, whereas the $\alpha_{a}$ data (shown for bulk PbTiO3 only) are
denoted by open black circles. The data for bulk PbTiO3 are identical in both
panels.

We can understand the trends shown in Figure 2 by considering individually the
contributions of elastic anharmonicity, thermal stress and linear elasticity
in Equation 13 to $\Delta\alpha_{c}$, the deviation in thermal expansion along
$c$ in strained PbTiO3 compared to bulk. Figures 3 and 4 show
$\Delta\alpha_{c}$ and its various contributions for PbTiO3 thin-films under
compressive and tensile strain for two different substrate thermal expansion
coefficients (the same as those considered in Figure 2). Both figures appear
to show that linear elasticity (the first term of Equation 13) makes the
dominant contribution to $\Delta\alpha_{c}$, with the red squares almost
overlaid on the black squares. In this scenario, the thermal expansion of the
$c$-axis as a function of temperature is controlled by Poisson-like linear
elastic coupling between the $a$ and $c$-axes. If we were to consider _only_
the linear elastic contribution to $\Delta\alpha_{c}$, the PbTiO3 films would
appear highly harmonic.

Figure 3: Contributions to $\Delta\alpha_{c}$, broken down according to
Equation 13, as a function of temperature for compressively strained films.
The misfit strain at 300 K ($\varepsilon_{a}^{300}$) is -0.75% while the
strain rate is varied from that of bulk PbTiO3. The linear elasticity
contribution refers to the first term of Equation 13, the anharmonic
elasticity contribution refers to the second and third terms of Equation 13
(grouped in parentheses), and the anharmonic thermal stress contribution
refers to the fourth term of Equation 13.

.

Figure 4: Contributions to $\Delta\alpha_{c}$, broken down according to
Equation 13, as a function of temperature for tensile strained films. The
misfit strain at 300 K ($\varepsilon_{a}^{300}$) is +0.75%, while the strain
rate is varied from that of bulk PbTiO3. The linear elasticity contribution
refers to the first term of Equation 13, the anharmonic elasticity
contribution refers to the second and third terms of Equation 13 (grouped in
parentheses), and the anharmonic thermal stress contribution refers to the
fourth term of Equation 13.

.

However, Figures 3 and 4 additionally show that there are also large
contributions to $\Delta\alpha_{c}$ from elastic anharmonicity and anharmonic
thermal stress – in the case of the PbTiO3 thin-film under compressive strain
and $\alpha_{a}$ = 2.5$\times$10-5, these terms are actually larger in
magnitude than the linear elastic contribution through nearly the entire
temperature range. Critically, the contributions to $\Delta\alpha_{c}$ from
elastic anharmonicity and anharmonic thermal stress are of _opposite sign_ ,
leading to a fortuitous near-cancellation (whereas the sign of the
contribution from elastic anharmonicity is always opposite that of the
anharmonic thermal stress, which one is positive and which negative depends on
the sign of the misfit strain at a given temperature; for elastic
anharmonicity the contribution tends to be the same sign as the misfit
strain). This fortuitous near-cancellation between different types of
anharmonicity, the origin of which will be discussed in the next section,
results in what looks like a linear elastic thermal response of strained
PbTiO3 thin-films.

###  IV.2 Dependence of elastic constants and Grüneisen parameters on
temperature and misfit strain

The near-cancellation of contributions from anharmonicity shown in Figures 3
and 4 can ultimately be attributed to changes in the elastic constants and
Grüneisen parameters with temperature and misfit strain. Figure 5 shows that
$C_{31}$, $C_{33}$ and $\gamma_{c}$ for the strained films can all differ
significantly from bulk PbTiO3, and also change significantly with temperature
depending on the misft strain and substrate thermal expansion coefficient.
Despite the fact that contributions to $\Delta\alpha_{c}$ from just the linear
elasticity term of Equation 13 closely match those calculated using the full
expression for $\Delta\alpha_{c}$, the large strain dependencies depicted in
Figure 5 corroborate the data shown in Figures 3 and 4 in illustrating that
strained PbTiO3 thin-films are highly anharmonic, especially with respect to
the elastic constants (note that even though the values of $\gamma_{c}$ tend
to be fairly small, they still experience large _changes_ with strain,
indicating the presence of anharmonicity). The details of how the elastic
constants and Grüneisen parameters change with temperature affect how elastic
anharmonicity and anharmonic thermal stress contribute to changes in the
$c$-axis of strained PbTiO3 films with temperature.

Figure 5: Comparison between strained films on substrates with two different
thermal expansion coefficients ($\alpha_{a}=1.0\times 10^{-5}$ and
$\alpha_{a}=2.5\times 10^{-5}$) and bulk PbTiO3 for $C_{33}$ (left), $C_{31}$
(middle) and $\gamma_{c}$ (right) from our first-principles calculations.

Considering first the compressively strained film, $C_{33}$ is softer than in
bulk PbTiO3 ($\Delta C_{33}<0$), whereas $C_{31}$ is stiffer ($\Delta
C_{31}>0$). The effect of these strain-induced changes to elasticity, which is
captured by the second term (grouped in square brackets) in Equation 13 and
the green squares in Figure 3, is to make a negative contribution to
$\Delta\alpha_{c}$. Films under tensile strain exhibit the opposite behavior
of the elastic constants compared to bulk – under tensile strain, elastic
anharmonicity makes a positive contribution to $\Delta\alpha_{c}$, as shown by
the green squares in Figure 4. The Supplementary Information contains a short
discussion of how the observed effect of changes in the elastic constants on
$\Delta\alpha_{c}$ can be predicted from a careful analysis of the sign of
each of the two terms contributing to anharmonic elasticity in Equation 13.

With respect to the Grüneisen parameters, Figure 5 shows that $\gamma_{c}$ is
larger in the compressively strained film than in bulk PbTiO3
($\Delta\gamma_{c}>0$), whereas it is smaller in the tensile strained film
($\Delta\gamma_{c}<0$). Referring back to Equation 13, compressive strain then
causes the anharmonic thermal stress (last) term of Equation 13 to make a
positive contribution to $\Delta\alpha_{c}$, as shown by the blue squares in
Figure 3. In contrast, the anharmonic thermal stress makes a negative
contribution to $\Delta\alpha_{c}$ in films under tensile strain, as shown by
the blue squares in Figure 4 (note that while $C_{33}$ appears in the
denominator of this term, it cannot change the sign of the term). In other
words, the contribution from anharmonic thermal stress to thermal expansion
along $c$ will be the same sign as $\Delta\gamma_{c}$, which is typically
opposite in sign to the misfit strain. The changes in contributions from the
anharmonic thermal stress with respect to misfit strain are opposite those of
anharmonic elasticity, as discussed above.

Hence, when PbTiO3 undergoes biaxial strain, there is a competition between
elastic anharmonicity (the second and third terms of Equation 13), and
anharmonic thermal stress (the last term of Equation 13), and a large portion
of their sum cancels. This near-cancellation is the result of large,
simultaneous changes in elastic and vibrational properties in the film, and
causes the response to appear nearly elastic, with a small departure from
linear elasticity due to the cancellation not being exact. The contribution to
$\alpha_{c}$ from anharmonic elasticity is proportional to $\Delta C_{33}$,
which is typically the same sign as the misfit strain, and opposes the
contribution from anharmonic thermal stress proportional to
$\Delta\gamma_{c}$, which is typically the opposite sign of the misfit strain.

###  IV.3 Understanding differences among strained films in terms of
contributions to $\Delta\alpha_{c}$

Finally, returning now to Figure 2, although all of the strained films we
studied exhibit suppressed negative thermal expansion along the $c$-axis
compared to bulk (except at low temperatures for the $\alpha_{a}$ =
2.5$\times$10-5 substrate), there are differences in thermal expansion among
the strained films, as noted earlier. If the behavior of the PbTiO3 thin-films
really was fully controlled by linear elasticity, the thermal expansion of the
$c$-axis would be completely determined by the substrate thermal expansion
coefficient $\alpha_{a}$ and the elastic properties of bulk PbTiO3. Hence, the
$c$-axis thermal expansion would be the same, regardless of the misfit strain.
Instead, the details of the thermal expansion of the films do depend on both
the misfit strain and the mismatch between the thermal expansion coefficients
of the substrate and bulk PbTiO3 ($\Delta\alpha_{a}$, the misfit strain
_rate_).

For example, compressively strained films on substrates with $\alpha_{a}$ =
1$\times$10-5 show the _least_ negative rate of thermal expansion along the
$c$-axis but when $\alpha_{a}$ increases to 2.5$\times$10-5, they show the
_most_ negative thermal expansion along $c$ among the thin-film systems. This
is because for compressive strain and $\alpha_{a}$ = 1$\times$10-5, the
positive contributions to $\Delta\alpha_{c}$ from anharmonic thermal stress
are larger in magnitude than the negative contributions from anharmonic
elasticity, pushing the thermal expansion of the $c$-axis to be slightly more
positive than for the other strained systems, as shown in the left panel of
Figure 3. When $\alpha_{a}$ = 2.5$\times$10-5, the negative contributions to
$\Delta\alpha_{c}$ from elastic anharmonicity are larger in magnitude than the
positive contributions from anharmonic thermal stress, pushing the thermal
expansion of $c$ to be slightly more negative than for the other strained
systems, as shown in the right panel of Figure 3. These trends are reversed
for films under tensile strain.

Figures 3 and 4 show that, for a given strain (compressive or tensile), the
anharmonic thermal stress is relatively insensitive to $\alpha_{a}$. In
addition, Equation 13 shows that neither $\alpha_{a}$ nor $\Delta\alpha_{a}$
appear in the (fourth) anharmonic stress term. This directs our attention to
the two terms associated with anharmonic elasticity. Figure 6 shows how the
average values of the two contributions to $\Delta\alpha_{c}$ from anharmonic
elasticity terms in Equation 13 change as a function of $\alpha_{a}$ for both
compressive and tensile strain. For films under compressive strain, it is the
first term,
$\left(2\frac{C^{\mathrm{bulk}}_{31}}{C^{\mathrm{bulk}}_{33}}\Delta\alpha_{a}-\alpha^{\mathrm{bulk}}_{c}\right)\Delta
C_{33}/(C^{\mathrm{bulk}}_{33}+\Delta C_{33})$, that has the most negative
contribution to $\Delta\alpha_{c}$ and the most negative slope with respect to
$\alpha_{a}$. This means that for compressively strained films, it is this
term that makes the contribution of elastic anharmonicity to
$\Delta\alpha_{c}$ negative and larger in magnitude than the positive
contribution from anharmonic thermal stress. Although the compressively
strained film initially shows the least negative rate of thermal expansion
along the $c$-axis among the strained films, as $\alpha_{a}$ increases the
magnitude of the negative contribution from elastic anharmonicity also
increases, leading the compressively strained film to exhibit the most
negative rate of thermal expansion along $c$ when $\alpha_{a}$ =
2.5$\times$10-5. For films under tensile strain, both anharmonic elastic terms
are positive and increase as $\alpha_{a}$ increases – the first term again
exhibits the largest contribution and slope with respect to $\alpha_{a}$, yet
this time both are positive. Hence, in this case the positive contributions to
$\Delta\alpha_{c}$ from anharmonic elasticity grow larger in magnitude than
the negative contributions from anharmonic thermal stress as $\alpha_{a}$
increases, leading to films under tensile strain having a slower rate of
thermal expansion along the $c$-axis than either bulk PbTiO3 or compressively
strained films. As discussed in Section IV.2 and the Supplementary
Information, the sign of this critical first term is determined by the sign of
$\Delta C_{33}$, which is the same as the sign of $\varepsilon_{a}$. The
behavior of this term, combined with the relative insensitivity of anharmonic
thermal stress on $\alpha_{a}$, explains the departure of $\Delta\alpha_{c}$
from what would be predicted by linear elasticity alone.

Figure 6: Variation in contribution to $\Delta\alpha_{c}$ of anharmonic
elasticty term of Equation 13 in terms of the second term of Equation 13
(blue) and the third term of Equation 13 (red). Data are shown as a function
of substrate thermal expansion coefficient ($\alpha_{a}$) for compressively
strained ($\epsilon_{a}^{300}$ = -0.75%) and tensile strained
($\epsilon_{a}^{300}$ = +0.75%) films. Each data point corresponds to the
average value of the indicated term from 0 K to 800 K. The lines are guides
for the eye.

###  IV.4 Explaining thermal expansion in experimentally synthesized PbTiO3
thin-films

We can use our findings above to explain the experimentally observed thermal
expansion behavior of PbTiO3 thin-films on LSAT,Chakoumakos  _et al._ (1998)
SrTiO3de Ligny and Richet (1996) and DyScO3Biegalski  _et al._ (2005)
substrates. Figure 7 shows the lattice parameters and thermal expansion
coefficients as functions of temperature for PbTiO3 thin-films on all three
substrates. All three substrates have similar thermal expansion coefficients –
from 0 – 800 K, $\alpha_{a}$ is fairly low and constant with temperature with
average values of $0.93\times 10^{-5}$ K-1, $1.10\times 10^{-5}$ K-1, and
$1.04\times 10^{-5}$ K-1, for LSAT, SrTiO3, and DyScO3, respectively. In each
case the thermal expansion along $c$ is less negative than in bulk PbTiO3,
however the $\alpha_{c}$ for each substrate is different, and the framework
established in this study explains why.

Figure 7: Comparison between lattice parameters (left) and thermal expansion
coefficients (right) of bulk PbTiO3 and strained PbTiO3 thin-films on LSAT,
DyScO3 and SrTiO3 substrates. Substantially similar data are reported in Ref.
Ritz and Benedek, 2020, however in this work we used a slightly more dense
grid for our QHA calculations. Open cirlces denote data for the $a$-axis
lattice parameters and $\alpha_{a}$ (experimental data for SrTiO3, DyScO3, and
LSATde Ligny and Richet (1996); Biegalski  _et al._ (2005); Chakoumakos  _et
al._ (1998)), whereas closed squares denote data for the $c$-axis lattice
parameters and $\alpha_{c}$ (predicted from our QHA calculations).

At 300 K, LSAT induces a large compressive misfit strain of -1.01%, SrTiO3 a
small strain of -0.01%, and DyScO3 induces a large tensile strain of +1.01%.
We showed above that under compressive strain and when $\alpha_{a}$ is small,
the negative magnitude of the contribution to $\Delta\alpha_{c}$ from elastic
anharmonicity is smaller than the positive contribution from anharmonic
thermal stress, leading to a smaller rate of negative thermal expansion along
$c$ compared to both bulk and tensile strained films. By this reasoning,
PbTiO3 on LSAT should show the least negative thermal expansion along $c$ and
PbTiO3 on DyScO3 should show the most negative thermal expansion along $c$;
this is indeed what is shown in Figure 7. Put another way, we showed earlier
that the sign of the contribution from anharmonic thermal stress to
$\Delta\alpha_{c}$ has the opposite sign of the misfit strain. Therefore, we
would expect $\Delta\alpha_{c}$ to be most positive in LSAT and most negative
in DyScO3, and this is indeed what we see.

##  V Summary and Conclusions

Our results reveal that although PbTiO3 thin-films are highly anharmonic, a
fortuitous near-cancellation between the thermal stress and anharmonic elastic
contributions to the temperature-dependent structural parameters means that
the thermal response of the $c$-axis of the films can be predicted very well
using only the misfit strain and bulk PbTiO3 properties. We note that the
$\Delta\gamma^{c}$ term, which contributes to the anharmonic thermal stress,
involves the _derivative_ of the Grüneisen parameter with respect to strain.
Since the Grüneisen parameter is a function of the third-order interatomic
force constants (and higher order, if sublattice displacements are present, as
is the case for ferroelectric PbTiO3 Horton and Maradudin (1974); Wallace
(1972)), $\Delta\gamma^{c}$ is a function of the _fourth_ and higher-order
force constants. Additionally, as elastic constants are related to the
curvature of the Helmholtz free energy surface, the large changes in the
elastic constants with strain are also driven by changes in both electronic
and vibrational free energy, and thus the third, fourth, and higher-order
force constants Dangić  _et al._ (2018). The low thermal conductivity observed
for bulk PbTiO3Tachibana  _et al._ (2008); Yoshida (1960) corroborates these
findings. In this light however, the approximately linear-elastic behavior of
the strained PbTiO3 thin-films is unexpected.

Our findings may also help explain why popular phenomenological models for
predicting the strain behavior of the polarization in ferroelectric PbTiO3
thin-films, such as those described in Refs. 38 and 39, are reasonably
effective at capturing the thermal expansion behavior of both bulk and
epitaxially strained films of PbTiO3 with temperature Highland  _et al._
(2014); Lichtensteiger  _et al._ (2005); Venkatesan  _et al._ (2008); Haun
_et al._ (1987), despite being functions of the elastic constants of the high-
symmetry cubic phase, and with no explicit elastic anharmonicity included
beyond its coupling to the polarization. If the structural behavior of
ferroelectric PbTiO3 films clamped to a substrate can be well-described by
linear elasticity with much of the higher-order anharmonicity ignored, it
stands to reason that the inclusion of higher-order elastic coupling terms
beyond the few judiciously chosen in the above models would be unnecessary to
accurately capture the structural evolution of the system with temperature.

Finally, our work shows that the thermal expansion behavior of PbTiO3 films
depends on the misfit strain throughout the _entire_ temperature range. This
means it depends on not only the misfit strain between film and substrate
lattice parameters at a single temperature, but on the misfit between their
thermal strain _rates_ , which determines how misfit strain accumulates with
temperature. These dependencies are entwined in a complex way – the misfit
strain between the lattice parameters of the substrate and the equilibrium
(bulk) lattice parameters of the film material can induce large changes in the
elastic properties and Grüniesen parameters of the film, which themselves
influence how the $c$-axis of the film responds to the misfit strain rate. Our
combined first-principles and phenomenological approach provides significant
insights into this behavior and offers a systematic framework for exploring
thermal expansion phenomena in both bulk and thin-film systems.

###### Acknowledgements.

This work was supported by the National Science Foundation. E. T. R. and N. A.
B. were supported by DMR-1550347. Computational resources were provided by the
Cornell Center for Advanced Computing and the Extreme Science and Engineering
Discovery Environment (XSEDE) through allocation DMR-160052.

## References

  * Tucker  _et al._ (2005) M. G. Tucker, A. L. Goodwin, M. T. Dove, D. A. Keen, S. A. Wells, and J. S. Evans, Phys. Rev. Lett. 95, 255501 (2005).
  * Mounet and Marzari (2005) N. Mounet and N. Marzari, Phys. Rev. B 71, 205214 (2005).
  * Kim  _et al._ (2018) D. S. Kim, O. Hellman, J. Herriman, H. Smith, J. Lin, N. Shulumba, J. Niedziela, C. Li, D. Abernathy, and B. Fultz, Proc. Natl. Acad. Sci. U.S.A. 115, 1992 (2018).
  * Ritz and Benedek (2018) E. T. Ritz and N. A. Benedek, Phys. Rev. Lett. 121, 255901 (2018).
  * Goodwin  _et al._ (2008a) A. L. Goodwin, D. A. Keen, and M. G. Tucker, Proc. Natl. Acad. Sci. U.S.A. 105, 18708 (2008a).
  * Goodwin  _et al._ (2008b) A. L. Goodwin, D. A. Keen, M. G. Tucker, M. T. Dove, L. Peters, and J. S. O. Evans, J. Am. Chem. Soc. 130, 9660 (2008b).
  * Goodwin  _et al._ (2009) A. L. Goodwin, B. J. Kennedy, and C. J. Kepert, J. Am. Chem. Soc. 131, 6334 (2009).
  * Lind (2012) C. Lind, Materials 5, 1125 (2012).
  * Miller  _et al._ (2009) W. Miller, C. Smith, D. Mackenzie, and K. Evans, J. Mater. Sci. 44, 5441 (2009).
  * Barrera  _et al._ (2005) G. D. Barrera, J. A. O. Bruno, T. H. K. Barron, and N. L. Allan, J. Phys. Condens. Matter 17, R217 (2005).
  * Yanase  _et al._ (2011) I. Yanase, T. Kojima, and H. Kobayashi, Solid State Commun. 151, 595 (2011).
  * Ablitt  _et al._ (2019) C. Ablitt, H. McCay, S. Craddock, L. Cooper, E. Reynolds, A. A. Mostofi, N. C. Bristowe, C. A. Murray, and M. S. Senn, Chem. Mater. (2019).
  * Chen  _et al._ (2005) J. Chen, X. Xing, R. Yu, and G. Liu, J. Am. Ceram. Soc. 88, 1356 (2005).
  * Shirane and Hoshino (1951) G. Shirane and S. Hoshino, J. Phys. Soc. Jpn. 6, 265 (1951).
  * Janolin  _et al._ (2007) P.-E. Janolin, F. Le Marrec, J. Chevreul, and B. Dkhil, Appl. Phys. Lett. 90, 192910 (2007).
  * Janolin (2009) P.-E. Janolin, Journal of Materials Science 44, 5025 (2009).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, _et al._ , J. Phys. Condens. Matter 21, 395502 (2009).
  * Wu and Cohen (2006) Z. Wu and R. E. Cohen, Phys. Rev. B 73, 235116 (2006).
  * Garrity  _et al._ (2014) K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Comput. Mater. Sci. 81, 446 (2014).
  * Ritz  _et al._ (2019) E. T. Ritz, S. J. Li, and N. A. Benedek, J. Appl. Phys. 126, 171102 (2019).
  * Ritz and Benedek (2020) E. T. Ritz and N. A. Benedek, Physical Review Materials 4, 084410 (2020).
  * (22) Our DFT simulations predict that the unit cell of bulk PbTiO3 exhibits unstable phonon modes at 0 K in its cubic phase, as well as for nearly-cubic unit cells with very low $c/a$ ratios. This results in difficulty when calculating QHA grid points very close to the tetragonal-cubic phase transition. The lowest value of $c$ in the set of unit cells used in our quasiharmonic grid is approximately 4.00 Å, meaning that for temperatures greater than 720 K, the point that minimizes the Helmholtz free energy surface fit is no longer bounded by that grid along the $c$ axis. While such an extrapolation of the derivatives of the vibrational free energy is within the spirit of the Grun̈eisen theory of thermal expansion, our simulation can be interpreted to conclude that the tetragonal-cubic phase transition occurs somewhere within the range of 720 K to 800 K.
  * Wallace (1972) D. C. Wallace,  _Thermodynamics of Crystals_ (Wiley, 1972).
  * Horton and Maradudin (1974) G. K. Horton and A. A. Maradudin,  _Dynamical Properties of Solids_ , Vol. 1 (Elsevier, 1974).
  * Choy  _et al._ (1984) C. Choy, S. Wong, and K. Young, Phys. Rev. B 29, 1741 (1984).
  * Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Holt, Rinehart and Winston, New York, 1976).
  * Munn (1972) R. Munn, J. Phys. C: Solid State Physics 5, 535 (1972).
  * Barron and Klein (1965) T. Barron and M. Klein, Proceedings of the Physical Society (1958-1967) 85, 523 (1965).
  * Wallace (1965) D. C. Wallace, Reviews of Modern Physics 37, 57 (1965).
  * Wang and Li (2012) H. Wang and M. Li, Physical Review B 85, 104103 (2012).
  * (31) The films under compressive and tensile strain at 300 K remain under compressive and tensile strain throughout nearly the entire temperature range studied (the only exception is the film constrained by the substrate with $\varepsilon_{a}^{300}=+0.75\%$ at $\alpha_{a}$ = 1$\times$10-5 for temperatures over 700 K).
  * Chakoumakos  _et al._ (1998) B. Chakoumakos, D. Schlom, M. Urbanik, and J. Luine, J. Appl. Phys. 83, 1979 (1998).
  * de Ligny and Richet (1996) D. de Ligny and P. Richet, Phys. Rev. B 53, 3013 (1996).
  * Biegalski  _et al._ (2005) M. Biegalski, J. Haeni, S. Trolier-McKinstry, D. Schlom, C. Brandle, and A. V. Graitis, J. Mater. Res. 20, 952 (2005).
  * Dangić  _et al._ (2018) D. Dangić, A. R. Murphy, É. D. Murray, S. Fahy, and I. Savić, Physical Review B 97, 224106 (2018).
  * Tachibana  _et al._ (2008) M. Tachibana, T. Kolodiazhnyi, and E. Takayama-Muromachi, Applied Physics Letters 93, 092902 (2008).
  * Yoshida (1960) I. Yoshida, Journal of the Physical Society of Japan 15, 2211 (1960).
  * Pertsev  _et al._ (1998) N. Pertsev, A. Zembilgotov, and A. Tagantsev, Phys. Rev. Lett. 80, 1988 (1998).
  * Zembilgotov  _et al._ (2002) A. G. Zembilgotov, N. A. Pertsev, H. Kohlstedt, and R. Waser, Journal of Applied Physics 91, 2247 (2002), cond-mat/0111218 .
  * Highland  _et al._ (2014) M. Highland, D. Fong, G. Stephenson, T. Fister, P. Fuoss, S. Streiffer, C. Thompson, M.-I. Richard, and J. Eastman, Appl. Phys. Lett. 104, 132901 (2014).
  * Lichtensteiger  _et al._ (2005) C. Lichtensteiger, J.-M. Triscone, J. Junquera, and P. Ghosez, Physical review letters 94, 047603 (2005).
  * Venkatesan  _et al._ (2008) S. Venkatesan, A. Vlooswijk, B. J. Kooi, A. Morelli, G. Palasantzas, J. T. De Hosson, and B. Noheda, Phys. Rev. B 78, 104112 (2008).
  * Haun  _et al._ (1987) M. J. Haun, E. Furman, S. Jang, H. McKinstry, and L. Cross, J. Appl. Phys. 62, 3331 (1987).