Fingerprints of hot-phonon physics in time-resolved correlated quantum lattice
dynamics

E. Cappelluti1*, D. Novko2,3

1 Istituto di Struttura della Materia, CNR (ISM-CNR), 34149 Trieste, Italy

2 Institute of Physics, 10000 Zagreb, Croatia

3 Donostia International Physics Center (DIPC), 20018 Donostia-San Sebastián,
Spain

*<EMAIL_ADDRESS>
## Abstract

The time dynamics of the energy flow from electronic to lattice degrees of
freedom in pump-probe setups could be strongly affected by the presence of a
hot-phonon bottleneck, which can sustain longer coherence of the optically
excited electronic states. Recently, hot-phonon physics has been
experimentally observed and theoretically described in MgB2, the electron-
phonon based superconductor with $T_{\rm c}\approx 39$ K. By employing a
combined ab-initio and a semiclassical approach and by taking MgB2 as an
example, here we propose a novel path for revealing the presence and
characterizing the properties of hot phonons through a direct analysis of the
information encoded in the lattice inter-atomic correlations. Such method
exploits the underlying symmetry of the $E_{2g}$ hot modes characterized by a
out-of-phase in-plane motion of the two boron atoms. Since hot phonons occur
typically at high-symmetry points of the Brillouin zone, with specific
symmetries of the lattice displacements, the present analysis is quite general
and it could aid in revealing the hot-phonon physics in other promising
materials, such as graphene, boron nitride, black phosphorus, or even
underdoped cuprates.

###### Contents

  1. 1 Introduction
  2. 2 Theoretical modelling
  3. 3 Time-resolved lattice dynamics
  4. 4 Wider perspective on other materials and experimental observations
  5. 5 Conclusions

##  1 Introduction

Ultrafast time-resolved measurements have proven in the past few decades to be
a powerful tool for investigating fundamental mechanisms of a variety of
physical phenomena and for revealing processes governing different states of
the matter, where several degrees of freedom (e.g., electrons, lattice and
magnetic modes) play a mutual role [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13]. In a typical pump-probe setup, energy is initially injected into the
electronic degrees of freedom by means of particle-hole excitations, giving
rise to a non-thermal electronic distribution. The time evolution of such non-
thermal excitations is governed by different scattering mechanisms that lead
to an energy transfer towards other degrees of freedom, in particular in the
lattice mode sector, until, on a long time scale, a new final steady state is
reached where all the different degrees of freedom are thermalized and at
equilibrium.

The possibility of sustaining non-thermal states for sufficiently long time
opens striking opportunities in the field of quantum information as long as
the scattering sources affecting such states can be kept under control. The
electron-phonon (el-ph) coupling is of major importance in this context since
it represents one of the primary channels responsible for the internal
thermalization of the electronic degrees of freedom, and because it leads
eventually, when it is not constrained, to a system heating detrimental for
any coherent process [14].

However, the thermalization process of some materials can be hindered by the
so-called phonon bottleneck[15, 16, 17, 18, 19, 20, 2, 3, 21, 22, 23, 24, 25,
26, 27, 28, 29, 30, 31, 32, 33, 34, 35] that quenches the energy transfer
between electron and lattice degrees of freedom. Such a physical phenomenon is
usually encountered in semiconductors and semimetals where the pump-driven
particle-hole excitations are restricted to few single points (valleys) in the
Brillouin zone. In such conditions only a few phonons with the right symmetry
and with selected momentum ${\bf q}$ connecting two valleys can be excited
efficiently. Energy from the electron sub-system is thus prevalently
transferred to these modes that get “hot”, meaning they acquire a phonon
population $b_{\bf q}$ much larger than other modes, giving rise to a non-
thermal phonon distribution.

Within this context, the possibility of a phonon bottleneck, and hence of hot
phonons, has been so far associated only to semiconductors and semimetals,
whereas conventional metals, with a large Fermi surface allowing scattering
with many ${\bf q}$-phonons, were thought to be incompatible with such a
scenario. At odds with this belief, a novel path for inducing hot phonons has
been recently proposed for unconventional metals characterized by a strong
anisotropy of the el-ph coupling [36, 37]. In the paradigmatic case of MgB2,
the el-ph coupling is concentrated in a few in-plane modes at the Brillouin
zone center possessing the $E_{2g}$ symmetry[38, 39, 40, 41, 42, 43, 44, 45],
with a strong resemblance to the relevant modes in single- and multi-layer
graphene. As a consequence of such remarkable anisotropy, the energy initially
pumped into the electron sector is efficiently and rapidly ($\sim 50$ fs)
transferred only to such few lattice degrees of freedom that get a much larger
population compared to the remaining lattice modes, whereas a final
thermalization among the whole phonon modes and among electron and lattice
degrees of freedom occurs on a much larger time scale ($\sim 0.4$ ps) [37]. On
the experimental ground, an indirect signature of such hot-phonon scenario has
been observed in the onset of unconventional features in the time-resolved
optical spectra [36]. More direct and precise fingerprints have been proposed
at the theoretical level in the analysis of time-resolved Raman spectroscopy,
both in spectral features as well as in integrated spectral weights [37]. An
accurate measurement of time-resolved Raman spectra in MgB2 needs to face with
the limitations of the uncertainty principle constraining time and energy
resolution [46]. In addition, it should be remarked that in both cases the
presence of a non-thermal population of the in-plane $E_{2g}$ modes is not
directly probed via the properties of the lattice degrees of freedom, but
indirectly through the el-ph driven many-body renormalization of the
electronic response, i.e. in the optical response [36] or in the phonon self-
energy [37, 47], which is strictly related to the electronic screening.

In the present work we suggest that a suitable evidence of the presence of hot
phonons in MgB2 can be attained in a more straightforward way by the analysis
of the time-resolved lattice dynamics, as obtained for instance by means of
ultrafast electron diffuse scattering [48, 49, 50, 11, 51, 12, 52, 53, 54, 13,
55, 35, 56, 57, 58]. In particular, we show that the inter-atomic pair-
distribution function, which encodes information about the correlated atomic
motion, is highly sensitive to the presence of a hot phonon and it provides a
useful tool for detecting them. As related to the analysis of time-independent
quantities, our proposal is not affected by the limitations of uncertainty
principle. Note that the present analysis, here applied to the paradigmatic
case of MgB2, does not rely on specific properties of this material but it can
be employed as well in other compounds, such as single-layer and multi-layer
graphene [59, 60], transition metal dichalcogenides [61], black phosphorous
[35, 55], hole-doped diamond [62, 63] and underdoped cuprates [64, 65], where
el-ph coupling plays a major role.

The paper is structured as follows: Sec. 2 introduces the current modelling
for hot phonons in MgB2 as previously discussed in Ref. [37]. Such basic
notions will be employed in Sec. 3 to investigate the effects of hot phonons
on the correlated lattice dynamics in MgB2. A wider perspective of the present
analysis in regards to other materials and possible experimental approaches is
finally provided in Sec. 4.

##  2 Theoretical modelling

Figure 1:  (a) Crystal structure of MgB2. Also displayed are the lattice
displacements of the in-plane $E_{2g}$ mode at ${\bf q}=0$. (b) Electronic
band dispersion of MgB2. Bands are labelled as $\sigma$\- and $\pi$-bands
according to their symmetry. (c) Phonon dispersions of MgB2. The color code
(color bar atop) represents the strength of the el-ph coupling
$\lambda_{\mathbf{q}\nu}$. (d) Corresponding phonon density of states
$F(\omega)$ (dashed line) and the total Eliashberg function
$\alpha^{2}F(\omega)$ (blue solid line). The solid red line shows the
contribution to the Eliashberg function associated with the hot $E_{2g}$
modes. Readaptation from Ref. [37].

The crystal structure of MgB2 is rather simple, with hexagonal graphene-like
planes of B atoms spaced vertically by Mg atoms located in the center of the
hexagons (see Fig. 1a) [66]. In order to reproduce an ab-initio ground state
properties of this material, we employ here the quantum espresso package [67].
Norm-conserving pseudopotentials were employed with the Perdew-Burke-Ernzerhof
exchange-correlation functional [68]. The in-plane lattice parameter and the
distance between two boron planes were set to 3.083 Å and 3.521 Å,
respectively. A $24\times 24\times 24$ Monkhorst-Pack grid in momentum space
and a plane-wave cutoff energy of 60 Ry were used for ground-state
calculations. The phonon dispersion $\omega_{{\bf q},\nu}$ was calculated on a
$12\times 12\times 12$ grid using density functional perturbation theory [69].
The el-ph coupling strength for each momentum ${\bf q}$ and each phonon branch
$\nu$ was computed by using the epw code [70] according the standard
expression:

$\displaystyle\lambda_{{\bf q},\nu}$ $\displaystyle=$
$\displaystyle\frac{2}{N(0)N}\sum_{{\bf k},n,m}\frac{|g_{{\bf k}n,{\bf
k+q}m}^{\nu}|^{2}}{\omega_{{\bf q},\nu}}\delta(\epsilon_{{\bf
k}n}-\mu)\delta(\epsilon_{{\bf k+q}m}-\mu),$ (1)

where $\mu$ is the electron chemical potential, $\epsilon_{{\bf k}n}$ are the
electron energies, $\omega_{{\bf q},\nu}$ are the phonon frequencies, $g_{{\bf
k}n,{\bf k+q}m}^{\nu}$ are the el-ph matrix elements, $N(\mu)=\sum_{{\bf
k},n}\delta(\epsilon_{{\bf k}n}-\mu)$ is the electronic density of states per
spin at the Fermi level, and $N$ is the total number of sampling points in the
Brillouin zone. The delta functions $\delta(x)$ are modeled by Gaussian
functions with finite broadenings. Electron energies, phonon frequencies, and
electron-phonon coupling matrix elements were interpolated using maximally-
localized Wannier functions [71]. The Eliashberg function, which shows a
strength of the el-ph coupling for a particular phonon energy, was obtained on
a $40\times 40\times 40$ grid of electron and phonon momenta.

