2022 171 08 August

Measurements of groomed charm-jet substructure

ALICE Collaboration ALICE Collaboration

Understanding the role of parton mass and Casimir colour factors in the
quantum chromodynamics parton shower represents an important step in
characterising the emission properties of heavy quarks. Recent experimental
advances in jet substructure techniques have provided the opportunity to
isolate and characterise gluon emissions from heavy quarks. In this work, the
first direct experimental constraint on the charm-quark splitting function is
presented, obtained via the measurement of the groomed shared momentum
fraction of the first splitting in charm jets, tagged by a reconstructed
$\textup{D}^{0}$ meson. The measurement is made in proton–proton collisions at
$\sqrt{s}=$13\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$, in the low jet
transverse-momentum interval of $15\leq p_{\textnormal{T}}^{\textnormal{jet
ch}}<$30\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$
where the emission properties are sensitive to parton mass effects. In
addition, the opening angle of the first perturbative emission of the charm
quark, as well as the number of perturbative emissions it undergoes, are
reported. Comparisons to measurements of an inclusive-jet sample show a
steeper splitting function for charm quarks compared to gluons and light
quarks. Charm quarks also undergo fewer perturbative emissions in the parton
shower, with a reduced probability of large-angle emissions.

In hadronic collisions, quarks and gluons (partons) with high transverse
momentum ($p_{\textnormal{T}}$) and/or large masses are produced in scattering
processes involving large momentum transfers. These energetic partons lose
their large initial virtuality by emitting partons in a cascade process known
as a parton shower, until they eventually form hadrons. Theoretical
descriptions of parton showers rely on determining the emission probability of
partons, which are described by splitting functions [1, 2]. Splitting
functions are universal properties of quantum chromodynamics (QCD) that
describe how the energy of a parton is shared as it fragments into further
partons. The splitting functions of light-flavour (up, down and strange)
quarks and gluons differ due to their different Casimir colour factors [2].
For heavy-flavour (charm and beauty with masses of $1.27\pm
0.02\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}{\mathrm{\text{$c$}}}^{2}}{\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}}$$
and $4.18^{+0.03}_{-0.02}$
$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}{\mathrm{\text{$c$}}}^{2}}{\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}}$,
respectively [3]) quarks the splitting functions are strongly influenced by
the large quark mass [4]. As a consequence, the radiation pattern of the
parton shower depends on the type of initiating parton.

Experimental access to the parton shower can be gained via measurements
involving jets, which are collimated bunches of hadrons resulting from the
fragmentation of partons scattered in the initial stages of collisions. Jets
represent a powerful tool for testing perturbative QCD, as their inner
structure (substructure) reflects the properties of the underlying parton
shower. Measurements of jet substructure have been performed in inclusive jets
(with no flavour tagging) to constrain the splitting functions of light quarks
and gluons [5, 6, 7, 8]. However, an experimentally clean separation of gluon-
initiated and quark-initiated jets has remained challenging, with measurements
constraining the splitting functions comprising an admixture of the two
subsets.

Recent experimental advances in techniques pertaining to the substructure of
heavy-flavour jets [9, 10] have allowed for the isolation and study of gluon
emissions from heavy quarks. This technique has been used to perform a first
direct measurement [11] of the dead-cone effect [12] in charm-quark jets. The
dead cone corresponds to an angular region along the flight direction of the
emitter, with an opening proportional to the ratio of mass to energy of the
emitter, within which the probability of emissions is suppressed. The
significant size of this region for low-energy heavy-flavour quarks gives rise
to mass-dependent effects in the parton shower. In addition to the sensitivity
to these mass effects in QCD, heavy-flavour jets represent an experimentally
enhanced quark-initiated jet sample, which can be further used to constrain
the flavour-dependent properties of QCD emissions arising from different
Casimir colour factors.

In this work, the first measurement directly constraining the charm-quark
splitting function is reported. This requires access to perturbative
splittings (the splitting scale is much larger than the QCD scale) of the
charm quark, where the shared momentum fraction of the splittings maps onto
the analytical splitting function [13, 14, 15]. Additionally, measurements of
the opening angle of these emissions, as well as the number of perturbative
emissions in charm jets, are reported. Measurements of these two observables
have also been reported for inclusive jets [16, 6, 7, 8, 17], with the opening
angle calculated at next-to-leading logarithmic accuracy [18]. The
measurements reported in this work are made in the low jet transverse-momentum
interval of $15\leq p_{\textnormal{T}}^{\textnormal{jet
ch}}<$30\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$,
where mass effects are expected to play a role in the fragmentation of the
charm quark and the inclusive sample is mainly comprised of gluon-initiated
jets.

Once jets have been identified from the sample of charged particles
(reconstructed as tracks) in pp collisions, using the
anti-$k_{\textnormal{T}}$ algorithm [19], track-based jet-substructure
observables can be constructed from the jet constituents. To reconstruct the
chain of emissions (splitting tree) inside a jet, the jet constituents can be
reclustered with the Cambridge–Aachen (C/A) algorithm [20], which is well
suited to jet substructure studies because it follows the angular ordering of
emissions in QCD [1] by clustering the jet constituents based on their angular
distance. By following a branch along the splitting tree returned by the
reclustering procedure, a sequence of emissions can be studied [21].

Charm jets are identified by the presence of a reconstructed $\textup{D}^{0}$
meson amongst their constituents. The full reconstruction of the
$\textup{D}^{0}\to\textup{K}^{-}\uppi^{+}$ decay allows for the replacement of
the $\textup{D}^{0}$ decay products with the four-momentum of the
$\textup{D}^{0}$ meson, prior to jet clustering. This ensures that the full
$\textup{D}^{0}$-meson momentum is contained within the jet cone. As the charm
flavour is conserved throughout the jet evolution, the charm quark can be
traced as it dynamically evolves in the showering process, by following the
branch containing the $\textup{D}^{0}$ meson through the splitting tree [9].

To reduce the contribution of non-perturbative effects and increase
sensitivity to perturbative emissions [22], the Soft Drop grooming procedure
[23] is applied, which only selects particular splittings of interest along
the followed branch. To do this, each splitting along this branch is tested
against the Soft Drop condition, as given by

$z\equiv\frac{\min(p_{\textnormal{T, }{1}},p_{\textnormal{T,
}{2}})}{p_{\textnormal{T, }{1}}+p_{\textnormal{T,
}{2}}}>z_{\textnormal{cut}}\left(\frac{\mathup{\Delta}R_{1,2}}{R}\right)^{\beta},$
(1)

where $p_{\textnormal{T, }{1}}$ and $p_{\textnormal{T, }{2}}$ are the
transverse momenta of the leading and subleading prongs of the splitting,
respectively, and $R$ is the jet resolution parameter. The grooming behaviour
is determined by the parameters $z_{\textnormal{cut}}$ and $\beta$, which
control the interplay between the shared momentum fraction, $z$, and the
aperture angle between the prongs,
$\mathup{\Delta}R_{1,2}\equiv\sqrt{(y_{1}-y_{2})^{2}+(\varphi_{1}-\varphi_{2})^{2}}$,
where $y$ and $\varphi$ are the rapidity and azimuth of prongs 1 and 2,
respectively. In this work, values of $z_{\textnormal{cut}}=0.1$ and $\beta=0$
are chosen, such that the Soft Drop condition is satisfied if the subleading
prong carries at least $10\text{\,}\mathrm{\char 37\relax}$ of the sum of the
transverse momenta of the prongs, enriching the selection with perturbative
splittings.

Two groomed substructure observables are constructed against the first
splitting that satisfies the Soft Drop condition: $z_{\textnormal{g}}=z$ and
$R_{\textnormal{g}}=\mathup{\Delta}R_{1,2}$, as given in Eq. 1. For the given
choices of grooming parameters, the shared momentum fraction,
$z_{\textnormal{g}}$, converges to the QCD splitting function at sufficiently
high jet energies [13]. As the emissions are angular ordered, the groomed
emission angle $R_{\textnormal{g}}$ characterises the widest emission in the
splitting tree which passes the Soft Drop condition and sets the geometrical
scale of the jet. This observable is expected to be sensitive to the dead cone
of the charm quark at small angles, whilst emissions from gluons, which have a
larger Casimir colour factor, are expected to dominate the large-angle region
[24, 25]. The number of emissions of the charm quark satisfying the Soft Drop
condition, $n_{\textnormal{SD}}$, is also measured by evaluating all
splittings along the branch containing the $\textup{D}^{0}$ meson. This
strongly correlates to the number of perturbative emissions of the charm
quark. The data were collected using the ALICE apparatus at the LHC, during
the years 2016, 2017, and 2018. Information about the detector configuration
and performance can be found in Refs. 26, 27. The main detector systems used
for this work were the central-barrel detectors, located at $|\eta|<0.9$
within a solenoidal magnet and used for charged-particle tracking and
identification. These include the Inner Tracking System, the Time Projection
Chamber, and the Time-Of-Flight detector. The V0 detector, consisting of two
scintillator arrays located at pseudorapidities of $2.8<\eta<5.1$ and
$-3.7<\eta<-1.7$, was used for the trigger and event selections. Minimum-bias
events were selected by requiring a signal above a given threshold in both V0
counters. The analysed data sample consisted of about
$1.7\text{\times}{10}^{9}$ pp collisions at a centre-of-mass energy of
$\sqrt{s}=$13\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$, corresponding to an
integrated luminosity of
$\mathcal{L}_{\text{int}}=$29\text{\,}\mathrm{n}\mathrm{b}^{-1}$$.

The $\textup{D}^{0}$ mesons were reconstructed in the $\textup{D}^{0}$
$p_{\textnormal{T}}$ interval of $5\leq
p_{\textnormal{T}}^{\textup{D}^{0}}<$30\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$,
through the $\textup{D}^{0}\to\textup{K}^{-}\uppi^{+}$ (and charge conjugate)
decay channel with a branching ratio of $3.95\pm 0.03\text{\,}\mathrm{\char
37\relax}$ [3]. Background candidates, built from pairs of tracks not
corresponding to a $\textup{D}^{0}$ decay, were suppressed by applying
geometrical selections on the displaced decay-vertex topology, as well as
particle identification on the decay-particle tracks [28]. The
$\textup{D}^{0}$ decay daughters were replaced in the event by the
$\textup{D}^{0}$-meson four-momentum, which was obtained by summing the four-
momenta of these decay daughters. Jet finding was then performed using the
FastJet package [29], with the anti-$k_{\textnormal{T}}$ algorithm, a
resolution parameter of $R=0.4$ and the $E$-scheme recombination. The inputs
to the jet finder consisted of tracks of charged particles with
$p_{\textnormal{T}}\geq$0.15\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$
and the reconstructed $\textup{D}^{0}$-meson candidate. The full set of
quality selections applied to the tracks, including the $\textup{D}^{0}$-meson
candidate decay tracks, is given in Ref. 30. The jets were required to be
fully contained within the central barrel acceptance, by imposing a
pseudorapidity selection of $|\eta|<0.5$ on the jet axis. For every
$\textup{D}^{0}$ meson candidate in an event, jet finding was performed
independently, with only the two tracks forming that particular candidate
replaced with the sum of their four-momenta prior to the jet finding pass.
Once jet finding had been performed, the $\textup{D}^{0}$-tagged jet
constituents were reclustered using the C/A algorithm and the
$z_{\textnormal{g}}$, $R_{\textnormal{g}}$, and $n_{\textnormal{SD}}$
observables were calculated for each jet. The contribution to the three
observable distributions from jets tagged by background $\textup{D}^{0}$
candidates that survived the selection was removed via a sideband subtraction
procedure, as detailed in Refs. 30, 31. This procedure provides the true
$\textup{D}^{0}$-tagged jet observable distributions in intervals of
$p_{\textnormal{T}}^{\textup{D}^{0}}$.

