# Thermal transport signatures of the excitonic transition and associated
phonon softening in the layered chalcogenide Ta2NiSe5

Yuan-Shan Zhang  Max Planck Institute for Solid State Research, 70569
Stuttgart, Germany  Jan A. N. Bruin<EMAIL_ADDRESS>Max Planck Institute
for Solid State Research, 70569 Stuttgart, Germany  Yosuke Matsumoto  Max
Planck Institute for Solid State Research, 70569 Stuttgart, Germany  Masahiko
Isobe  Max Planck Institute for Solid State Research, 70569 Stuttgart, Germany
Hidenori Takagi<EMAIL_ADDRESS>Max Planck Institute for Solid State
Research, 70569 Stuttgart, Germany  Institute for Functional Matter and
Quantum Technologies, University of Stuttgart, 70569 Stuttgart, Germany
Department of Physics, University of Tokyo, 113-0033 Tokyo, Japan

###### Abstract

The layered compound $\mathrm{Ta_{2}NiSe_{5}}$ is a quasi-one-dimensional and
narrow-gap semiconductor, which is proposed to undergo a transition to an
excitonic insulator at $T_{\mathrm{c}}=326$ K. We found a clear anomaly at
$T_{\mathrm{c}}$ in the in-plane thermal conductivities both parallel
($\parallel$ $a$) and perpendicular ($\parallel$ $c$) to the one-dimensional
chains, $\kappa_{\mathrm{a}}$ and $\kappa_{\mathrm{c}}$. While
$\kappa_{\mathrm{a}}$ shows a rapid decrease below $T_{\mathrm{c}}$,
$\kappa_{\mathrm{c}}$ shows a pronounced V-shaped suppression centered at
$T_{\mathrm{c}}$. We argue that the decrease of $\kappa_{\mathrm{a}}$
represents the suppression of the quasiparticle contribution below
$T_{\mathrm{c}}$ due to excitonic condensation. On the other hand, the
V-shaped suppression of $\kappa_{\mathrm{c}}$ comes from the enhanced phonon
scattering by soft phonons associated with the monoclinic transition with
momentum $\mathbf{q}\parallel c$. The continued suppression of
$\kappa_{\mathrm{c}}$ up to an extremely high temperature above
$T_{\mathrm{c}}$ suggests the persistence of phonon softening likely coupled
to electronic, presumably excitonic, fluctuations.

The excitonic insulator (EI), in which bound electron-hole pairs form a
condensate below a transition temperature $T_{\mathrm{c}}$, has been proposed
theoretically more than 50 years ago [1], but its realization in existing
materials remains an active subject of investigation. The excitonic transition
is anticipated in a nearly-zero-gap semiconductor when the exciton binding
energy $E_{\mathrm{b}}$ exceeds the one-electron gap $E_{\mathrm{g}}$, or in a
semimetal when the attractive interaction between electrons and holes is not
screened out well. The realization in bulk materials has been discussed, for
example, in $\mathrm{TmSe_{0.45}Te_{0.55}}$ [2] and $1T$-$\mathrm{TiSe_{2}}$
[3], but the evidence for excitonic condensation in these compounds has been
far from convincing. The former is a $4f^{13}$ Mott-like system, not a band
semiconductor/semimetal [4], and the latter remains metallic and may suffer
from strong screening of electron-hole interactions even below the proposed
excitonic transition temperature [5].

Very recently, $\mathrm{Ta_{2}NiSe_{5}}$ has emerged as a leading candidate
for an excitonic insulator [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18].
It is an almost-zero-energy- and direct-gap semiconductor with Ta $5d$
conduction bands and Ni $3d$ valence bands [8]. It shows a transition to an
insulator at $T_{\mathrm{c}}=326$ K, which was proposed to be an excitonic
transition. The crystal structure of $\mathrm{Ta_{2}NiSe_{5}}$ is layered, and
as shown in Fig. 1(a), each layer consists of an array of Ta and Ni chains
running along the $a$ axis, which renders the system electronically quasi-one-
dimensional and gives rise to spatially isolated electron and hole chains. The
one dimensionality is known to enhance the exciton binding energy, which makes
$\mathrm{Ta_{2}NiSe_{5}}$ an ideal platform for the search for an excitonic
insulator. Well below $T_{\mathrm{c}}$, an optical gap of 0.16 eV develops,
which is close to the exciton binding energy of 0.25 eV experimentally
observed in the optical conductivity spectrum of the sister compound
$\mathrm{Ta_{2}NiS_{5}}$. $\mathrm{Ta_{2}NiS_{5}}$ has a much larger gap of
0.6 eV and exhibits no signatures of a phase transition [8, 9]. This strongly
suggests the excitonic origin of the insulating gap in
$\mathrm{Ta_{2}NiSe_{5}}$ below $T_{\mathrm{c}}$. As a function of the one-
electron gap $E_{\mathrm{g}}$ controlled by pressure and chemical
substitution, the transition temperature $T_{\mathrm{c}}$ shows a dome-shaped
behavior centered at $E_{\mathrm{g}}\approx 0$, corresponding to pure
$\mathrm{Ta_{2}NiSe_{5}}$ at ambient pressure [8]. This is fully consistent
with what is expected for the canonical excitonic insulator [1].

The aforementioned experimental results clearly point to the excitonic
character of the transition. Nevertheless, it has become increasingly clear
that the transition cannot be purely excitonic and that the actual situation
is more complicated [19, 20, 21], which makes the physics of the putative
excitonic transition even more challenging and attractive. The putative
excitonic transition at $T_{\mathrm{c}}=326$ K is accompanied by a structural
phase transition from an orthorhombic to a monoclinic structure [22, 23]. In
the monoclinic phase, the $\mathrm{TaSe_{2}}$ chain-ladder unit distorts, such
that the Ta chain and the surrounding Se ladders move along the $a$ axis, but
in opposite directions [24]. The movement of two neighboring
$\mathrm{TaSe_{2}}$ chain-ladder units surrounding the Ni chain is asymmetric
with respect to the Ni chain. This shear-type distortion with
antiferroelectric-type displacement of the Ta and the Se chains does not bring
about any superstructure (i.e., $\mathbf{q}=0$) but allows for the
hybridization of Ta conduction bands and Ni valence bands around the $\Gamma$
point. The resultant hybridization gap could account for at least a part of
the insulating gap and should be one of the important ingredients of the
transition at 326 K. A symmetry analysis demonstrated that the excitonic order
parameter must be accompanied by the monoclinic distortion and hence the
admixture of the excitonic gap and the hybridization gap is unavoidable [25].

In this study, we measure the in-plane anisotropic thermal conductivities
$\kappa_{\mathrm{a}}(T)$ along the chains ($\parallel$ $a$) and
$\kappa_{\mathrm{c}}(T)$ perpendicular to the chains ($\parallel$ $c$) of
$\mathrm{Ta_{2}NiSe_{5}}$ single crystals, in order to unveil the
quasiparticle thermal conductivity of the putative excitonic insulator and the
phonon transport around the monoclinic transition. In excitonic insulators,
quasiparticle transport analogous to Bardeen-Cooper-Schrieffer (BCS)
superconductors may be expected [1, 26]. Thermal conductivity is an especially
useful probe of the quasiparticle transport in an excitonic insulator as the
cancelation between the electron and hole quasiparticles does not take place.
In addition, thermal conductivity has a large phonon contribution and is
sensitive to soft phonons associated with the monoclinic transition, which may
give us hints to decipher the role of the lattice. We show that
$\kappa_{\mathrm{a}}(T)$ along the chains is rapidly suppressed below
$T_{\mathrm{c}}$ and that it comes from the decrease in number of
quasiparticles (electrons and holes) below $T_{\mathrm{c}}$, which is expected
for a canonical excitonic insulator [27]. In contrast, $\kappa_{\mathrm{c}}$
perpendicular to the chains shows a V-shaped suppression due to enhanced
scattering by soft phonons associated with the orthorhombic-to-monoclinic
transition with momentum $\mathbf{q}$ perpendicular to the chains ($\parallel$
$c$). The strong suppression of $\kappa_{\mathrm{c}}$ persists up to more than
100 K above $T_{\mathrm{c}}$, which suggests a robust phonon softening likely
coupled to excitonic (and hybridization-gap) fluctuations.