The well-known band structure is shown in Fig. 1b. Due to the symmetry of the
systems, the in-plane $\sigma$-bands retain a strong two-dimensional character
with a very weak amount of Mg, very similar to the $\sigma$ bands of graphene,
whereas the $p_{z}$ B-orbitals strongly hybridize with Mg giving rise to
$\pi$-bands with a strong three-dimensional character. As a consequence of
such strong anisotropy in the electronic degrees of freedom, the el-ph
coupling results to be highly anisotropic, as depicted in Fig. 1c. The el-ph
coupling among states in the $\pi$ bands and among $\sigma$ and $\pi$ states,
involving $\pi$ bands with a good metallic character, is quite weak. On the
other hand, the $\sigma$ bands appear to be only slightly doped, resulting in
a poor screening of the in-plane B-B lattice modes, and in a corresponding
large el-ph coupling between states in the $\sigma$ bands. A further
consequence of the small hole-doping of the two-dimensional $\sigma$ bands is
the intrinsic restriction, obeying to the Lindhard model, of a sizable el-ph
coupling in a small momentum region $|{\bf q}_{||}|\leq 2k_{\rm F}$,[43, 72,
73, 74, 75] where ${\bf q}_{||}=(q_{x},q_{y})$ and $k_{\rm F}$ is the in-plane
Fermi momentum of the $\sigma$ bands which is only weakly dependent on
$k_{z}$.

Within the above framework, a predominant role is played by the in-plane
$E_{2g}$ phonon modes at small ${\bf q}_{||}$, characterized by a strong el-ph
coupling and which corresponds to out-of-phase displacements of the two boron
atoms within the unit cell, as depicted in Fig. 1a. The overall strongly-
coupled modes are thus limited, in a schematic way, to only two phonon
branches within the energy window $\omega_{{\bf q},\nu}\in[60:75]$ meV, with
roughly $|{\bf q}_{||}|\leq 2k_{\rm F}$ and weak dispersion along $q_{z}$,
amounting roughly to 5 % of the total phonon modes, and corresponding thus to
a very small specific heat capacity. The crucial role of the few $E_{2g}$
modes in the total coupling can be clearly pointed out in the Eliashberg
function (Fig. 1d) which shows a remarkable peak in the range $\omega_{{\bf
q},\nu}\in[60:75]$ meV, in spite of the absence of any particular feature in
the corresponding phonon density of states.

Sophisticated and powerful tools for evaluating the non-thermal time-evolution
of the coupled electron-lattice system have been developed in the recent years
[34, 35, 61, 58, 57]. Such techniques appear particularly useful when in
investigating the time range when effective electron and/or lattice
temperatures cannot be properly defined. The hot-phonon physics we discuss
here involves however macroscopical quantities integrated over momenta, where
such level of accuracy is not needed and effective-temperature models can
provide as well a reliable description. Within this context, following Ref.
[37], the predominance of the hot-phonon modes in the total coupling can thus
be captured by splitting the total Eliashberg function in a hot and cold
component, $\alpha^{2}F(\omega)=\alpha^{2}F_{\rm hot}(\omega)+\alpha^{2}F_{\rm
cold}(\omega)$, where $\alpha^{2}F_{\rm hot}(\omega)$ contains the
contribution of the hot $E_{2g}$ modes close to the
$\mathrm{\overline{\Gamma}-\overline{A}}$ path in the relevant energy range
$\omega\in[60:75]$ meV (the phonon modes colored in red in Fig. 1c), while
$\alpha^{2}F_{\rm cold}(\omega)$ takes into account the other weakly coupled
cold modes. A similar splitting can be performed for the phonon density of
states $F(\omega)=F_{\rm hot}(\omega)+F_{\rm cold}(\omega)$.

Equipped with such first-principles input, we can describe the time-dynamics
of the energy transfer between the different degrees of freedom in terms of
three temperatures [76, 1, 77, 78, 9, 37, 59, 60, 79], i.e., an electron one
$T_{\rm e}$, a hot-phonon temperature $T_{\rm hot}$, which governs the
population of the $E_{2g}$-like strongly-coupled hot-phonon modes, and the
cold-lattice temperature $T_{\rm cold}$ that describes the temperature of the
weakly coupled cold modes. The time evolution of the these characteristic
temperatures upon a pump pulse can be described by the set of coupled
equations

$\displaystyle C_{\rm e}\frac{\partial T_{\rm e}}{\partial t}$
$\displaystyle=$ $\displaystyle S(z,t)+\nabla_{z}(\kappa\nabla_{z}T_{\rm
e})-G_{\rm hot}(T_{\rm e}-T_{\rm hot})-G_{\rm cold}(T_{\rm e}-T_{\rm cold}),$
(2) $\displaystyle C_{\rm hot}\frac{\partial T_{\rm hot}}{\partial t}$
$\displaystyle=$ $\displaystyle G_{\rm hot}(T_{\rm e}-T_{\rm hot})-C_{\rm
hot}\frac{T_{\rm hot}-T_{\rm cold}}{\tau_{0}},$ (3) $\displaystyle C_{\rm
cold}\frac{\partial T_{\rm cold}}{\partial t}$ $\displaystyle=$ $\displaystyle
G_{\rm cold}(T_{\rm e}-T_{\rm cold})+C_{\rm hot}\frac{T_{\rm hot}-T_{\rm
cold}}{\tau_{0}}.$ (4)

Here $C_{\rm e}$, $C_{\rm hot}$, and $C_{\rm cold}$ are the specific heat
capacities for the electron, hot-phonon, and cold-phonon states, respectively,
and $G_{\rm hot}$/$G_{\rm cold}$ are the electron-phonon relaxation rates
between electronic states and hot/cold phonons modes. They can be computed as

$\displaystyle C_{\rm e}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}d\varepsilon
N(\varepsilon)\varepsilon\frac{\partial f(\varepsilon-\mu;T_{\rm e})}{\partial
T_{\rm e}},$ (5) $\displaystyle C_{\rm hot}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\omega F_{E_{\rm
hot}}(\omega)\omega\frac{\partial b(\omega;T_{E_{2g}})}{\partial T_{E_{2g}}},$
(6) $\displaystyle C_{\rm cold}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\omega F_{\rm
cold}(\omega)\omega\frac{\partial b(\omega;T_{\rm cold})}{\partial T_{\rm
cold}},$ (7) $\displaystyle G_{\rm hot}$ $\displaystyle=$
$\displaystyle\frac{2\pi k_{B}}{\hbar}N(\mu)\int
d\Omega\Omega\alpha^{2}F_{E_{2g}}(\Omega),$ (8) $\displaystyle G_{\rm cold}$
$\displaystyle=$ $\displaystyle\frac{2\pi k_{B}}{\hbar}N(\mu)\int
d\Omega\Omega\alpha^{2}F_{\rm cold}(\Omega).$ (9)

Here $N(\varepsilon)$ is the electronic density of states, $\mu$ the
electronic chemical potential, and $f(x;T)=1/[\exp(x/T)+1]$,
$b(x;T)=1/[\exp(x/T)-1]$ are the Fermi-Dirac and the Bose-Einstein
distribution functions, respectively. The values obtained from the first-
principles calculations are $C_{\mathrm{e}}=90\,\mathrm{J/m^{3}K^{2}}\times
T_{\mathrm{e}}$, $C_{\rm hot}=0.13\,\mathrm{J/m^{3}K}$, and
$C_{\mathrm{cold}}=4.1\,\mathrm{J/m^{3}K}$, $G_{\rm hot}=2.8\times
10^{18}\,\mathrm{W/m^{3}K}$ and $G_{\rm cold}=3.6\times
10^{18}\,\mathrm{W/m^{3}K}$. Furthermore, modelling a pump-probe setup, the
term $S(z,t)=I(t)e^{-z/\delta}/\delta$ takes into account the energy absoption
from the laser pulse into the electronic degrees of freedom, where $I(t)$ is
the intensity of the absorbed fraction of the laser pulse (with a Gaussian
profile) and $\delta$ is the penetration depth. Finally, the relaxation time
$\tau_{0}$ takes into account the scattering between hot and cold modes driven
by the anharmonic phonon-phonon coupling, for which we take an estimate of 400
fs [80].

Figure 2: (a) Time evolution of the characteristic effective temperatures
$T_{\mathrm{e}}$, $T_{\rm hot}$, $T_{\rm cold}$ for the three-temperature
model suitable for hot phonons. The grey dashed line shows the pulse profile,
with the pulse duration of 45 fs and an absorbed fluence of 12 J/m2. (b)
Similar as in panel (a) within the assumption of thermal distribution (two-
temperature model) of the lattice modes as described with Eqs. (10) and (11).
Partial readaptation from Ref. [37].

The time evolution of the three characteristic temperatures $T_{\rm e}$,
$T_{\rm hot}$, $T_{\rm cold}$ (Fig. 2a) shows a remarkable increase of the
effective temperature of the hot modes $T_{\rm hot}$, which reaches a peak
$T_{\rm hot}^{\rm max}\approx 1250$ K at $t^{*}\approx 50$ fs with a small
delay compared to the time behavior of the electronic temperature $T_{\rm e}$.
The remnant cold phonon modes follow, on the other hand, a completely
different behavior with a slow monotonic increase towards the final
equilibrium state at $t\geq 0.3-0.4$ ps where all the degrees of freedom are
thermalized with each other.

Such a behavior, showing a clear hot-phonon scenario, can be compared with the
isotropic case where the energy from the electronic degrees of freedom is
transferred in a equal way to all the lattice modes, without a preferential
channel, and one can thus define a unique phonon temperature $T_{\rm ph}$ for
all the modes [76]. The coupled equations of the two temperature model can be
obtained by setting $T_{\rm hot}=T_{\rm cold}\equiv T_{\rm ph}$ in Eqs.
(2)-(4):

$\displaystyle C_{\rm e}\frac{\partial T_{\rm e}}{\partial t}$
$\displaystyle=$ $\displaystyle S(z,t)+\nabla_{z}(\kappa\nabla_{z}T_{\rm
e})-G_{\rm ph}(T_{\rm e}-T_{\rm ph}),$ (10) $\displaystyle C_{\rm
ph}\frac{\partial T_{\rm ph}}{\partial t}$ $\displaystyle=$ $\displaystyle
G_{\rm ph}(T_{\rm e}-T_{\rm ph}),,$ (11)

where $C_{\rm ph}=C_{\rm hot}+C_{\rm cold}$ and $G_{\rm ph}=G_{\rm hot}+G_{\rm
cold}$. The time evolution of $T_{\rm e}$, $T_{\rm ph}$ is displayed in Fig.
2b. In the absence of the hot-phonon bottleneck, a much faster thermalization
between the electrons and the lattice degrees of freedom is obtained, within a
time-scale of $t\sim 0.08$ ps, while the average phonon temperature $T_{\rm
ph}$ reaches the maximum of about 420 K.

##  3 Time-resolved lattice dynamics