Each of these background-subtracted distributions was scaled by the selection
and reconstruction efficiency of prompt (charm-showering)
$\textup{D}^{0}$-tagged jets, estimated in
$p_{\textnormal{T}}^{\textup{D}^{0}}$ intervals, from Monte Carlo (MC)
simulations using PYTHIA 8.243 [32, 33, 34] and GEANT 3 [35]. The efficiency-
corrected distributions were summed over the full
$p_{\textnormal{T}}^{\textup{D}^{0}}$ range considered for the analysis. In
order to study emissions in the charm-quark shower, the contribution from non-
prompt $\textup{D}^{0}$ mesons, originating from the decay of beauty hadrons,
was estimated using POWHEG [36] simulations with PYTHIA 6.425 [33] showering
and decays implemented through the EvtGen package [37]. It was then subtracted
from the measured distributions. The magnitude of this non-prompt contribution
ranged from $20\text{\,}\mathrm{\char 37\relax}$ to $40\text{\,}\mathrm{\char
37\relax}$ in most of the $z_{\textnormal{g}}$, $R_{\textnormal{g}}$, and
$n_{\textnormal{SD}}$ intervals considered in the analysis. The
$p_{\textnormal{T}}^{\textnormal{jet ch}}$ distributions and each prompt
$\textup{D}^{0}$-tagged jet substructure observable were simultaneously
corrected for detector effects via a two-dimensional unfolding, using the
iterative Bayesian unfolding algorithm with four iterations [38].

The $z_{\textnormal{g}}$, $R_{\textnormal{g}}$, and $n_{\textnormal{SD}}$ jet
observables were also measured for inclusive jets. To allow for direct
comparisons with the heavy-flavour results, the constraint on the $Q^{2}$
value of the parton scattering implied by the selection of $\textup{D}^{0}$
mesons with
$p_{\textnormal{T}}^{\textup{D}^{0}}\geq$5\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$
was mimicked for inclusive jets, by requiring a leading track with
$p_{\textnormal{T}}\geq$5.33\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$.
This corresponds to the $p_{\textnormal{T}}$ of a charged pion with the same
transverse mass as a $\textup{D}^{0}$ meson with
$p_{\textnormal{T}}^{\textup{D}^{0}}=$5\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$.
In addition, during the grooming procedure, the hardest branch of the
splitting tree was followed. In the heavy-flavour sample, following the
hardest branch was found to coincide with following the branch containing the
$\textup{D}^{0}$, in over $99\text{\,}\mathrm{\char 37\relax}$ of cases. The
inclusive-jet distributions were then corrected for detector effects via the
same unfolding procedure as in the $\textup{D}^{0}$-tagged jet case.

Systematic uncertainties related to the reconstruction and identification of
prompt $\textup{D}^{0}$-tagged jets were studied. The stability of the
selections applied to the $\textup{D}^{0}$-meson candidates were estimated by
loosening and tightening the topological and PID criteria used to identify the
$\textup{D}^{0}$ decays. The uncertainty pertaining to the signal extraction
procedure was estimated by changing the fitting parameters on the invariant-
mass fits, as well as the widths of the signal and sideband regions used for
the sideband subtraction procedure. The uncertainty on the $\textup{D}^{0}$
feed-down estimation was obtained by varying the renormalisation and
factorisation scales, as well as the mass of the beauty quark, in the POWHEG
simulations [39]. The uncertainty on the jet energy resolution was estimated
by artificially reducing the tracking efficiency used to generate the
unfolding response matrix by $4\%$, in line with the tracking efficiency
uncertainty of the ALICE detector. Finally, the uncertainty associated with
the unfolding procedure was obtained by varying the number of iterations, the
choice of prior, and the ranges of the response matrix. Additionally,
variations of the fragmentation properties of the MC simulation used to
construct the response matrix, obtained by selecting subsamples of the jets
based on jet-angularity observables [40], were also included.

The dominant contributions in most of the intervals resulted from variations
of the topological selections, signal extraction, and feed-down subtraction.
For the feed-down subtraction uncertainty, the maximum negative and the
maximum positive deviations were taken, respectively. For each of the other
variation categories, the root mean square of deviations from the central
values was calculated and assigned as a systematic uncertainty, which was then
symmetrised around the central value. The uncertainties from all categories
were then added in quadrature to obtain the full systematic uncertainty. The
total systematic uncertainties range from $7\text{\,}\mathrm{\char 37\relax}$
to $77\text{\,}\mathrm{\char 37\relax}$, from $9\text{\,}\mathrm{\char
37\relax}$ to $12\text{\,}\mathrm{\char 37\relax}$, and from
$5\text{\,}\mathrm{\char 37\relax}$ to $59\text{\,}\mathrm{\char 37\relax}$
for the $\textup{D}^{0}$-tagged jet $z_{\textnormal{g}}$,
$R_{\textnormal{g}}$, and $n_{\textnormal{SD}}$ distributions, respectively.

The systematic uncertainties for the inclusive jet distributions were
estimated by performing variations of the tracking efficiency and unfolding,
in the same way as in the $\textup{D}^{0}$-tagged jet case. In addition, the
$p_{\textnormal{T}}$ selection on the leading track of the jet was also varied
in line with the track $p_{\textnormal{T}}$ resolution of the
$\textup{D}^{0}$-decay tracks. The total systematic uncertainties range from
$1.4\text{\,}\mathrm{\char 37\relax}$ to $4.3\text{\,}\mathrm{\char
37\relax}$, from $2.2\text{\,}\mathrm{\char 37\relax}$ to
$5.7\text{\,}\mathrm{\char 37\relax}$, and from $3.3\text{\,}\mathrm{\char
37\relax}$ to $8.1\text{\,}\mathrm{\char 37\relax}$ for the inclusive jet
$z_{\textnormal{g}}$, $R_{\textnormal{g}}$, and $n_{\textnormal{SD}}$
distributions, respectively.

The first measurement mapping on to the charm-quark splitting function is
reported via the $z_{\textnormal{g}}$ distribution in Fig. 1, for charm jets
tagged by a prompt $\textup{D}^{0}$ meson. The $R_{\textnormal{g}}$ and
$n_{\textnormal{SD}}$ distributions are shown in the left and right panels of
Fig. 2, respectively. These distributions are fully corrected for detector
effects and are reported in the jet transverse-momentum interval of $15\leq
p_{\textnormal{T}}^{\textnormal{jet
ch}}<$30\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$.
Each charm-tagged jet measurement is accompanied by a measurement of an
inclusive-jet sample in the same kinematic interval, also fully corrected for
detector effects. Both the charm-tagged jet and inclusive jet distributions
are normalised to the total number of jets in each respective category,
irrespective of whether they had a splittings passing the Soft Drop condition
or not. In this way, the distributions are sensitive to the proportion of jets
that do not have any splittings passing the Soft Drop condition, denoted as
the SD-untagged fraction. These fractions are $22\text{\,}\mathrm{\char
37\relax}$ and $2.8\text{\,}\mathrm{\char 37\relax}$, in the charm-tagged jet
and inclusive jet samples, respectively. Both samples are compared to PYTHIA 8
simulations and also to POWHEG + PYTHIA 6 simulations for the charm-tagged
jets. The ratios of predicted distributions from parton-shower models to the
measured distributions are also shown in the bottom panels of Figs. 1 and 2.

The $z_{\textnormal{g}}$ distributions show that charm-tagged jets have
significantly fewer symmetric splittings (large $z_{\textnormal{g}}$ values),
compared to inclusive jets. This is consistent with theoretical predictions
[10] of the role of mass effects in the QCD splitting function. The
$R_{\textnormal{g}}$ distribution for charm quarks shows a reduction at large-
angles compared to inclusive jets. This can be due to the differences between
quark and gluon fragmentation, with the gluon-dominated inclusive sample
expected to feature larger-angle perturbative emissions. The charm and
inclusive jet distributions are consistent at small $R_{\textnormal{g}}$. This
could be due to an interplay between the dead-cone effect, which suppresses
small-angle emissions from the charm quark, and the more abundant emissions
from quarks compared to gluons at small angles, as gluons lose more of their
initial virtuality through larger-angle emissions. The $n_{\textnormal{SD}}$
distribution shows a significant shift to smaller values for the charm-tagged
jets ($\langle n_{\textnormal{SD}}\rangle=1.34\pm 0.06\text{ (stat.)}\pm
0.01\text{ (syst.)}$) compared to inclusive jets ($\langle
n_{\textnormal{SD}}\rangle=2.31\pm 0.01\text{ (stat.)}\pm 0.01\text{
(syst.)}$). This indicates that, compared to light and massless partons, charm
quarks on average emit fewer gluons with a large enough $p_{\textnormal{T}}$
to pass the Soft Drop condition in the showering process. This is consistent
with the expectation from the presence of a dead cone for charm quarks, which
results in a harder fragmentation of charm quarks compared to light quarks and
gluons.

For all three observables, both the PYTHIA 8 and POWHEG + PYTHIA 6 predictions
for charm-quark jets describe the measurements within uncertainties. For the
more precise inclusive measurements, PYTHIA 8 exhibits some tension,
particularly for the $n_{\textnormal{SD}}$ observable. This may be partially
due to the poorly constrained quark and gluon fractions in the inclusive
sample at low $p_{\textnormal{T}}^{\textnormal{jet ch}}$. In contrast, as the
$\textup{D}^{0}$-tagged jets represent a quark-enriched sample, the comparison
of the shower generators to data are less reliant on this fraction and can
instead be used to more accurately study the mechanisms underpinning the
shower properties.

Figure 1: The $z_{\textnormal{g}}$ distribution of prompt
$\textup{D}^{0}$-tagged jets compared to that of inclusive jets for $15\leq
p_{\textnormal{T}}^{\textnormal{jet
ch}}<$30\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$
in pp collisions at $\sqrt{s}=$13\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$,
normalised to the total number of jets . Model/data ratios are shown in the
bottom panels for PYTHIA 8 [32, 33, 34] and POWHEG [36] \+ PYTHIA 6 [33]
simulations.

Figure 2: The $R_{\textnormal{g}}$ (left) and $n_{\textnormal{SD}}$ (right)
distributions of prompt $\textup{D}^{0}$-tagged jets compared to those of
inclusive jets for $15\leq p_{\textnormal{T}}^{\textnormal{jet
ch}}<$30\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$
in pp collisions at $\sqrt{s}=$13\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$.
Model/data ratios are shown in the bottom panels for PYTHIA 8 [32, 33, 34] and
POWHEG [36] \+ PYTHIA 6 [33] simulations.

The measurement of the groomed shared momentum fraction, $z_{\textnormal{g}}$,
of charm-tagged jets, reported in this letter, represents the first direct
experimental constraint of the splitting function of heavy-flavour quarks. The
$z_{\textnormal{g}}$ distribution appears steeper than that of light quarks
and gluons, with a suppressed probability of symmetric emissions. Measurements
of the number of splittings passing the Soft Drop condition as well as the
emission angle of the first of these splittings, show that heavy-flavour
quarks on average have fewer perturbative emissions compared to light quarks
and gluons, with a lower probability of these emissions occurring at large
angles. These different characteristics for heavy-quark emissions, compared to
light quarks and gluons, constrain the roles of quark mass and Casimir colour
factors in the parton shower, which are different between the two samples. The
upgraded ALICE detector in the LHC Run 3 will extend this measurement to jets
tagged with a fully reconstructed beauty meson [41]. Direct comparison of the
two quark-enriched samples of charm-tagged and beauty-tagged jets will enable
the isolation of mass effects from the effects due to Casimir colour factors.
The larger projected integrated luminosity in Run 3 will also allow us to
extend the charm-tagged to inclusive jet comparisons to higher
$p_{\textnormal{T}}^{\textnormal{jet ch}}$, where mass effects become
negligible and the remaining impact of the Casimir colour factors can be
cleanly studied.