Needle-shaped single crystals of $\mathrm{Ta_{2}NiSe_{5}}$ with typical size
of $\mathrm{4~{}mm\times 1~{}mm\times 0.1~{}mm}$ were grown by chemical vapor
transport [8, 28], with their longest dimension along the Ta/Ni chains ($a$
axis) and their shortest dimension along the interlayer $b$ axis. Rectangular-
shaped samples with their long axis along the $a$ and $c$ axes were cut out
from the single crystals. The electrical resistivity of the samples was
measured by a conventional four-probe technique, while thermal conductivity
measurements along the same axes were performed via a “heat-pipe” method [29,
30]. The temperature gradients were detected with a differential type-E
thermocouple. A transition was observed at $T_{\mathrm{c}}=326$ K as a clear
kink in the electrical resistivity along the chain direction
($\rho_{\mathrm{a}}$), as reported previously [8, 28].

Figure 1: (a) Crystal structure of $\mathrm{Ta_{2}NiSe_{5}}$ showing the stack
of layers along the $b$ axis (left) and a single layer in the $a$-$c$ plane
consisting of an array of two chains of Ta octahedrally coordinated with Se
and one chain of Ni in between, tetrahedrally coordinated with Se (right)
[31]. (b) Thermal conductivities $\kappa_{\mathrm{a}}$ with heat current
parallel to the chains ($\parallel$ $a$; red) and $\kappa_{\mathrm{c}}$ with
heat current perpendicular to the chains ($\parallel$ $c$; blue) as a function
of temperature. The dashed black line indicates $T_{\mathrm{c}}=326$ K. Inset:
enlarged view of $\kappa_{\mathrm{a}}$ and $\kappa_{\mathrm{c}}$ near
$T_{\mathrm{c}}$.

As shown in Fig. 1(b), $\kappa_{\mathrm{a}}(T)$ along the chains ($\parallel$
$a$) is a factor of 3-4 larger than $\kappa_{\mathrm{c}}$ perpendicular to the
chains ($\parallel$ $c$). With increasing temperature from 10 K, both
$\kappa_{\mathrm{a}}$ and $\kappa_{\mathrm{c}}$ increase and show a peak
around 30-40 K, followed by a gradual decrease that continues above room
temperature. This behavior is common for a crystalline solid with phonon-
dominated thermal conductivity [32, 33].

At the transition temperature $T_{\mathrm{c}}=326$ K, we observe a clear
anomaly both in $\kappa_{\mathrm{a}}(T)$ and $\kappa_{\mathrm{c}}(T)$. The
anomaly at $T_{\mathrm{c}}$, however, shows up in a contrasted manner between
the two directions $a$ and $c$, which can be seen more clearly in the enlarged
plot around $T_{\mathrm{c}}$ in the inset of Fig. 1(b).
$\kappa_{\mathrm{a}}(T)$ parallel to the chains shows a shoulder structure
with a well-defined kink at $T_{\mathrm{c}}$, followed by an abrupt drop below
$T_{\mathrm{c}}$. Upon further lowering of the temperature,
$\kappa_{\mathrm{a}}(T)$ recovers and begins increasing below $\sim$315 K. In
contrast, $\kappa_{\mathrm{c}}(T)$ shows a discontinuous change of slope at
$T_{\mathrm{c}}$ from a weak increase to a very rapid increase with decreasing
temperature, which can be viewed as a V-shaped suppression from the $1/T$-like
increase expected for phonon-dominant, high-temperature thermal transport
[32]. The contrast is even clearer in the temperature derivatives,
$\partial\kappa_{\mathrm{a}}/\partial T$ and
$\partial\kappa_{\mathrm{c}}/\partial T$ [28]. We argue that the contrasted
behaviors of the anomaly at $T_{\mathrm{c}}$ between $\kappa_{\mathrm{a}}(T)$
and $\kappa_{\mathrm{c}}(T)$ come from the predominant electronic and phonon
contributions in the temperature dependence of $\kappa_{\mathrm{a}}(T)$ and
$\kappa_{\mathrm{c}}(T)$, respectively, near $T_{\mathrm{c}}$, which allows us
to discuss the critical behaviors of quasiparticles (electrons and holes) and
phonons around $T_{\mathrm{c}}$ separately.

The measured thermal conductivity $\kappa$ comprises of contributions from
electrons, $\kappa^{(\mathrm{e})}$, and from phonons,
$\kappa^{(\mathrm{ph})}$. Around $T_{\mathrm{c}}$, $\kappa^{(\mathrm{e})}$ can
be estimated from the electrical conductivity $\sigma$ using the Wiedemann-
Franz (W-F) law. Strictly speaking, the W-F law is valid only at $T=0$ K.
Nevertheless, it is also a realistic approximation at finite temperatures,
where large-angle scattering processes dominate [32, 33].
$\kappa^{(\mathrm{ph})}$ can be estimated as the remnant
$\kappa-\kappa^{(\mathrm{e})}$. As shown in Fig. 2(a), the electrical
conductivity along the chains, $\sigma_{\mathrm{a}}$, is an order of magnitude
larger than that perpendicular to the chains, $\sigma_{\mathrm{c}}$,
reflecting the quasi-one-dimensional band structure of
$\mathrm{Ta_{2}NiSe_{5}}$. In addition, the anomaly at $T_{\mathrm{c}}$ is
much more pronounced in $\sigma_{\mathrm{a}}(T)$ as compared with
$\sigma_{\mathrm{c}}(T)$. A much larger electronic contribution to the
temperature dependence of $\kappa_{\mathrm{a}}(T)$ than
$\kappa_{\mathrm{c}}(T)$ is thus expected.

As shown in Fig. 2(b), the electronic thermal conductivity along the chains
$\kappa_{\mathrm{a}}^{(\mathrm{e})}(T)$, estimated from $\sigma_{\mathrm{a}}$,
is only $\sim$5$\%$ of the total $\kappa_{\mathrm{a}}(T)$, but fully accounts
for the shoulder structure in $\kappa_{\mathrm{a}}(T)$ with the rapid decrease
below $T_{\mathrm{c}}$. On the other hand, the phonon thermal conductivity
$\kappa_{\mathrm{a}}^{(\mathrm{ph})}(T)$, estimated by subtracting
$\kappa_{\mathrm{a}}^{(\mathrm{e})}(T)$ from $\kappa_{\mathrm{a}}(T)$, is
surprisingly featureless and shows a monotonic increase through
$T_{\mathrm{c}}$, indicating that the anomaly at $T_{\mathrm{c}}$ in
$\kappa_{\mathrm{a}}(T)$ is purely electronic in origin. The rapid decrease
below $T_{\mathrm{c}}$ therefore reflects the decrease in number of
quasiparticles, electrons and holes, in analogy with the prototypical behavior
of thermal conductivity in a superconductor [34, 35]. We therefore conclude
that the quasiparticle transport of a putative excitonic insulator is captured
by thermal transport.

Figure 2: Separation of electron and phonon contributions to the thermal
conductivities $\kappa_{\mathrm{a}}$ with heat current parallel to the chains
($\parallel$ $a$; red) and $\kappa_{\mathrm{c}}$ with heat current
perpendicular to the chains ($\parallel$ $c$; blue) in Fig. 1, using the
Wiedemann-Franz law. (a) The temperature-dependent electrical conductivities,
$\sigma_{\mathrm{a}}(T)$ (red) parallel to the chains and
$\sigma_{\mathrm{c}}(T)$ (blue) perpendicular to the chains. (b) The total,
electron, and phonon thermal conductivities along the $a$ axis, parallel to
the chains, denoted by $\kappa_{\mathrm{a}}$ (circles),
$\kappa_{\mathrm{a}}^{(\mathrm{e})}$ (broken), and
$\kappa_{\mathrm{a}}^{(\mathrm{ph})}$ (solid), respectively. (c) The total,
electron, and phonon thermal conductivities along the $c$ axis, perpendicular
to the chains, denoted $\kappa_{\mathrm{c}}$ (circles),
$\kappa_{\mathrm{c}}^{(\mathrm{e})}$ (broken), and
$\kappa_{\mathrm{c}}^{(\mathrm{ph})}$ (solid), respectively.