As discussed above, the hot-phonon scenario is characterized by the regime
$T_{\rm hot}\gg T_{\rm cold}$, where the population of the hot-lattice modes
is singled out and is significantly larger than the population of the other
cold modes. Equipped with the detailed input from the _ab initio_
calculations, we discuss now how this scenario in MgB2 can be revealed
directly in the dynamical lattice properties, i.e., in the amplitude of the
mean square lattice displacements, resolved for atomic species and for
different directions; and more strikingly in the degree of correlation of the
atomic motion. The peculiar properties of hot phonons in MgB2 rely on the
$E_{2g}$ symmetry of the strongly coupled hot modes, in particular: ($i$) the
pure in-plane B character of the $E_{2g}$ modes; ($ii$) the counter-phase
motion of the two B atoms per cell that move in opposite direction, as
depicted in Fig. 1a.

The evaluation of the mean square lattice displacements based on the force-
constant-model calculations has been detailed in Refs. [81, 82]. Here, we
extend such computation by using fully _ab initio_ methodology and, more
importantly, by introducing the hot-phonon scenario and specifying a different
temperatures for each relevant modes. More explicitly, we can define the
projected mean-square lattice displacement $\sigma^{2}(i_{\alpha})$ as

$\sigma^{2}(i_{\alpha})=\langle[{\bf u}_{i}\cdot\hat{\bf
r}_{\alpha}]^{2}\rangle,$ (12)

where ${\bf u}_{i}$ is the lattice displacement of atom $i$ from its average
position and $\hat{\bf r}_{\alpha}$ is the unit vector pointing along the
direction $\alpha=x,y,z$. On the microscopic ground we can write

$\sigma^{2}(i_{\alpha})=\frac{\hbar}{N}\sum_{{\bf q},\mu}\frac{|\epsilon_{{\bf
q},\mu}^{i\alpha}|^{2}}{M_{i}\omega_{{\bf
q},\mu}}\left[\frac{1}{2}+b\left(\frac{\omega_{{\bf q},\mu}}{T_{{\bf
q},\mu}}\right)\right],$ (13)

where $N$ is the total number of ${\bf q}$-points considered in the phonon
Brillouin zone, $M_{i}$ is the atomic mass of atom $i$, $\omega_{{\bf q},\mu}$
is the frequency of a phonon mode of branch $\mu$ with momentum ${\rm q}$,
$T_{{\bf q},\mu}=T_{\rm hot},T_{\rm cold}$ is the appropriate temperature for
such mode, and $\epsilon_{{\bf q},\mu}^{i\alpha}$ is the component of the
corresponding eigenvector $\hat{\epsilon}_{{\bf q},\mu}$ describing the
displacement of the $i$ atom along the $\alpha$ direction.

Figure 3: Time evolution of the mean square lattice displacements $\sigma^{2}$
for each atom and for different axis directions in the case of the hot-phonon
scenario, described by Eqs. (2)-(4) (full lines). The corresponding results
for the thermal phonon distribution [described by Eqs. (10)-(11)] are shown
with the dashed lines.

The behavior of the mean square lattice displacements resolved for each atom
and for different axis directions in the case of hot non-thermal phonon
distribution is shown in Fig. 3 (solid lines). As expected, compared with the
case of a thermal phonon distribution described by Eqs. (10)-(11) (dashed
lines), we do not see any remarkable difference for the phonon modes with Mg
character and with out-of-plane lattice displacements, apart from a slower
thermalization in the case of hot phonons, due to the longer storage of energy
in the hot-phonon modes. A slight effect of the presence of hot phonons can be
detected in the behavior of the mean square in-plane B lattice displacements
which shows a non-monotonic behavior with a peak value larger than the final
equilibrium state, at $t^{*}=50$ fs, in fair accordance with the peak of
$T_{\rm hot}$.

A more striking signature of the hot-phonon scenario can be detected from the
analysis of the correlated lattice dynamics. More explicitly, we focus on the
mean-square relative displacement of atomic pairs projected onto the vector
joining the atom pairs [81, 82]:

$\sigma^{2}_{ij}=\langle[({\bf u}_{i}-{\bf u}_{j})\cdot\hat{\bf
r}_{ij}]^{2}\rangle,$ (14)

where ${\bf u}_{i}$ and ${\bf u}_{j}$ are the lattice displacements of atoms
$i$ and $j$ from their average positions, $\hat{\bf r}_{ij}$ is the unit
vector connecting atoms $i$ and $j$.

This physical quantity can be evaluated microscopically as:

$\displaystyle\sigma^{2}_{ij}$ $\displaystyle=$
$\displaystyle\frac{\hbar}{N}\sum_{\bf
q,\mu}\left[\frac{1}{2}+b\left(\frac{\omega_{{\bf q},\mu}}{T_{{\bf
q},\mu}}\right)\right]\Bigg{\\{}\frac{|\hat{\epsilon}^{i}_{{\bf
q},\mu}\cdot\hat{\bf r}_{ij}|^{2}}{M_{i}\omega_{\bf
q,\mu}}+\frac{|\hat{\epsilon}^{j}_{{\bf q},\mu}\cdot\hat{\bf
r}_{ij}|^{2}}{M_{j}\omega_{\bf q,\mu}}$ (15)
$\displaystyle-\frac{2\mbox{Re}\left[(\hat{\epsilon}^{i}_{{\bf
q},\mu}\cdot\hat{\bf r}_{ij})(\hat{\epsilon}^{j*}_{{\bf q},\mu}\cdot\hat{\bf
r}_{ij})\mbox{e}^{i{\bf q}\cdot{\bf r}_{ij}}\right]}{\omega_{\bf
q,\mu}\sqrt{M_{i}M_{j}}}\Bigg{\\}}.$

It is useful to quantify the degree of correlation by introducing a convenient
dimensionless correlation factor $\rho_{ij}$ defined as

$\displaystyle\sigma^{2}_{ij}$ $\displaystyle=$
$\displaystyle\sigma^{2}(i_{j})+\sigma^{2}(j_{i})-2\sigma(i_{j})\sigma(j_{i})\rho_{ij},$
(16)

where