## References

  * [1] Y. L. Dokshitzer, V. A. Khoze, A. H. Mueller, and S. I. Troian, Basics of perturbative QCD.  Ed. Frontieres, 1991.
  * [2] G. Altarelli and G. Parisi, “Asymptotic freedom in parton language”, Nuclear Physics B 126 no. 2, (1977) 298–318.
  * [3] Particle Data Group Collaboration, P. A. Zyla et al., “Review of Particle Physics”, PTEP 2020 no. 8, (2020) 083C01. And 2021 update.
  * [4] E. G. de Oliveira, A. D. Martin, M. G. Ryskin, and A. G. Shuvaev, “Treatment of heavy quarks in QCD”, Eur. Phys. J. C 73 no. 10, (2013) 2616, arXiv:1307.3508 [hep-ph].
  * [5] CMS Collaboration, A. M. Sirunyan et al., “Measurement of the Splitting Function in pp and Pb–Pb Collisions at $\sqrt{s_{\textup{NN}}}=$5.02\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$”, Phys. Rev. Lett. 120 no. 14, (2018) 142302, arXiv:1708.09429 [nucl-ex].
  * [6] ALICE Collaboration, S. Acharya et al., “Measurement of the groomed jet radius and momentum splitting fraction in pp and Pb–Pb collisions at $\sqrt{s_{\textup{NN}}}=$5.02\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$”, Phys. Rev. Lett. 128 no. 10, (2022) 102001, arXiv:2107.12984 [nucl-ex].
  * [7] ATLAS Collaboration, G. Aad et al., “Measurement of soft-drop jet observables in pp collisions with the ATLAS detector at $\sqrt{s}=$13\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$”, Phys. Rev. D 101 no. 5, (2020) 052007, arXiv:1912.09837 [hep-ex].
  * [8] STAR Collaboration, J. Adam et al., “Measurement of groomed jet substructure observables in p+p collisions at $\sqrt{s}=$200\text{\,}$gigaelectronvolt{\mathrm{GeV}}$$$ with STAR”, Phys. Lett. B 811 (2020) 135846, arXiv:2003.02114 [hep-ex].
  * [9] L. Cunqueiro and M. Płoskoń, “Searching for the dead cone effects with iterative declustering of heavy-flavor jets”, Phys. Rev. D 99 no. 7, (2019) 074027, arXiv:1812.00102 [hep-ph].
  * [10] P. Ilten, N. L. Rodd, J. Thaler, and M. Williams, “Disentangling heavy flavor at colliders”, Phys. Rev. D 96 no. 5, (2017) 054019, arXiv:1702.02947 [hep-ph].
  * [11] ALICE Collaboration, S. Acharya et al., “Direct observation of the dead-cone effect in QCD”, arXiv:2106.05713 [nucl-ex].
  * [12] Y. L. Dokshitzer, V. A. Khoze, and S. I. Troian, “On specific QCD properties of heavy quark fragmentation (‘dead cone’)”, J. Phys. G 17 (1991) 1602–1604.
  * [13] A. Larkoski, S. Marzani, J. Thaler, A. Tripathee, and W. Xue, “Exposing the QCD Splitting Function with CMS Open Data”, Phys. Rev. Lett. 119 no. 13, (2017) 132003, arXiv:1704.05066 [hep-ph].
  * [14] A. J. Larkoski, S. Marzani, and J. Thaler, “Sudakov safety in perturbative QCD”, Phys. Rev. D 91 no. 11, (2015) 111501, arXiv:1502.01719 [hep-ph].
  * [15] P. Cal, K. Lee, F. Ringer, and W. J. Waalewijn, “The soft drop momentum sharing fraction $z_{\textnormal{g}}$ beyond leading-logarithmic accuracy”, arXiv:2106.04589 [hep-ph].
  * [16] ALICE Collaboration, S. Acharya et al., “Exploration of jet substructure using iterative declustering in pp and Pb–Pb collisions at LHC energies”, Phys. Lett. B 802 (2020) 135227, arXiv:1905.02512 [nucl-ex].
  * [17] ALICE Collaboration, S. Acharya et al., “Measurements of the groomed jet radius and momentum splitting fraction with the soft drop and dynamical grooming algorithms in pp collisions at $\sqrt{s}=5.02\ \mathrm{TeV}$”, arXiv:2204.10246 [nucl-ex].
  * [18] Z.-B. Kang, K. Lee, X. Liu, D. Neill, and F. Ringer, “The soft drop groomed jet radius at NLL”, JHEP 02 (2020) 054, arXiv:1908.01783 [hep-ph].
  * [19] M. Cacciari, G. P. Salam, and G. Soyez, “The anti-$k_{\textnormal{t}}$ jet clustering algorithm”, JHEP 04 (2008) 063, arXiv:0802.1189 [hep-ph].
  * [20] Y. L. Dokshitzer, G. D. Leder, S. Moretti, and B. R. Webber, “Better jet clustering algorithms”, JHEP 08 (1997) 001, arXiv:hep-ph/9707323.
  * [21] F. A. Dreyer, G. P. Salam, and G. Soyez, “The Lund Jet Plane”, JHEP 12 (2018) 064, arXiv:1807.04758 [hep-ph].
  * [22] ALICE Collaboration, S. Acharya et al., “Measurements of the groomed and ungroomed jet angularities in pp collisions at $\sqrt{s}=$5.02\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$”, arXiv:2107.11303 [nucl-ex].
  * [23] A. J. Larkoski, S. Marzani, G. Soyez, and J. Thaler, “Soft drop”, JHEP 05 (2014) 146, arXiv:1402.2657 [hep-ph].
  * [24] OPAL Collaboration, P. D. Acton et al., “A Study of differences between quark and gluon jets using vertex tagging of quark jets”, Z. Phys. C 58 (1993) 387–404.
  * [25] CMS Collaboration, A. Tumasyan et al., “Study of quark and gluon jet substructure in Z+jet and dijet events from pp collisions”, JHEP 01 (2022) 188, arXiv:2109.03340 [hep-ex].
  * [26] ALICE Collaboration, K. Aamodt et al., “The ALICE experiment at the CERN LHC”, JINST 3 (2008) S08002.
  * [27] ALICE Collaboration, B. Abelev et al., “Performance of the ALICE Experiment at the CERN LHC”, Int. J. Mod. Phys. A 29 (2014) 1430044, arXiv:1402.4476 [nucl-ex].
  * [28] ALICE Collaboration, S. Acharya et al., “Measurement of the production of charm jets tagged with D0 mesons in pp collisions at $\sqrt{s}=$7\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$”, JHEP 08 (2019) 133, arXiv:1905.02510 [nucl-ex].
  * [29] M. Cacciari, G. P. Salam, and G. Soyez, “FastJet user manual”, Eur. Phys. J. C 72 (2012) 1896, arXiv:1111.6097 [hep-ph].
  * [30] ALICE Collaboration, S. Acharya et al., “Measurement of the production of charm jets tagged with $\textup{D}^{0}$ mesons in pp collisions at $\sqrt{s}$ = 5.02 and 13 TeV”, arXiv:2204.10167 [nucl-ex].
  * [31] ALICE Collaboration, S. Acharya et al., “Groomed jet substructure measurements of charm jets tagged with $\textup{D}^{0}$ mesons in pp collisions at $\sqrt{s}=$13\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$”, 2020.  https://cds.cern.ch/record/2719005. ALICE-PUBLIC-2020-002.
  * [32] P. Skands, S. Carrazza, and J. Rojo, “Tuning PYTHIA 8.1: the Monash 2013 tune”, Eur. Phys. J. C 74 no. 8, (2014) 3024, arXiv:1404.5630 [hep-ph].
  * [33] T. Sjostrand, S. Mrenna, and P. Z. Skands, “PYTHIA 6.4 physics and manual”, JHEP 05 (2006) 026, arXiv:hep-ph/0603175.
  * [34] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, “An introduction to PYTHIA 8.2”, Comput. Phys. Commun. 191 (2015) 159–177, arXiv:1410.3012 [hep-ph].
  * [35] R. Brun, F. Bruyant, M. Maire, A. C. McPherson, and P. Zanarini, GEANT 3 : user’s guide Geant 3.10, Geant 3.11; rev. version.  CERN, Geneva, 1987.  https://cds.cern.ch/record/1119728.
  * [36] S. Alioli, P. Nason, C. Oleari, and E. Re, “A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX”, JHEP 06 (2010) 043, arXiv:1002.2581 [hep-ph].
  * [37] D. J. Lange, “The EvtGen particle decay simulation package”, Nucl. Instrum. Meth. A 462 (2001) 152–155.
  * [38] G. D’Agostini, “A multidimensional unfolding method based on Bayes’ theorem”, Nucl. Instrum. Meth. A 362 (1995) 487–498.
  * [39] M. Cacciari, P. Nason, and R. Vogt, “QCD predictions for charm and bottom production at RHIC”, Phys. Rev. Lett. 95 (2005) 122001, arXiv:hep-ph/0502203.
  * [40] A. J. Larkoski, J. Thaler, and W. J. Waalewijn, “Gaining (Mutual) Information about Quark/Gluon Discrimination”, JHEP 11 (2014) 129, arXiv:1408.3122 [hep-ph].
  * [41] ALICE Collaboration, S. Acharya et al., “Future high-energy pp programme with ALICE”,. https://cds.cern.ch/record/2724925. ALICE-PUBLIC-2020-005.
  * [42] H. T. Li and I. Vitev, “Inverting the mass hierarchy of jet quenching effects with prompt b-jet substructure”, Phys. Lett. B 793 (2019) 259–264, arXiv:1801.00008 [hep-ph].
  * [43] H. T. Li, Z. L. Liu, and I. Vitev, “Heavy flavor jet production and substructure in electron–nucleus collisions”, Phys. Lett. B 827 (2022) 137007, arXiv:2108.07809 [hep-ph].

##  Appendix A Supplementary material

The measured $z_{\textnormal{g}}$ distribution is also compared to analytical
calculations performed in the soft-collinear effective theory (SCET) framework
with resummation to modified leading-logarithmic accuracy (MLL) [42, 43],
presented in Fig. 3, which are consistent with the data. The calculations do
not include hadronisation effects and are performed utilising a leading-order
evaluation of the charm cross section, which experimentally corresponds to
jets with both charged and neutral information. It is expected that the
low-$z_{\textnormal{g}}$ region of the calculations are affected by non-
perturbative effects. For comparison with the SCET predictions, the
measurement is normalised to the number of jets passing the Soft Drop
condition in each sample, as the calculations are provided for
$z_{\textnormal{g}}\geq 0.1$ and do not include an untagged fraction.

Figure 3: The $z_{\textnormal{g}}$ distribution of prompt
$\textup{D}^{0}$-tagged jets compared to that of inclusive jets for $15\leq
p_{\textnormal{T}}^{\textnormal{jet
ch}}<$30\text{\,}$\mathrm{missing}\mathrm{revert_as}\mathrm{=}{presentation}\mathrm{]}{\text{/}\mathrm{GeV}\mathrm{\text{c}}}{\mathrm{GeV}\text{/}\mathrm{\text{c}}}$$$
in pp collisions at $\sqrt{s}=$13\text{\,}$teraelectronvolt{\mathrm{TeV}}$$$,
normalised to the number of jets passing the Soft Drop condition. Model/data
ratios are shown in the bottom panels for soft-collinear effective theory
calculations [42, 43].


# Calculating polaron mobility in halide perovskites

Jarvist Moore Frost  Department of Materials, Imperial College London,
Exhibition Road, London SW7 2AZ, UK<EMAIL_ADDRESS>Centre for
Sustainable Chemical Technologies and Department of Chemistry, University of
Bath, Claverton Down, Bath BA2 7AY, UK

###### Abstract