In superconductors, the thermal conductivity often does not show a rapid
decrease right below $T_{\mathrm{c}}$, the prototypical behavior, but instead
shows an increase. This can be ascribed to the suppression of scattering of
electrons and phonons by electrons due to the opening of an electronic gap
below $T_{\mathrm{c}}$. In certain strongly correlated superconductors with
anisotropic pairing, such as high-$T_{\mathrm{c}}$ cuprates [36], the
enhancement of quasiparticle lifetime due to the suppression of strong
electron-electron scattering overcomes the decrease of quasiparticle density
and thereby results in an increase of $\kappa$ right below $T_{\mathrm{c}}$.
In certain conventional electron-phonon superconductors, the increase of
phonon mean free path due to the suppression of strong scattering of phonons
by electrons may also result in the enhancement of phonon thermal conductivity
below $T_{\mathrm{c}}$, which brings about the increase of $\kappa$ right
below $T_{\mathrm{c}}$ [37]. The observation of a rapid decrease in
$\kappa_{\mathrm{a}}(T)$ below $T_{\mathrm{c}}$ in $\mathrm{Ta_{2}NiSe_{5}}$
thus indicates that the scattering of phonons and electrons with momentum
parallel to the chains by electrons should be reasonably weak as compared with
other dominant scattering mechanisms, such as phonon-phonon scattering for
phonon transport.

As $\sigma_{\mathrm{c}}$ is one order of magnitude smaller and exhibits much
weaker structure at $T_{\mathrm{c}}$ than $\sigma_{\mathrm{a}}$, the
electronic thermal conductivity $\kappa_{\mathrm{c}}^{\mathrm{(e)}}(T)$
perpendicular to the chain direction is as small as $\sim$1$\%$ of the total
$\kappa_{\mathrm{c}}(T)$ in magnitude, and, more importantly, does not
contribute appreciably to the kink-like anomaly at $T_{\mathrm{c}}$ in the
temperature dependence of $\kappa_{\mathrm{c}}(T)$, as seen in Fig. 2(c). In
contrast to $\kappa_{\mathrm{a}}(T)$ along the chains,
$\kappa_{\mathrm{c}}(T)$ around $T_{\mathrm{c}}$ is dominated by
$\kappa_{\mathrm{c}}^{(\mathrm{ph})}(T)$ and almost purely phononic in origin.
The comparison of $\kappa_{\mathrm{c}}^{(\mathrm{ph})}(T)$ and
$\kappa_{\mathrm{a}}^{(\mathrm{ph})}(T)$ reveals a substantial anisotropy not
only in electron transport but also in phonon transport between the directions
parallel and perpendicular to the chains. As clearly seen in Fig. 3, while a
$1/T$-like monotonous decrease with featureless behavior around
$T_{\mathrm{c}}$ is observed in $\kappa_{\mathrm{a}}^{(\mathrm{ph})}(T)$, a
clear V-shaped suppression from the $1/T$-like decrease centered at
$T_{\mathrm{c}}$ is seen in $\kappa_{\mathrm{c}}^{(\mathrm{ph})}(T)$. The
suppression is broad in temperature range and appears to be asymmetric above
and below $T_{\mathrm{c}}$. While $\kappa_{\mathrm{c}}^{(\mathrm{ph})}(T)$
shows a rapid increase and crossover to $1/T$-like behavior with decreasing
temperatures right below $T_{\mathrm{c}}$, it is almost flat up to at least
400 K with increasing temperature from $T_{\mathrm{c}}$, implying the
persistence of strong suppression over a wide temperature range above
$T_{\mathrm{c}}$.

Phonon thermal conductivity is generally expressed as
$\kappa^{\mathrm{(ph)}}=\frac{1}{3}C^{\mathrm{(ph)}}v^{\mathrm{(ph)}}l^{(\mathrm{ph})}$,
where $C^{\mathrm{(ph)}}$, $v^{\mathrm{(ph)}}$, and $l^{(\mathrm{ph})}$ are
the specific heat, the velocity, and the mean free path of phonons [32]. The
temperature dependence of the total heat capacity $C$ for
$\mathrm{Ta_{2}NiSe_{5}}$ around $T_{\mathrm{c}}=326$ K is rather small and
the small anomaly at $T_{\mathrm{c}}$ is highly likely to be purely electronic
in origin [8]. The temperature dependence of $C^{\mathrm{(ph)}}$ can thus be
ignored in discussing the origin of the V-shaped suppression of
$\kappa_{\mathrm{c}}^{\mathrm{(ph)}}(T)$. The suppression of
$\kappa_{\mathrm{c}}^{\mathrm{(ph)}}(T)$ should originate instead from the
reduction of $v^{\mathrm{(ph)}}$ and/or $l^{(\mathrm{ph})}$ around
$T_{\mathrm{c}}$. As it occurs only for phonons with momentum
$\mathbf{q}\parallel c$, it is natural to ascribe its origin to the reported
softening of the acoustic shear-mode phonon with $\mathbf{q}\parallel c$,
which corresponds to the monoclinic distortion of lattice frozen at the
excitonic transition $T_{\mathrm{c}}=326$ K [24]. As pronounced softening
occurs only for one specific acoustic phonon mode, the reduction of
$v^{\mathrm{(ph)}}$ is unlikely to account for the large suppression of
$\kappa_{\mathrm{c}}^{\mathrm{(ph)}}(T)$ at $T_{\mathrm{c}}$ [Fig. 2 (c)]. The
primary origin of the suppression of $\kappa_{\mathrm{c}}^{\mathrm{(ph)}}(T)$
is highly likely the increased scattering of phonons, namely the reduction of
$l^{(\mathrm{ph})}$, by the soft acoustic shear-mode phonon. The strong
suppression up to at least 400 K above $T_{\mathrm{c}}$ should then imply the
robust softening of the shear mode acoustic phonon above $T_{\mathrm{c}}$,
which is indeed quite consistent with the reported phonon dispersion
relationship of $\mathrm{Ta_{2}NiSe_{5}}$ measured by inelastic X-ray
scattering up to 400 K [24].

Figure 3: Comparison of temperature-dependent thermal conductivities
$\kappa(T)$ for $\mathrm{Ta_{2}NiSe_{5}}$ ($\parallel$ $a$ red; $\parallel$
$c$ blue) with typical order-disorder-type ferroelectrics, $\mathrm{NH_{4}Cl}$
(reproduced from [38]) and $\mathrm{KH_{2}PO_{4}}$ (reproduced from [39]). The
solid black line is a guide line for the $1/T$ power law.

We should point out that the asymmetric V-shaped suppression of
$\kappa_{\mathrm{c}}^{\mathrm{(ph)}}(T)$, more robust above $T_{\mathrm{c}}$,
is analogous to those observed around the order-disorder-type ferroelectric
transition in $\mathrm{KH_{2}PO_{4}}$ and $\mathrm{NH_{4}Cl}$, as shown in
Fig. 3. The order-disorder-type ferroelectrics have extended ferroelectric
fluctuations above $T_{\mathrm{c}}$, where well-developed but fluctuating
local dipoles give rise to soft local phonons and enhanced phonon scattering
over a wide temperature range, and hence suppress $\kappa(T)$ strongly over a
wide temperature range above $T_{\mathrm{c}}$ [40, 39, 38, 41]. It is tempting
to infer that, in analogy with order-disorder type ferroelectrics, large and
local excitonic and/or hybridization gap ($\mathbf{q}=0$ charge density wave)
fluctuations may be present above $T_{\mathrm{c}}$ in
$\mathrm{Ta_{2}NiSe_{5}}$. The electronic fluctuations are inherently and
strongly coupled to the shear mode phonon, which may give rise to the extended
phonon softening above $T_{c}$ and may suppress
$\kappa_{\mathrm{c}}^{\mathrm{(ph)}}(T)$ over a wide temperature range. Such
giant excitonic and/or hybridization gap fluctuations may be inferred from the
temperature dependence of the bulk magnetic susceptibility $\chi(T)$ of
$\mathrm{Ta_{2}NiSe_{5}}$ [22], where $\chi(T)$ at temperatures well above
$T_{\mathrm{c}}$ shows a pronounced downturn upon lower temperatures towards
$T_{\mathrm{c}}$, indicative of such fluctuations.