$\displaystyle\sigma^{2}(i_{j})$ $\displaystyle=$ $\displaystyle\langle[({\bf
u}_{i}\cdot\hat{{\bf r}}_{ij}]^{2}\rangle$ (17) $\displaystyle=$
$\displaystyle\frac{\hbar}{N}\sum_{\bf
q,\mu}\left[\frac{1}{2}+b\left(\frac{\omega_{{\bf q},\mu}}{T_{{\bf
q},\mu}}\right)\right]\frac{|\hat{\epsilon}^{i}_{{\bf q},\mu}\cdot\hat{\bf
r}_{ij}|^{2}}{M_{i}\omega_{\bf q,\mu}}.$

The inter-atomic correlation factor can be thus calculated as[82]

$\displaystyle\rho_{ij}$ $\displaystyle=$
$\displaystyle\frac{\sigma^{2}(i_{j})+\sigma^{2}(j_{i})-\sigma^{2}_{ij}}{2\sigma(i_{j})\sigma(j_{i})}.$
(18)

Positive values of correlation factor $\rho_{ij}>0$ describe a situation where
the couple of atoms $i$, $j$ move in phase, so that the resulting value of
$\sigma^{2}_{ij}$ is smaller than for the uncorrelated case. On the other
hand, a predominance of counter-phase atomic vibrations is expected to result
in a negative correlation factor $\rho_{ij}<0$.

Equations (13) and (18) provide fundamental tools for investigating the
effects of the hot-phonon scenario on the time evolution of the correlated
lattice dynamics. At the initial thermal equilibrium at room temperature we
find positive values $\rho_{\rm B-B}\approx 0.24$ [83], $\rho_{\rm
B-Mg}\approx 0.22$, reflecting a slight dominance of the acoustic (in-phase)
modes in contributing to the lattice dynamics.

Figure 4: (a) Time evolution of the mean-square relative lattice displacement
$\sigma_{ij}^{2}$ as defined by Eq. (15), for the in-plane nearest neighbor
B-B bond, and for the nearest neighbor B-Mg bond. Solid lines show the hot-
phonon scenario, whereas dashed lines represent the time behavior assuming a a
thermal phonon distribution as modelled by Eqs. (10)-(11). (b) Corresponding
correlation functions $\rho_{ij}$.

The time dependence of the mean-square relative displacement for the in-plane
boron-boron pair, $\sigma_{\rm B-B}$, and for the interplane boron-magnesium
one $\sigma_{\rm B-Mg}$ is reported in Fig. 4a. The inter-atomic boron-
magnesium lattice displacement (blue full line) does not show any significant
feature, just a monotonic increase between the two equilibrium states, as in
the absence of hot phonons (blue dashed line).

A definitive fingerprint of the crucial role of the hot-phonon is provided by
the time evolution of the mean-square relative in-plane boron-boron
displacements $\sigma_{\rm B-B}$ as shown in Fig. 4a. For the assumption of a
thermal phonon distribution (red dashed line), $\sigma_{\rm B-B}^{2}$
increases monotonically obeying to a thermal heating of the lattice (similar
as $\sigma_{\rm B-Mg}^{2}$). Within the hot-phonon scenario, however, the in-
plane nearest neighbor mean-square relative lattice displacement $\sigma_{\rm
B-B}^{2}$ presents a remarkable peak at the highest $T_{\rm hot}$ achievable
at $t^{*}=50$ fs, with a following decrease towards a new equilibrium
conditions (red full line). Such a peak is a direct consequence of a dominance
of a larger population of the in-plane boron modes with the $E_{2g}$ symmetry.

The relevance of hot-phonon physics is revealed in a even more striking way in
the analysis of the lattice displacement correlation factor $\rho_{ij}$. As
shown in Fig. 4b, the onset of $E_{2g}$ hot phonons in MgB2 is reflected in a
direct way in a remarkable anomaly in the time evolution of the relative
boron-boron in-plane lattice correlation, with a sharp decrease of $\rho_{\rm
B-B}$ in the first 50 fs. The presence of hot-phonon modes with pure in-plane
boron character is visible also in the interplane boron-magnesium correlation
factor $\rho_{\rm B-Mg}$, although in a weaker way. Such anomalies are
peculiar of the hot-phonon scenario and are not predicted to show when the
lattice modes are assumed to obey a thermal distribution (dashed lines).

Our analysis provides thus a simple and straightforward way to probe the
presence of hot phonons in MgB2 directly in the analysis of the lattice
dynamics properties, as it can be revealed in time-resolved pump-probe
experiments. More precisely, we predict a marked anomaly in the time behavior
of lattice correlation factor $\rho_{\rm B-B}$ occuring at the delay time
$t^{*}$ when the hot-phonon lattice modes reach their largest population. The
time scale over which such anomaly disappears prompts also a possible direct
procedure for evaluating experimentally the time scale when all the lattice
degrees of freedom reach their final thermalization.

##  4 Wider perspective on other materials and experimental observations

The approach described in the present work suggests a possible way to detect
hot phonons not in an indirect way through the effects on the electronic
properties but directly in the lattice dynamics, exploiting the consequences
in the correlation of inter-atomic motion. In the specific case of MgB2, the
candidate modes for hot phonons are the ${\bf q}=0$ optical $E_{2g}$ ones that
correspond to out-of-phase in-plane motion of the two B atoms. The large non-
thermal population of these modes gives rise thus to a marked anomaly in the
in-plane correlation factor $\rho_{\rm B-B}$ with respect to the steady
conditions, as depicted in Fig. 4.

As we stressed above, such approach is not specific of MgB2 but it is quite
general and it can be usefully applied to other materials. The overall scheme,
for a given material candidate to hot-phonon physics, can be summarized as
follows ($i$) identifying the lattice modes with particularly large el-ph
coupling, likely thus to sustain hot-phonon physics; ($ii$) analyzing how the
eigenvectors of such modes affect the inter-atomic correlated motion; ($iii$)
thus investigating and predicting specific anomalies in the appropriate
correlation factors, validating such prediction with the comparison with a
overall thermal heating of the lattice with no hot phonons included. The most
straightforward generalization is probably graphene, where hot phonons have
been assessed in many ways, and which shares many similarities with MgB2. The
$E_{2g}$ modes of the hexagonal B plane MgB2 are the same $E_{2g}$ modes of
the hexagonal C plane that get hot in graphene. We can hence predict that a
similar anomaly can be observable also for graphene.

A similar approach is expected to work also for more complex materials with
several atoms in unit cell. A suitable class of compounds, along this
perspective, is the family of superconducting cuprates, as for instance the
representative system La2-xSrxCuO4 where the occurence of hot-phonon physics
has been already discussed in literature [65]. The lattice modes that get hot
in this compound are predicted to belong mainly to the $A_{1g}$ symmetry,
involving out-of-plane displacements either of the La atoms or of the apical
oxygen atoms [often labelled as O(3)] [64]. A residual el-ph coupling also
involves $E_{g}$ modes of the O(3) oxygens. A sketch of the relevant modes is
provided in Fig. 5. A time-resolved assessement of the correlated inter-atomic
motion can be a poweful tool also in this case since it could distinguish the
role of the different modes in the La-La or in the O($i$)-O($i$) inter-atomic
correlations ($i=1,2,3$). Note that the different contribution of the in-plane
or apical oxygens in the inter-atomic lattice correlations, as probed for
instance by the inter-atomic probability distribution function, could be
easily reveled due to the different bond lengths of the different oxygens [84,
85].

Figure 5: Main optical lattice modes of La2-xSrxCuO4 displaying sizable el-ph
coupling and which are thought to support hot-phonon physics. Sr is a random
substitutional atom for La and it is not shown here. La2-xSrxCuO4 has an
layered perovskite structure with apical-oxygen distance larger than the in-
plane oxygens.

Other more complex materials with several atoms for cell, sustaining hot-
phonon physics, could be in principle devised. In such a framework, on the
theoretical side, one cannot exclude that in-phase motions can be favoured,
resulting in a positive anomaly in a specific inter-atomic correlation
function, when the out-of-phase motion of a subset of atoms is accompanied by
an in-phase motion of another subset of atoms.

From the experimental point of view, it should be noticed that the most
accurate probes of the inter-atomic correlated motion under equilibrium
conditions are presently provided by neutron diffraction measurements that are
probably not best suited for time-resolved femtosecond resolution. Our
theoretical proposal does not however rely on a specific experimental setup.
It should be noticed indeed that, even in the steady neutron diffraction
measurements, raw data are collected in the momentum space, and later
translated in the real-space by Fourier transform. In this perspective
ultrafast electron diffuse scattering [86, 11, 53], as well as time-resolved
resonant inelastic X-ray scattering [87, 88], represent two promising
alternative techniques which are rapidly ramping up and constantly expanding
their accuracy and limits.

##  5 Conclusions

In this paper, using first-principle calculations, we have analyzed the
presence of the hot-phonon physics in MgB2 and its effect on the time
evolution of lattice dynamics in a typical pump-probe experiment. We have
shown that hot phonons in MgB2, with their characteristic $E_{2g}$ symmetry
that implies an in-plane counter-phase motion of the boron atoms, can be
directly traced down in the analysis of the atom-resolved mean-square lattice
displacements, and more evidently in the analysis of the mean square relative
lattice displacements $\sigma_{ij}^{2}$ and lattice correlation factors
$\rho_{ij}$. Even though we apply our investigation to the specific case of
MgB2, the analysis here described is quite generic and it can be extended to a
wide range of materials (semiconductors, metals or semimetals) sustaining hot
phonons, since they are commonly established at high-symmetry points of the
Brillouin zone, with well-defined and well-known symmetry and atomic contents.
Promising compounds might be for instance single-layer and multi-layer
graphene [59, 60], transition metal dichalcogenides [61], black phosphorous
[35, 55], hole-doped diamond [62, 63] and underdoped cuprates [64, 65]. Our
work provides thus a guideline for future direct experimental verification and
characterization of hot-phonon mechanisms in several families of compounds.

## Acknowledgements

D.N. acknowledges financial support from the Croatian Science Foundation
(Grant no. UIP-2019-04-6869) and from the European Regional Development Fund
for the “Center of Excellence for Advanced Materials and Sensing Devices”
(Grant No. KK.01.1.1.01.0001).

## References

  * [1] L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, H. Eisaki and M. Wolf,  _Ultrafast Electron Relaxation in Superconducting ${\mathrm{Bi}}_{2}{\mathrm{Sr}}_{2}{\mathrm{CaCu}}_{2}{\mathrm{O}}_{8+\delta}$ by Time-Resolved Photoelectron Spectroscopy_,  Phys. Rev. Lett. 99, 197001 (2007),  10.1103/PhysRevLett.99.197001.
  * [2] K. Ishioka, M. Hase, M. Kitajima, L. Wirtz, A. Rubio and H. Petek,  _Ultrafast electron-phonon decoupling in graphite_ ,  Phys. Rev. B 77, 121402 (2008),  10.1103/PhysRevB.77.121402.
  * [3] H. Yan, D. Song, K. F. Mak, I. Chatzakis, J. Maultzsch and T. F. Heinz,  _Time-resolved Raman spectroscopy of optical phonons in graphite: Phonon anharmonic coupling and anomalous stiffening_ ,  Phys. Rev. B 80, 121403 (2009),  10.1103/PhysRevB.80.121403.
  * [4] C. Gadermaier, A. S. Alexandrov, V. V. Kabanov, P. Kusar, T. Mertelj, X. Yao, C. Manzoni, D. Brida, G. Cerullo and D. Mihailovic,  _Electron-Phonon Coupling in High-Temperature Cuprate Superconductors Determined from Electron Relaxation Rates_ ,  Phys. Rev. Lett. 105, 257001 (2010),  10.1103/PhysRevLett.105.257001.
  * [5] R. Ulbricht, E. Hendry, J. Shan, T. F. Heinz and M. Bonn,  _Carrier dynamics in semiconductors studied with time-resolved terahertz spectroscopy_ ,  Rev. Mod. Phys. 83, 543 (2011),  10.1103/RevModPhys.83.543.
  * [6] R. Cortés, L. Rettig, Y. Yoshida, H. Eisaki, M. Wolf and U. Bovensiepen,  _Momentum-Resolved Ultrafast Electron Dynamics in Superconducting ${\mathrm{Bi}}_{2}{\mathrm{Sr}}_{2}{\mathrm{CaCu}}_{2}{\mathbf{O}}_{8+\delta}$_,  Phys. Rev. Lett. 107, 097002 (2011),  10.1103/PhysRevLett.107.097002.
  * [7] D. Rudolf, C. La-O-Vorakiat, M. Battiato, R. Adam, J. M. Shaw, E. Turgut, P. Maldonado, S. Mathias, P. Grychtol, H. T. Nembach, T. J. Silva, M. Aeschlimann _et al._ ,  _Ultrafast magnetization enhancement in metallic multilayers driven by superdiffusive spin current_ ,  Nature Communications 3(1) (2012),  10.1038/ncomms2029.
  * [8] B. Arnaud and Y. Giret,  _Electron Cooling and Debye-Waller Effect in Photoexcited Bismuth_ ,  Phys. Rev. Lett. 110, 016405 (2013),  10.1103/PhysRevLett.110.016405.
  * [9] J. C. Johannsen, S. Ulstrup, F. Cilento, A. Crepaldi, M. Zacchigna, C. Cacho, I. C. E. Turcu, E. Springate, F. Fromm, C. Raidel, T. Seyller, F. Parmigiani _et al._ ,  _Direct View of Hot Carrier Dynamics in Graphene_ ,  Phys. Rev. Lett. 111, 027403 (2013),  10.1103/PhysRevLett.111.027403.
  * [10] S. A. Jensen, Z. Mics, I. Ivanov, H. S. Varol, D. Turchinovich, F. H. L. Koppens, M. Bonn and K. J. Tielrooij,  _Competing Ultrafast Energy Relaxation Pathways in Photoexcited Graphene_ ,  Nano Letters 14(10), 5839 (2014),  10.1021/nl502740g.
  * [11] T. Chase, M. Trigo, A. H. Reid, R. Li, T. Vecchione, X. Shen, S. Weathersby, R. Coffee, N. Hartmann, D. A. Reis, X. J. Wang and H. A. Dürr,  _Ultrafast electron diffraction from non-equilibrium phonons in femtosecond laser heated Au films_ ,  Applied Physics Letters 108(4), 041909 (2016),  10.1063/1.4940981.
  * [12] L. Waldecker, R. Bertoni, R. Ernstorfer and J. Vorberger,  _Electron-Phonon Coupling and Energy Flow in a Simple Metal beyond the Two-Temperature Approximation_ ,  Phys. Rev. X 6, 021003 (2016),  10.1103/PhysRevX.6.021003.
  * [13] L. P. René de Cotret, J.-H. Pöhls, M. J. Stern, M. R. Otto, M. Sutton and B. J. Siwick,  _Time- and momentum-resolved phonon population dynamics with ultrafast electron diffuse scattering_ ,  Phys. Rev. B 100, 214115 (2019),  10.1103/PhysRevB.100.214115.
  * [14] M. Sentef, A. F. Kemper, B. Moritz, J. K. Freericks, Z.-X. Shen and T. P. Devereaux,  _Examining Electron-Boson Coupling Using Time-Resolved Spectroscopy_ ,  Phys. Rev. X 3, 041033 (2013),  10.1103/PhysRevX.3.041033.
  * [15] P. J. Price,  _Hot phonon effects in heterolayers_ ,  Superlattices and Microstructures 1(3), 255 (1985),  https://doi.org/10.1016/0749-6036(85)90013-8.
  * [16] P. Kocevar,  _Hot phonon dynamics_ ,  Physica B+C 134(1), 155 (1985),  https://doi.org/10.1016/0378-4363(85)90336-5.
  * [17] W. Pötz,  _Hot-phonon effects in bulk GaAs_ ,  Phys. Rev. B 36, 5016 (1987),  10.1103/PhysRevB.36.5016.
  * [18] R. P. Joshi and D. K. Ferry,  _Hot-phonon effects and interband relaxation processes in photoexcited GaAs quantum wells_ ,  Phys. Rev. B 39, 1180 (1989),  10.1103/PhysRevB.39.1180.
  * [19] P. Langot, N. Del Fatti, D. Christofilos, R. Tommasi and F. Vallée,  _Femtosecond investigation of the hot-phonon effect in GaAs at room temperature_ ,  Phys. Rev. B 54, 14487 (1996),  10.1103/PhysRevB.54.14487.
  * [20] S. Butscher, F. Milde, M. Hirtschulz, E. Malić and A. Knorr,  _Hot electron relaxation and phonon dynamics in graphene_ ,  Applied Physics Letters 91(20), 203103 (2007),  10.1063/1.2809413.
  * [21] K. Kang, D. Abdula, D. G. Cahill and M. Shim,  _Lifetimes of optical phonons in graphene and graphite by time-resolved incoherent anti-Stokes Raman scattering_ ,  Phys. Rev. B 81, 165405 (2010),  10.1103/PhysRevB.81.165405.
  * [22] S. Berciaud, M. Y. Han, K. F. Mak, L. E. Brus, P. Kim and T. F. Heinz,  _Electron and Optical Phonon Temperatures in Electrically Biased Graphene_ ,  Phys. Rev. Lett. 104, 227401 (2010),  10.1103/PhysRevLett.104.227401.
  * [23] C. H. Lui, K. F. Mak, J. Shan and T. F. Heinz,  _Ultrafast Photoluminescence from Graphene_ ,  Phys. Rev. Lett. 105, 127404 (2010),  10.1103/PhysRevLett.105.127404.
  * [24] L. Huang, B. Gao, G. Hartland, M. Kelly and H. Xing,  _Ultrafast relaxation of hot optical phonons in monolayer and multilayer graphene on different substrates_ ,  Surface Science 605(17), 1657 (2011),  https://doi.org/10.1016/j.susc.2010.12.009.
  * [25] M. Scheuch, T. Kampfrath, M. Wolf, K. von Volkmann, C. Frischkorn and L. Perfetti,  _Temperature dependence of ultrafast phonon dynamics in graphite_ ,  Applied Physics Letters 99(21), 211908 (2011),  10.1063/1.3663867.
  * [26] S. Wu, W.-T. Liu, X. Liang, P. J. Schuck, F. Wang, Y. R. Shen and M. Salmeron,  _Hot Phonon Dynamics in Graphene_ ,  Nano Letters 12(11), 5495 (2012),  10.1021/nl301997r.
  * [27] D. C. Hannah, K. E. Brown, R. M. Young, M. R. Wasielewski, G. C. Schatz, D. T. Co and R. D. Schaller,  _Direct Measurement of Lattice Dynamics and Optical Phonon Excitation in Semiconductor Nanocrystals Using Femtosecond Stimulated Raman Spectroscopy_ ,  Phys. Rev. Lett. 111, 107401 (2013),  10.1103/PhysRevLett.111.107401.
  * [28] D. Golla, A. Brasington, B. J. LeRoy and A. Sandhu,  _Ultrafast relaxation of hot phonons in graphene-hBN heterostructures_ ,  APL Materials 5(5), 056101 (2017),  10.1063/1.4982738.
  * [29] Y. Yang, D. Ostrowski, R. France, K. Zhu, J. van de Lagemaat, J. Luther and M. Beard,  _Observation of a hot-phonon bottleneck in lead-iodide perovskites_ ,  Nat. Photon. 10, 53 (2016),  doi.org/10.1038/nphoton.2015.213.
  * [30] J. Koivistoinen, P. Myllyperkiö and M. Pettersson,  _Time-Resolved Coherent Anti-Stokes Raman Scattering of Graphene: Dephasing Dynamics of Optical Phonon_ ,  The Journal of Physical Chemistry Letters 8(17), 4108 (2017),  10.1021/acs.jpclett.7b01711.
  * [31] E. M. Hamham, J. M. Iglesias, E. Pascual, M. J. Martín and R. Rengel,  _Impact of the hot phonon effect on electronic transport in monolayer silicene_ ,  Journal of Physics D: Applied Physics 51(41), 415102 (2018),  10.1088/1361-6463/aad94c.
  * [32] C. C. S. Chan, K. Fan, H. Wang, Z. Huang, D. Novko, K. Yan, J. Xu, W. C. H. Choy, I. Lončarić and K. S. Wong,  _Uncovering the Electron-Phonon Interplay and Dynamical Energy-Dissipation Mechanisms of Hot Carriers in Hybrid Lead Halide Perovskites_ ,  Advanced Energy Materials 11(9), 2003071 (2021),  https://doi.org/10.1002/aenm.202003071.
  * [33] F. Sekiguchi, H. Hirori, G. Yumoto, A. Shimazaki, T. Nakamura, A. Wakamiya and Y. Kanemitsu,  _Enhancing the Hot-Phonon Bottleneck Effect in a Metal Halide Perovskite by Terahertz Phonon Excitation_ ,  Phys. Rev. Lett. 126, 077401 (2021),  10.1103/PhysRevLett.126.077401.
  * [34] X. Tong and M. Bernardi,  _Toward precise simulations of the coupled ultrafast dynamics of electrons and atomic vibrations in materials_ ,  Phys. Rev. Res. 3, 023072 (2021),  10.1103/PhysRevResearch.3.023072.
  * [35] H. Seiler, D. Zahn, M. Zacharias, P.-N. Hildebrandt, T. Vasileiadis, Y. W. Windsor, Y. Qi, C. Carbogno, C. Draxl, R. Ernstorfer and F. Caruso,  _Accessing the Anisotropic Nonthermal Phonon Populations in Black Phosphorus_ ,  Nano Letters 0(0), null (0),  10.1021/acs.nanolett.1c01786.
  * [36] E. Baldini, A. Mann, L. Benfatto, E. Cappelluti, A. Acocella, V. M. Silkin, S. V. Eremeev, A. B. Kuzmenko, S. Borroni, T. Tan, X. X. Xi, F. Zerbetto _et al._ ,  _Real-Time Observation of Phonon-Mediated $\sigma\text{$-$}\pi$ Interband Scattering in ${\mathrm{MgB}}_{2}$_,  Phys. Rev. Lett. 119, 097002 (2017),  10.1103/PhysRevLett.119.097002.
  * [37] D. Novko, F. Caruso, C. Draxl and E. Cappelluti,  _Ultrafast Hot Phonon Dynamics in ${\mathrm{MgB}}_{2}$ Driven by Anisotropic Electron-Phonon Coupling_,  Phys. Rev. Lett. 124, 077001 (2020),  10.1103/PhysRevLett.124.077001.
  * [38] A. Y. Liu, I. I. Mazin and J. Kortus,  _Beyond Eliashberg Superconductivity in ${\mathrm{MgB}}_{2}$: Anharmonicity, Two-Phonon Scattering, and Multiple Gaps_,  Phys. Rev. Lett. 87, 087005 (2001),  10.1103/PhysRevLett.87.087005.
  * [39] J. M. An and W. E. Pickett,  _Superconductivity of MgB 2: Covalent Bonds Driven Metallic_,  Phys. Rev. Lett. 86, 4366 (2001),  10.1103/PhysRevLett.86.4366.
  * [40] T. Yildirim, O. Gülseren, J. Lynn, C. Brown, T. Udovic, Q. Huang, N. Rogado, K. Regan, M. Hayward, J. Slusky, T. He, M. Haas _et al._ ,  _Giant Anharmonicity and Nonlinear Electron-Phonon Coupling in MgB 2: A Combined First-Principles Calculation and Neutron Scattering Study_,  Phys. Rev. Lett. 87, 037001 (2001),  10.1103/PhysRevLett.87.037001.
  * [41] H. J. Choi, D. Roundy, H. Sun, M. L. Cohen and S. G. Louie,  _First-principles calculation of the superconducting transition in ${\mathrm{MgB}}_{2}$ within the anisotropic Eliashberg formalism_,  Phys. Rev. B 66, 020513 (2002),  10.1103/PhysRevB.66.020513.
  * [42] H. J. Choi, D. Roundy, H. Sun, M. L. Cohen and S. G. Louie,  _The origin of the anomalous superconducting properties of MgB 2_,  Nature 418, 758 (2002),  10.1038/nature00898.
  * [43] Y. Kong, O. Dolgov, O. Jepsen and O. Andersen,  _Electron-phonon interaction in the normal and superconducting states of MgB 2_,  Phys. Rev. B 64, 020501 (2001),  10.1103/PhysRevB.64.020501.
  * [44] A. Golubov, J. Kortus, O. Dolgov, O. Jepsen, Y. Kong, O. Andersen, B. Gibson, K. Ahn and R. Kremer,  _Specific heat of MgB 2 in a one- and a two-band model from first-principles calculations_,  J. Phys.: Condens. Matter 14, 1353 (2002),  10.1088/0953-8984/14/6/320.
  * [45] R. Gonnelli, D. Daghero, G. Ummarino, V. Stepanov, J. Jun, S. Kazakov and J. Karpinski,  _Direct Evidence for Two-Band Superconductivity in MgB 2 Single Crystals from Directional Point-Contact Spectroscopy in Magnetic Fields_,  Phys. Rev. Lett. 89, 247004 (2002),  10.1103/PhysRevLett.89.247004.
  * [46] R. B. Versteeg, J. Zhu, P. Padmanabhan, C. Boguschewski, R. German, M. Goedecke, P. Becker and P. H. M. van Loosdrecht,  _A tunable time-resolved spontaneous Raman spectroscopy setup for probing ultrafast collective excitation and quasiparticle dynamics in quantum materials_ ,  Structural Dynamics 5(4), 044301 (2018),  10.1063/1.5037784.
  * [47] D. Novko,  _Nonadiabatic coupling effects in MgB 2 reexamined_,  Phys. Rev. B 98, 041112 (2018),  10.1103/PhysRevB.98.041112.
  * [48] F. Carbone, D.-S. Yang, E. Giannini and A. H. Zewail,  _Direct role of structural dynamics in electron-lattice coupling of superconducting cuprates_ ,  Proc. Nat. Ac. of Sci. 105, 20161 (2008),  10.1073/pnas.0811335106.
  * [49] F. Carbone, N. Gedik, J. Lorenzana and A. H. Zewail,  _Real-Time Observation of Cuprates Structural Dynamics by Ultrafast Electron Crystallography_ ,  Adv. Condens. Matter Phys. 2010, 958618 (2010),  doi.org/10.1155/2010/958618.
  * [50] R. P. Chatelain, V. R. Morrison, B. L. M. Klarenaar and B. J. Siwick,  _Coherent and Incoherent Electron-Phonon Coupling in Graphite Observed with Radio-Frequency Compressed Ultrafast Electron Diffraction_ ,  Phys. Rev. Lett. 113, 235502 (2014),  10.1103/PhysRevLett.113.235502.
  * [51] M. Harb, H. Enquist, A. Jurgilaitis, F. T. Tuyakova, A. N. Obraztsov and J. Larsson,  _Phonon-phonon interactions in photoexcited graphite studied by ultrafast electron diffraction_ ,  Phys. Rev. B 93, 104104 (2016),  10.1103/PhysRevB.93.104104.
  * [52] L. Waldecker, R. Bertoni, H. Hübener, T. Brumme, T. Vasileiadis, D. Zahn, A. Rubio and R. Ernstorfer,  _Momentum-Resolved View of Electron-Phonon Coupling in Multilayer ${\mathrm{WSe}}_{2}$_,  Phys. Rev. Lett. 119, 036803 (2017),  10.1103/PhysRevLett.119.036803.
  * [53] M. J. Stern, L. P. René de Cotret, M. R. Otto, R. P. Chatelain, J.-P. Boisvert, M. Sutton and B. J. Siwick,  _Mapping momentum-dependent electron-phonon coupling and nonequilibrium phonon dynamics with ultrafast electron diffuse scattering_ ,  Phys. Rev. B 97, 165416 (2018),  10.1103/PhysRevB.97.165416.
  * [54] T. Konstantinova, J. D. Rameau, A. H. Reid, O. Abdurazakov, L. Wu, R. Li, X. Shen, G. Gu, Y. Huang, L. Rettig, I. Avigo, M. Ligges _et al._ ,  _Nonequilibrium electron and lattice dynamics of strongly correlated Bi 2Sr2CaCu2O8+δ single crystals_,  Sci. Adv. 4, eaap7427 (2018),  10.1126/sciadv.aap7427.
  * [55] D. Zahn, P.-N. Hildebrandt, T. Vasileiadis, Y. W. Windsor, Y. Qi, H. Seiler and R. Ernstorfer,  _Anisotropic Nonequilibrium Lattice Dynamics of Black Phosphorus_ ,  Nano Letters 20(5), 3728 (2020),  10.1021/acs.nanolett.0c00734.
  * [56] M. R. Otto, J.-H. Pöhls, L. P. R. de Cotret, M. J. Stern, M. Sutton and B. J. Siwick,  _Mechanisms of electron-phonon coupling unraveled in momentum and time: The case of soft phonons in TiSe 2_,  Sci. Adv. 7(20), eabf2810 (2021),  10.1126/sciadv.abf2810.
  * [57] M. Zacharias, H. Seiler, F. Caruso, D. Zahn, F. Giustino, P. C. Kelires and R. Ernstorfer,  _Multiphonon diffuse scattering in solids from first principles: Application to layered crystals and two-dimensional materials_ ,  Phys. Rev. B 104, 205109 (2021),  10.1103/PhysRevB.104.205109.
  * [58] M. Zacharias, H. Seiler, F. Caruso, D. Zahn, F. Giustino, P. C. Kelires and R. Ernstorfer,  _Efficient First-Principles Methodology for the Calculation of the All-Phonon Inelastic Scattering in Solids_ ,  Phys. Rev. Lett. 127, 207401 (2021),  10.1103/PhysRevLett.127.207401.
  * [59] D. Novko and M. Kralj,  _Phonon-assisted processes in the ultraviolet-transient optical response of graphene_ ,  npj 2D Materials and Applications 3(1) (2019),  10.1038/s41699-019-0131-5.
  * [60] F. Caruso, D. Novko and C. Draxl,  _Photoemission signatures of nonequilibrium carrier dynamics from first principles_ ,  Phys. Rev. B 101, 035128 (2020),  10.1103/PhysRevB.101.035128.
  * [61] F. Caruso,  _Nonequilibrium Lattice Dynamics in Monolayer MoS 2_,  The Journal of Physical Chemistry Letters 12(6), 1734 (2021),  10.1021/acs.jpclett.0c03616.
  * [62] L. Boeri, J. Kortus and O. K. Andersen,  _Three-Dimensional ${\mathrm{M}\mathrm{g}\mathrm{B}}_{2}$-Type Superconductivity in Hole-Doped Diamond_,  Phys. Rev. Lett. 93, 237002 (2004),  10.1103/PhysRevLett.93.237002.
  * [63] F. Giustino, J. R. Yates, I. Souza, M. L. Cohen and S. G. Louie,  _Electron-Phonon Interaction via Electronic and Lattice Wannier Functions: Superconductivity in Boron-Doped Diamond Reexamined_ ,  Phys. Rev. Lett. 98, 047005 (2007),  10.1103/PhysRevLett.98.047005.
  * [64] B. Mansart, M. J. G. Cottet, G. F. Mancini, T. Jarlborg, S. B. Dugdale, S. L. Johnson, S. O. Mariager, C. J. Milne, P. Beaud, S. Grübel, J. A. Johnson, T. Kubacka _et al._ ,  _Temperature-dependent electron-phonon coupling in La 2-xSrxCuO4 probed by femtosecond x-ray diffraction_,  Phys. Rev. B 88, 054507 (2013),  10.1103/PhysRevB.88.054507.
  * [65] S. L. Johnson, M. Savoini, P. Beaud, G. Ingold, U. Staub, F. Carbone, L. Castiglioni, M. Hengsberger and J. Osterwalder,  _Watching ultrafast responses of structure and magnetism in condensed matter with momentum-resolved probes_ ,  Struct. Dyn. 4, 061506 (2017),  10.1063/1.4996176.
  * [66] J. Kortus, I. I. Mazin, K. D. Belashchenko, V. P. Antropov and L. L. Boyer,  _Superconductivity of Metallic Boron in ${\mathrm{MgB}}_{2}$_,  Phys. Rev. Lett. 86, 4656 (2001),  10.1103/PhysRevLett.86.4656.
  * [67] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo and et al.,  _QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials_ ,  Journal of Physics: Condensed Matter 21(39), 395502 (2009),  10.1088/0953-8984/21/39/395502.
  * [68] J. P. Perdew, K. Burke and M. Ernzerhof,  _Generalized Gradient Approximation Made Simple_ ,  Phys. Rev. Lett. 77, 3865 (1996),  10.1103/PhysRevLett.77.3865.
  * [69] S. Baroni, S. de Gironcoli, A. Dal Corso and P. Giannozzi,  _Phonons and related crystal properties from density-functional perturbation theory_ ,  Rev. Mod. Phys. 73, 515 (2001),  10.1103/RevModPhys.73.515.
  * [70] S. Poncé, E. Margine, C. Verdi and F. Giustino,  _EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions_ ,  Computer Physics Communications 209, 116 (2016),  https://doi.org/10.1016/j.cpc.2016.07.028.
  * [71] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza and D. Vanderbilt,  _Maximally localized Wannier functions: Theory and applications_ ,  Rev. Mod. Phys. 84, 1419 (2012),  10.1103/RevModPhys.84.1419.
  * [72] K.-P. Bohnen, R. Heid and B. Renker,  _Phonon Dispersion and Electron-Phonon Coupling in MgB 2 and AlB2_,  Phys. Rev. Lett. 86, 5771 (2001),  10.1103/PhysRevLett.86.5771.
  * [73] B. Renker, K. Bohnen, R. Heid, D. Ernst, H. Schober, M. Koza, P. Adelmann, P. Schweiss and T. Wolf,  _Strong Renormalization of Phonon Frequencies in Mg 1-xAlxB2 _,  Phys. Rev. Lett. 88, 067001 (2002),  10.1103/PhysRevLett.88.067001.
  * [74] A. Shukla, M. Calandra, M. d’Astuto, M. Lazzeri, F. Mauri, C. Bellin, M. Krisch, J. Karpinski, S. M. Kazakov, J. Jun, D. Daghero and K. Parlinski,  _Phonon Dispersion and Lifetimes in ${\mathrm{M}\mathrm{g}\mathrm{B}}_{2}$_,  Phys. Rev. Lett. 90, 095506 (2003),  10.1103/PhysRevLett.90.095506.
  * [75] M. d’Astuto, M. Calandra, S. Reich, A. Shukla, M. Lazzeri, F. Mauri, J. Karpinski, N. D. Zhigadlo, A. Bossak and M. Krisch,  _Weak anharmonic effects in $\mathrm{Mg}{\mathrm{B}}_{2}$: A comparative inelastic x-ray scattering and Raman study_,  Phys. Rev. B 75, 174508 (2007),  10.1103/PhysRevB.75.174508.
  * [76] P. B. Allen,  _Theory of thermal relaxation of electrons in metals_ ,  Phys. Rev. Lett. 59, 1460 (1987),  10.1103/PhysRevLett.59.1460.
  * [77] C. Lui, K. Mak, J. Shan and T. Heinz,  _Ultrafast Photoluminescence from Graphene_ ,  Phys. Rev. Lett. 105, 127404 (2010),  10.1103/PhysRevLett.105.127404.
  * [78] S. Dal Conte, C. Giannetti, G. Coslovich, F. Cilento, D. Bossini, T. Abebaw, F. Banfi, G. Ferrini, H. Eisaki, M. Greven, A. Damascelli, D. van der Marel _et al._ ,  _Disentangling the Electronic and Phononic Glue in a High-Tc Superconductor_ ,  Science 335(6076), 1600 (2012),  10.1126/science.1216765.
  * [79] D. Novko,  _First-principles study of ultrafast dynamics of Dirac plasmon in graphene_ ,  New Journal of Physics 23(4), 043023 (2021),  10.1088/1367-2630/abec0f.
  * [80] A. Shukla, M. Calandra, M. d’Astuto, M. Lazzeri, F. Mauri, C. Bellin, M. Krisch, J. Karpinski, S. M. Kazakov, J. Jun, D. Daghero and K. Parlinski,  _Phonon Dispersion and Lifetimes in ${\mathrm{M}\mathrm{g}\mathrm{B}}_{2}$_,  Phys. Rev. Lett. 90, 095506 (2003),  10.1103/PhysRevLett.90.095506.
  * [81] I.-K. Jeong, R. H. Heffner, M. J. Graf and S. J. L. Billinge,  _Lattice dynamics and correlated atomic motion from the atomic pair distribution function_ ,  Phys. Rev. B 67, 104301 (2003),  10.1103/PhysRevB.67.104301.
  * [82] G. Campi, E. Cappelluti, T. Proffen, X. Qiu, E. Bozin, S. Billinge, S. Agrestini, N. Saini and A. Bianconi,  _Study of temperature dependent atomic correlations in MgB 2_,  Eur. Phys. J. B 52, 15 (2006),  doi.org/10.1140/epjb/e2006-00269-7.
  * [83] The values of $\rho_{\rm B-B}$ at steady conditions computed here appears larger than the values experimentally estimated in Ref. [82]. One should consider however that we evaluate here only the contribution of the lattice dynamics to the lattice correlations $\sigma_{\rm B-B}$. The effects of static disorder was modelled in Ref. [82] by means of a constant correction term $\sigma^{2}_{0}(i_{\alpha})$ in the lattice mean-square displacements $\sigma^{2}(i_{\alpha})$, but, in the absence of a well-accepted modelling of the effects of disorder of the interatomic correlation, no effect of disorder was included in the analysis of $\sigma_{\rm B-B}$. We note that a constant upwards offshift of the theoretical $\sigma_{\rm B-B}$ would improve significantly the agreement between theoretical and experimental estimates of $\sigma_{\rm B-B}$, and bring thus the experimental estimation of $\rho_{\rm B-B}$ closer to the theoretical predictions $\rho_{\rm B-B}\approx 0.24$. Also, the presence of nonadiabatic effects and the possible breaking of the harmonic approximation might account for such discrepancy.
  * [84] R. E. Cohen, W. E. Pickett and H. Krakauer,  _First-Principles Phonon Calculations for ${\mathrm{La}}_{2}$Cu${\mathrm{O}}_{4}$_,  Phys. Rev. Lett. 62, 831 (1989),  10.1103/PhysRevLett.62.831.
  * [85] W. E. Pickett,  _Electronic structure of the high-temperature oxide superconductors_ ,  Rev. Mod. Phys. 61, 433 (1989),  10.1103/RevModPhys.61.433.
  * [86] M. Trigo, J. Chen, V. H. Vishwanath, Y. M. Sheu, T. Graber, R. Henning and D. A. Reis,  _Imaging nonequilibrium atomic vibrations with x-ray diffuse scattering_ ,  Phys. Rev. B 82, 235205 (2010),  10.1103/PhysRevB.82.235205.
  * [87] L. J. P. Ament, M. van Veenendaal, T. P. Devereaux, J. P. Hill and J. van den Brink,  _Resonant inelastic x-ray scattering studies of elementary excitations_ ,  Rev. Mod. Phys. 83, 705 (2011),  10.1103/RevModPhys.83.705.
  * [88] M. Mitrano and Y. Wang,  _Probing light-driven quantum materials with ultrafast resonant inelastic X-ray scattering_ ,  Comm. Phys. 3 (2020),  10.1038/s42005-020-00447-6.


# fplyr: the split-apply-combine strategy for big data in R

by Federico Marotta

###### Abstract

We present fplyr, a new package for the R language to deal with big files. It
allows users to easily implement the split-apply-combine strategy for files
that are too big to fit into the available memory, without relying on data
bases nor introducing non-native R classes. A custom function can be applied
independently to each group of observations, and the results may be either
returned or directly printed to one or more output files.

## Introduction

Many fields of science and industry are witnessing an expansion in the volume
of data collected by their statistical experiments, and the size of the files
to be analysed is growing accordingly. On the other hand, the hardware
resources needed to scale up the analyses are not available to every
practitioner, and sometimes are even beyond the current technological reach.
In particular, it can happen that the size of a file exceeds the total RAM on
the machine where the analysis is to be performed. In such cases the R (R Core
Team, 2020) programmer has several tools at his or her disposal, each with its
own advantages and limitations. The existing tools can be broadly classified
in three families: database-backed, file-backed, and MapReduce. In the next
few paragraphs, these three approaches will be briefly reviewed; the purpose,
however, is not to compare them, as each is suited to different situations.

To the first family, database-backed, belong packages such as sqldf
(Grothendieck, 2017) and Rdbi, which set up a database behind the scenes and
perform SQL queries on it without loading the whole file into memory. While
this approach can leverage the speed and efficiency of database routines, it
is sometimes necessary to perform on the data particular operations that are
not entailed by the database software. For example, while it is possible to
use the groupby SQL keyword to aggregate the results of a query according to
the value of one field, only a small number of simple ‘aggregating operations’
are supported; the exact number and type of operations can differ according to
the database backend, but in general they are not much more sophisticated than
summary statistics like AVG, MAX, or STDEV.

A filebacking strategy is employed in the bigmemory (Kane et al., 2013) and ff
(Adler et al., 2014) packages, which rely on low-level operating system
functions to create a virtual map of the file, once again without loading it
into memory. These packages are highly optimised, but their efficiency comes
with tradeoffs; for instance, bigmemory is restricted to numerical matrices,
and with ff only a (wide but) limited set of preimplemented operations are
available. The ffbase (De Jonge et al., 2014) considerably enlarges the scope
of action of ff. One potential disadvantage that remains, however, is that ff
and related packages implement their own set of classes, rather than using the
native ones offered by R.

The iotools package (Arnold et al., 2015) offers both an efficient way to read
a file from disk and a set of functions to parse the file chunk by chunk.
These latter functions, such as chunk.apply() and chunk.map(), are reminescent
of the MapReduce paradigm, which gives the name to the last family of our
classification. Here, the large file is read piece by piece in such a way
that, at any given moment, only a limited number of rows of the file are
present in the RAM; at the same time, an arbitrary function can be applied to
each ‘chunk’. The results of the processing are then combined and returned.

In this paper we propose a new approach, implemented in the fplyr pacakge, to
deal with big files; it aims to be a simple and user-friendly solution that
integrates well with the existing R functions. Nonetheless, as all of the
approaches described above, it also has a restricted applicability, as it
deals only with a particular class of big files: those that are amenable to
the split-apply-combine approach (Wickham, 2011). For example, the famous iris
data set (Anderson, 1936; Fisher, 1936) would be a good candidate because the
‘Species’ field defines a partition of the data into independent groups of
observations, so that a particular function can be applied separately to each
species. Another example is that of gene expression files, which are often
organised as matrices where each row refers to a gene, and each column reports
the expression of that gene in a different individual. In this case, one may
apply a function separately to each gene. More specifically, the files on
which fplyr operates should be formatted in such a way that consecutive rows
contain all the measurements relative to the same subject, and the first field
contains the subject IDs. We refer to each group of consecutive rows
pertaining to the same subject as a block.

Much as apply() applies a function to each row or column of a matrix, the
functions in fplyr apply custom functions to each block of data, independently
of all the other blocks. Thus, at its core, this package enables one to mimick
the behaviour of the by() function, without requiring that the whole file be
loaded into memory. Indeed, the file being read block by block, fplyr’s
functions run with an $O(1)$ space complexity.

As an illustration of the possible usage of the package, suppose that the path
to a big file to be processed is stored as a character string in the variable
f. Then, the following code computes and returns the summary() of each block:

flply(f, summary)

## Comparison with existing packages

The approach we presented may appear similar to iotools’ chunk-wise
processing, but it differs from it in several respects. Most importantly, one
of iotools’ chunks can contain measurements about many subjects, while an
fplyr’s block contains each and every measurement about one subject only. In
other words, one of iotools’ chunks may contain many of fplyr’s blocks.
Moreover, in fplyr, each block is treated as independent of all the others;
the aim is to obtain a list of values, one for each block, whereas with
iotools the typical aim is to obtain a single value for the whole file.

Similar differences distinguish fplyr by ffbase’s way of performing operations
by chunks, a function called ffdfdply(). As with iotools, here a chunk can
contain measurements from many subjects, and the task of further separating
them is left to the user. Moreover, ffdfdply() can only return an "ffdf", the
equivalent of a "data.frame" for ff, whereas with fplyr it is possible to
return any R object, or even to directly write something to an output file
block by block, during the processing.

As previously stated, aggregate operations with database-backed packages are
limited to simple operations. If the operation to be applied is more complex,
it is still possible to use sqldf, but the analysis must be performed in at
least two steps: first the relevant group of observations must be selected,
and only then can the custom function be applied. In order to replicate the
behaviour of the groupby keyword with arbitrary aggregating functions, a loop
should be manually set up where at each iteration a different group of
observations is retrieved and analysed; furthermore, all the possible values
of the ‘groupby’ field must be known in advance. fplyr offers an effortless
way of doing the same thing.

## Implementation

fplyr is mainly built upon two other R packages: iotools and data.table, the
former being used to read the file chunk by chunk, the latter providing some
efficient ways to split the chunk into its constituent blocks, applying the
user-specified function to each, combining the results, and (when applicable)
writing back the result to the disk (Dowle and Srinivasan, 2019). From this
point onwards, for the sake of distinguishing the user-specified function to
be applied to each block from other functions, we shall refer to it as FUN.
fplyr’s algorithm is, in essence, as follows: one chunk is read with the help
of the iotools package, then it is split into its constituent blocks using the
by() function, and FUN is applied to each. Once the whole file has been
processed, it is returned a list where each element corresponds to a block.
This algorithm, with additional technicalities, is implemented in the flply()
function. However, if the output of FUN is a "data.frame", the data.table
package allows us to replace the base R by() function with a faster
alternative; this second algorithm is implemented in the ftply() function (see
The names of the functions for an interpretation of the names of the functions
in this package). Since the user knows _a priori_ the type of output that FUN
returns, he or she can choose which function to use accordingly. Note,
however, that the second approach is implemented only as a shortcut, because
the same result could be achieved using by() followed by rbind(), albeit it
would be achieved much more slowly. In the following paragraphs we shall
discuss several aspects of the implementation.

One of the reasons why iotools is faster than base R functions in reading a
file from disk is that it reads the file in bulk as a "raw" vector and only
afterwards does it convert the data into a suitable structure, like a
"data.frame" or a "matrix". In fact, for instance, the mstrsplit() function
takes a "raw" or "character" vector and arranges it into a matrix, and
dstrsplit() returns a "data.frame" instead. Functions like mstrsplit() and
dstrsplit() are called formatters in iotools jargon (Arnold et al., 2015). In
fplyr, we defined a new formatter which takes a "raw" vector and returns a
"data.table". In particular, it first converts the "raw" vector to a
"character" one, and then feeds the "character" vector to data.table’s
fread(). (Indeed, not only can fread() take the path to a file as input, but
it can also handle character strings directly and cast them into
"data.table"s.) By relying on it for the low-level reading of the files, fplyr
inherits some of iotools’ features, such as the ability to read compressed
files. Incidentally, benchmarks show that our data.table-based formatter,
dtstrsplit(), is even faster than iotools’ builtin formatter (Figure 1).

Figure 1: Benchmark of the two formatters: fplyr’s, called dtstrsplit, and
iotools’, called dstrsplit. The two were used to read the same file (253316
lines, 9.2 MB) for 200 times without performing any other operation. As a
reference, the time taken by data.table’s fread() is also reported. The
benchmark was performed on an Acer Aspire 5750G, 2.4GHz Intel Core i5 with 4
cores, 4GB DDR3 RAM. The code to reproduce the benchmarks is available on the
GitHub repository of the package at https://github.com/fmarotta/fplyr.

Once a chunk has been read and formatted to a "data.table", it must be split
into its blocks and FUN must be applied to each of them. According to the type
of FUN’s output, as we mentioned, two different strategies are employed. If
FUN returns a "data.table", then we use the following construct:

output <\- chunk[, FUN(.SD, .BY, ...), by = blocks]

In this case, FUN must take at least two arguments: at evaluation time, when
FUN is called, a "data.table" containing one block (except the first field) is
passed to FUN as the first argument, and a "character" vector with the name of
the block is passed as the second argument. Additional arguments can be
specified by the user, much as in the apply() family. It is up to the user to
write FUN satisfying these specifications.

On the other hand, if FUN returns something different than a "data.table",
then by() is used and the result is returned as a list where each element
corresponds to a block. In this case, FUN must take at least one argument: at
evaluation time, the whole block, including the first field, is passed to FUN.
Once again, it is up to the user to write or find a function that acts on each
block and satisfies this specification, and additional arguments can be passed
by the user like in the previous case.

Due to the mutual independence of the blocks, the processing of a file can be
easily parallelised111The parallelisation is only possible on *nix operating
systems.; indeed, all the functions in fplyr support the parallel argument to
specify the number of workers to be initialised. In particular, we adopted
iotools’ pipeline parallelism (Arnold et al., 2015), where the master process
reads from and writes to the disk sequentially, but each chunk is then
relegated to the workers, which process it in parallel as described
previously; while the workers are busy, the next chunk is pre-allocated by the
master. The sequential reading and writing avoids conflicts—for instance if
multiple processes try to write to the same file at the same time—but can
result in bottlenecks if the processing of each chunk takes an amount of time
that differs too much from the time needed for the reading and writing.
Furthermore, some workers can become idle if the amount of time required to
process different chunks is different.

For convenience, we also implemented two additional functions which may be
useful when even the result of the processing of the file would not fit into
the RAM. In these cases, it is not possible to simply return everything at the
end. If the output of FUN is a "data.table", one solution could be to print
the resulting data for each block to an output file as soon as it is ready,
proceeding to append the results of all the blocks to the same file. This
strategy is implemented in the ffply() function. The last function, fmply(),
combines the behaviours of ffply() and flply(): for each block, it is possible
to write to one or more files and, at the same time, to return an arbitrary R
object.

## Examples

To illustrate the usefulness of the package, we shall now discuss some
examples. Throughout this section, we assume that the path to the file of
interest is stored in a variable called fin, while the path to an output file
is stored inside fout.

Tabular data with multiple measurements can be represented in two ways: the
long and the wide format (Wickham, 2007). These two representations are
equivalent in that data can always be converted from one to the other, but
when the data are in long format, they are ‘tidy’ and easier to analyse
(Wickham, 2014), and some R functions require their input to be in this
format. However, if some big file is only available in wide format, reshaping
it could become a problem due to the scarcity of RAM. Indeed, even if the
wide-formatted file fits into the memory, the long one may not. With fplyr the
reshaping can be performed block by block:

ffply(fin, fout, function(d, by) melt(d, measure.vars = names(d))

Here the by argument, which contains the subject ID, was ignored, but in
principle it can be used inside an if condition to select only some of the
blocks. The d argument contains one whole block _without the first column_ ,
and it is reshaped using data.table’s melt() function. The first column,
containing the IDs, will be automatically added after FUN has been applied.

Sometimes a group of observations must be analysed in many ways, producing
different outputs. For instance, when performing a linear regression by
blocks, it may be convenient to print the coefficients to an output file,
while at the same time returning the full "lm" object for future in-depth
inspection. The fmply() function takes as arguments the path to the input
file, a vector of paths to the (possibly many) output files, and a function
that returns a list of objects. If there are, say, d output files, the first d
elements of the list must be "data.table"s (or "data.frame"s) and are printed
to the corresponding output files; optionally, FUN can return d+1 elements, in
which case the last one is returned by fmply(). In this example, the
coefficients are written to fout, and the "lm" objects themselves are returned
at the end. In particular, l will be a list where each element is an "lm"
object corresponding to one block.

l <\- fmply(fin, fout, function(d) { lm.fit <\- lm(Y ~ ., data = d[, -1]) #
Add the name of the block as the first field lm.coef <\-
as.data.table(cbind(d[1, 1], t(coef(lm.fit)))) # The coefficients will be
printed, the fitted object will be returned return(list(lm.coef, lm.fit))})

Additional examples can be found in the package vignette.

## The names of the functions

Despite the names of this package is reminescent of other packages belonging
to the tidyverse (Wickham et al., 2019), fplyr bears no relation with them. We
did, however, follow the same conventions of one of Hadley Wickham’s packages
for naming the functions (Wickham, 2011). All the names consist of two letters
followed by ‘ply’: the first letter represent the type of input, whereas the
second letter characterises the type of output, and the final ‘ply’ clinches
the relation with the existing ‘apply’ family of functions. The first letter
is usually ‘f’, because the input is the path to a file. The second letter is
‘l’ if the output is a list, as in flply(), it is ‘t’ if the output is a
"data.table", ‘f’ if it is another file, and ‘m’ if it can be multiple things.

## Discussion

We believe that fplyr fills a gap in the landscape of the existing tools to
process big files in R. It adresses a problem that in principle could be
solved by other packages as well, but only with workarounds. Furthermore, our
implementation combines the strengths of two other packages, iotools and
data.table, and is therefore reasonably efficient. A variety of features, such
as the transparent parallelisation, the ability to read compressed files, and
the possibility to specify the maximum number of blocks, make it also user-
friendly.

There are, however, also some limitations. First and foremost, the file on
which fplyr operates must contain observations which can be assigned to
several independent ‘subjects’, and the subject IDs can only be in the first
field, the reason being that this is also how iotools works. Although it is
possible to pre-process the file with *nix command-line tools such as awk and
sort to ensure that the IDs are in the first column and that the rows
referring to the same subject are adjacent, possible future work could extend
the package in order to support a custom field for the subject IDs. A related
extension could allow the blocks to be defined by the combined values of two
or more columns.

Another possible weakness is the pipeline parallelism algorithm, which,
besides not being available for Windows-based operating systems, can cause
bottlenecks if the times required for the operation on the block on the one
hand, and for reading the block on the other, are on very different scales.
Nevertheless, this algorithm has been proven to be more efficient in a variety
of common situations (Arnold et al., 2015).

In summary, the main strength of the package is that it simplifies the task of
splitting, applying and combining for files too big to fit in the available
memory.

## URLs

  * •

GitHub repository

  * •

CRAN page of the package

  * •

Package vignette

## Aknowledgements

The author would like to thank Prof. Paolo Provero (Università degli Studi di
Torino) for his insightful comments and suggestions.

## References

  * Adler et al. (2014) D. Adler, C. Gläser, O. Nenadic, J. Oehlschlägel, and W. Zucchini.  _ff: Memory-efficient Storage of Large Data on Disk and Fast Access Functions_ , 2014.  URL https://CRAN.R-project.org/package=ff.  R package version 2.2-14.
  * Anderson (1936) E. Anderson.  The species problem in iris.  _Annals of the Missouri Botanical Garden_ , 23(3):457–509, 1936.
  * Arnold et al. (2015) T. Arnold, M. Kane, and S. Urbanek.  iotools: High-performance i/o tools for r.  _arXiv preprint arXiv:1510.00041_ , 2015.
  * De Jonge et al. (2014) E. De Jonge, J. Wijffels, and J. van der Laan.  _ffbase: Basic Statistical Functions for Package ff. R package version 0.11. 3_ , 2014.  URL https://CRAN.R-project.org/package=ffbase.  R package version 0.12.7.
  * Dowle and Srinivasan (2019) M. Dowle and A. Srinivasan.  _data.table: Extension of ‘data.frame‘_ , 2019.  URL https://CRAN.R-project.org/package=data.table.  R package version 1.12.8.
  * Fisher (1936) R. A. Fisher.  The use of multiple measurements in taxonomic problems.  _Annals of eugenics_ , 7(2):179–188, 1936.
  * Grothendieck (2017) G. Grothendieck.  _sqldf: Manipulate R Data Frames Using SQL_ , 2017.  URL https://CRAN.R-project.org/package=sqldf.  R package version 0.4-11.
  * Kane et al. (2013) M. J. Kane, J. Emerson, and S. Weston.  Scalable strategies for computing with massive data.  _Journal of Statistical Software_ , 55(14):1–19, 2013.  URL http://www.jstatsoft.org/v55/i14/.
  * R Core Team (2020) R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2020.  URL https://www.R-project.org/.
  * Wickham (2007) H. Wickham.  Reshaping data with the reshape package.  _Journal of statistical software_ , 21(12):1–20, 2007.
  * Wickham (2011) H. Wickham.  The split-apply-combine strategy for data analysis.  _Journal of Statistical Software_ , 40(1):1–29, 2011.
  * Wickham (2014) H. Wickham.  Tidy data.  _Journal of Statistical Software_ , 59(10):1–23, 2014.
  * Wickham et al. (2019) H. Wickham, M. Averick, J. Bryan, W. Chang, L. McGowan, R. François, G. Grolemund, A. Hayes, L. Henry, J. Hester, et al.  Welcome to the tidyverse.  _Journal of Open Source Software_ , 4(43):1686, 2019.

_Federico Marotta
Università degli Studi di Torino
I-10126, Torino (TO), Italy
ORCiD: 0000-0002-0174-3901
<EMAIL_ADDRESS>