Lead halide perovskite semiconductors are soft, polar, materials. The strong
driving force for polaron formation (the dielectric electron-phonon coupling)
is balanced by the light band effective-masses, leading to a strongly-
interacting large-polaron. A first-principles prediction of mobility would
help understand the fundamental mobility limits. Theories of mobility need to
consider the polaron (rather than free-carrier) state due to the strong
interactions. In this material we expect that at room temperature polar-
optical phonon mode scattering will dominate, and so limit mobility. We
calculate the temperature-dependent polaron mobility of hybrid halide
perovskites by variationally solving the Feynman polaron model with the
finite-temperature free-energies of Ōsaka. This model considers a simplified
effective-mass band-structure interacting with a continuum dielectric of
characteristic response frequency. We parametrise the model fully from
electronic-structure calculations. In methylammonium lead iodide at
$300\text{\,}\mathrm{K}$ we predict electron and hole mobilities of $133$ and
$94\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
respectively. These are in acceptable agreement with single-crystal
measurements, suggesting that the intrinsic limit of the polaron charge
carrier state has been reached. Repercussions for hot-electron photo-excited
states are discussed. As well as mobility, the model also exposes the dynamic
structure of the polaron. This can be used to interpret impedance measurements
of the charge-carrier state. We provide the phonon-drag mass-renormalisation,
and scattering time constants. These could be used as parameters for larger-
scale device models and band-structure dependent mobility simulations.

##  I Introduction

Hybrid halide perovskites are a new class of solution processed
semiconductorsManser  _et al._ (2016). They are of considerable practical
interest, mainly due to the high photovoltaic action and potential low cost.
One important figure of merit for a semiconductor is the charge-carrier
mobility.

The material is polar and soft. This leads to a large dielectric electron-
phonon coupling. Bare charge-carriers are dressed in phonon excitations,
forming a polaron.

In this work we calculate the temperature-dependent polaron mobility. We write
custom codes to variationally solve the Feynman polaron modelFeynman (1955)
for finite temperatureOsaka (1959). The inputs to this model are fully
specified by electronic-structure calculations, with no free or empirical
parameters. We then solve for low-field temperature-dependent mobility with
the FHIPFeynman  _et al._ (1962), KadanoffKadanoff (1963) and Hellwarth et
al.Hellwarth and Biaggio (1999) methods.

In a photovoltaic device architecture, the mobility limits the thickness of
the active layer. In a laser diode, carrier-mobility suggests whether
population-inversion (and thus lasing) is feasible. Different measures of
mobility have different systematic errors, probe different charge carrier
densities, are only sometimes selective to hole and electron mobilities, and
pose general experimental challenges that limit temperature-dependent
(cryostat) measurements. To measure intrinsic mobilities, considerable
technical effort must be made in sample preparation, to create sufficient pure
and crystalline materials.

Methylammonium lead iodide perovskite (MAPI) is the most well studied of the
halide perovskites. Time-resolved microwave-conductivity (TRMC) is a
contactless method of measuring mobility, suitable for temperature dependent
measurements. Measures on polycrystalline MAPI give a room-temperature
mobility-yield product of
$35\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$Milot
_et al._ (2015). A more recent measure of a single-crystal gives
$115\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$Semonin
_et al._ (2016). It is not known a priori whether these are the maximum
achievable, or which scattering process limits the mobility.

A fully ab-initio method to predict mobilities for infinite perfect crystals
would help understand the implicit limitations of a new material.

From an electronic structure calculation, one can fit the dispersion relation
(energy vs crystal momentum) for the band extrema to a quadratic function. The
curvature of this quadratic is the bare-electron band effective-mass. This
describes a quasi-particle for an singly-charged excitation in the interacting
electron system, with a mass renormalised relative to that of the free
electron. In lead halide perovskites proper inclusion of spin-orbit-coupling
leads to a complex and non-quadratic dispersion relationBrivio  _et al._
(2014). Nevertheless effective-masses of of 0.12 $m_{e}$ and 0.15 $m_{e}$ can
be ascribedBrivio  _et al._ (2014) to electrons and holes respectively
($m_{e}$ is the bare-electron mass).

Continuous thermal disorder interacting with spin-orbit coupling further
complicates this picture. The band-structure is a dynamic and spin-dependent
object, responding to local fluctuations in electric fieldAzarhoosh  _et al._
(2016). The lattice dynamically responds to the presence of a charge-carrier,
perturbing the electronic structure of the orbitals in which the charge
carrier sits. The full treatment of charge-carrier mobility in these materials
is clearly complex, and well beyond what has been attempted so far. Given the
dynamic nature of the electronic structure, it is not at all clear whether
there is anything to be gained by including the rigid-band structure of a
particular stochastic realisation of a disordered material, in a mobility
calculation.

The bare-electron rigid-band effective-mass is an inertia term. Irrespective
of the effective-mass, in the absence of scattering, mobilities are infinite.
Yet a small effective-mass (highly disperse bands) is often used as a
predictor of a high mobility. This correlation assumes that the scattering
processes are independent of (or at least weakly correlated with) the
effective-mass. This is the constant scattering time approximation. Small
effective-masses therefore correlate with the charge-carrier accelerating to a
higher velocity in the time before relaxation events. This scattering time has
to be inserted into the theory, and is often used as a free parameter. For
hybrid halide perovskites, in analogy with the similar dispersion relation and
effective mass present in CdTe, and assuming typical covalent semiconductor
scattering times, a mobility as high as
$1000\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
would be inferred. This is a considerable overestimate compared to what has
been observed.

As the material is soft, the Debye temperature for low-energy optical modes in
this system is $\approx$ $110\text{\,}\mathrm{K}$ ($2.25\text{\,}\mathrm{THz}$
Brivio  _et al._ (2015)). Charge carriers above this temperature can scatter
inelastically by emitting optical phononsRidley (2013), a highly dissipative
process. Impurity and acoustic-phonon scattering are elastic, so less
dissipative. They would be expected to dominate mobility below this threshold
temperature.

The same processes are at play in a standard covalent semiconductorRidley
(2013), but the relevant temperatures are shifted in in proportion to the
stiffening of the phonon modes. For example, in GaP the mobility above
$400\text{\,}\mathrm{K}$ is limited by optical-mode scatteringRode (1972). The
Debye temperature is $580\text{\,}\mathrm{K}$ for the
$12.06\text{\,}\mathrm{THz}$ LO mode of GaPPödör (1983)).

There is evidence that anti-bonding frontier orbitals (equivalently,
’inverted’ band structure) in the material classyi Huang and Lambrecht (2013)
leads to a ’defect-tolerance’, with defect states repelled away from the
gapBrandt  _et al._ (2015). With fewer electrically-active impurities there is
less impurity scattering. The large static dielectric constant also shields
charged impurities, reducing the cross-section for scattering.

Though crystalline, halide perovskites are polar and soft. This leads to a
large zero-frequency (static) dielectric constant. The lattice deforms around
a charge carrier, localising it. This provides a large driving force for
stabilising the polaron.

The halide perovskites are unusual in that they are highly polar yet possess
light effective masses. There is a large dielectric electron-phonon
interaction, yet the kinetic energy of the electron is sufficient to keep the
polaron from collapsing into a localised (small-polaron) state. A correct
transport theory in this material must incorporate the strongly interacting
nature of the polaronic charge-carriers. Previous theoretical studies of
mobility in this materialMotta  _et al._ (2015); Zhao  _et al._ (2016);
Filippetti  _et al._ (2016) have generally solved the Boltzmann equation in
the relaxation time approximation, with no explicit treatment of the polaron
state.

In this work we return to direct theories of polaron mobility. We develop
general codes to apply these methods to arbitrary polar systems, and present
results for halide perovskite materials. We predict the temperature dependence
and absolute value of polaron mobility in halide perovskites, without
empirical parameters. The system we consider is highly idealised. The
electronic band-structure is present only as an effective mass. The physical
response of the lattice is parametrised by optical and static dielectric
constants, and an effective dielectric-response phonon frequency. We only
consider the polaron state, and its scattering with this characteristic
phonon. As such, the calculated mobilities are an upper bound for a perfect
single crystal.

We predict the same temperature dependence of mobility as that shown by time-
resolved microwave conductivity of a polycrystalline sampleMilot  _et al._
(2015). Good absolute agreement is found with (room-temperature only) single
crystal terahertz conductivity measurementsSaidaminov  _et al._ (2015);
Semonin  _et al._ (2016).

In our model we ignore the perovskite phase transitions present in the system,
assuming a constant effective mass and phonon response. We do not consider
additional scattering from dissipative rearrangement of the ions, which may be
particularly relevant for the (dynamic) cubic phase.

Our calculations explain the published data—the model has predictive power.
This suggests that the intrinsic performance limits of the material have been
realised, and that polaron optical-phonon scattering dominates room
temperature mobility.

As well as predicting the mobility (a phenomenological quantity), we
characterise the nature of the polaron state, calculating phonon-drag mass
renormalisation and scattering rates, and a polaron size. These can complement
experiments to characterise the charge carrier state. In particular, Hendry et
al.Hendry  _et al._ (2004) applied the same models we use here to polaron
mobility, effective mass, and scattering time in TiO2.

Further measurements of temperature-dependent mobility in the halide
perovskites will help understand the charge-carrier state and scattering
processes in these materials.

##  II Methods

###  II.1 Electron phonon coupling

The interaction between charge-carriers and lattice vibrations is mediated by
the electron-phonon coupling. This is challenging to calculate by first
principleGiustino (2017). A common method is to use density-functional
perturbation-theory. This should continue the greater contribution for a hard
(covalent) system, where the atomic motion is small. For a soft system such as
hybrid perovskites, this lattice distortion contribution is large. Non-
perturbative methods may be necessaryWhalley  _et al._ (2016); Saidi  _et al._
(2016). Within perturbation theory, there is uncertainty about which diagrams
(types and orders of perturbative interaction) to include in the summation.
The spiky nature of the numerical integration across the double (phonon and
electron) Brillouin zones makes the calculations heavy and convergence
difficult.

For polar systems, the major contribution to electron-phonon coupling comes
from the long-ranged electric fields generated by the atomic displacements.
This is specified by the dielectric response of the material. This
considerably simplifies the problem.

FröhlichFrohlich (1952) first constructed a Hamiltonian for a system of
independent (i.e. low density) electrons interacting with harmonic (non-
interacting) polar optical phonons. The dimensionless Fröhlich parameter
$\alpha$ of dielectric electron-phonon coupling is

$\alpha=\frac{1}{4\pi\epsilon_{0}}\frac{1}{2}\left(\frac{1}{\epsilon_{\infty}}-\frac{1}{\epsilon_{S}}\right)\frac{e^{2}}{\hbar\Omega}\left(\frac{2m_{b}\Omega}{\hbar}\right)^{1/2}.$
(1)

This is fully defined by the material specific properties of the optical
($\epsilon_{\infty}$) and static ($\epsilon_{S}$) dielectric constants, the
bare-electron band effective-mass ($m_{b}$), and a characteristic phonon
angular frequency ($\Omega$). As usual, $2\pi\hbar$ is Planck’s constant,
$\epsilon_{0}$ the permittivity of free space, and $e$ the electron charge.

The optical and static dielectric constants form a pre-factor for this
electron-phonon interaction. This is the response of the lattice, modelled as
a continuum dielectric, beyond the boundary of the polaron. Within the
polaron, only the electronic excitations (optical dielectric constant) can
keep up with the motion of the electron. Outside the boundary of the polaron,
the polar lattice excitations can also respond.

###  II.2 Multiple phonon branches

The static dielectric constant can be calculated by a summation over the
Brillouin zone centre (gamma point) harmonic phonon modes. Only these modes
contribute, as they have no phase factor which otherwise leads to a zero
contribution when integrating over real-space. These phonon mode eigenvectors
are used to project the Born effective charges, to give an effective dipole.
This is also the infrared activity of the mode. Integrating through these
individual Lorentzian responses leads to the dielectric functionGonze and Lee
(1997),