In summary, we have discovered an extremely anisotropic behavior of the
thermal conductivity around the transition from an almost-zero gap
semiconductor to an insulator in a putative excitonic insulator, the layered
chalcogenide $\mathrm{Ta_{2}NiSe_{5}}$. The distinct temperature-dependent in-
plane $\kappa_{\mathrm{a}}(T)$ and $\kappa_{\mathrm{c}}(T)$, parallel and
perpendicular to the Ta and Ni chains, respectively, allow us to unveil the
behaviors of quasiparticles (electrons and holes) and phonons around the
transition temperature $T_{\mathrm{c}}=326$ K. The rapid suppression in the
number of quasiparticles below $T_{\mathrm{c}}$, analogous to superconductors,
is captured in $\kappa_{\mathrm{a}}(T)$, due to the quasi-one-dimensional
electronic structure and the $\mathbf{q}\parallel a$ phonons being insensitive
to the transition. The soft shear-mode acoustic phonons with momentum
$\mathbf{q}\parallel c$, which corresponds to the monoclinic structural
distortion below $T_{\mathrm{c}}$, manifest themselves only in
$\kappa_{\mathrm{c}}(T)$ perpendicular to the chains as a pronounced
suppression of dominant phonon thermal transport. The suppression of
$\kappa_{\mathrm{c}}$ extends up to high temperatures well above
$T_{\mathrm{c}}$, indicating the persistence of phonon softening over a wide
temperature range. This may suggest the presence of giant excitonic and/or
hybridization gap fluctuations above $T_{\mathrm{c}}$ and the strong electron-
lattice coupling for the monoclinic distortion in $\mathrm{Ta_{2}NiSe_{5}}$.

###### Acknowledgements.

This work is supported in part by the Alexander von Humboldt foundation. We
thank L. Dorner-Finkbeiner and C. Bush for experimental support in crystal
growth and characterization, A. Bangura for suggesting the “heat-pipe”
technique, and D. Huang for useful discussions.

## References

  * Jérome  _et al._ [1967] D. Jérome, T. M. Rice, and W. Kohn, Phys. Rev. 158, 462 (1967).
  * Wachter  _et al._ [2004] P. Wachter, B. Bucher, and J. Malar, Phys. Rev. B 69, 094502 (2004).
  * Kogar  _et al._ [2017] A. Kogar, M. S. Rak, S. Vig, A. A. H., F. Flicker, Y. I. Joe, L. Venema, G. J. MacDougall, T. C. Chiang, and E. Fradkin, Science 358, 1314 (2017).
  * Wachter [2018] P. Wachter, Adv. Mater. Phys. Chem. 8, 120 (2018).
  * Di Salvo  _et al._ [1976] F. J. Di Salvo, D. E. Moncton, and J. V. Waszczak, Phys. Rev. B 14, 4321 (1976).
  * Wakisaka  _et al._ [2009] Y. Wakisaka, T. Sudayama, K. Takubo, T. Mizokawa, M. Arita, H. Namatame, M. Taniguchi, N. Katayama, M. Nohara, and H. Takagi, Phys. Rev. Lett. 103, 026402 (2009).
  * Kim  _et al._ [2016] S. Kim, Y. Kim, C. Kang, E. An, H. K. Kim, M. J. Eom, M. Lee, C. Park, T. Kim, H. Choi, B. Min, and J. Kim, ACS Nano 10, 8888 (2016).
  * Lu  _et al._ [2017] Y. F. Lu, H. Kono, T. I. Larkin, A. W. Rost, T. Takayama, A. V. Boris, B. Keimer, and H. Takagi, Nat. Commun. 8, 14408 (2017).
  * Larkin  _et al._ [2017] T. I. Larkin, A. N. Yaresko, D. Pröpper, K. A. Kikoin, Y. F. Lu, T. Takayama, Y.-L. Mathis, A. W. Rost, H. Takagi, B. Keimer, and A. V. Boris, Phys. Rev. B 95, 195144 (2017).
  * Mor  _et al._ [2017] S. Mor, M. Herzog, D. Golež, P. Werner, M. Eckstein, N. Katayama, M. Nohara, H. Takagi, T. Mizokawa, C. Monney, and J. Stähler, Phys. Rev. Lett. 119, 086401 (2017).
  * Werdehausen  _et al._ [2018] D. Werdehausen, T. Takayama, M. Höppner, G. Albrecht, A. W. Rost, Y. Lu, D. Manske, H. Takagi, and S. Kaiser, Sci. Adv. 4, eaap8652 (2018).
  * Seo  _et al._ [2018] Y. Seo, M. Eom, J. Kim, C. Kang, B. I. Min, and J. Hwang, Sci. Rep. 8, 11961 (2018).
  * Okazaki  _et al._ [2018] K. Okazaki, Y. Ogawa, T. Suzuki, T. Yamamoto, T. Someya, S. Michimae, M. Watanabe, Y. Lu, M. Nohara, H. Takagi, N. Katayama, H. Sawa, M. Fujisawa, T. Kanai, N. Ishii, J. Itatani, T. Mizokawa, and S. Shin, Nat. Commun. 9, 4322 (2018).
  * Fukutani  _et al._ [2019] K. Fukutani, R. Stania, J. Jung, E. F. Schwier, K. Shimada, C. I. Kwon, J. S. Kim, and H. W. Yeom, Phys. Rev. Lett. 123, 206401 (2019).
  * Chen  _et al._ [2020] L. Chen, T. T. Han, C. Cai, Z. G. Wang, Y. D. Wang, Z. M. Xin, and Y. Zhang, Phys. Rev. B 102, 161116(R) (2020).
  * [16] P. A. Volkov, M. Ye, H. Lohani, I. Feldman, A. Kanigel, K. Haule, and G. Blumberg,  arXiv:2007.07344 [cond-mat.str-el] .
  * [17] K. Kim, H. Kim, J. Kim, C. Kwon, J. Kim, and B. Kim,  arXiv:2007.08212 [cond-mat.str-el] .
  * [18] Q. He, X. Que, L. Zhou, M. Isobe, D. Huang, and H. Takagi, arXiv:2012.08395 [cond-mat.str-el] .
  * Watson  _et al._ [2020] M. D. Watson, I. Marković, E. A. Morales, P. Le Fèvre, M. Merz, A. A. Haghighirad, and P. D. C. King, Phys. Rev. Res. 2, 013236 (2020).
  * Subedi [2020] A. Subedi, Phys. Rev. Mater. 4, 083601 (2020).
  * [21] E. Baldini, A. Zong, D. Choi, C. Lee, M. H. Michael, L. Windgaetter, I. I. Mazin, S. Latini, D. Azoury, B. Lv, A. Kogar, Y. Wang, Y. Lu, T. Takayama, H. Takagi, A. J. Millis, A. Rubio, E. Demler, and N. Gedik, arXiv:2007.02909 [cond-mat.str-el] .
  * Di Salvo  _et al._ [1986] F. J. Di Salvo, C. H. Chen, R. M. Fleming, J. V. Waszczak, R. G. Dunn, S. A. Sunshine, and J. A. Ibers, J. Less Common Met. 116, 51 (1986).
  * Nakano  _et al._ [2018a] A. Nakano, K. Sugawara, S. Tamura, N. Katayama, K. Matsubayashi, T. Okada, Y. Uwatoko, K. Munakata, A. Nakao, and H. Sagayama, IUCrJ 5, 158 (2018a).
  * Nakano  _et al._ [2018b] A. Nakano, T. Hasegawa, S. Tamura, N. Katayama, S. Tsutsui, and H. Sawa, Phys. Rev. B 98, 045139 (2018b).
  * Mazza  _et al._ [2020] G. Mazza, M. Rösner, L. Windgätter, S. Latini, H. Hübener, A. J. Millis, A. Rubio, and A. Georges, Phys. Rev. Lett. 124, 197601 (2020).
  * Tinkham [2004] M. Tinkham,  _Introduction to Superconductivity_ (Courier Corporation, 2004).
  * Zittartz [1968] J. Zittartz, Phys. Rev. 165, 612 (1968).
  * [28] See Supplemental Material at [URL will be inserted by publisher] for further information regarding crystal preparation, electrical resistivity, temperature derivatives of thermal conductivities, phonon heat capacity and phonon thermal diffusivity.
  * Allen  _et al._ [1994] P. B. Allen, X. Du, L. Mihaly, and L. Forro, Phys. Rev. B 49, 9073 (1994).
  * Wakeham  _et al._ [2011] N. Wakeham, A. F. Bangura, X. Xu, J.-F. Mercure, M. Greenblatt, and N. E. Hussey, Nat. Commun. 2, 1 (2011).
  * Momma and Izumi [2011] K. Momma and F. Izumi, J. App. Cry. 44, 1272 (2011).
  * Berman [1976] R. Berman,  _Thermal conduction in solids_ (Clarendon Press, 1976).
  * Tritt [2005] T. M. Tritt,  _Thermal conductivity: theory, properties, and applications_ (Springer Science & Business Media, 2005).
  * Satterthwaite [1962] C. B. Satterthwaite, Phys. Rev. 125, 873 (1962).
  * Ambegaokar and Woo [1965] V. Ambegaokar and J. Woo, Phys. Rev. 139, A1818 (1965).
  * Yu  _et al._ [1992] R. C. Yu, M. B. Salamon, J. P. Lu, and W. C. Lee, Phys. Rev. Lett. 69, 1431 (1992).
  * Carlson and Satterthwaite [1970] J. R. Carlson and C. B. Satterthwaite, Phys. Rev. Lett. 24, 461 (1970).
  * Bausch and Waidelich [1969] W. Bausch and W. Waidelich, Phys. Lett. A 30, 190 (1969).
  * Suemune [1967] Y. Suemune, J. Phys. Soc. Jap. 22, 735 (1967).
  * Lines and Glass [2001] M. E. Lines and A. M. Glass,  _Principles and applications of ferroelectrics and related materials_ (Oxford University Press, 2001).
  * Strukov and Belov [1994] B. A. Strukov and A. A. Belov, Phase Transit. 51, 175 (1994).