$\epsilon_{S}=\epsilon_{\infty}+\frac{4\pi}{V}\sum_{i=0}^{N}\frac{(Z^{*}U^{*}_{i})(Z^{*}U_{i})}{\Omega_{i}^{2}}.$
(2)

Here, the summation is over the $N$ phonon modes, normalised by $V$ the unit-
cell volume. $Z^{*}$ are the Born effective charges, $U_{i}$ a specific phonon
eigenmode, $\Omega_{i}$ the phonon frequency.

The factor of $\Omega_{i}^{-2}$ means the summation is dominated by the lowest
energy (frequency) modes which are infrared-active (polar). In the case of a
two-atom unit cell such as most covalent semiconductors, there is one
infrared-active mode, the linear-optical (LO) mode. Most of the polaron
literature is constructed in this framework of a single polar response
frequency.

Hellwarth et al.Hellwarth and Biaggio (1999) provide a prescription to reduce
multiple infrared-active phonon branches to a single equivalent dielectric-
response phonon frequency.

###  II.3 Feynman polaron model

Feynman solved the Frölich Hamiltonian in an innovative mannerFeynman (1955),
where the electron interacting with a cloud of independent (harmonic) phonon
excitations is path-integrated over the (phonon) quantum field. The electron
interacts with the disturbance it has previous created in passing through the
lattice, which exponentially dies out in time. This, Coulomb-like, interaction
would be expected to depend on inverse-distance. Instead, a harmonic model is
constructed. This can be analytically path-integrated. An exponential
dampening factor $w$ adds a degree of freedom for anharmonicity, and the
interaction path-integral is scaled by a coupling strength $C$.

These model parameters are then varied to minimise the (athermal) ground state
energy, for a given coupling strength $\alpha$ and phonon response frequency
$\Omega$.

The resulting model is a single particle system where the electron interacts
with a time-retarded potential. The Hamiltonian in the centre-of-mass frame
can be expressed as an electron coupled to a finite mass ($M$, expressed in
units of the charge carrier effective mass) with a (harmonic) spring-constant
($k$). These give rise to an angular-frequency for the oscillation of the
energy between the electron and phonons of $w=\sqrt{\frac{k}{M}}$.

It is convenient to work in this angular-frequency variable and make a further
substitution of $v^{2}=w^{2}+\frac{4C}{w}$. Here $C$ is the electron-phonon
coupling coefficient of the model Lagrangian, related back to the spring-
constant by $k=\frac{4C}{w}$Feynman (1955).

Feynman’s model is non-perturbative and so correct to all orders in the
coupling constant $\alpha$. Up until this point, theories had either been
based on an assumption of weak coupling ($\alpha<<1$) and perturbative, or
assumed strong coupling ($\alpha>10$). Many systems of experimental interest
have an intermediate coupling. We will show later that inorganic halide
perovskites have a characteristic $\alpha$ coupling of 1 to 2, whereas the
additional dielectric response of the molecular cation gives the hybrid halide
perovskites an $\alpha$ of 2 to 3.

ŌsakaOsaka (1959) extended Feynman’s athermal variational solution by
providing finite-temperature free-energies of the coupled electron-phonon
system. The parameters can then be varied to minimise the total polaron free-
energy (at specific temperature). The model thus becomes explicitly
temperature dependent, rather than using the athermal variational solution
with a temperature-dependent mobility theory. Here we use a more recent
presentation of the Ōsaka free-energies by Hellwarth et al.Hellwarth and
Biaggio (1999), which have been made more amenable to numeric computation.

###  II.4 Polaron mobility

By itself, the variational (finite-temperature extended) Feynman model
provides an effective-mass of the polaron quasi-particle representing the
phonon-drag term, and a spring coupling constant for the exchange of energy
between the electron and its coupled phonon cloud. As with a bare-electron
band effective-mass, there is no dissipative (frictional) term. To calculate a
polaron mobility, we must analyse the response of this system to a
perturbation. As our main concern is with application to photovoltaic
materials, we are interested in the small-field (direct current, DC) limit.

The models used herein consider the electronic band structure only in an
effective-mass approximation.

The initial mobility work of FHIPFeynman  _et al._ (1962) directly used the
response of the polaron centre-of-mass coordinate to an applied field,
containing interactions with a dynamically maintained steady state of
thermally excited phonons. A general expression for the impedance is given
therein (equations 46 and 47). A low temperature approximation is made (a
power-series expansion in the internal parameter $b$, small for low-
temperature), leading to a mobility generally known as $\mu_{FHIP}$,

$\mu_{FHIP}=\left(\frac{w}{v}\right)^{3}\frac{3e}{4m_{b}}\frac{exp\left(\beta\right)}{\Omega\alpha\beta}exp\left(\frac{v^{2}-w^{2}}{w^{2}v}\right).$
(3)

Here $w$ and $v$ are the (variational) parameters specifying the polaron
model, $m_{b}$ is the bare-electron band effective-mass, $\Omega$ the phonon
angular frequency, and $\beta=\hbar\Omega/k_{B}T$ is a reduced thermodynamic
temperature in units of the phonon energy.

This formalism was quickly realised to be pathological for high
temperatures—at the thermal energy ($k_{b}T$) equaling the phonon energy
($\hbar\Omega$) there is a resonance minimising the mobility, but the mobility
thereafter increases as a function of temperature. Soon afterwards,
KadanoffKadanoff (1963) provided a similar mobility identity based on a
solution to the Boltzmann equation, applicable to all $\alpha$. This Boltzmann
solution implicitly assumes independent scattering events, still formally
limiting the model to low temperature. We refer to this as $\mu_{K}$,

$\mu_{K}=\left(\frac{w}{v}\right)^{3}\frac{e}{2m_{b}}\frac{exp\left(\beta\right)}{\Omega\alpha}exp\left(\frac{v^{2}-w^{2}}{w^{2}v}\right).$
(4)

Though derived from different assumptions, compared to the FHIP there is a
additional factor of $3/2\beta$. This is now understoodPeeters and Devreese
(1983) to formally arise from the order in which the limits are taken in the
asymptotic approximations. Physically, the FHIP considers stimulated-emission
and absorption of phonons, which saturates and balances, once there is a
thermal population (Bose-Einstein statistics). The Kadanoff identity adds
spontaneous emission of phonons. Empirically this corrects the pathological
behaviour of FHIP for temperatures above the Debye temperature, the Kadanoff
mobility asymptotically approaching a constant.

The relaxation-time approximation (independent scattering) used by Kadanoff
allows a direct evaluation of a scattering rate, based on the thermal
population of phonons $\bar{N}=exp(\beta)$. The rate ($\Gamma$) is evaluated
for small momentum exchanges (low-fields) as

$\Gamma_{0}=2\alpha\bar{N}\sqrt{(M+1)}exp(-M/v).$ (5)

Here $M$ is the phonon-drag effective-mass (in units of the band effective-
mass, $m_{b}$).

This rate $\Gamma_{0}$ is expressed in reduced phonon units. Multiplying with
the phonon frequency $\frac{\Omega}{2\pi}\Gamma_{0}$ gives a real-time rate.
This scattering rate can be directly related to $\mu_{K}$ (as it is assumed to
follow a Boltzmann equation) by

$\mu_{K}=\frac{e}{\Omega m_{b}(M+1)\Gamma_{0}}.$ (6)

More recently, Hellwarth et al.Biaggio  _et al._ (1997); Hellwarth and Biaggio
(1999) returned to the general response theory of FHIP (Equations 46 and 47 in
Feynman  _et al._ (1962)). Rather than taking a low-temperature limit and
perform a power-expansion, they contour integrate for the self-energy of the
perturbed polaron.

$a^{2}=(\beta/2)^{2}+\left(\frac{v^{2}-w^{2}}{w^{2}v}\right)\beta coth(\beta
v/2),$ (7)

$b=\left(\frac{v^{2}-w^{2}}{w^{2}v}\right)\frac{\beta}{sinh(\beta v/2)},$ (8)

$K=\int_{0}^{\infty}du\left(u^{2}+a^{2}-bcos(vu)\right)^{-\frac{3}{2}}cos(u).$
(9)

This integral gives the polaron response to a first order change in the
driving force, providing an analogous role to the scattering rate in
Kadanoff’s Boltzmann construction, to give,

$\mu_{H}=\left(\frac{w}{v}\right)^{3}\frac{3e}{m_{b}}\frac{\sqrt{\pi}sinh(\beta/2)}{\Omega\alpha\beta^{\frac{5}{2}}}K^{-1}.$
(10)

This is as in reference Hellwarth and Biaggio (1999), transformed into S.I.
units and rearranged for easier comparison with $\mu_{FHIP}$ and $\mu_{K}$.

Both Biaggio et al.Biaggio  _et al._ (1997) and Hellwarth et al.Hellwarth and
Biaggio (1999), approximate $b=0$, allowing for an analytic solution of $K$
with modified Bessel functions. Here we do the full integration. We note that
both papers have a typographic error in the formula for $b$, possessing a
spurious term of $b$ on the right hand side. The form here is as given by
FHIPFeynman  _et al._ (1962).

In custom codes we reimplement the Hellwarth et al.Hellwarth and Biaggio
(1999) posing of Ōsaka’sOsaka (1959) finite-temperature variational solution
to Feynman’sFeynman (1955) model. Integration is numeric using an adaptive
Gauss-Kronrod quadrature algorithm. The Ōsaka free-energies (which include
numeric integrals) are automatically differentiated in a forward-mode to
produce gradients. These gradients are used by a BFGS algorithm to find the
optimal (temperature-dependent) $v$ and $w$. These model parameters directly
give the FHIPFeynman  _et al._ (1962) and KadanoffKadanoff (1963) mobility. By
performing the contour integration for the polaron self-energy numerically, we
calculate an Hellwarth et al.Hellwarth and Biaggio (1999) mobility (with a
slight refinement by considering $b\neq 0$). Codes are providedFrost (2017) to
encourage the application of these methods to other systems of interest.

##  III Results

###  III.1 Methylammonium lead iodide perovskite

Methlyammonium lead halide is the most well studied of the hybrid halide
perovskites. We use the previously mentioned QSGW effective masses of 0.12
(electron) and 0.15 (hole)Brivio  _et al._ (2014). We take $4.5$ and
$24.1$Frost  _et al._ (2014) as the optical (QSGW) and static (harmonic
phonon, DFT) dielectric constants.

We note that this static dielectric constant only includes the harmonic
response of the phonons, additional (slower) terms may come from the
(anharmonic) realignment of the polar methylammonium. We can estimate this
contribution from a generalised form of Onsager theory for polar
liquidsKirkwood (1939). Ignoring the response of the local environment and
assuming the ions totally free to realign (i.e. most applicable to the cubic
phase), this is $\epsilon_{p}=\frac{4}{3}\pi N\frac{\mu^{2}}{3kT}$. With the
methlyammonium-dipole $\mu=2.2$D Frost  _et al._ (2014), and $N$ the number
density (one methylammonium per $6.2\text{\,}\mathrm{\SIUnitSymbolAngstrom}$
cubic unit cell), this gives $\epsilon_{p}=+8.9$.

It is well known that the increasing frustration of motion of these dipoles as
the material passes through the second-order tetragonal phase transition
towards the fully hindered orthorhombic ground state gives rise to a divergent
dielectric responseOnoda-Yamamuro  _et al._ (1992). It is not clear how to
integrate this (slow, dissipative) response into a formal polaron transport
theory.

One might assume that the additional response in the cubic phase, and in the
divergence approaching the orthorhombic phase, would provide an additional
dissipation of electron energy, and so lower mobility. We ignore these terms
in the present work. Disorder in the hybrid material will generate further
localisation pressure, perhaps even forming small polaronsMa and Wang (2015);
Neukirch  _et al._ (2016).

Hellwarth et al.Hellwarth and Biaggio (1999) provides two approximation
schemes for reducing multiple phonon branches in the polaron problem to an
single dielectric-response frequency. In this work we use the more simple
athermal ’B’ scheme. The temperature-dependent ’A’ scheme may offer more
temperature-dependent phonon structure in the polaron model, and therefore the
resulting temperature-dependent mobilities. This will be the subject of future
work.

We take our modes from ab-initio phonon and infrared-activity
calculationsBrivio  _et al._ (2015). Applying the ’B’ scheme to all modes
(including the 18 high frequency molecular modes) gives a seemingly unphysical
characteristic response of $10.0\text{\,}\mathrm{THz}$. It is not clear
whether this is due to errors in our estimates of the infrared activity of the
molecular modes, or whether the Hellwarth et al. summation is only well
defined for long range consummate lattice distortions.

Applied to the first 15 (non-intramolecular) modes, the effective dielectric-
response (’LO’) phonon frequency is $2.25\text{\,}\mathrm{THz}$. We use this
value in all work presented herein.

Combined, these data specify the Fröhlich electron-phonon parameter as
$\alpha=2.40$ (electron) and $\alpha=2.68$ (hole).

Bokdam et al.Bokdam  _et al._ (2016) provide a direct evaluation of the total
dielectric function. From a density functional perturbation theory calculation
(harmonic athermal response), the imaginary dielectric function has a peak at
$8\text{\,}\mathrm{meV}$ ($1.9\text{\,}\mathrm{THz}$). This corresponds well
with the estimate from the Hellwarth et al. B scheme applied to the lattice
modes. Including some anharmonic contributions via ab-initio molecular
dynamics, they find a softening of the modes to $4\text{\,}\mathrm{meV}$ at
$300\text{\,}\mathrm{K}$. This implies that polaron effects are actively
strengthened at higher temperature, as the anharmonic softening of the lattice
increases the effective dielectric electron-phonon coupling.

Figure 1:  Predicted Hellwarth hole (cross, red) and electron (circle, blue)
polaron mobility versus temperature, co-plotted against temperature-dependent
time-resolved-microwave-conductivity measurements on polycrystalline films
from Milot et al.Milot  _et al._ (2015) (green, saltire), room temperature
space-charge-limited-current hole-mobility estimate on a single crystal from
Saidaminov et al.Saidaminov  _et al._ (2015) (diamond, brown), and room
temperature time-resolved-microwave-conductivity on a single crystal Semonin
et al.Semonin  _et al._ (2016) (square, purple). As the TRMC data is a
mobility-yield product, the low $8\text{\,}\mathrm{K}$ data point may be
related to the stability of the exciton (none-unitary free carrier yield) at
those temperatures, rather than the mobility actually reducing.

At 300 K we predict a Kadanoff polaron mobility of
$\mu_{e}=$$197\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$,
$\mu_{h}=$$136\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$;
and a Hellwarth mobility of
$\mu_{e}=$$133\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$,
$\mu_{h}=$$94\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$.
The phonon mass renormalisation is $0.37$ for the electron polaron and $0.43$
for the hole polaron. This agrees well with the athermal perturbative (small
$\alpha$) estimate $m_{p}^{*}=\frac{\alpha}{6}+\frac{\alpha^{2}}{40}$ Feynman
(1955) which gives 0.54 and 0.63 respectively.

Table 1: Parameters of the Feynman polaron model as used in this work. Relative high frequency ($\epsilon_{\infty}$) and static ($\epsilon_{S}$) dielectric constants are given in units of the permittivity of free space ($\epsilon_{0}$). Frequency (f) is in $\mathrm{THz}$. Effective mass ($m^{*}$) is in units of the bare electron mass.  Material | $\epsilon_{\infty}$ | $\epsilon_{S}$ | $f$ | $m^{*}$
---|---|---|---|---
MAPbI3-e | 4.5 | 24.1 | 2.25 | 0.12
MAPbI3-h |  |  |  | 0.15
CsPbI3 | 6.1 | 18.1 | 2.57 | 0.12
MAPbI3 111These parameters specify the model as in Sendner et al.Sendner  _et al._ (2016). Effective mode frequency, reduced by a Hellwarth et al. scheme, is from prviate communication. | 5.0 | 33.5 | 3.38 | 0.104
MAPbBr3 | 4.7 | 32.3 | 4.47 | 0.117
MAPbCl3 | 4.0 | 29.8 | 6.42 | 0.2
CsSnI3 222Parameters for holes in cesium tin halides are taken from Huang et al.yi Huang and Lambrecht (2013). Dielectric constants are from Table VII therein, effective masses from Table VI, phonon frequency from table VII. | 6.05 | 48.2 | 4.56 | 0.069
CsSnBr3 | 5.35 | 32.4 | 5.48 | 0.082
CsSnCl3 | 4.80 | 29.4 | 7.28 | 0.140
Table 2: Dielectric electron-phonon coupling ($\alpha$, athermal). Predicted ($300\text{\,}\mathrm{K}$) mobilities (Kadanoff, $\mu_{K}$; Hellwarth, $\mu_{H}$), polaron effective mass renormalisation ($m_{h}^{*}$), Feynman-model variational parameters ($v$ and $w$), Kadanoff polaron relaxation time ($\tau$, $\mathrm{ps}$), and Schultz polaron radius ($r_{f}$, $\mathrm{\SIUnitSymbolAngstrom}$).  Material | $\alpha$ | $\mu_{K}$ | $\mu_{H}$ | $m_{h}^{*}$ | $v$ | $w$ | $\tau$ | $r_{f}$
---|---|---|---|---|---|---|---|---
MAPbI3-e | 2.39 | 197 | 136 | +0.37 | 19.9 | 17.0 | 0.12 | 26.8
MAPbI3-h | 2.68 | 133 | 94 | +0.43 | 20.1 | 16.8 | 0.10 | 25.3
CsPbI3 | 1.35 | 389 | 258 | +0.21 | 16.5 | 15.1 | 0.20 | 39.6
MAPbI3 | 1.71 | 272 | 195 | +0.31 | 13.2 | 11.5 | 0.13 | 43.1
MAPbBr3 | 1.69 | 212 | 157 | +0.36 | 10.1 | 8.7 | 0.12 | 53.6
MAPbCl3 | 2.19 | 73 | 57 | +0.62 | 7.6 | 6.0 | 0.08 | 61.4
CsSnI3 | 1.02 | 703 | 487 | +0.20 | 9.64 | 8.81 | 0.21 | 70.0
CsSnBr3 | 1.09 | 511 | 356 | +0.24 | 8.14 | 7.32 | 0.19 | 77.7
CsSnCl3 | 1.39 | 212 | 147 | +0.36 | 6.51 | 5.59 | 0.14 | 85.4

The variational parameters $v$ and $w$ specify the polaron state. These come
from minimising the finite-temperature ŌsakaOsaka (1959) free energies (Figure
2). These parameters can be mapped to a phonon-drag mass term $M$, and spring
constant $k$, for the effective single-particle coupled harmonic system
(Figure 3).

Following SchultzSchultz (1959), our variational parameters at
$300\text{\,}\mathrm{K}$ give a polaron radius (for holes and electrons) of
$\approx$ $26\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. (See Table 2 for
further values.) This radius is defined as the standard deviation of the
Gaussian wavefunction which would form in the harmonic confining potential of
the polaron. Though the definition is fairly arbitrary, it is of interest in
understanding how commensurate the polaron is to real-space fluctuations in
electrostatic potential.

As with Hellwarth et al.Hellwarth and Biaggio (1999), we find that setting b=0
in our calculation of $\mu_{H}$ makes only $<0.2\%$ difference to the
mobility. This difference increases with temperature.

As a function of temperature (Figure 1), the mobility decays from an infinite
quantity at zero temperature. The Kadanoff mobility asymptotically approaches
$190\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
at high temperature. This is associated with the variational fit approaching a
phonon effective mass of $0.37m_{e}$, while the spring constant $k$ increases
linearly with temperature (Figure 3). By comparison, the (low-temperature)
FHIP mobility has a minimum at $\hbar\Omega=k_{B}T$, and increases linearly
with temperature. The Hellwarth mobility has weak positive temperature
dependence at high temperatures.

The high temperature behaviour of polaron mobility is important for
photovoltaic device operation. Hybrid halide perovskites are predicted to have
extremely low thermal conductivityWhalley  _et al._ (2016), and cooling of
photoexcited states is known to be unusually slowPrice  _et al._ (2015); Yang
_et al._ (2015). As such, it is likely that the initial electron temperature
of a photo-excited state is extremely high. Optical phonons emitted by
electron thermalisation and carrier scattering may remain in the region of the
polaron, reheating the electron state. It is therefore useful and interesting
to know that even with a high temperature ’hot carrier’, the polaron mobility
is finite, and sufficient for photovoltaic device operation.

Figure 2:  Finite-temperature free-energies used in the variational fitting of
the $v$ and $w$ parameters for the coupled electron-phonon system. The upper
bound of the free energy (cross, red) is $F=-(A+B+C)$. Here A (circle, blue)
is total thermodynamic energy of the system, via the log of the partition
function of the density matrix. B (saltire, green) and C (squares, purple) are
the expectation values of the action for the electron-phonon system and the
trial action respectively. See reference Osaka (1959). Data for the electron
(0.12) effective mass in MAPI.

Figure 3:  (Top) Temperature-dependent phonon-drag effective mass (cross,red)
in units of the bare-electron band effective-mass. Polaron scattering-time
(circle, blue) in units of $\mathrm{ps}$. (Bottom) Temperature-dependent
variational parameters $v$ (cross, red) and $w$ (circle, blue) for the coupled
electron phonon system. Data for the electron (0.12) effective mass in MAPI.

The Kadanoff and FHIP mobilities posess at their core a calculation of the
rate of emission and absorption of optical phonons $\Gamma$, taken to the low-
momentum limit $\Gamma_{0}$. KadanoffKadanoff (1963) directly relates this
rate to a relaxation-time approximation of the Boltzmann equation,
$\tau=1/\Gamma$. We thereby access the scattering time for the polaron as a
function of temperature, finding for the electron-polaron
$\tau=$$0.12\text{\,}\mathrm{ps}$ at a temperature of
$300\text{\,}\mathrm{K}$. Combined with the polaron mass renormalisation (the
phonon drag), these values parametrise a temperature-dependent relaxation-time
approximation Boltzmann equation for the polarons in this material, and may be
of use in larger-scale device models. Further, these parameters can be probed
experimentally. Previous workHendry  _et al._ (2004) used Terahertz
spectroscopy to probe the scattering in TiO2, observing a counter-intuitive
relationship between scattering rate and mobility, due to the complementary
relationship between effective-mass and scattering-rate $\mu=e/(m_{p}\Gamma)$.

###  III.2 Powerlaw temperature dependence of mobility

Inference of the dominant electron scattering mechanism in a material often
comes from the circumstantial behaviour of the mobility as a function of
temperature. This temperature-dependent mobility (a phenomenological quantity)
often shows a power-law scaling. For the halide perovskites, early data
suggests a $T^{-\frac{3}{2}}$ exponent, consistent with textbook descriptions
of acoustic-mode scattering dominated mobilitiesRidley (2013). We previous
suggested that activation of multiple soft optical phonon modes could mimic
such (large exponent) behaviourLeguy  _et al._ (2016).