# Exploring the difficulty of estimating win probability: a simulation study

Ryan S. Brill and Abraham J. Wyner Graduate Group in Applied Mathematics and
Computational Science, University of Pennsylvania. Correspondence to:
<EMAIL_ADDRESS>of Statistics and Data Science, The Wharton
School, University of Pennsylvania

###### Abstract

Estimating win probability is one of the classic modeling tasks of sports
analytics. Many widely used win probability estimators are statistical win
probability models, which fit the relationship between a binary win/loss
outcome variable and certain game-state variables using data-driven regression
or machine learning approaches. To illustrate just how difficult it is to
accurately fit a statistical win probability model from noisy and highly
correlated observational data, in this paper we conduct a simulation study. We
create a simplified random walk version of football in which true win
probability at each game-state is known, and we see how well a model recovers
it. We find that the dependence structure of observational play-by-play data
substantially inflates the bias and variance of estimators and lowers the
effective sample size. This makes it essential to quantify uncertainty in win
probability estimates, but typical bootstrapped confidence intervals are too
narrow and don’t achieve nominal coverage. Hence, we introduce a novel method,
the fractional bootstrap, to calibrate these intervals to achieve adequate
coverage.

##  1 Introduction

Win probability (${\mathsf{WP}}$) as a function of game-state is a canonical
value function in sports analytics. In-game win probability is the crux of
strategic decision making – make the decision that maximizes win probability –
and is central to live betting on game outcomes. Fourth-down decision making
is a prime example: choose between a conversion attempt, field goal attempt,
and punt according to win probability (Brill et al.,, 2024). Win probability,
however, is not an observable quantity. Rather, it is defined by a model that
must be estimated from data. Estimating win probability is one of the classic
modeling tasks of sports analytics.

Win probability estimates arise broadly from one of two classes of models,
probabilistic state-space models or statistical models. On one hand, state-
space models simplify the game of football into a series of transitions
between game-states. Transition probabilities are estimated from play-level
data and are then propagated into win probability by simulating games. When
implemented correctly, these models are sensible ways to estimate
${\mathsf{WP}}$. However, they are difficult in practice, as they require: a
careful encoding of the convoluted rules of a sport into a set of states and
the actions between those states, careful estimation of transition
probabilities, and enough computing power to run enough simulated games to
achieve desired granularity. Each of these are nontrivial.

On the other hand, statistical models are fit entirely from historical data.
Given the results of a set of observed plays, statistical models fit the
relationship between a binary win/loss outcome variable and certain game-state
variables using data-driven regression or machine learning approaches.
Notably, these models are widely used today in football analytics because of
publicly available play-by-play data (e.g., nflFastR (Carl and Baldwin,,
2022)) and accessible off-the-shelf machine learning models (e.g.,
${\mathsf{XGBoost}}$ (Chen and Guestrin,, 2016)). Additionally, due to a
perceived abundance of data, flexible machine learning models are viewed as
more “trustworthy” than previous mathematical models that make more
simplifying assumptions. For these reasons, the open source win probability
models used today for fourth-down decision making in American football are
statistical / machine learning models, which we focus on in this paper.

Statistical win probability models are fit from a historical play-by-play
dataset whose outcome variable is binary win/loss, indicating whether the team
with possession won the game, and whose covariates encode the game-state.
Analysts then typically use blackbox machine learning models like a random
forest (e.g., Lock and Nettleton, (2014)) or ${\mathsf{XGBoost}}$ (e.g.,
Baldwin, (2021)) to fit the model. The outcome variable, however, is noisy and
features a strong dependence structure: each play in the same game shares the
same draw of the win/loss outcome. Accordingly, Brill et al., (2024) found
that these estimators are subject to substantial uncertainty and produce wide
confidence intervals, even when fit from a large dataset featuring $229,635$
first-down plays across $4,101$ games and $16$ years. Nonetheless, such models
continue to be widely used across sports analytics, often without quantifying
uncertainty in their estimates.

To illustrate just how difficult it is to accurately fit a statistical win
probability model from noisy and highly correlated observational data, in this
paper we conduct a simulation study. We create a simplified random walk
version of football in which true win probability at each game-state is known.
Then, we see how well a statistical ${\mathsf{WP}}$ model recovers true win
probability. We find that the dependence structure of observational play-by-
play data inflates both the bias and variance of these estimators. We also
calculate the effective sample size of the observational play-by-play dataset.
An auto-correlated observational play-by-play dataset consisting of $4,101$
games produces a ${\mathsf{WP}}$ estimator that has the same accuracy as one
fit from a dataset having half as many games but consisting entirely of
independent outcomes. In other words, due to the dependence structure of
observational data, we have half as much data as we think. Finally, we explore
the efficacy of bootstrapped confidence intervals in quantifying uncertainty
in win probability estimates. Naive bootstrapped confidence intervals don’t
achieve nominal marginal coverage. Hence, we introduce the fractional
bootstrap, which can be tuned to produce high enough coverage. We find that to
cover true win probability $90\%$ of the time, confidence intervals need to be
substantially wide (i.e., a mean width of $6.3\%$ ${\mathsf{WP}}$).

The remainder of this paper is organized as follows. In Section 2.1 we specify
random walk football, including the rules of the game, how we generate
historical play-by-play datasets for that sport, and our statistical win
probability model. In Section 2.2 we compute the bias-variance decomposition
of win probability estimators fit from various versions of observational
datasets. In Section 2.3 we use this bias-variance decomposition to compute
the effective sample size of a dataset that mimics the historical dataset of
real American football plays. In Section 2.4 we show that naive bootstrapped
win probability confidence intervals don’t achieve nominal coverage and are
too narrow. Finally, in Section 2.5 we introduce the fractional bootstrap,
which produces calibrated bootstrapped confidence intervals that achieve
adequate coverage.

##  2 The simulation study

###  2.1 Introducing random walk football

We begin by describing the rules of random walk football. Random walk football
begins at midfield. Each play, the ball moves left or right by one yardline
with equal probability. If the ball reaches the left (right) end of the field,
team one (two) scores a touchdown, worth $+1$ ($-1$) point. The ball resets to
midfield after each touchdown. After $T$ plays, the game ends. If the game is
still tied after $T$ plays, a fair coin is flipped to determine the winner. We
include the formal mathematical specification of the game in Appendix A. We
also explicitly compute true win probability as a function of time, field
position, and score differential using dynamic programming in Appendix A. In
this study, we let $T=56$, the average number of first-down plays per game in
the dataset of American football plays from Brill et al., (2024). We use $L=4$
yardlines so that the average number of plays between each score is similar to
the average number of first-down plays in a game of American football.