Do polaron mobility theories follow a power-law scaling, and if so, what is
the critical exponent? Displaying the Hellwarth mobility on a log-log scale
(Figure 4 reveals straight-line (power-law) behaviour, for temperatures above
the critical phonon emission energy. Fitting these data to a power-law
($\mu=aT^{k}$) by least-squares minimisation of data above
$100\text{\,}\mathrm{K}$ finds an exponent of $k=-0.46$. This is close to
$-\frac{1}{2}$, which we propose as our best estimate of the power-law
behaviour of polaron optical-phonon scattering dominated mobility, once above
the phonon emission threshold. (This will hold even in the case where there
are multiple distinct optical phonon thresholds, as all modes will produce a
$k=-\frac{1}{2}$ scaling once above their emission threshold.)

This exponent is considerly less than the $k=-1.33$ extracted from a full fit
to the TRMC data of Milot et al.Milot  _et al._ (2015). However, if we
similarly constrain the experimental fit to above $100\text{\,}\mathrm{K}$
(motivated by the uncertainity of measures in the low temperature phase, and
to avoid complications due to multiple phonon scattering routes), we find an
exponent of $k=-0.95$. Further temperature-dependent measures of mobility will
provide stronger statistical evidence from which to infer the nature of
mobility-limit scattering in these materials.

Figure 4:  Hellwarth electron mobility for MAPI (down-triangles, blue)
coplotted on a log-log axis with TRMC temperature dependent data of (hexagons,
green) of Milot et al.Milot  _et al._ (2015). Fit lines are by a least-squares
fit to a power-law ($\mu=aT^{k}$), in a linear space. Fitting to above
$100\text{\,}\mathrm{K}$, the predicted Hellwarth mobility (this work) has an
exponent (orange trend line) of $-0.46\approx-\frac{1}{2}$. Fitting all but
the spurious first point of the Milot et al. data produces an exponent (purple
trend line) of $k=-0.95$, whereas fitting just above $100\text{\,}\mathrm{K}$
finds a lower exponent (brown trend line) of $k=-0.95$.

###  III.3 Comparison to Sendner et al.

Recently, Sendner et al.Sendner  _et al._ (2016) approached the polaron
mobility from an experimental point of view. They fitted four critical points
to infrared transmissivity measures of thin-films. These provide a simplified
phonon spectrum. This is then used with a generalisation of the Lyddane-Sachs-
Teller relation to derive dielectric constants. This assumes that there are
well defined linear and transverse optical modes as would be found in a
crystal of cubic symmetry. These data were then reduced with a Hellwarth
scheme to an effective dielectric-response phonon frequency (for MAPI) of
$3.38\text{\,}\mathrm{THz}$ (prviate communication). Combined with a DFT
effective mass (0.104), this was then used with a Hellwarth model to derive a
room-temperature mobility of $200\pm
30\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
for MAPI. See Table 1 for a statement of their model, and Table 2 for a cross-
validation with the codes developed in this work. With their parameters and
our custom codes, we can predict the temperature-dependent mobility (Figure
LABEL:figs/Rob-comparison-COLOUR).

Figure 5:  Calculation of temperature-dependent solutions to the Sendner et
al.Sendner  _et al._ (2016) (data therein and by private communication)
parameters, compared to this work. Sendner’s parameters for iodine (satire,
green), bromine (square, purple) and chlorine (cross, orange). Data from this
work are iodide electrons (cross, red), and holes (circle, blue). Approximate
agreement is found in mobility trends, but due to quite different underlying
parameters for the polaron state.

By comparison to the present work, they extrapolate a greater dielectric
response than we calculate from harmonic phonons. This is to be expected as
their measurements will include all anharmonicity (including the molecular
rotation contribution), whereas we consider only the harmonic response.
However it is not immediately clear that the contribution from this
dissipative molecular response should be integrated into polaron theories
built on the perfectly elastic (non dissipative) response of harmonic phonons.

Though the mobility results of Sendner et al. and this work are in broad
agreement, the underlying parameters of the polaron system disagree. Here we
suggest a more polaronic state, with a larger dielectric coupling resulting in
a greater electron-mass renormalisation. Sendner et al. suggest a less
polaronic charge-carrier state, but with a greater scattering strength due to
the larger dielectric constant. Overall this leads to a very similar Kadanoff
relaxation time (see Table 2).

###  III.4 Inorganic lead halide perovskite

To understand the role of the organic cation in polaron formation and mobility
limits of halide perovskites, we make a comparison with the fully inorganic
analogue, cesium lead iodide. We calculate the dielectric constants for cubic
cesium lead iodide (r=$6.78\text{\,}\mathrm{\SIUnitSymbolAngstrom}$), by
density functional perturbation theory. We use the generalised gradient
approximation, in the plane-wave VASP code with a $9\times 9\times 9$
Monkhorst-pack k-space integration, $700\text{\,}\mathrm{eV}$ plane wave cut-
off, and the PBESol functional. We find optical and static dielectric
constants of $6.1$ and $18.1$. The Hellwarth (B-scheme) dielectric effective
phonon frequency is slightly stiffened to $2.57\text{\,}\mathrm{THz}$.
Assuming the same electron effective mass as electrons in MAPI ($0.12$), this
gives a Kadanoff mobility of
$389\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
and a Hellwarth mobility of
$258\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
at $300\text{\,}\mathrm{K}$ (see Figure 6 for temperature-dependence).

This suggests that the additional (harmonic) dielectric response of the
organic ion in MAPI reduces the mobility by a factor of two, as a result of
doubling the dielectric electron-phonon coupling, leading to a doubling in the
polaron effective mass and a halving of the scattering time.

Figure 6:  Predicted mobility for CsPbI3, assuming an effective mass for the
carriers of 0.12. The theories used are Hellwarth (cross, red), Kadanoff
(circle, blue) and FHIP (saltire, green).

###  III.5 Inorganic tin halide perovskite

The cesium tin halide perovskites have been well characterisedyi Huang and
Lambrecht (2013) with the QSGW method for effective masses, and lattice
dynamics with the LDA density functional. Taking values from this publication,
we can solve for the polaron model and mobilities. The calculated hole
mobilities for CsSnI3 at $300\text{\,}\mathrm{K}$ are
$703\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
(Kadanoff) and
$487\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
(Hellwarth). This compares well to measurements of
$585\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
(Hall-effect) and
$400\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
(by ’transport property’)Chung  _et al._ (2012). Temperature dependence is
presented in Figure LABEL:figs/CsSn-halide-mobility-COLOUR.

Figure 7:  Predicted Hellwarth mobility for cesium tin halide perovskites
(Iodine: cross, red; Bromine: circle, blue; Chlorine: saltire, green),
assuming an effective mass for the carriers of 0.12.

In comparison to the lead halide perovskites, the material is less polaronic.
The driving force for polaron stabilisation is greater with the larger
dielectric constants, but the effective masses are smaller, the phonon
spectrum is stiffer. As with MAPI, the agreement between predicted and
measured mobility suggests that polaron mobility describes the charge carrier
state and scattering processes well.

##  IV Discussion

###  IV.1 Prior theoretical studies

Most prior theoretical studies on mobility in hybrid halide perovskitesMotta
_et al._ (2015); Zhao  _et al._ (2016); Filippetti  _et al._ (2016) have
solved the Boltzmann equation in the relaxation time approximation using the
BoltzTraP codesMadsen and Singh (2006).

An effective relaxation-time (scattering-constant) $\tau$ is inserted into a
theory otherwise parametrised from electronic structure calculations. These
works include the band-structure explicitly, whereas we use an effective mass
approximation throughout. In hybrid halide perovskites, the materials are new
and unusual, so there is no equivalent well characterised system (such as in
covalent semiconductors) on which to base this relaxation time.

The first work by Motta et al.Motta  _et al._ (2015) used an empirically
motivated scattering time constant of $1\text{\,}\mathrm{ps}$. An analogy is
made by relating this to the time constant of motion of the methylammonium
ion. However there is no direct correspondence between the timescale of motion
of the ion and the timescale of charge carrier scattering. This led to a
predicted mobility of
$5\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}12\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
for holes and
$2.5\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}10\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
for electrons. At the time these values were comparable to the highest
measured mobilities.

Later, Zhao et al.Zhao  _et al._ (2016) found a relaxation time of
$0.1\text{\,}\mathrm{ps}$ from acoustic-phonon scattering in a golden-rule
formalism with deformation-potential matrix elements. The mobility predicted
from this is
$1000\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$.
Considering charged-impurity scattering with a static dielectric constant of
6.5, and assuming large non-compensated charge defect densities of
$1\text{\times}{10}^{18}\text{\,}{\mathrm{cm}}^{-3}$, these mobilities reduce
to
$100\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$.
The density of charged impurities and static dielectric constant used in this
work seems unphysical.