A simulated observational dataset of random walk football plays consists of
$G$ games, each with $T$ plays per game. Such a dataset has the form

${\mathscr{D}}=\\{(t,X_{gt},S_{gt},y_{gt}):t=1,...,T\text{ and }g=1,...,G\\}.$
(1)

For each play of game $g$, we record the timestep $t$, the field position
$X_{gt}$, the score differential $S_{gt}$, and a binary variable $y_{gt}$
indicating whether the team with possession wins the game. We also know the
true win probability ${\mathsf{WP}}_{gt}$ at each play, which we use to
evaluate win probability models but not to train them.

We want to assess the impact of the dependence structure of observational
football data on the accuracy of a statistical win probability estimator. To
do so, we compare the accuracy of win probability estimators fit from datasets
of varying degress of dependence. We introduce a parameter $K$ that controls
the degree of dependence of the win/loss outcome variable in a generated
dataset: we keep a random subsample of $K$ plays per generated game. $K=1$
reflects independence and $K=T$ reflects complete dependence within each game.
For $K=T$, we keep all $T$ plays in each of the $G$ simulated games (recall
$T=56$). When $K=T$, the outcome variable $y_{gt}$ for each play $t$ in game
$g$ reflects the same draw of the win/loss outcome of the game. For $K=1$, we
keep just $1$ randomly sampled play in each of the $G$ simulated games. When
$K=1$, the outcome variable $y_{gt}$ reflects an independent draw of the
win/loss outcome of the game since each play is filtered from a separate game.
Integer values of $K$ between $1$ and $T$ reflect intermediate degrees of
dependence.

Generating a dataset with $G$ games and $T=56$ plays per game, keeping just
$K$ randomly sampled plays per game, has a nominal sample size (number of
rows) of $G\cdot K$ plays. We want to compare the performance of win
probability estimators fit from datasets of varying degrees of dependence $K$
that have the same nominal sample size. To do so, given $K$ we generate
$G={\mathrm{round}}(\zeta\cdot T/K)$ games and keep $K$ plays per game,
yielding a dataset consisting of (approximately) $\zeta\cdot T$ plays. Here,
$\zeta$ parameterizes the sample size.

Finally, continuing the tradition of Baldwin, (2021), in this paper we
estimate win probability using ${\mathsf{XGBoost}}$. The covariates are
${\mathbf{x}}=(t,x,s)$, where $t$ denotes time, $x$ denotes field position,
and $s$ denotes score differential. The outcome variable is binary win/loss
$y$. We use half of the games from the training set as a validation set to
tune the ${\mathsf{XGBoost}}$ model.

###  2.2 Bias-variance decomposition

In this section, we compute the bias-variance decomposition of win probability
estimators fit from datasets having the same sample size $\zeta\cdot T$ but
generated with varying degrees of dependence $K$. Given a combination of data
generating parameters, we generate $M=100$ training datasets and fit a win
probability estimator from each dataset,
$\\{\widehat{{\mathsf{WP}}}^{(m)}\\}_{m=1}^{M}$. We also generate $M=100$ out-
of-sample testing datasets $\\{{\mathscr{D}}_{\text{test}}^{(m)}\\}_{m=1}^{M}$
using $(G=10,000,T=56,K=1)$. Then, we calculate the squared bias of the
$m^{th}$ estimator by

${\mathrm{bias}}^{2}_{m}=\frac{1}{|{\mathscr{D}}_{\text{test}}^{(m)}|}\sum_{{\mathbf{x}}\in{\mathscr{D}}_{\text{test}}^{(m)}}\big{(}{\mathsf{WP}}({\mathbf{x}})-\widehat{{\mathsf{WP}}}^{(m)}({\mathbf{x}})\big{)}^{2}$
(2)

and the variance of the $m^{th}$ estimator by

${\mathrm{var}}_{m}=\frac{1}{|{\mathscr{D}}_{\text{test}}^{(m)}|}\sum_{{\mathbf{x}}\in{\mathscr{D}}_{\text{test}}^{(m)}}\big{(}\widehat{{\mathsf{WP}}}^{(m)}({\mathbf{x}})-\frac{1}{M}\sum_{j=1}^{M}\widehat{{\mathsf{WP}}}^{(j)}({\mathbf{x}})\big{)}^{2}.$
(3)

The root mean squared error of the $m^{th}$ estimator is
${\mathrm{RMSE}}_{m}=\sqrt{{\mathrm{bias}}^{2}_{m}+{\mathrm{var}}_{m}}$. We
then calculate the average squared bias, variance, and ${\mathrm{RMSE}}$
across the $M$ simulations and their standard errors.

First, let $\zeta=4,101$ to mimic the dataset of real American football plays.
We compare the accuracy of an estimator fit from a
$(G={\mathrm{round}}(\zeta\cdot T/K),\ T=56,\ K)$ dataset as the degree of
dependence $K$ varies. The sample size in each of these datasets is
(approximately) the same, $G\cdot K=\zeta\cdot T$. We visualize this bias-
variance decomposition as $K$ varies in Figure 1. As the strength $K$ of the
correlation increases, accuracy decreases linearly. In other words, fixing the
number of observations in the dataset but increasing the degree of dependence
across outcomes, which decreases the effective sample eize, reduces model
accuracy. We compute the explicit reduction in effective sample size in the
next Section 2.3.

Figure 1:  Squared bias (left), variance (middle), and ${\mathrm{RMSE}}$ of a
win probability estimator fit from a $(G={\mathrm{round}}(\zeta\cdot T/K),\
T=56,\ K)$ dataset as a function of $K$, where $\zeta=4,101$. The dots denote
the average values across $M=100$ simulations and the bars denote plus/minus
twice the standard errors. The gray line is the regression line.

Next, as a function of sample size $\zeta\cdot T$ (with $T=56$), we compare
the accuracy of estimators fit from three types of datasets. First, we
consider a $(G=\zeta,\ K=T)$ dataset, which involves keeping each generated
play per game. This dataset mimics the historical dataset of real American
football plays. Then, we consider a $(G=\zeta,\ K=1)$ dataset, which keeps
just one randomly sampled play per game. This dataset consists entirely of
independent outcomes, and can be formed wholly from a $(G=\zeta,\ K=T)$
dataset, but its sample size is much smaller ($\zeta$ rather than $\zeta\cdot
T$). Finally, we consider a $(G=\zeta\cdot T,\ K=1)$ dataset, which has the
same sample size (number of rows) $\zeta\cdot T$ as the first dataset, but
consists entirely of independent outcomes.

In Figure 2 we visualize the bias-variance decomposition of win probability
estimators fit from these three types of datasets as a function of $\zeta$.
The $x$-axis is $\log_{4}(\zeta)$ because $4^{6}=4,096\approx 4,101$ is the
sample size of the historical dataset of American football plays. We see that
it is much better to use all plays per game rather than one independent play
per game. Despite the strong dependence structure, keeping all the plays
provides information about the structure of the covariate space. We also see
that it would be much better if the plays had independent outcomes. This
suggests that the dependence structure reduces the effective sample size of
the dataset. We explore the extent of this reduction in effective sample size
in the next section.

(a)

(b)

Figure 2:  Squared bias (left), variance (middle), and ${\mathrm{RMSE}}$ of a
win probability estimator fit from three datasets as a function of $\zeta$.
The $(G=\zeta,\ K=T)$ dataset (blue) involves keeping each generated play per
game, which mimics the historical dataset of real American football plays. The
$(G=\zeta,\ K=1)$ dataset (red) is formed by keeping just 1 play per game. The
$(G=\zeta\cdot T,\ K=1)$ dataset (orange) has the same sample size (number of
rows) $\zeta\cdot T$ as the first dataset but consists entirely of independent
outcomes.

###  2.3 Effective sample size