Filippetti et al.Filippetti  _et al._ (2016) combined contributions from
impurity scattering, acoustic-phonon scattering (approximated by the
deformation potential), and a relaxation time approximation recasting of the
Fröhlich optical electron-phonon scattering element. They use a static
dielectric constant of $60$. Overall they predict a bare-charge electron
scattering time to $9.2\text{\,}\mathrm{fs}$. They find a low charge-density
mobility of
$57\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
(electron) and
$40\text{\,}{\mathrm{cm}}^{2}\text{\,}{\mathrm{V}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
(hole).

Very recently Zhang et al.Zhang  _et al._ (2017) have provided a semi-
classical model of charge transport in hybrid halide perovskites, based on
polarons scattering from acoustic modes. They point out that, by definition,
acoustic scattering produces the experimentally observed $T^{-\frac{3}{2}}$
dependence. However they do not consider inelastic scattering processes, such
as emission of optical phonons.

###  IV.2 This work

This work provides the first prediction of (temperature-dependent) polaron
mobility in hybrid halide perovskites. There are no free parameters.

The bare-electron band effective-masses are taken from QSGW electronic
structure calculationsBrivio  _et al._ (2014). The parameters used in the
calculation are the optical and static dielectric constants, and an effective
phonon frequency. These values we take from prior density-functional-theory
calculations on the harmonic (phonon) response of the materialBrivio  _et al._
(2015).

This fully specifies a model Hamiltonian, for which the finite-temperature
Feynman variational solution for the polaron is made. Numerical solution of
the full DC-mobility theory (without making a Boltzmann transport equation
approximation) provides a temperature-dependent mobility.

These mobilities are an upper bound, as they consider only one scattering
contribution (that of the polaron scattering with optical-phonons) in an
otherwise perfect effective-mass system.

Both the Kadanoff and Hellwarth polaron mobility as a function of temperature
show a trend that agrees well with available data. As the Kadanoff mobility is
based on assuming a Boltzmann equation, this suggests that the underlying
assumption of independent scattering events is approximately correct. The
Hellwarth mobility is constructed by an explicit contour integration for the
self-energy of the perturbed polaron. We expect it to be more generally
applicable, and more accurate for high temperatures.

The agreement with measured mobilities suggests that we have captured the
essential physics of the system.

This suggests that impurity scattering is a relatively minor process. This
adds circumstantial evidence that these semiconductors are ’defect tolerant’,
at least as far as mobility is concerned. Scattering by acoustic phonon modes
is neglected, as it is low-energy and elastic (due to the thermal population
of these low-energy modes). Where we predict a divergence to infinite mobility
at zero temperature, acoustic phonon and impurity scattering will come to
dominate.

In these theories we have assumed a single characteristic dielectric phonon
frequency. Though Hellwarth et al.Hellwarth and Biaggio (1999) shows that
there are good physical reasons for the approximation of the action of
multiple phonon modes with a single effective mode, this still loses structure
in the temperature-dependent response. In the case of lead-halide perovskites,
there are infrared active phonon branches on
$1\text{\,}\mathrm{THz}1.5\text{\,}\mathrm{THz}$ and
$2\text{\,}\mathrm{THz}2.5\text{\,}\mathrm{THz}$. Phonon perturbation theory
calculations show that phonons in this system are highly anharmonicWhalley
_et al._ (2016), with a broad range of energy and small lifetime. This may
invalidate the non-interacting (independent) phonon approximation.

We have shown that the Hellwarth mobility follows an approximately $T^{-0.5}$
power-law dependence, when above the phonon emission threshold. More
temperature-dependent mobility data will help understand the microscopic
processes responsible for charge carrier scattering in these materials.

With the robust codes developed for this work, temperature-dependent large-
polaron mobility calculations are made simple. The material-specific
parameters required from electronic structure calculations are well defined.
Modest computational resources are then required to solve the polaron and
mobility theories. Unlike other models of halide perovskite mobility presented
in the literature, there is no relaxation-time or other empirical parameter,
inferred by analogy to other systems.

We have shown that models predict the experimentally observed temperature-
dependent trend in hybrid halide perovskites, and that the quantitative
agreement suggests that polaron optical-phonon scattering dominates room-
temperature mobility.

###  IV.3 Future work

Several classes of novel electronic materials (transparent conducting oxides,
oxide-based thermoelectrics and organic semiconductors) are expected to be
dominated by polaronic transport. The calculation of temperature-dependent
mobility described in this paper is useful for both design of new materials,
and the interpretation of measurements and characterisation of existing
materials.

Considering the evident utility, it is therefore perhaps surprising that there
has not been more application of these methods. Certainly some of this is due
to the underlying theoretical work being couched in the rather esoteric
language of path integrals. A more practical issue is the lack of an available
computational implementation. Here we developed custom codes in the Julia high
level mathematical languageBezanson  _et al._ (2017). This offered automatic-
differentiation (to calculate the gradients of the objective function of the
Ōsaka free-energies, stabilising the optimisation procedure) and strong
control of numeric errors. The resulting codes appear numerically stable and
general purpose. We provide our codesFrost (2017) to encourage community reuse
of these methods.

With a validated set of codes, we can now easily apply these polaron mobility
theories to polar materials as a natural component of materials design. Though
many idealisations are made in the polaron mobility theories, they offer a
first-principles method to arrive at a phenomenological value of extreme
technical interest. The electronic structure calculations necessary for inputs
to the model are now a standard part of computational materials discovery and
characterisation. The methods and codes are well matched to high-throughput
calculation of mobilities.

Within hybrid perovskites, we plan to combine these models with experimental
measures to characterise the internal structure of the polaron state. As well
as the zero-frequency susceptibility (and mobility), the frequency-dependent
mobility can be evaluated by an additional (numeric) integrationFeynman  _et
al._ (1962), and the optical response of the polaron calculatedDevreese  _et
al._ (1972).

It will be interesting to see whether a more sophisticated method of
calculating electron-phonon matrix elements from a band-structureGiustino
(2017) can complement the use of the simple Fröhlich dielectric $\alpha$
parameter.

The Hellwarth et al.Hellwarth and Biaggio (1999) method of deriving a single
effective polar mode could be supplemented by an explicit calculation of
contributions per mode. This would allow for some of the structure in the low-
temperature mobility (smoothed-out by the choice of a single effective mode)
to return.

###### Acknowledgements.

We thank Piers Barnes for discussion and a critical reading of a draft of this
manuscript. We thank Robert Lovrincic for discussion and bringing to our
attention his prior workSendner  _et al._ (2016). We acknowledge membership of
the UK’s HPC Materials Chemistry Consortium, which is funded by EPSRC grant
EP/F067496. J.M.F. is funded by EPSRC grant EP/K016288/1.

##  V Supplementary Information

Computer codes to replicate all data and figures in this publication are
shared via a repository on GitHubFrost (2017),
https://github.com/jarvist/PolaronMobility.jl.

## References

  * Manser  _et al._ (2016) J. S. Manser, J. A. Christians, and P. V. Kamat, Chemical Reviews 116, 12956 (2016).
  * Feynman (1955) R. P. Feynman, Physical Review 97, 660 (1955).
  * Osaka (1959) Y. Osaka, Progress of Theoretical Physics 22, 437 (1959).
  * Feynman  _et al._ (1962) R. P. Feynman, R. W. Hellwarth, C. K. Iddings, and P. M. Platzman, Physical Review 127, 1004 (1962).
  * Kadanoff (1963) L. P. Kadanoff, Physical Review 130, 1364 (1963).
  * Hellwarth and Biaggio (1999) R. W. Hellwarth and I. Biaggio, Physical Review B 60, 299 (1999).
  * Milot  _et al._ (2015) R. L. Milot, G. E. Eperon, H. J. Snaith, M. B. Johnston, and L. M. Herz, Advanced Functional Materials 25, 6218 (2015).
  * Semonin  _et al._ (2016) O. E. Semonin, G. A. Elbaz, D. B. Straus, T. D. Hull, D. W. Paley, A. M. van der Zande, J. C. Hone, I. Kymissis, C. R. Kagan, X. Roy, and J. S. Owen, The Journal of Physical Chemistry Letters 7, 3510 (2016).
  * Brivio  _et al._ (2014) F. Brivio, K. T. Butler, A. Walsh, and M. van Schilfgaarde, Physical Review B 89 (2014), 10.1103/physrevb.89.155204.
  * Azarhoosh  _et al._ (2016) P. Azarhoosh, S. McKechnie, J. M. Frost, A. Walsh, and M. van Schilfgaarde, APL Materials 4, 091501 (2016).
  * Brivio  _et al._ (2015) F. Brivio, J. M. Frost, J. M. Skelton, A. J. Jackson, O. J. Weber, M. T. Weller, A. R. Goñi, A. M. A. Leguy, P. R. F. Barnes, and A. Walsh, Physical Review B 92 (2015), 10.1103/physrevb.92.144308.
  * Ridley (2013) B. Ridley,  _Quantum Processes in Semiconductors_ , 5th ed. (Oxford University Press, 2013).
  * Rode (1972) D. L. Rode, Physica Status Solidi (b) 53, 245 (1972).
  * Pödör (1983) B. Pödör, physica status solidi (b) 120, 207 (1983).
  * yi Huang and Lambrecht (2013) L. yi Huang and W. R. L. Lambrecht, Physical Review B 88 (2013), 10.1103/physrevb.88.165203.
  * Brandt  _et al._ (2015) R. E. Brandt, V. Stevanović, D. S. Ginley, and T. Buonassisi, MRS Communications 5, 265 (2015).
  * Motta  _et al._ (2015) C. Motta, F. El-Mellouhi, and S. Sanvito, Scientific Reports 5 (2015), 10.1038/srep12746.
  * Zhao  _et al._ (2016) T. Zhao, W. Shi, J. Xi, D. Wang, and Z. Shuai, Scientific Reports 7, 19968 (2016).
  * Filippetti  _et al._ (2016) A. Filippetti, A. Mattoni, C. Caddeo, M. I. Saba, and P. Delugas, Phys. Chem. Chem. Phys. 18, 15352 (2016).
  * Saidaminov  _et al._ (2015) M. I. Saidaminov, A. L. Abdelhady, B. Murali, E. Alarousu, V. M. Burlakov, W. Peng, I. Dursun, L. Wang, Y. He, G. Maculan, A. Goriely, T. Wu, O. F. Mohammed, and O. M. Bakr, Nature Communications 6, 7586 (2015).
  * Hendry  _et al._ (2004) E. Hendry, F. Wang, J. Shan, T. F. Heinz, and M. Bonn, Physical Review B 69 (2004), 10.1103/physrevb.69.081101.
  * Giustino (2017) F. Giustino, Reviews of Modern Physics 89 (2017), 10.1103/revmodphys.89.015003.
  * Whalley  _et al._ (2016) L. D. Whalley, J. M. Skelton, J. M. Frost, and A. Walsh, Physical Review B 94 (2016), 10.1103/physrevb.94.220301.
  * Saidi  _et al._ (2016) W. A. Saidi, S. Poncé, and B. Monserrat, The Journal of Physical Chemistry Letters 7, 5247 (2016).
  * Frohlich (1952) H. Frohlich, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 215, 291 (1952).
  * Gonze and Lee (1997) X. Gonze and C. Lee, Physical Review B 55, 10355 (1997).
  * Peeters and Devreese (1983) F. M. Peeters and J. T. Devreese, physica status solidi (b) 115, 539 (1983).
  * Biaggio  _et al._ (1997) I. Biaggio, R. W. Hellwarth, and J. P. Partanen, Physical Review Letters 78, 891 (1997).
  * Frost (2017) J. M. Frost, https://github.com/jarvist/PolaronMobility.jl (2017).
  * Frost  _et al._ (2014) J. M. Frost, K. T. Butler, F. Brivio, C. H. Hendon, M. van Schilfgaarde, and A. Walsh, Nano Letters 14, 2584 (2014).
  * Kirkwood (1939) J. G. Kirkwood, The Journal of Chemical Physics 7, 911 (1939).
  * Onoda-Yamamuro  _et al._ (1992) N. Onoda-Yamamuro, T. Matsuo, and H. Suga, Journal of Physics and Chemistry of Solids 53, 935 (1992).
  * Ma and Wang (2015) J. Ma and L.-W. Wang, Nano Letters 15, 248 (2015).
  * Neukirch  _et al._ (2016) A. J. Neukirch, W. Nie, J.-C. Blancon, K. Appavoo, H. Tsai, M. Y. Sfeir, C. Katan, L. Pedesseau, J. Even, J. J. Crochet, G. Gupta, A. D. Mohite, and S. Tretiak, Nano Letters 16, 3809 (2016).
  * Bokdam  _et al._ (2016) M. Bokdam, T. Sander, A. Stroppa, S. Picozzi, D. D. Sarma, C. Franchini, and G. Kresse, Scientific Reports 6 (2016), 10.1038/srep28618.
  * Sendner  _et al._ (2016) M. Sendner, P. K. Nayak, D. A. Egger, S. Beck, C. Müller, B. Epding, W. Kowalsky, L. Kronik, H. J. Snaith, A. Pucci, and R. Lovrinčić, Mater. Horiz. 3, 613 (2016).
  * Schultz (1959) T. D. Schultz, Physical Review 116, 526 (1959).
  * Price  _et al._ (2015) M. B. Price, J. Butkus, T. C. Jellicoe, A. Sadhanala, A. Briane, J. E. Halpert, K. Broch, J. M. Hodgkiss, R. H. Friend, and F. Deschler, Nature Communications 6, 8420 (2015).
  * Yang  _et al._ (2015) Y. Yang, D. P. Ostrowski, R. M. France, K. Zhu, J. van de Lagemaat, J. M. Luther, and M. C. Beard, Nature Photonics 10, 53 (2015).
  * Leguy  _et al._ (2016) A. M. A. Leguy, A. R. Goñi, J. M. Frost, J. Skelton, F. Brivio, X. Rodríguez-Martínez, O. J. Weber, A. Pallipurath, M. I. Alonso, M. Campoy-Quiles, M. T. Weller, J. Nelson, A. Walsh, and P. R. F. Barnes, Phys. Chem. Chem. Phys. 18, 27051 (2016).
  * Chung  _et al._ (2012) I. Chung, J.-H. Song, J. Im, J. Androulakis, C. D. Malliakas, H. Li, A. J. Freeman, J. T. Kenney, and M. G. Kanatzidis, Journal of the American Chemical Society 134, 8579 (2012).
  * Madsen and Singh (2006) G. K. Madsen and D. J. Singh, Computer Physics Communications 175, 67 (2006).
  * Zhang  _et al._ (2017) M. Zhang, X. Zhang, L.-Y. Huang, H.-Q. Lin, and G. Lu, “Charge transport in hybrid halide perovskites,” (2017), arXiv:1708.09469 .
  * Bezanson  _et al._ (2017) J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, SIAM Review 59, 65 (2017).
  * Devreese  _et al._ (1972) J. Devreese, J. D. Sitter, and M. Goovaerts, Physical Review B 5, 2367 (1972).