As discussed in the previous section, we can calculate the the accuracy of a
win probability estimator fit from a $(G=\zeta,\ T=56,\ K)$ dataset, denoted
${\mathrm{RMSE}}(\zeta,K)$. The sample size (number of rows) of such a dataset
with $\zeta=4,101$ and $K=T$, which mimics the historical dataset of American
football plays, is $\zeta\cdot T$. We saw that the dependence structure of
this dataset reduces the accuracy of our estimator, but we’d like to
understand the extent of this reduction. In particular, we are interested in
the effective sample size (${\mathrm{ESS}}$) of that dataset. The
${\mathrm{ESS}}$ is the sample size $\zeta^{\prime}\cdot T$ of a
$(G=\zeta^{\prime}\cdot T,\ T=56,\ K=1)$ dataset consisting of independent
outcomes, which produces an estimator having the same accuracy as one fit from
the original dataset. For brevity, we refer to the sample size as $\zeta$ and
the effective sample size as $\zeta^{\prime}$, dropping the $T$ since we use
$T=56$ throughout this study.

To estimate this effective sample size, we begin by fitting the $K=1$ and
$K=T$ accuracy curves $\zeta\mapsto{\mathrm{RMSE}}(\zeta,K)$ from Figure 2b.
For each curve, we fit a biexponential model using nonlinear least squares.
Then, as a function of $\zeta$, the effective sample size is the value
$\zeta^{\prime}$ satisfying
${\mathrm{RMSE}}(\zeta^{\prime},K=1)={\mathrm{RMSE}}(\zeta,K=T)$. In Figure 3
we visualize the effective sample size $\zeta^{\prime}$ as a function of
$\zeta$. The effective sample size of a $(\zeta=4,101,K=T)$ dataset is
$\zeta^{\prime}=2,291$, or $56\%$ of the nominal sample size. This result is
striking: we estimate that the historical dataset of American football plays
(where $\zeta=4,101$), consists of about half as much data as suggested by the
number of plays. In other words, we are effectively fitting win probability
from models $8$ years worth of independent win/loss outcomes, not $16$ years
worth. Additionally, if we halved the size of our $K=T$ dataset, fitting a win
probability model from $\zeta=2,050$ games ($8$ seasons) in order to, say, use
just the most recent data, we estimate the effective sample size is
$\zeta^{\prime}=645$, or just $31\%$ of the nominal sample size. If we doubled
the size of our $K=T$ dataset, fitting a win probability model from
$\zeta=8,202$ games ($32$ seasons), we estimate the effective sample size is
$\zeta^{\prime}=6,911$, or $84\%$ of the nominal sample size.

Figure 3:  The effective sample size $\zeta^{\prime}$ of a $(G=\zeta,\ T=56,\
K=T)$ dataset as a function of $\zeta$. The red dot denotes $\zeta=4,101$, the
number of games in the historical dataset of real American football plays.

###  2.4 Coverage of bootstrapped confidence intervals

We have seen that machine learning win probability models fit from noisy and
highly correlated observational data are high variance estimators. Due to the
dependence structure of historical football data, the effective sample size is
much smaller than the nominal sample size. Therefore, it is critical to
quantify uncertainty in win probability point estimates, as the point
estimates alone may not be trustworthy. The bootstrap is a natural choice to
capture such uncertainty since it is non-parametric and doesn’t make strong
assumptions. Hence, in this section we explore the efficacy of bootstrapped
win probability confidence intervals.

We begin with the standard (i.i.d.) bootstrap, which assumes each row (play)
of the dataset is independently drawn. In the standard bootstrap, each of $B$
bootstrapped datasets are formed by re-sampling $G\cdot T$ plays uniformly
with replacement (recall $G$ is the number of games, $T$ is the number of
plays per game, and $G\cdot T$ is the total number of plays in a random walk
football observational dataset). The standard bootstrap doesn’t match the
dependence structure of our observational dataset. Hence, we also try the
cluster bootstrap, in which each of $B$ bootstrapped datasets are formed by
re-sampling $G$ games uniformly with replacement, keeping each observed play
within each re-sampled game. Finally, in the randomized cluster bootstrap,
each of $B$ bootstrapped datasets are formed by re-sampling $G$ games
uniformly with replacement, and within each game re-sampling $T$ plays
uniformly with replacement.

Each type of bootstrap produces $B$ bootstrapped datasets
$\\{{\mathscr{D}}_{\text{train}}^{(b)}\\}_{b=1}^{B}$ from the training dataset
${\mathscr{D}}_{\text{train}}$. We then fit a win probability model to each
bootstrapped dataset using ${\mathsf{XGBoost}}$,
$\\{\widehat{\mathsf{WP}}_{b}\\}_{b=1}^{B}$. From these, we form a $90\%$
confidence interval for ${\mathsf{WP}}({\mathbf{x}})$ at game-state
${\mathbf{x}}$ by the $5^{th}$ and $95^{th}$ quantiles of
$\\{\widehat{\mathsf{WP}}_{b}({\mathbf{x}})\\}_{b=1}^{B}$. Letting $B=101$ in
this section, we form a $90\%$ confidence interval by
$[\widehat{\mathsf{WP}}_{(6)}({\mathbf{x}}),\
\widehat{\mathsf{WP}}_{(96)}({\mathbf{x}})]$. To avoid substantially low
coverage near the extremes (${\mathsf{WP}}({\mathbf{x}})\approx 0$ or
${\mathsf{WP}}({\mathbf{x}})\approx 1$), we widen our confidence intervals
when $\widehat{{\mathsf{WP}}}({\mathbf{x}})<0.025$ to have a lower bound of 0
and when $\widehat{{\mathsf{WP}}}({\mathbf{x}})>0.975$ to have an upper bound
of 1.

We evaluate the efficacy of these bootstrapped confidence intervals by their
coverage and width. For each type of bootstrap (standard bootstrap, cluster
bootstrap, and randomized cluster bootstrap) and each simulation
$m\in\\{1,...,M=100\\}$, we compute the pointwise marginal coverage of
bootstrapped confidence intervals,

${\mathrm{coverage}}_{m}=\frac{1}{|{\mathscr{D}}_{\text{test}}^{(m)}|}\sum_{{\mathbf{x}}\in{\mathscr{D}}_{\text{test}}^{(m)}}\mathbb{I}\left({\mathsf{WP}}({\mathbf{x}})\in{\mathsf{CI}}({\mathbf{x}})\right),$
(4)

and the mean width,

$0pt_{m}=\frac{1}{|{\mathscr{D}}_{\text{test}}^{(m)}|}\sum_{{\mathbf{x}}\in{\mathscr{D}}_{\text{test}}^{(m)}}0pt({\mathsf{CI}}({\mathbf{x}})).$
(5)

We report the average of these values and the standard error in Table 1. To
mimic the historical dataset of American football plays, each simulated
dataset here consists of $G=4,101$ games, $T=56$ plays per game, and $K=T$
plays per game that share the same outcome.

$90\%$ ${\mathsf{CI}}$ method | ${\mathrm{coverage}}$ | $0pt$
---|---|---
standard bootstrap | $0.60\pm 0.01$ | $0.027\pm 0.0005$
cluster bootstrap | $0.71\pm 0.01$ | $0.036\pm 0.0004$
randomized cluster bootstrap | $0.76\pm 0.01$ | $0.042\pm 0.0003$
Table 1:  Pointwise marginal coverage and mean width of nominally $90\%$
confidence intervals formed from each type of bootstrap with $B=101$
bootstrapped re-samples. We report these values averaged over the $M=100$
simulations plus/minus twice their standard errors. Each simulation uses
$G=4,101$, $T=56$, and $K=T$.

Even in our simplified setting of random walk football, each of these
bootstrapped win probability confidence intervals are undercovered. The naive
standard bootstrap in particular achieves dismally low marginal coverage. Even
the randomized cluster bootstrap that accounts for the dependence structure
doesn’t achieve high enough coverage. We suspect this coverage issue would be
even worse for real American football, which is exponentially more complex
than random walk football.

###  2.5 Calibrating bootstrapped confidence intervals

Since these naive bootstrapped confidence intervals are not wide enough, we
wish to tune the bootstrap so that its nominally $90\%$ intervals indeed
achieve $90\%$ marginal coverage. The traditional method of tuning non-
parametric bootstrapped confidence intervals is to calibrate the bootstrapped
quantiles (DiCiccio and Efron,, 1996). For instance, instead of using the
$\alpha/2^{th}$ and $(1-\alpha/2)^{th}$ quantiles of
$\\{\widehat{\mathsf{WP}}_{b}({\mathbf{x}})\\}_{b=1}^{B}$ to form a $1-\alpha$
confidence interval, use the $\beta/2^{th}$ and $(1-\beta/2)^{th}$ quantiles
for some $\beta<\alpha$.

In order for this traditional calibration method to work, $B$ would have to be
much larger than $101$, likely an order of magnitude larger (e.g., $B=1001$).
We, however, prefer to use lower values of $B$ (e.g., closer to $101$) than
larger ones for several reasons. First, in order to develop applications that
utilize bootstrapped win probability confidence intervals in real-time, such
as fourth-down decision bots à la Brill et al., (2024), it is critical to keep
$B$ small. Their fourth-down recommendations for one play take about 15
seconds when $B=101$ and about $2.5$ minutes when $B=1001$. The former can be
run before a fourth-down play begins and the latter takes far too long.
Additionally, storing $1001$ machine learning models is much more expensive
than storing $101$ of them. For these reasons, in this study we stray away
from the traditional bootstrap calibration method.

Instead, we introduce a novel method to calibrate bootstrapped confidence
intervals, the fractional bootstrap, that has the same time and storage
complexity as traditional bootstrap methods. Specifically, we introduce a
tuning parameter $\phi\in(0,1]$ representing the fraction of data to be re-
sampled in generating a bootstrapped dataset. By re-sampling less data than in
the original dataset, we widen bootstrapped confidence intervals and increase
coverage. In the fractional standard bootstrap, we re-sample $T\cdot
G\cdot\phi$ plays (rows) uniformly with replacement. In the fractional cluster
bootstrap, we re-sample $G\cdot\phi$ games uniformly with replacement, keeping
each observed play within each re-sampled game. Finally, in the fractional
randomized cluster bootstrap, we re-sample $G\cdot\phi$ games uniformly with
replacement, and within each game re-sample $T$ plays uniformly with
replacement.

In Table 2 we report the results of our simulation study in which we use
$\phi$ to tune the randomized cluster bootstrap. To mimic the historical
dataset of American football plays, each simulated dataset consists of
$G=4,101$ games, $T=56$ plays per game, and $K=T$ plays per game that share
the same outcome. As expected, lowering $\phi$ widens the confidence intervals
and increases marginal coverage. In order for fractional randomized cluster
bootstrapped confidence intervals to achieve nominal $90\%$ marginal coverage,
$\phi$ needs to be as small as $0.35$. Those intervals have a mean width of
$6.3\%$. This result is striking: in our simplified setting of random walk
football, win probability confidence intervals need to be extremely wide to
achieve nominal coverage. This exemplifies the difficulty of accurately
estimating win probability by fitting a machine learning model from noisy and
highly correlated football game outcomes. These are high variance estimators
subject to large uncertainty.

$\phi$ | ${\mathrm{coverage}}$ | $0pt$
---|---|---
1 | $0.76\pm 0.01$ | $0.042\pm 0.0003$
0.75 | $0.80\pm 0.01$ | $0.047\pm 0.0003$
0.5 | $0.85\pm 0.01$ | $0.055\pm 0.0004$
0.35 | $0.90\pm 0.01$ | $0.063\pm 0.0004$
Table 2:  Pointwise marginal coverage and mean width of $90\%$ confidence
intervals formed from the fractional randomized cluster bootstrap with $B=101$
bootstrapped re-samples. We report these values averaged over the $M=100$
simulations plus/minus twice their standard errors. Each simulation uses
$T=56$, $K=56$, and $G=4,101$.

Marginal coverage is a sufficient condition for confidence intervals to be
“good,” but it is not a necessary condition. Even with $90\%$ marginal
coverage, it could be that ${\mathsf{CI}}({\mathbf{x}})$ always covers
${\mathsf{WP}}({\mathbf{x}})$ for $90\%$ of game-states ${\mathbf{x}}$ and
never covers ${\mathsf{WP}}({\mathbf{x}})$ for the other $10\%$ of game-states
${\mathbf{x}}$. In that case, it may be disastrous to make decisions at the
game-states for which intervals never cover. To check that our fractional
randomized cluster bootstrapped confidence intervals achieve reasonable
coverage across the space of game-states, we bin game-states ${\mathbf{x}}$ by
their true win probability ${\mathsf{WP}}({\mathbf{x}})$ and consider coverage
in each bin. In Figure 4 we visualize coverage and its standard error across
such bins. Although the intervals don’t achieve the nominal $90\%$ coverage in
each bin, coverage is high enough in each bin for us to comfortably make
decisions using these intervals.

Figure 4:  Coverage of $90\%$ bootstrapped confidence intervals, via the
fractional randomized cluster bootstrap for $\phi=0.35$, across the space of
game-states ${\mathbf{x}}$ binned by ${\mathsf{WP}}({\mathbf{x}})$.

## References

  * Baldwin, (2021) Baldwin, B. (2021).  NFL win probability from scratch using xgboost in R.
  * Brill et al., (2024) Brill, R. S., Yurko, R., and Wyner, A. J. (2024).  Analytics, have some humility: a statistical view of fourth-down decision making.
  * Carl and Baldwin, (2022) Carl, S. and Baldwin, B. (2022).  nflfastR: Functions to Efficiently Access NFL Play by Play Data.  https://www.nflfastr.com/.
  * Chen and Guestrin, (2016) Chen, T. and Guestrin, C. (2016).  XGBoost: A scalable tree boosting system.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, pages 785–794, New York, NY, USA. ACM.
  * DiCiccio and Efron, (1996) DiCiccio, T. J. and Efron, B. (1996).  Bootstrap confidence intervals.  Statistical Science, 11(3):189–212.
  * Lock and Nettleton, (2014) Lock, D. and Nettleton, D. (2014).  Using random forests to estimate win probability before each play of an nfl game.  Journal of Quantitative Analysis in Sports, 10.

SUPPLEMENTARY MATERIAL

##  Appendix A Simulation study details

Formally, the outcome of the $t^{th}$ play of the $g^{th}$ game is

$\xi_{gt}\overset{iid}{\sim}\pm 1.$ (6)

The game starts at midfield, $X_{g0}=L/2$, and the game begins tied,
$S_{g0}=0$. The field position at the start of play $t$ is

$X_{g,t+1}:=\begin{cases}X_{gt}+\xi_{gt}&\text{ if }0<X_{gt}+\xi_{gt}<L\text{
(not a TD)}\\\ L/2&\text{ else,}\end{cases}$ (7)

and the score differential at the start of play $t$ is

$S_{g,t+1}:=\begin{cases}S_{gt}+1&\text{ if }X_{gt}+\xi_{gt}=0\text{ (TD)}\\\
S_{gt}-1&\text{ if }X_{gt}+\xi_{gt}=L\text{ (opp. TD)}\\\ S_{gt}&\text{
else.}\end{cases}$ (8)

The response column win is

$y_{gt}\equiv y_{g,T+1}:=\begin{cases}1&\text{ if }S_{g,T+1}>0\\\ 0&\text{ if
}S_{g,T+1}<0\\\ \text{Bernoulli}(1/2)&\text{ else (overtime). }\end{cases}$
(9)

As in our dataset of real football plays, this response column is highly
correlated – plays from the same game share the same draw of the winner of the
game.

The true win probability

${\mathsf{WP}}(t,x,s):={\mathbb{P}}(S_{g,T+1}>0|X_{gt}=x,S_{gt}=s)$ (10)

of random walk football is computed explicitly using dynamic programming,

${\mathsf{WP}}(T+1,x,s)=\begin{cases}1&\text{if }s>0\\\ 1/2&\text{if }s=0\\\
0&\text{if }s<0,\\\ \end{cases}$ (11)

and

${\mathsf{WP}}(t-1,x,s)=\begin{cases}\frac{1}{2}{\mathsf{WP}}(t,\frac{L}{2},s+1)+\frac{1}{2}{\mathsf{WP}}(t,x+1,s)&\text{if
}x=1\\\
\frac{1}{2}{\mathsf{WP}}(t,x-1,s)+\frac{1}{2}{\mathsf{WP}}(t,\frac{L}{2},s-1)&\text{if
}x=L-1\\\
\frac{1}{2}{\mathsf{WP}}(t,x-1,s)+\frac{1}{2}{\mathsf{WP}}(t,x+1,s)&\text{else.}\end{cases}$
(12)