# Ice phonon spectra and Bayes inference: a gateway to a new understanding of
terahertz sound propagation in water

Alessio De Francesco<EMAIL_ADDRESS>CNR-IOM & INSIDE@ILL c/o Operative Group
in Grenoble (OGG), F-38042 and Institut Laue Langevin, Grenoble, France  Luisa
Scaccia  Dipartimento di Economia e Diritto, Università di Macerata, Via
Crescimbeni 20, 62100 Macerata, Italy  Ferdinando Formisano  CNR-IOM &
INSIDE@ILL c/o Operative Group in Grenoble (OGG), F-38042 and Institut Laue
Langevin, Grenoble, France  Eleonora Guarini  Dipartimento di Fisica e
Astronomia, Università di Firenze, via G. Sansone 1, I-50019 Sesto Fiorentino,
Italy  Ubaldo Bafile  Consiglio Nazionale delle Ricerche, Istituto di Fisica
Applicata ”Nello Carrara”, via Madonna del Piano 10, I-50019 Sesto Fiorentino,
Italy  Ahmet Alatas  Argonne National Laboratory, Advanced Photon Source, P.O.
Box 5000 Upton, 11973 NY, USA  Scott T. Lynch  Department of Physics,
University of Wisconsin at Madison, 1150 University Avenue, Madison, WI, USA
Alessandro Cunsolo  Department of Physics, University of Wisconsin at Madison,
1150 University Avenue, Madison, WI, USA

###### Abstract

Understanding how molecules engage in collective motions in a liquid where a
network of bonds exists has both fundamental and applied relevance. On the one
hand, it can elucidate the “ordering” role of long-range correlations in an
otherwise strongly dissipative system; on the other hand, it can inspire new
avenues to control such order to implement sound manipulation. Water
represents an ideal investigation case to unfold these general aspects and,
across the decades, it has been the focus of thorough scrutiny. Despite this
investigative effort, the spectrum of terahertz density fluctuations of water
largely remains a puzzle for Condensed Matter physicists. To unravel it, we
compare previous scattering measurements of water spectra with new ones on
ice. Thanks to the unique asset of Bayesian inference, we draw a more detailed
portrayal of the phonon response of ice. The comparison with the one of liquid
water challenges the current understanding of density fluctuations in water,
or more in general, of any networked liquid.

##  I Introduction

It is a matter of common experience that liquids and solids oppose a different
resistance to mechanical attempts to change their macroscopic shape. A
possible way to test such a different rigidity is to induce a macroscopic
mechanical perturbation generating the propagation of a density wave. An
alternative pathway followed in inelastic scattering measurements consists in
stimulating the propagation of density fluctuations at mesoscopic scales by
irradiating a material with a beam of particle waves, such as neutrons or
photons. When dealing with a liquid system, such a density wave propagates
with a strong damping and at macroscopic scales, uniquely in a direction
parallel to the triggering force. In a structurally ordered solid such as a
crystal, applied stresses propagate for a longer time and in directions both
parallel and orthogonal to the applied force, in the form of longitudinal (LA)
or transverse acoustic (TA) phonons, respectively. However, this clear-cut
distinction becomes more elusive over distances and time lapses respectively
approaching the size of a single atom first neighbors cage and the period of
such an atom’s “in cage” bounces.

Indeed, abundant experimental evidence endorses the conclusion that density
waves in liquids bear strong evidence for terahertz viscoelasticity [1], i.e.,
they combine liquid-like and solid-like aspects [2, 3, 4, 5]. For instance,
the transition between the viscous and the elastic regime in liquid water at
about 280 K occurs at a few nm-1, according to the combined results of Refs.
[6] and [7]. Relevant points of the broad picture emerging from the intensive
scrutiny of terahertz viscoelasticity in fluids can be summarized as follows:

1) The first spectacular manifestation of mesoscale viscoelasticity is the
increase of sound velocity when probing distances and times roughly as small
as those featuring first-neighbor atomic interactions. This trend reflects the
enhanced rigidity of the system (elastic regime) at these so-called mesoscopic
scales, typically much shorter than those covered by dissipative - diffusion
or relaxation - processes in liquids.

2) The transition to such an elastic regime is paralleled by an increased
ability of the liquid to support the propagation of transverse acoustic waves.

3) As opposed to a long-standing belief, this viscoelasticity persists in
thermodynamic domains extending well above the critical point [8, 9], although
their actual delimitation raised some controversy [10, 11, 12, 13].

4) Aside from its unsettled dependence on thermodynamic conditions, high-
frequency viscoelasticity also has a still poorly understood link with the
nature of microscopic interaction; for instance, it is unclear why it is
substantially more pronounced in specific relatively low-viscosity systems as
water [14, 15, 3], carbon dioxide [16], deuterium [17], and noble gases [18,
19, 20], while it is often barely detectable in liquid metals [21, 22, 23, 24,
25].

Far from being a mere academic endeavor, gaining insight on these topics could
deepen our understanding of subjects as fundamental as the very nature of
liquid aggregation. Also, it will likewise disclose new avenues in the
emerging domain of sound propagation design and engineering, a field whose
enormous practical interest has been recognized [26], yet not fully explored.

On the experimental side, terahertz phonon propagation in condensed materials
can be directly probed by inelastic scattering methods, such as Inelastic
X-Ray (IXS [27, 28]) and Neutron Scattering (INS, [29, 30]). Conceptually, IXS
or INS spectroscopic probes resemble large microscopes “pointed on the
dynamics”, which can be zoomed in to focus on dynamic events occurring over
various scales. This goal is ordinarily achieved by suitable tuning of the two
main variables of the scattering event, namely energy and momentum exchanged
between the probe wave particles (here assumed to be photons), respectively
referred to as $E=\hbar\omega$, and $\bf{p}=\hbar\bf{Q}$. Here, $\hbar$ is the
reduced Planck constant, while $\bf{Q}$ ($\omega$) represents the exchanged
wave-vector (frequency), i.e., the difference between before-scattering and
after-scattering photon wave-vector (frequency). For small enough values of
$Q=|\bf{Q}|$ and $E$, the target system appears as a continuum whose dynamic
response is probed as an average over many microscopic events. Conversely, at
extremely large $(Q,E)$’s, distances and times spanned become so small that
the only event observed is the free recoil of the single struck atom after the
collision with the photon and before subsequent first-neighbor interactions.
Clearly, all dynamic processes characterizing the interatomic dynamics of the
sample can be captured, somewhere between these two limits, by adequately
tuning the spectroscopic probe. Not surprisingly, the developments of the
scattering techniques directly mapping this broad dynamic domain, INS and IXS,
have dramatically improved the current understanding of high-frequency
viscoelasticity phenomena in an impressively varied family of disordered
systems. Although water is one of the disordered systems most thoroughly
studied across various decades, its deep scrutiny has often evidenced
seemingly anomalous and poorly understood behaviors, albeit sometimes later
recognized as “normal effects” that water shares with many other liquids.
Notable examples include the high propagation speed of the terahertz sound
mode, once called “fast sound” [6, 7], and the onset of a shear mode
propagation [31, 32, 33, 34, 35]. Although similar phenomena are being
reported for an increasing number of disordered systems, the case of water
found a straightforward and broadly accepted explanation in terms of hydrogen
bond (HB) dynamics. According to this interpretative scheme, [36, 37] high
frequency transverse and longitudinal acoustic waves in water respectively
couple with the bending of two HBs linking triplets of oxygen atoms, or to the
stretching of the HBs connecting two adjacent oxygens. Despite the remarkable
insight achieved over the decades, further experimental effort is still needed
to clarify a few more subtle aspects of water dynamics. These include the
microscopic mechanism leading shear acoustic modes of water to exhibit a
nearly flat $Q$ dependence at the crossover between quasi-macroscopic and
mesoscopic distances, while gradually becoming predominant over the
longitudinal acoustic modes.

From a merely experimental perspective, a natural pathway to understand better
the collective movements of water molecules is to observe how they compare
with those in the solid. Indeed, the study of the acoustic response of
(polycrystalline) ice presents significant simplifications: 1) both atomic
diffusion and HB network relaxations are frozen, thus yielding no visible
contribution to sound propagation; 2) the lifetime of acoustic modes is much
longer, and the related spectral features correspondingly sharper.
Nonetheless, as illustrated in the remainder of this paper, the phonon
response of ice is far from being trivially interpreted and characterized. In
this paper, we discuss an experimental attempt to elucidate similarities and
distinctive behaviors of the terahertz dynamics of ice and water by comparing
their IXS spectra. Specifically, we jointly investigate collective excitations
in water and phonon modes in polycrystalline hexagonal ice (Ih). On a rigorous
ground, inelastic excitations in the spectrum of a polycrystal should not be
referred to as phonons owing to their projection along multiple
crystallographic directions. However, in the following we will resort to this
broadly used nomenclature for consistency with existing literature.

The use of a Bayesian inference-based modeling of IXS lineshapes reveals an
unprecedentedly complex phonon behavior of ice, while rectifying and
complementing the current understanding of terahertz acoustic excitations in
liquid water.

##  II The measurement

Measurements were executed at the Sector 30 beamline [38, 39] of the Advanced
Photon Source at Argonne National Laboratory. The instrument was operated
using the $\approx$ 23.7 keV harmonic of the undulator source, which
corresponds to the Si(12 12 12) backscattering reflection from both the
monochromator and energy analyzers. Spectral acquisitions covered the 3 $\div$
21 nm-1 $Q$ range with a 2 nm-1 step. The instrumental resolution function was
measured through the IXS signal from a Plexiglas sample at the $Q$ of its
first sharp diffraction maximum, i.e., about 10 nm-1. The resulting spectrum
had a 0.8 meV broad (half width at half maximum, HWHM) nearly Lorentzian
profile, sufficient to properly resolve the relevant spectral features
discussed in the remainder of this paper. All measurements were executed at a
temperature of 225 K. The polycrystalline nature of the sample makes its
global orientation irrelevant. As a consequence, the measurement probed the
constant-$Q$ phonon response of ice as an average over the various
crystallographic directions. Fig. 1 illustrates the emergence of phonon
excitations in the IXS spectra of ice by comparing raw and best-fitting model
spectra for selected $Q$ values; the semilogarithmic plot enables a full
appreciation of relevant lineshape features across various intensity decades.
Model lineshapes, obtained as discussed in Sec. I of Supplemental Materials,
consist of a $\delta$-function describing the elastic scattering and
individual DHO terms accounting for inelastic excitations in the spectrum. The
most plausible number of DHO components was determined a posteriori, based on
the outcome of the Bayesian inferential analysis of measured data, as
described in our previous works [40, 41, 42, 43, 44]. Individual DHO
components of the best-fitting total spectral distribution are also included
in Fig. 1 for reference. The assignment of each spectral contribution is
summarized in the caption of Fig. 1 and illustrated in Table I.

Table 1: Polycrystalline Ice modes assignment Mode | Label | Assignment |  | Symbol color
---|---|---|---|---
1 | TA | Transverse Acoustic |  | blue
2 | OL | Lower Optical mode |  | grey
3 | OH | Higher Optical mode |  | orange
4 | LA | Longitudinal acoustic mode of center of mass |  | magenta
5 | HFM | highest frequency steeply dispersive mode |  | wine

Notice that the same color code is maintained consistent throughout all graphs
included in this work.

Figure 1: Inelastic X-Ray Scattering (IXS) spectra of polycrystalline ice and
their model lineshapes. IXS spectra measured at representative $Q$’s and
$T=225$ K discussed in this work (open circles) compared with corresponding
best-fit lineshapes (red lines through data) along with their inelastic DHO
model components (see text): the low frequency transverse mode (blue line),
the two intermediate frequency optic-like modes (grey and orange lines), the
longitudinal acoustic mode (magenta) and a higher frequency mode (wine line).
The green line represents the instrumental energy resolution. Error bars are
estimated as the square root of the scattering intensity counts.

The IXS spectra in Fig. 1 are reported in the region where they have been
modeled. This region, varying in extent and symmetry with respect to the
central peak position, essentially coincided with the energy range covered by
the corresponding energy resolution measurement.

One can readily notice that spectral shapes displayed in Fig. 1 are dominated
by multiple inelastic features, which, in some cases, partially overlap giving
rise to inelastic peaks visibly broader than the resolution function, as
particularly evident in the $Q$ = 21 nm-1 spectrum. Overall, Fig. 1 provides a
comprehensive rendering of the complexity of the phonon response of
polycrystalline ice down to mesoscopic distances. The $Q$-evolution of various
phonon modes will be discussed in further detail in the next Section, here we
only stress that these modes display quite distinctive $Q$ dependencies
revealing their either acoustical or optical origin.

##  III Discussion of results

The good consistency between the model and measured lineshapes clearly emerges
from the Fig. 1. Selected IXS spectra are compared with their best-fitting
model lineshapes discussed in the Sec. I of Supplemental Materials. Overall,
phonon spectra presented in this paper are qualitatively similar to those
reported in the IXS measurement described in Ref. [2]. However, the latter was
performed on a polycrystalline sample seemingly more ordered than the present
one, as suggested by the weaker elastic contribution, yet with a resolution
significantly coarser, as to be expected considering the still pioneering
stage of high-resolution IXS. Once assessed the ability of the model to
reproduce the detail of the spectral shape, we can now focus on best-fit
values of phonon frequencies, whose $Q$ dependence and assignment are included
in Fig. 2.

Figure 2: Phonon branches of polycrystalline ice and corresponding sound
dispersions in liquid water. Panel a: The plot displays the following
dispersion branches: the transverse acoustic (TA) mode (blue squares); the
longitudinal acoustic (LA) mode (magenta dots); the high-frequency mode (HFM,
wine dots); the lower and higher energy optical modes OL and OH (grey and
orange lozenges, respectively. Error bars in the dispersion curves are
estimated through the standard deviations of the corresponding posterior
distribution functions. The lines refer to literature results and are colored
in black or in red for ice and liquid water respectively. Specifically, they
represent the longitudinal mode of liquid water at 263 K and 2 kbar [34]
(solid red line); the transverse mode of room temperature D2O [31] (dashed red
line) and the LA mode of polycrystalline ice Ref. [45] (black dash-dotted
line). The horizontal red arrow points to the optical frequency measured by
optical spectroscopy [46]. Panel b: The phonon Density of States measured by
neutron scattering on single hexagonal ice ([47], blue line). Horizontal
arrows are guides to the eye showing the rough correspondence between DoS
features and dispersion curves extrema.

As apparent from the dispersion curves displayed in Fig. 2, the low-energy
($<$ 25 meV) phonons of ice probed by this measurement show some rather
peculiar features, previously either undetected [2, 48] or only partially
reported by studies in the literature [33, 45]. Specifically:

1) The longitudinal acoustic phonon of ice, i.e., the one characterized by the
sharpest (linear) low $Q$ growth in Fig. 2a, actually splits into two
branches, consistent with what was reported in a previous IXS work by
Cimatoribus and collaborators [33], and previous lattice calculations [49].
Most importantly, this trend might suggest reconsidering the surprising line
broadening of the longitudinal acoustic mode of both water and polycrystalline
ice reported in Ref. 50. The resolution limitation of such a pioneering IXS
measurement is the most likely explanation for the mode-splitting appearing as
a line-broadening. At low enough $Q$’s, the two split branches attain similar
values following a parallel $Q$-dependence while they spread apart at larger
$Q$’s. To estimate the low wavevector propagation speed of these two modes,
say ${c_{1}}$ and ${c_{2}}$, we looked at the slope of the three lowest $Q$
points of the corresponding phonon branches (Fig. 2). We found that the ratio
of the sound speeds ${c_{1}}/{c_{2}}$ = 1.09 $\pm$ 0.07, a value consistent
within the error to that of $\sqrt{M_{CM}/M_{O}}=1.06$ where $M_{CM}$ and
$M_{O}$ are the masses of the water molecule and the oxygen atom,
respectively. This finding might suggest that, in the linear dispersion
regime, the modes reported in Fig. 2 as LA and HFM, might be connected to
collective vibrations involving either the whole H2O molecule or the oxygen
atoms only. This is the main rationale behind the identification of the lower
frequency component of the doublet to the ordinary longitudinal acoustic (LA)
mode.

2) From moderate to high $Q$ values, the LA branch acquires a negative $Q$
slope, bending downward to some minimum at about 17 nm-1, in agreement with
what was previously observed by the molecular dynamic simulation by Criado and
collaborators [45]. Conversely, at about the same $Q$ value, the HFM reaches
its maximum at about 25 meV. Noticeably the plot suggests that the HFM branch
is the one most closely resembling the longitudinal acoustic dispersion of
liquid water, which is also reported in Fig. 2a, as measured by IXS at 263 K
and 2 kbar [34].

4) Noticeably, two additional weakly dispersing modes show up in the spectrum
inside the 7-12 meV, $Q>$ 10 nm-1 window: a lower and a higher optical mode
here labeled as OL and OH, respectively. We assign an optic origin to these
branches owing to their mild $Q$ dependence and seemingly non-vanishing low
$Q$ trend, which distinguishes them from the transverse acoustic mode, sitting
at an even lower frequency. The existence of multiple phonon branches at
similar energies was reported in a computational study by Criado and
collaborators on polycrystalline ice [45] and by previous works in the single
crystal at ambient [48], or higher pressures [51]. Also, present results are
consistent with Density of State measurements by incoherent INS [52, 53] and
optical spectroscopy [46]. In Fig. 2a, a horizontal arrow points to the $Q$ =
0 (energy axis) value corresponding to one of the optical mode energies
reported by the latter work. Revisited a posteriori in the light of current
results, the optical-like behavior of medium $Q$ phonons of ice is also
consistent with the findings of coherent scattering measurements by INS [48,
31] and IXS [2, 54].

Fig. 2a suggests that a proper discernment of two distinct modes, LA and HFM,
may be challenging at some $Q$ values. However, the inference of a mode
doublet with such a narrow energy separation is pondered against a precise
estimation of a joint probability distribution. The involved energies are so
close that the modes partially overlap, giving origin to a single peak yet
still significantly broader than the instrumental energy resolution and
somewhat flat atop as expected for an excitation composed by multiple spectral
modes of similar amplitude. According to the probability rating of the
Bayesian algorithm, the most plausible number of DHO modes participating to
such an excitation is two. It would be tempting to compare the inelastic modes
of the considered polycrystal to its single-crystal counterpart thoroughly
studied in the literature [48, 55]. However, this comparison would be
problematic since, as mentioned, excitations of our sample are averaged over
the various crystallographic directions as opposed to the single-crystal case.
For instance, the HFM dispersion curve might be contaminated by higher energy
optical modes especially in the apex region. In this perspective, the single-
crystal phonon Density of States (DoS), being a wavevector-averaged (non-
directional) property, provides a better quantity to compare with. For
reference, we included this spectral function in Fig. 2b as measured by
neutron scattering on single crystal Ih [47]. The horizontal arrows serve as
guides to the eye, highlighting the overall consistency between relevant DoS
features and the extremes of dispersion curves derived in the current
measurement.

Looking at the sound dispersions drawn in Fig. 2a, one may notice the
disappearance of the transverse acoustic (TA) branch for $Q>$ 14 nm-1. We tend
to rule out any physical rationale behind this trend, instead blaming the
tendency of several inelastic features to cram into a narrow energy window,
thus challenging the algorithm’s ability to distinguish them, especially those
having a faint spectral fingerprint. In fact at $Q=$21 nm-1 there might be a
hint for the presence of the transverse mode (see Fig. S1 and S2 of the SM)
but, in our opinion, not with enough experimental evidence. Very likely at
this $Q$ this mode is still overwhelmed by the near more intense optic modes
which makes the detection uncertain.

###  III.1 Comparing polycrystalline ice and liquid water

Fig. 3 compares, in a more restricted intensity and energy window, the phonon
lineshapes of ice with the low energy mode we measured in a previous work on
liquid water [56]. This comparison becomes especially informative when
considered in combination with results in Fig. 2.

Before digging further into this aspect, we recall here that the emergence of
the low energy mode in the spectrum of water has been the focus of intensive
experimental [31, 32, 33, 34, 35], and computational [35, 57] studies, which,
despite some controversy [58, 59], endorsed the now broadly accepted
assignment to a transverse acoustic excitation. However, results in Fig. 2 and
Fig. 3 call for a reconsideration of this assignment, which, ironically
enough, was originally inspired right by the comparison with the ice spectrum
[7, 2]. Indeed, Fig. 2 suggests that the transverse mode of water does
coincide with the TA branch of ice at the lowest $Q$’s only. Conversely, upon
$Q$-increase, it seems to mimic the lower energy optical phonon of ice, while
acquiring, in turn, a weaker $Q$ dependence. This conclusion is further
endorsed by Fig. 3, which shows that the low-energy spectral feature of water
closely resembles its transverse acoustic counterpart in ice for $Q\leq$ 9
nm-1 only. Conversely, at larger $Q$’s, it rather parallels the optical
phonons dominating the ice spectrum, even though the latter, owing to the lack
of damping, have visibly higher inelastic shift. These evidences urge us to
infer that collective modes in water correspondingly acquire a dominating
optic-like character in this $Q$ region. On a more general perspective, the
above arguments might seem to question the very nature of a collective mode in
a fluid, as they imply that such a mode can actually combine together
independent phonon vibrations. Moreover, the damping enhancement stemming from
diffusive and relaxation processes can make even fuzzier the physical
interpretation of the inelastic mode of a liquid. However, the scenario is not
as bewildering for water, at least when its spectrum is compared to that of
its solid-state, ordered, counterpart. In fact, although vibrational modes of
different nature seem to participate in the low energy mode of water, they
dominate complementary $Q$ windows: optical modes become preponderant upon $Q$
increase while (transverse) acoustic ones follow just the opposite trend. From
this perspective, one is authorized to conclude that, upon $Q$ increase, the
low frequency bump in the water spectrum changes its dominant character from
primarily acoustic to mainly optical.

Figure 3: Comparison of low and intermediate frequency phonon modes in
polycrystalline ice and in pure water. Selected spectra measured in this work
at the indicated $Q$s are here displayed in expanded energy transfer and
intensity windows. Experimental spectra in ice (black circles) are compared
with best-fitting model profiles (red line), the transverse low energy
component (blue line), and the two optical-like ones (grey and orange lines),
whenever present. The thick olive green line represents the low energy mode in
pure water at 300 K as obtained through the lineshape analysis discussed in
Ref. [56].

We believe that this finding is of high relevance since, thus far, no
experimental or computational evidence has ever challenged the consensual
assignment of such a spectral feature to an acoustic phonon-like (transverse)
mode. One of the most visible signs of this acoustic-to-optic transformation
is the gradually vanishing group velocity $v_{g}=\partial\Omega(Q)/\partial Q$
(with $\Omega(Q)$ being the phonon frequency), which reveals an underlying
mode localization. This conclusion rests on the assumption that similar
vibration modes are shared by water and ice, although only in ice they can be
properly resolved due to their sharp profile and relatively tiny energy
offset. We finally notice that, even in ice, in some $Q$ window, low energy
modes tend to overcrowd the quasielastic region of the spectrum, thus
challenging the ability of the Bayesian algorithm to properly identify their
presence. This “mode overcrowding” combined with the overwhelming intensity of
adjacent optical modes, is likely, as mentioned, the reason of an apparent
disappearance at a certain $Q$ value of the TA phonon (see Fig. 2a).

Figure 4: High-frequency dynamics in polycrystalline ice and in water. IXS
spectra (black circles) measured at two selected $Q$ values are reported along
with corresponding best-fitting curves (red lines) and both the HFM (wine
line) and the LA (magenta line) energy model (DHO) components. With the orange
and grey lines the optical-like modes are reported as in the previous figures.
The corresponding IXS spectra of pure water measured in Ref. [56] are also
shown (blue circles) after arbitrary vertical shift. The two vertical arrows
point at the position of the transverse mode energy in pure water.

Let us now comment on the comparison between current measurements on ice and
those on liquid water reported in Ref. [56] and also included in Fig. 4 after
arbitrary vertical shift. Both the HFM mode of ice and its counterpart in
water weaken upon $Q$-increase, down to almost vanish at high $Q$’s. In ice
the two optical branches, OL and OH, essentially take it over in the inelastic
wings. Notice that the naked eye’s perception of this mode in the 21 nm-1
spectrum of ice largely owes to the semilogarithmic representation. This
strong suppression of the HFM is also witnessed by the rapid increase of the
damping parameter with increasing $Q$ (Fig. S4 in the Supplemental Materials).
In water, such an effect is even more pronounced, due to the enhanced damping
which makes this mode barely distinguishable (if at all) from the large
inelastic wings of the dominating transverse mode, whose position is roughly
indicated by the two vertical arrows (see also Fig. 2 of Ref. [56]). Although
the naked eye can barely discern the HFM mode in the 21 nm-1 spectrum of Fig.
4, the Bayesian algorithm endorsed without ambiguity its presence. One can
quantify the statistical significance of this excitation through the high
probability rating of a lineshape model explicitly containing it (66 $\%$).
Nonetheless, the posterior of the related excitation frequency (see Fig. S2 of
the Supplemental Information) is quite broad and weak, thus revealing the
uncertain location of this otherwise genuine phonon excitation. Of course, the
significance of a poorly discernible spectral feature should be pondered
against all diagnostic factors available, mainly two in the present case: the
probabilistic support of the Bayesian inference and the continuity with the
trend followed by neighboring Q-points. Both indicators persuaded us about the
presence of a loose HFM feature at about 23 meV in the $Q$ = 21 nm-1 spectrum.

###  III.2 The key role of Bayesian analysis

At this stage, it is crucial to understand how the reported result fit into
the broader context of available literature data, and, in particular, why the
outcome of this study brings to the surface a more complex phonon behavior
than previously observed by individual works. Aside from some improvement in
the statistical accuracy in the currently measured spectra, the use of
Bayesian analysis here provides the critical asset of a minimally biased
modelling of measured spectra. In this approach, the model itself is entrusted
to establish on a probabilistic basis the most plausible number of phonon
modes in the measured spectrum. This provides a natural protection against two
somewhat opposite risks: the confirmation bias, holding off the discovery of
new spectral excitations, and the model over parametrization, as would be the
inclusion of too many independent excitations in the lineshape model. Bayesian
inference endows researchers with these antidotes also thanks to the intrinsic
implementation of the Occam razor principle, or ”lex parsimoniae”, prescribing
that, among equally plausible competing models, the simplest, i.e., the one
with a smaller number of free parameters, is always to be privileged [60].
Obviously, this aspect can become a game changer when facing inherent
ambiguities of a measured spectral line, perhaps exacerbated by a less than
adequate statistical accuracy. As an example we can consider the observation
of an anomalously wide spectral line, to be alternatively interpreted as a
genuine line broadening, secondary to an increased excitation lifetime, or the
emergence of a new line, indicating the onset of a mode-splitting. Of course,
Bayesian analysis only identifies the probabilistically most grounded
hypothesis based on the measurement obtained, yet it cannot substitute
investigators in difficult assessment of its physical plausibility. Further
details on the Bayesian inference method, model choice and results can be
found in the Supplemental Materials , specifically in Figs. S1, S2, S3 and S5.

##  IV Conclusion

We have presented an Inelastic X-Ray Scattering study of the terahertz dynamic
response of water, in which new measurements in the solid phase are compared
to our previous results in the liquid. In this work, a Bayesian inference
based analysis of measured spectra represented a pivotal tool to unveil the
full complexity of the phonon response of ice, thereby improving, via direct
comparison, our understanding of the water’s one.

In summary, our data show that, on approaching microscopic scales, high energy
modes becomes increasingly impeded. This hindrance is reflected, on one side,
by the rapid increase of the longitudinal acoustic damping of the highest
frequency mode, and, on the other, by a relative attenuation compared to optic
modes which gradually become dominant in adjacent spectral windows. Indeed, by
observing how measured spectra and dispersion curves derived from them compare
in liquid and solid phases, we are urged to conclude that optic terahertz
modes dominate the spectrum of density fluctuations of both water and ice at
sufficiently short distances. This finding is especially noteworthy as clearly
at variance with our long-lasting assumption of a dominating acoustic
character of collective modes of water. Despite this evidence, our results
also stimulate the scientific community with new interpretative challenges.
Likewise, the most urgent one is whether the observed behavior is distinctive
or shared with other networked fluids. A dedicated experimental and
computational effort is planned to elucidate this aspect hopefully improving
our understanding of the high-frequency response of disordered systems.

###### Acknowledgements.

This research was funded by Ministero dell’Istruzione dell’Università e della
Ricerca Italiano (Grant No. PRIN2017-2017Z55KCW).

## References

  * Christensen [2012] R. Christensen,  _Theory of viscoelasticity: an introduction_ (Elsevier, 2012).
  * Ruocco  _et al._ [1996a] G. Ruocco, F. Sette, U. Bergmann, M. Krisch, C. Masciovecchio, V. Mazzacurati, G. Signorelli, and R. Verbeni, Equivalence of the sound velocity in water and ice at mesoscopic wavelengths, Nature 379, 521 (1996a).
  * Monaco  _et al._ [1999] G. Monaco, A. Cunsolo, G. Ruocco, and F. Sette, Viscoelastic behavior of water in the terahertz-frequency range: an inelastic x-ray scattering study, Physical Review E 60, 5505 (1999).
  * Scopigno  _et al._ [2000] T. Scopigno, U. Balucani, G. Ruocco, and F. Sette, Evidence of two viscous relaxation processes in the collective dynamics of liquid lithium, Physical Review Letters 85, 4076 (2000).
  * Boon and Yip [1991] J. P. Boon and S. Yip,  _Molecular hydrodynamics_ (Courier Corporation, 1991).
  * Santucci  _et al._ [2006] S. Santucci, D. Fioretto, L. Comez, A. Gessini, and C. Masciovecchio, Is there any fast sound in water?, Physical Review Letters 97, 225701 (2006).
  * Sette  _et al._ [1996] F. Sette, G. Ruocco, M. Krisch, C. Masciovecchio, R. Verbeni, and U. Bergmann, Transition from normal to fast sound in liquid water, Physical Review Letters 77, 83 (1996).
  * Gorelli  _et al._ [2006] F. Gorelli, M. Santoro, T. Scopigno, M. Krisch, and G. Ruocco, Liquidlike behavior of supercritical fluids, Physical Review Letters 97, 245702 (2006).
  * Simeoni  _et al._ [2010] G. Simeoni, T. Bryk, F. Gorelli, M. Krisch, G. Ruocco, M. Santoro, and T. Scopigno, The Widom line as the crossover between liquid-like and gas-like behaviour in supercritical fluids, Nature Physics 6, 503 (2010).
  * Bryk  _et al._ [2017] T. Bryk, F. Gorelli, I. Mryglod, G. Ruocco, M. Santoro, and T. Scopigno, Behavior of supercritical fluids across the “frenkel line”, The Journal of Physical Chemistry Letters 8, 4995 (2017).
  * Bolmatov  _et al._ [2013] D. Bolmatov, V. Brazhkin, and K. Trachenko, Thermodynamic behaviour of supercritical matter, Nature Communications 4, 1 (2013).
  * Yang  _et al._ [2017] C. Yang, M. T. Dove, V. V. Brazhkin, and K. Trachenko, Emergence and evolution of the $k$ gap in spectra of liquid and supercritical states, Physical Review Letters 118, 215502 (2017).
  * Bolmatov  _et al._ [2017] D. Bolmatov, M. Zhernenkov, L. Sharpnack, D. M. Agra-Kooijman, S. Kumar, A. Suvorov, R. Pindak, Y. Q. Cai, and A. Cunsolo, Emergent optical phononic modes upon nanoscale mesogenic phase transitions, Nano Letters 17, 3870 (2017).
  * Cunsolo [2015] A. Cunsolo, The THz spectrum of density fluctuations of water: The viscoelastic regime, Advances in Condensed Matter Physics 2015 (2015).
  * Cunsolo  _et al._ [1999] A. Cunsolo, G. Ruocco, F. Sette, C. Masciovecchio, A. Mermet, G. Monaco, M. Sampoli, and R. Verbeni, Experimental determination of the structural relaxation in liquid water, Physical Review Letters 82, 775 (1999).
  * Sampoli  _et al._ [2009] M. Sampoli, U. Bafile, E. Guarini, and F. Barocchi, Collective dynamics and molecular interactions in liquid CO2 by inelastic neutron scattering and computer simulations, Physical Review B 79, 214203 (2009).
  * Guarini  _et al._ [2021] E. Guarini, F. Barocchi, A. De Francesco, F. Formisano, A. Laloni, U. Bafile, M. Celli, D. Colognesi, R. Magli, A. Cunsolo, and M. Neumann, Collective dynamics of liquid deuterium: Neutron scattering and approximate quantum simulation methods, Physical Review B 104, 174204 (2021).
  * de Schepper  _et al._ [1983] I. M. de Schepper, P. Verkerk, A. A. Van Well, and L. A. de Graaf, Short-wavelength sound modes in liquid argon, Physical Review Letters 50, 974 (1983).
  * Van Well  _et al._ [1985] A. A. Van Well, P. Verkerk, L. A. de Graaf, J.-B. Suck, and J. R. D. Copley, Density fluctuations in liquid argon: coherent dynamic structure factor along the 120-K isotherm obtained by neutron scattering, Physical Review A 31, 3391 (1985).
  * Cunsolo  _et al._ [2001] A. Cunsolo, G. Pratesi, R. Verbeni, D. Colognesi, C. Masciovecchio, G. Monaco, G. Ruocco, and F. Sette, Microscopic relaxation in supercritical and liquid neon, The Journal of Chemical Physics 114, 2259 (2001).
  * Balucani and Zoppi [1994] U. Balucani and M. Zoppi,  _Dynamics of the Liquid State_ (Clarendon, Oxford, 1994).
  * Balucani  _et al._ [1993] U. Balucani, A. Torcini, and R. Vallauri, Liquid alkali metals at the melting point: structural and dynamical properties, Physical Review B 47, 3011 (1993).
  * Guarini  _et al._ [2013] E. Guarini, U. Bafile, F. Barocchi, A. De Francesco, E. Farhi, A. Laloni, A. Orecchini, A. Polidori, M. Puglini, and F. Sacchetti, Dynamics of liquid Au from neutron Brillouin scattering and ab initio simulations: Analogies in the behavior of metallic and insulating liquids, Physical Review B 88, 104201 (2013).
  * Guarini  _et al._ [2020] E. Guarini, A. De Francesco, U. Bafile, A. Laloni, B. G. del Rio, D. J. Gonzalez, L. E. Gonzalez, F. Barocchi, and F. Formisano, Neutron Brillouin scattering and ab initio simulation study of the collective dynamics of liquid silver, Physical Review B 102, 054210 (2020).
  * De Francesco  _et al._ [2021] A. De Francesco, U. Bafile, A. Cunsolo, L. Scaccia, and E. Guarini, Searching for a second excitation in the inelastic neutron scattering spectrum of a liquid metal: a Bayesian analysis, Scientific Reports 11, 13974 (2021).
  * Maldovan [2013] M. Maldovan, Sound and heat revolutions in phononics, Nature 503, 209 (2013).
  * Cunsolo [2021] A. Cunsolo,  _The THz dynamics of liquids probed by inelastic X-ray scattering_ (World Scientific, 2021).
  * Sinha [2001] S. K. Sinha, Theory of inelastic x-ray scattering from condensed matter, Journal of Physics: Condensed Matter 13, 7511 (2001).
  * Squires [1996] G. L. Squires,  _Introduction to the theory of thermal neutron scattering_ (Courier Corporation, 1996).
  * Lovesey [1984] S. W. Lovesey,  _Theory of Neutron Scattering from Condensed Matter. Vol. 1: Nuclear Scattering_ (Oxford, 1984).
  * Cunsolo  _et al._ [2012] A. Cunsolo, C. Kodituwakku, F. Bencivenga, M. Frontzek, B. Leu, and A. Said, Transverse dynamics of water across the melting point: A parallel neutron and x-ray inelastic scattering study, Physical Review B 85, 174305 (2012).
  * Cunsolo  _et al._ [2016] A. Cunsolo, A. Suvorov, and Y. Q. Cai, The onset of shear modes in the high frequency spectrum of simple disordered systems: current knowledge and perspectives, Philosophical Magazine 96, 732 (2016).
  * Cimatoribus  _et al._ [2010] A. Cimatoribus, S. Saccani, F. Bencivenga, A. Gessini, M. Izzo, and C. Masciovecchio, The mixed longitudinal–transverse nature of collective modes in water, New Journal of Physics 12, 053008 (2010).
  * Pontecorvo  _et al._ [2005] E. Pontecorvo, M. Krisch, A. Cunsolo, G. Monaco, A. Mermet, R. Verbeni, F. Sette, and G. Ruocco, High-frequency longitudinal and transverse dynamics in water, Physical Review E 71, 011501 (2005).
  * Sampoli  _et al._ [1997] M. Sampoli, G. Ruocco, and F. Sette, Mixing of longitudinal and transverse dynamics in liquid water, Physical Review Letters 79, 1678 (1997).
  * Walrafen  _et al._ [1996] G. Walrafen, Y. Chu, and G. Piermarini, Low-frequency Raman scattering from water at high pressures and high temperatures, The Journal of Physical Chemistry 100, 10363 (1996).
  * Walrafen [1964] G. Walrafen, Raman spectral studies of water structure, The Journal of Chemical Physics 40, 3249 (1964).
  * Said  _et al._ [2011] A. H. Said, H. Sinn, and R. Divan, New developments in fabrication of high-energy-resolution analyzers for inelastic x-ray spectroscopy, Journal of Synchrotron Radiation 18, 492 (2011).
  * Toellner  _et al._ [2011] T. Toellner, A. Alatas, and A. Said, Six-reflection meV-monochromator for synchrotron radiation, Journal of Synchrotron Radiation 18, 605 (2011).
  * De Francesco  _et al._ [2016] A. De Francesco, E. Guarini, U. Bafile, F. Formisano, and L. Scaccia, Bayesian approach to the analysis of neutron Brillouin scattering data on liquid metals, Physical Review E 94, 023305 (2016).
  * De Francesco  _et al._ [2018] A. De Francesco, L. Scaccia, M. Maccarini, F. Formisano, Y. Zhang, O. Gang, D. Nykypanchuk, A. H. Said, B. M. Leu, A. Alatas, _et al._ , Damping off terahertz sound modes of a liquid upon immersion of nanoparticles, ACS Nano 12, 8867 (2018).
  * De Francesco  _et al._ [2019a] A. De Francesco, L. Scaccia, R. B. Lennox, E. Guarini, U. Bafile, P. Falus, and M. Maccarini, Model-free description of polymer-coated gold nanoparticle dynamics in aqueous solutions obtained by Bayesian analysis of neutron spin echo data, Physical Review E 99, 052504 (2019a).
  * De Francesco  _et al._ [2020a] A. De Francesco, L. Scaccia, F. Formisano, M. Maccarini, F. De Luca, A. Parmentier, A. Alatas, A. Suvorov, Y. Q. Cai, R. Li, _et al._ , Shaping the terahertz sound propagation in water under highly directional confinement, Physical Review B 101, 054306 (2020a).
  * De Francesco  _et al._ [2020b] A. De Francesco, L. Scaccia, F. Formisano, E. Guarini, U. Bafile, M. Maccarini, A. Alatas, Y. Q. Cai, D. Nykypanchuk, and A. Cunsolo, Onset of interfacial waves in the terahertz spectrum of a nanoparticle suspension, Physical Review E 102, 022601 (2020b).
  * Criado  _et al._ [1993] A. Criado, F. Bermejo, M. Garcia-Hernandez, and J. Martinez, Phonon dispersion in polycrystalline ice: Implications for the collective behavior of liquid water, Physical Review E 47, 3516 (1993).
  * Mazzacurati  _et al._ [1981] V. Mazzacurati, C. Pona, G. Signorelli, G. Briganti, M. Ricci, E. Mazzega, M. Nardone, A. De Santis, and M. Sampoli, Interaction induced light scattering: The translational spectrum of ice Ih single crystals, Molecular Physics 44, 1163 (1981).
  * del Rosso  _et al._ [2021] L. del Rosso, M. Celli, D. Colognesi, S. Rudic, N. J. English, and L. Ulivi, Density of phonon states in cubic ice Ic, The Journal of Physical Chemistry C 125, 23533 (2021).
  * Renker [1969] B. Renker, Phonon dispersion in D2O-ice, Physics Letters A 30, 493 (1969).
  * Johansson [1996] P. Johansson, Theory of inelastic x-ray scattering by phonons in ice, Physical Review B 54, 2988 (1996).
  * Ruocco  _et al._ [1996b] G. Ruocco, F. Sette, M. Krisch, U. Bergmann, C. Masciovecchio, and R. Verbeni, Line broadening in the collective dynamics of liquid and solid water, Physical Review B 54, 14892 (1996b).
  * Strässle  _et al._ [2004] T. Strässle, A. Saitta, S. Klotz, and M. Braden, Phonon dispersion of ice under pressure, Physical Review Letters 93, 225901 (2004).
  * Li  _et al._ [1991] J.-C. Li, J. Londono, D. Ross, J. Finney, J. Tomkinson, and W. Sherman, An inelastic incoherent neutron scattering study of ice II, IX, V, and VI—in the range from 2 to 140 meV, The Journal of Chemical Physics 94, 6770 (1991).
  * Klug  _et al._ [1991] D. D. Klug, E. Whalley, E. C. Svensson, J. H. Root, and V. F. Sears, Densities of vibrational states and heat capacities of crystalline and amorphous H2O ice determined by neutron scattering, Physical Review B 44, 841 (1991).
  * Schober  _et al._ [2000] H. Schober, M. M. Koza, A. Tölle, C. Masciovecchio, F. Sette, and F. Fujara, Crystal-like high frequency phonons in the amorphous phases of solid water, Physical Review Letters 85, 4100 (2000).
  * Wehinger  _et al._ [2014] B. Wehinger, D. Chernyshov, M. Krisch, S. Bulat, V. Ezhov, and A. Bosak, Diffuse scattering in ih ice, Journal of Physics: Condensed Matter 26, 265401 (2014).
  * De Francesco  _et al._ [2020c] A. De Francesco, L. Scaccia, F. Formisano, E. Guarini, U. Bafile, M. Maccarini, A. Alatas, Y. Q. Cai, and A. Cunsolo, The terahertz dynamics of an aqueous nanoparticle suspension: An inelastic x-ray scattering study, Nanomaterials 10, 860 (2020c).
  * Balucani  _et al._ [1996] U. Balucani, J. P. Brodholt, and R. Vallauri, Dynamical properties of liquid water, Journal of Physics: Condensed Matter 8, 9269 (1996).
  * Sacchetti  _et al._ [2004] F. Sacchetti, J.-B. Suck, C. Petrillo, and B. Dorner, Brillouin neutron scattering in heavy water: Evidence for two-mode collective dynamics, Physical Review E 69, 061203 (2004).
  * Petrillo  _et al._ [2000] C. Petrillo, F. Sacchetti, B. Dorner, and J.-B. Suck, High-resolution neutron scattering measurement of the dynamic structure factor of heavy water, Physical Review E 62, 3611 (2000).
  * MacKay [2003] D. MacKay,  _Information Theory, Inference and Learning algorithms_ (Cambridge University press, 2003).
  * Gilks  _et al._ [1996] W. R. Gilks, S. Richardson, and D. J. Spiegelhalter,  _Markov chain Monte Carlo in practice_ (Chapman & Hall/CRC, 1996).
  * Green [1995] P. J. Green, Reversible jump Markov chain Monte Carlo computation and Bayesian model determination, Biometrika 82, 711 (1995).
  * De Francesco  _et al._ [2020d] A. De Francesco, A. Cunsolo, and L. Scaccia, Bayesian approach for x-ray and neutron scattering spectroscopy, in  _Inelastic X-Ray Scattering and X-Ray Powder Diffraction Applications_ , edited by A. Cunsolo, M. K. K. D. Franco, and F. Yokaichiya (IntechOpen, 2020) Chap. 2, p. 26.
  * De Francesco  _et al._ [2019b] A. De Francesco, L. Scaccia, M. Maccarini, F. Formisano, E. Guarini, U. Bafile, and A. Cunsolo, Interpreting the terahertz spectrum of complex materials: the unique contribution of the Bayesian analysis, Materials 12, 2914 (2019b).

Supplemental Materials: Ice phonon spectra and Bayes inference: a gateway to a
new understanding of terahertz sound propagation in water

##  I Experimental details and Bayesian approach

To adequately reproduce the highly structured spectrum of density fluctuations
of ice, i.e., the dynamic structure factor $S(Q,E)$, we considered minimally
invasive hypothesis on its shape, which was assumed to include the
superposition of an unknown number of elementary excitations. Namely:

$\displaystyle S(Q,E)$ $\displaystyle=$ $\displaystyle
A_{e}(Q)\delta(E)+[n(E)+1]\frac{E}{k_{B}T}\Bigg{\\{}\sum_{j=1}^{k}A_{j}(Q){\rm
DHO}_{j}(Q,E)\Bigg{\\}},$ (S1)

here $\delta(E)$ is the Dirac delta function of $E$ describing the elastic
response and having integral $A_{e}$ and $n(E)=(e^{E/k_{B}T}-1)^{-1}$ is the
Bose statistics factor accounting for the detailed balance condition. The $k$
phonon excitations contributing to the spectrum are approximated by Damped
Harmonic Oscillator (${\rm DHO}_{j}(Q,E)$) multiplied by intensity factors
$A_{j}(Q)$ where:

$\displaystyle\noindent{\rm
DHO}_{j}(Q,E)=\frac{2}{\pi}\frac{\Omega_{j}^{2}(Q)\Gamma_{j}(Q)}{(E^{2}-\Omega_{j}^{2}(Q))^{2}+4[E\Gamma_{j}(Q)]^{2}},$
(S2)

and where $\Omega_{j}(Q)$ and $\Gamma_{j}(Q)$ are the undamped energy and the
damping coefficient of the _j_ th DHO excitation. We stress again that the
number $k$ of the ${\rm DHO}_{j}(Q,E)$ excitations and their shape
coefficients are equally treated as adjustable parameters of the model. To
describe accurately the measured spectrum, the model function in Eq. (S1) must
be convoluted with the instrument resolution function $R(Q,E)$ and a spectral
background should also be added to such a convolution. Explicitly:

$\tilde{S}(Q,E)=R(Q,E)\otimes S(Q,E)+B(E),$ (S3)

where $B(E)$ is a typically mildly $E$-dependent background intensity. A
Bayesian inferential procedure can be applied to determine, based on the
measurements and on our prior knowledge of the physical problem under
investigation, the most plausible shape, intensity and background parameters
entering in Eq. (S3) through Eqs. (S1) and (S2).

Best-fitting model lineshapes were determined using a Bayesian inferential
analysis implemented through a Markov Chain Monte Carlo (MCMC) [61] routine
with reversible jump (RJ) [62] steps. This approach can be used to
probabilistically infer, based on a given measurement, the joint posterior
probability distribution, or, in short, the posterior, of each parameter of
the model. Notice that in the present case a prior uniform distribution for
the number $k$ of inelastic modes has been deliberately chosen. In this
respect the Bayesian procedure can be used to probabilistically “rate” a guess
on the number of excitations in the sample scattering signal. In general, the
knowledge of the entire posterior distribution of each parameter also
authorizes the interpretation of the maximum of such a posterior as the
optimal (most plausible, or best-fitting) value of such a parameter. This
assignment is possible provided the posterior, albeit not necessarily
symmetric, is sharply peaked, well-shaped and unimodal. The shape itself of
the posterior carries important information on the precision and the
likelihood of a given best-fitting value. Foundations and working principle of
the performed Bayesian analysis are discussed in great detail in Ref. [40, 42,
63, 64].

##  II Comparing measured and model spectral shapes

The posterior distributions delivered by the MCMC-RJ algorithm (see main text)
for energy shifts and damping coefficients of the Damped Harmonic Oscillator
(DHO) model components are displayed in Figs. S1 and S2, respectively, as an
example at $Q$= 21 nm-1. The plots demonstrate how the Bayesian approach can
draw the entire probability distribution for all model parameters, rather than
providing best-fitting values only, as in a more traditional
$\chi^{2}$-minimization approach. From Fig. S1, it readily appears that the
majority of posteriors drawn for the DHO energy shifts are very sharp,
unimodal, and symmetric about their maxima, their mode thus yielding a
straightforward and unambiguous estimate of the best-fitting parameter value.
However, as illustrated in Fig. S2, the scenario is not always as rosy for the
damping coefficients, i.e., the width of the DHO model profiles; this
partially stems from the circumstance that is usually more challenging to
determine reliably the widths than the corresponding peak positions. In
general, the width of the posterior distribution provides a direct measure of
the uncertainty affecting the estimate of an optimal parameter value. This
uncertainty can be secondary to the poor statistical accuracy of data, or can
reveal a trend physically significant. For instance, a large posterior width
for a mode frequency can indicate a gradual approach to a mode overdamping
regime.

To provide a meaningful example, we can consider the $Q$ increasing damping of
the HFM, illustrated in Fig. S4. As discussed in the main manuscript in
further detail, such a high-energy phonon gradually damps off as $Q$
increases. This trend is also mirrored by the shape of the posterior
distribution of both frequency and damping of such mode, as clearly emerges
from considering the posteriors drawn for the highest energy phonon mode in
Figs. S1 and S2. These relatively unstructured posteriors hamper a precise
determination of the optimal parameter value. In the specific case, this
difficulty is exacerbated by the finite energy window of the available
experimental data, which in some cases prevented the measurement to fully
capture the position of the tails of the highest energy phonon mode.

Figure S1: Posterior distribution functions for the phonon mode frequencies in
polycrystalline ice.  The plot displays the posterior distribution functions
conditional on the experimental data $y$ for the phonon mode frequency (phonon
peak position) in the spectrum of polycrystalline ice at 225 K and wavevector
transfer $Q=$ 21 nm-1 as obtained by the MCMC-RJ algorithm (see main text).
Color code is as in the main text. Figure S2: Posterior distribution function
for the collective excitations dampings in pure polycrystalline ice. In the
four panels we report the posterior distribution drawn based on experimental
data $y$ for the damping coefficient of the five phonon modes identified in
the spectrum of pure ice at $Q=$ 21 nm-1. Figure S3: Posterior distribution
function for the excitation frequencies modes in ice at $T$ = 225 K. In the
left column are reported the posterior distribution functions for the
excitation frequencies modes revealed on ice at four selected $Q$ values. The
blue, grey, orange, magenta and wine colors indicate the distributions for the
transverse mode, the low frequency optical mode, the high frequency optical
mode, the longitudinal center of mass acoustic mode and the oxygen high
frequency mode, respectively. In the right column the corresponding traceplots
are reported with the same color code. Figure S4: Damping of the highest
energy mode in pure ice at $T$ = 225 K. Damping of the HFM associated to the
collective excitations of oxygen atoms as a function of the wavevector
transfer $Q$.

In Fig. S3 we give a graphical demonstration of how, depending on the value of
the wavevector $Q$, the number of modes that it is possible to recognize in
the IXS spectra changes as well as the greater or lesser “easiness” with which
they are solved. According to their dispersion curve and intensity, some modes
may come out from the resolution, be clearly visible, succumb to more intense
ones, or, finally, move out the investigated energy window. As in Fig. S1, in
the left column of the figure we report the fairly well-shaped posteriors of
the excitation frequencies of the different modes. In the right column, we
show the corresponding traceplots for the revealed mode frequencies. These
traceplots are one of the possible graphical (and diagnostic) methods MCMC
users adopt to decide with sound theoretical support when the algorithm
reaches convergence. The appearance of traceplots provides a clear
visualization of the efficiency of the Markov Chain in exploring the parameter
space. Bayesian statisticians usually call the correct look of a traceplot a
hairy caterpillar.

To provide an example of the uniqueness and the soundness of the model
description within the Bayes inference framework, in Fig. S5, we report the
posterior probability distribution of the number of excitations, conditional
to the measurements obtained at different $Q$’s. This distribution is
represented by the histograms drawn by the Bayesian algorithm. The
distribution mode, i.e., the position of the tallest histogram rectangle,
identifies the optimal parameter value. Adopting a model with such an optimal
value might appear a non-unequivocal choice if competing suboptimal model
options have comparable probability; however, such a choice is the one
probabilistically most grounded. Notice that suboptimal solutions having a too
small probability are not readily visible in the histograms, for this reason,
all probabilities are also reported in the numerical table here below.

Figure S5: Posterior distribution function for the number of modes in ice at selected momentum transfer values. Posterior distribution for the number of inelastic component $k$ in the spectra of polycrystalline ice at the indicated momentum transfer Q conditional to the experimental data $y$. The values of such probabilities are reported also in the Table here below, in fact too small probabilities are not detectable in the histograms. Table S1:  Posterior probabilities for the number of inelastic component in the spectra at the indicated different $Q$ values conditional to the experimental data $y$ | $k$
---|---
Q (nm-1) | 1 | 2 | 3 | 4 | 5
3 | 0 | 0 | 0.9995 | 0.0005 | 0
5 | 0 | 0 | 0.6883 | 0.3095 | 0.0022
7 | 0 | 0 | 0.6413 | 0.3401 | 0.0186
9 | 0 | 0.9981 | 0.0019 | 0 | 0
11 | 0 | 0 | 0.9225 | 0.0775 |
13 | 0 | 0 | 0 | 0.2741 | 0.7259
15 | 0 | 0 | 0.9832 | 0.0168 | 0
17 | 0 | 0 | 0.9382 | 0.0612 | 0.0006
19 | 0 | 0 | 0 | 1 | 0
21 | 0 | 0 | 0 | 0.3381 | 0.6619


# Determination of the proton spectral function of [12][]C from
$(e,e^{\prime}p)$ data

Artur M. Ankowski  Institute of Theoretical Physics, University of Wrocław,
50-204 Wrocław, Poland  Omar Benhar  INFN, Sezione di Roma , I-00185 Roma,
Italy  Makoto Sakuda  Physics Department, Okayama University, Okayama
700-8530, Japan

###### Abstract

The determination of the nuclear spectral function from the measured cross
section of the electron-nucleus scattering process $e+A\to e^{\prime}+p+(A-1)$
is discussed, and illustrated for the case of a carbon target. The theoretical
model based on the local density approximation, previously employed to derive
the spectral function from a combination of accurate theoretical calculations
and experimental data, has been developed further by including additional
information obtained from measurements performed with high missing energy
resolution. The implications for the analysis of $\gamma$-ray emission
associated with nuclear deexcitation are considered.

##  I Introduction

The electron-induced proton knockout reaction, in which the scattered electron
and the outgoing nucleon are detected in coincidence, has long been recognised
as a natural and powerful tool for the experimental determination of the
spectral function describing the energy and momentum distribution of protons
in the target nucleus; for a short historical account of the developments of
this field of research, see Ref. NPN .

Early studies of the $(e,e^{\prime}p)$ cross sections provided clear cut
evidence of the validity of the mean field approximation underlying the
nuclear shell model, exposing at the same time its conspicuous limitations
frullani_mougey . The spectroscopic lines corresponding to knockout of protons
occupying the shell model orbitals were, in fact, clearly observed in the
measured missing energy spectra, but the corresponding integrated strengths,
providing information on the normalisation of the proton wave functions,
turned out to be significantly lower than expected, regardless of the nuclear
target.

In the 1970s, the spectral functions of nuclei ranging from [12][]C to
[58][]Ni have been obtained from measurements carried out by Mougey et al.
using the 600 MeV electron beam available at CEN Saclay saclay . These
data—spanning the kinematic region corresponding to missing momentum
$p_{m}~{}<~{}320$ MeV and missing energy $E_{m}~{}<~{}80$ MeV—have been later
combined with the results of accurate nuclear matter calculations PkE:nm , to
model the full $p_{m}$ and $E_{m}$ dependence of the spectral functions within
the framework of the Local Density Approximation (LDA) LDA .

The spectral function $P({\bf k},E)$—trivially related to the two-point
Green’s function Green —is an intrinsic property of the target nucleus,
describing the probability of removing a particle of momentum ${\bf k}$ from
its ground state leaving the residual system with energy $E+E_{0}-m$, with
$E_{0}$ and $m$ being the ground-state energy and the mass of the knocked out
nucleon, respectively. As a consequence, it is a fundamental tool for the
analysis of a variety of nuclear scattering processes in the kinematic region
in which factorisation of the transition amplitude is expected to be
applicable. Notable applications include the electron- and neutrino-nucleus
scattering cross sections; see, e.g., Refs. LDA ; Benhar:2006wy ;
Benhar:2005dj ; Benhar:2006nr ; Ankowski:2014yfa ; gamma:16O . It should be
pointed out, however, that different reaction mechanisms and different
kinematic setups are sensitive to different features of the spectral function.
Theoretical analyses of inclusive $(e,e^{\prime})$ cross sections, involving
integration over the phase space of the knocked out nucleon, turn out to be
largely unaffected by the details of the energy dependence of $P({\bf k},E)$.
On the other hand, accurate calculations of processes involving transitions
between energy levels predicted by the nuclear shell model require a detailed
description of the removal energy distributions, allowing to pin down the
relevant spectroscopic strengths.

The authors of Ref. gamma:16O  applied the spectral function formalism in a
study of $\gamma$-ray emission from deexcitation of the residual nucleus
produced in [16][8]O$(\nu,\nu^{\prime}N)$ reactions, with $N$ denoting either
a proton or a neutron, which provides a signature of neutral-current neutrino-
nucleus interactions in water-Cherenkov detectors. Their analysis employed the
oxygen spectral function of Benhar et al. Benhar:2005dj , obtained from the
LDA-based approach of Ref. LDA  using the Saclay data of Bernheim et al.
Bernheim:1981si . This spectral function model allows the determination of the
spectroscopic strengths of the oxygen $p$-states, the values of which turn out
to be in remarkable agreement with those obtained from the analysis of the
high-resolution missing energy spectra of the [16][8]O$(e,e^{\prime},p)$
process performed by Leuschner et al. nikhef:oxygen .

In this article, we discuss the derivation of an improved model of the carbon
spectral function of Ref. LDA , designed to describe the missing energy
dependence of the [12][6]C$(e,e^{\prime},p)$ cross section measured at NIKHEF
using a 500 MeV high-duty-cycle electron beam and a pair of high-resolution
spectrometers nikhef . The improved spectral function model makes it possible
to pin down the spectroscopic strengths of the shell model states of the
carbon ground state, which are needed to extend the analysis of Ref. gamma:16O
to the case of neutral-current neutrino interactions in liquid scintillator
detectors.

The reminder of this article is organised as follows. The decomposition into
pole and continuum contributions and the phenomenological procedure employed
to obtain the nuclear spectral function within LDA are discussed in Sections
II.1 and II.2, respectively, whereas Section III describes the derivation of
an improved model of the carbon spectral function, including the information
provided by the high-resolution missing energy spectrum reported in Ref.
nikhef . Finally, in Section IV we summarise the results of our work and
outline their applications.

##  II The Nuclear Spectral function

When a nucleon is knocked out from a nucleus, the residual $(A-1)$-particle
system may be left in either a bound or a continuum state, in which at least
one of the spectators is excited to a positive energy level. These two
reactions mechanisms give rise to distinctively different missing energy
distributions.

###  II.1 Analytic properties of the missing energy distribution

In $(e,e^{\prime}p)$ processes, the missing energy is defined as

$\displaystyle E_{m}=\omega-T_{p}-T_{A-1}\approx\omega-T_{p}\ ,$ (1)

where $\omega=E_{e}-E_{e^{\prime}}$ denotes the electron energy loss, while
$T_{p}$ and $T_{A-1}$, with typically $T_{A-1}\ll T_{p}$, are the kinetic
energies of the outgoing proton and the residual nucleus, respectively.
Equation (1) shows that, to the extent to which Final State Interactions (FSI)
between the struck nucleon and the recoiling system can be neglected or
reliably described by an optical potential, the measured value of $E_{m}$ can
be identified with the removal energy of the knocked out particle. It follows
that transitions to bound ($A-1$)-nucleon states—associated with one-
particle–one-hole excitations of the target nucleus—lead to the appearance of
a missing energy distribution consisting of narrow peaks corresponding to the
poles of the spectral function at $E_{m}\approx|E_{\alpha}|$, with
$E_{\alpha}<0$ being the energy of the shell model state $\alpha$. On the
other hand, a smooth distribution, extending into the region of missing
energies above the two-nucleon emission threshold, signals the occurrence of
$n$-particle–$n$-hole excitations, with $n\geq 2$, arising from strong
correlations among the nucleons in the target ground state.

Based on the above considerations, the proton spectral function can be written
in model-independent form as

$\displaystyle P({\bf k},E)=P_{\rm MF}({\bf k},E)+P_{\rm corr}({\bf k},E)\ ,$
(2)

where $P_{\rm MF}$ and $P_{\rm corr}$ denote the pole, or mean-field (MF), and
continuum, or correlation (corr), contributions, respectively. Integration
over the full energy and momentum ranges yields the normalisation

$\displaystyle\int d^{3}kdE~{}P({\bf k},E)=Z\ ,$ (3)

with $Z$ being the nuclear charge number.

According to the shell model, in which correlations are disregarded
altogether, $P_{\rm corr}({\bf k},E)=0$, and the spectral function reduces to

$\displaystyle P_{\rm SM}({\bf
k},E)=\sum_{\alpha\in\\{F\\}}|\phi_{\alpha}({\bf
k})|^{2}\delta(E-|E_{\alpha}|)\ .$ (4)

In the above equation, $\alpha\equiv\\{n_{\alpha},\ell_{\alpha},j_{\alpha}\\}$
denotes the set of quantum numbers specifying the proton state $\alpha$,
$\phi_{\alpha}({\bf k})$ is the corresponding unit-normalised wave function in
momentum space, and the sum is extended to all states belonging to the Fermi
sea $\\{F\\}$.

The right-hand side of Eq. (4) is substantially modified by the effects of
nuclear dynamics beyond the shell model. Ground-state correlations lead to the
appearance of a non-vanishing continuum contribution, which turns out to
account for $\sim$20% of the normalisation PkE:nm . Furthermore, the pole
contribution—the normalisation of which is obviously reduced by the same
amount with respect to the shell-model prediction—is modified by the coupling
between one-hole and two-hole–one-particle states of the recoiling nucleus,
giving rise to a broadening of the energy distribution of the shell model
states Z . The resulting expression can be conveniently written in the form

$\displaystyle P_{\rm MF}({\bf k},E)$
$\displaystyle=\sum_{\alpha\in\\{F\\}}|\phi_{\alpha}({\bf
k})|^{2}f_{\alpha}(E)\ ,$ (5)

with the normalisation of the finite-width distribution $f_{\alpha}(E)$ being
given by

$\displaystyle\int dE\ f_{\alpha}(E)=Z_{\alpha}\ .$ (6)

The decomposition of the spectral function into discrete and continuum
contributions, which follows immediately from the Källen-Léhman representation
of the two-point Green’s function FW , is a powerful method for the model-
independent identification of correlation effects, the size of which is gauged
by the difference $[1-(\sum_{\alpha}Z_{\alpha})/Z]\geq 0$. A detailed
discussion of this topic can be found in Ref. Z .

It should be kept in mind that, although Eqs. (2) and (5) provide a useful and
physically motivated scheme for the parametrisation of the full spectral
function, the width of the distributions $f_{\alpha}(E)$, which is known to be
vanishingly small for valence states in the vicinity of the Fermi surface,
becomes very large for $E_{\alpha}\ll E_{F}$, with $E_{F}$ being the proton
Fermi energy. As a consequence, deeply bound states turn out to be very short-
lived, and their interpretation in terms of shell model orbits must be
carefully analysed.

Theoretical calculations of the nuclear spectral function based on advanced
microscopic models of nuclear dynamics and accurate computational techniques
have been carried out for the three-nucleon system Dieperink:1976wy ;
CiofidegliAtti:1980dbw ; Meier-Hajduk:1983pii , the doubly closed shell
nucleus [16][8]O Geurts:1996zza ; Rijsdijk:1996zz ; Polls:1996hg ; noemi:1 ,
and isospin-symmetric nuclear matter PkE:nm ; Ramos:1989hqs . Calculations
involving a somewhat simplified treatment of the continuum component have been
performed for [4][]He CiofidegliAtti:1990dh ; Morita:1991ka ; Benhar:1993ja .
More recently, the spectral functions of nuclei as heavy as [40][20]Ca or
having significant neutron excess, such as [40][18]Ar, have been also studied
using the Green’s function Monte Carlo method noemi:2 .

###  II.2 Phenmomenological determination of the spectral function

In the absence of FSI, the cross section of the process $e+A\to
e^{\prime}+p+(A-1)$, can be written in the simple and transparent form

$\displaystyle\frac{d\sigma_{A}}{dE_{e^{\prime}}d\Omega_{e^{\prime}}dE_{p^{\prime}}d\Omega_{p^{\prime}}}=E_{{p^{\prime}}}|{\bf
p}^{\prime}|\ \sigma_{ep}\ P(-{\bf p}_{m},E_{m})\ .$ (7)

Here $\Omega_{e^{\prime}}$ and $\Omega_{p^{\prime}}$ are the solid angles
specifying the directions of the scattered electron and the outgoing proton,
whose energy and momentum are denoted $E_{{p^{\prime}}}$ and ${\bf
p}^{\prime}$, respectively. The cross section $\sigma_{ep}$ describes elastic
electron scattering off a bound proton of momentum $-{\bf p}_{m}$ and removal
energy $E_{m}$. The missing momentum is defined as ${\bf p}_{m}={\bf q}-{\bf
p}^{\prime}$, with ${\bf q}$ being the momentum transfer.

The large data set of precisely measured $(e,e^{\prime}p)$ cross sections
collected at Saclay and NIKHEF from the 1970s to the 1980s covers mainly the
kinematic region where mean-field dynamics dominates. As a consequence, Eq.
(7) implies that, as long as the effects of FSI can be reliably taken into
account, these data have the potential to allow an accurate determination of
the spectroscopic strengths $S_{\alpha}$, obtained from integration of the
measured missing energy spectra over the regions
$|E_{\alpha}|-\Gamma_{\alpha}\lesssim
E_{m}\lesssim|E_{\alpha}|+\Gamma_{\alpha}$, with $\Gamma_{\alpha}>0$ being the
width of the peak of $f_{\alpha}(E)$ at $E\approx|E_{\alpha}|$. Note that, by
definition, $S_{\alpha}\leq Z_{\alpha}$, because the integration of Eq. (6) is
extended to all values of $E$. Therefore, $Z_{\alpha}$ includes contributions
associated with the occurrence of two-hole–one-particle intermediate states of
the residual nucleus, giving rise to the appearance of appreciable tails in
the distribution $f_{\alpha}(E)$ Z .

In experimental analyses, the corrections originating from FSI are usually
treated within the well established framework known as Distorted Wave Impulse
Approximation (DWIA), based on the use of complex proton-nucleus optical
potentials; for a detailed discussion of the formalism of DWIA, see, e.g.,
Refs. DWIA:1 ; DWIA:2 .

Experimental studies aimed at pinning down the correlation contribution
$P_{\rm corr}({\bf k},E)$ require a careful choice of the kimematic setup,
needed to minimise the background due to the occurrence of inelastic
processes, as well as a more refined treatment of FSI daniela_FSI . The data
collected at Jefferson Lab using a carbon target, while not providing a
detailed description of the full energy and momentum dependence, show a
remarkable consistency between the measured correlation strength and the
missing strength in the mean-field sector observed by previous experiments
daniela .

The formalism based on LDA, developed by the authors of Ref. LDA  in the
1990s, allows to determine the nuclear spectral function by exploiting the
decomposition of Eq. (2) and the availability of accurate measurements of the
mean-field contribution, Eq. (5). The main tenet underlying this
approach—strongly supported by the results of both theoretical and
experimental studies Alvioli:2016wwp —is that short-range nuclear dynamics,
driving the appearance of correlations among the nucleons, are largely
unaffected by surface and shell effects. As a consequence, correlation effects
can be reliably evaluated in uniform nuclear matter.

Within LDA, the continuum component of the spectral function of a nucleus of
mass number $A$ is approximated by

$\displaystyle P_{\rm corr}^{\rm LDA}({\bf k},E)=\int d^{3}r\ \varrho_{\rm
A}({\bf r})~{}P_{\rm corr}^{\rm NM}[{\bf k},E;\varrho=\varrho_{A}({\bf r})]\
,$ (8)

where $\varrho_{A}({\bf r})$ is the unit-normalised nuclear density
distribution, and $P_{\rm corr}^{\rm NM}({\bf k},E;\varrho)$ denotes the
continuum component of the spectral function of nuclear matter at density
$\varrho$. A detailed derivation of Eq. (8) can be found in Ref. LDA .

Accurate calculations of $P_{\rm corr}^{\rm NM}$ in the density range
$0.25\leq(\varrho/\varrho_{0})\leq 1$, with $\varrho_{0}\approx 0.16\ {\rm
fm}^{-3}$ being the equilibrium density of isospin-symmetric matter, have been
carried out using an advanced many-body approach and a realistic model of the
nuclear Hamiltonian PkE:nm ; LDA .

##  III The proton spectral function of carbon

###  III.1 LDA model based on Saclay data

Mougey et al. obtained the mean-field contribution to the proton spectral
function of carbon from the $(e,e^{\prime}p)$ data reported in Ref. saclay .
Based on the combined analysis of the missing momentum and missing energy
dependence of the measured cross section they were able to single out the
contribution of proton knockout from different shell-model states. The peak
observed at $15\leq E_{m}\leq 21.5$ MeV, associated with a maximum of the
corresponding momentum distribution at $p_{m}\approx 100$ MeV, was attributed
to the valence $p$-state, whereas the strength at $30\leq E_{m}\leq 50$ MeV,
whose $p_{m}$-dependence features a maximum at $p_{m}=0$, was interpreted as
an $s$-state, that is, $\ell=0$, contribution. This decomposition allowed to
extract the spectroscopic strengths of the $p$\- and $s$-states from the data,
corrected to take into account the effects of FSI.

The authors of Ref. LDA  have combined the empirical information reported in
Ref. saclay  with the nuclear matter results of Ref. PkE:nm  to obtain the
proton spectral function of carbon within LDA. The corresponding energy
distribution, defined as

$\displaystyle f_{\rm LDA}(E)=\int d^{3}k~{}\big{[}P_{\rm MF}({\bf
k},E)+P_{\rm corr}^{\rm LDA}({\bf k},E)\big{]}\ ,$ (9)

with $P_{\rm MF}$ and $P_{\rm corr}^{\rm LDA}$ given by Eqs. (5) and (8),
respectively, is displayed in Fig. 1. Note that $f(E)$ is nonzero in the
region $E\geq E_{\rm thr}$, with the proton knockout threshold $E_{\rm thr}$
being defined as

$\displaystyle E_{\rm
thr}=|M(\isotope[12][]{C})-M(\isotope[11][]{B})-m_{p}-m_{e}|\ ,$ (10)

where $m_{p}=938.27$ MeV and $m_{e}=0.51$ MeV denote the proton and electron
mass, respectively PDG . Using $M(\isotope[12][]{C})~{}=~{}11177.93$ MeV and
$M(\isotope[11][]{B})~{}=~{}10255.11$ MeV NDT , Eq. (10) yields $E_{\rm
thr}=15.96$ MeV.

Integration of the distribution of Fig. 1 yields the spectroscopic factors
$S_{p}=2.8$ and $S_{s}=1.2$, corresponding to occupation probabilities
$S_{\alpha}/(2j_{\alpha}+1)$ = 0.7 and 0.6, respectively. Note that the shell-
model predicts $S_{p}=4$ and $S_{s}=2$.

Figure 1:  Energy distribution obtained from momentum integration of the LDA
spectral function of Ref. LDA .

###  III.2 Parametrisation of NIKHEF data

In Ref. nikhef , Van der Steenhoven et al. reported the missing energy
spectrum of the [12][]C$(e,e^{\prime}p)$ reaction in the range corresponding
to knock out of a proton occupying the valence $p$-states of carbon, measured
at NIKHEF with 150 keV resolution (FWHM). The resulting distribution,
$f(E_{X})$, displayed in Fig. 2, exhibits three prominent peaks at excitation
energy $E_{X}=E_{m}-E_{\rm thr}=$ 0—that is, the process in which the
recoiling nucleus of [11][]B is left in its ground state—a $p_{1/2}$ state at
$E_{X}=2.125$ MeV, and a $p_{3/2}$ state at $E_{X}=5.020$ MeV. The primary aim
of this work is improving the spectral function model or Ref. LDA , based on
Saclay data, by inclusion of the additional information provided by the data
of Ref. nikhef .

As a first step, the distribution of Fig. 2 has been modeled assuming
gaussian-shaped peaks whose position and width, as well as the corresponding
spectroscopic strengths, have been determined using a code based on the
Marquardt fitting algorithm ingo_code .

Figure 2:  Comparison between the energy distributions extacted from NIKHEF
data nikhef  and the results of the gaussian fit employed in this work,
represented by the solid lines.

In order to combine NIKHEF data with the LDA spectral function based on Saclay
data, the fit shown by the solid line of Fig. 2, hereafter $f_{\rm fit}$, has
been renormalised in such a way as to match the $p$-state normalisation
resulting from integration in the LDA spectral function of Ref. LDA  in the
range $15\leq E_{\rm m}\leq 21.5$ MeV, which will be denoted $\Delta E$. After
renormalisation, the fitted energy distribution fulfils the condition

$\displaystyle\int_{\Delta E}dE_{\rm m}f_{\rm fit}(E_{\rm m})=\int_{\Delta
E}dE_{\rm m}f_{\rm LDA}(E_{\rm m})=S_{p}\ ,$ (11)

with $S_{p}=2.8$, by construction.

In Table 1, the spectroscopic strengths defined as

$\displaystyle S_{\alpha}=\int_{\Delta{E_{\alpha}}}dE_{\rm m}f_{\rm
fit}(E_{\rm m})\ ,$ (12)

as well as the values of the occupation probabilities of the shell-model
orbits, $Z_{\alpha}=S_{\alpha}/(2j_{\alpha}+1)$, are compared to the
experimental results reported in Ref. nikhef . Here the index $\alpha$ runs
over the states at $E_{m}=$ 15.96, 18.08, and 20.97 MeV, and
$\Delta{E_{\alpha}}$ is a narrow integration region designed to pin down the
associated pole strength.

$\alpha$ |  $p_{3/2}$ |  $p_{1/2}$ |  $p_{3/2}$
---|---|---|---
$S_{\alpha}^{\rm expt}$ | 1.72 | 0.26 | 0.20
$Z_{\alpha}^{\rm expt}$ | 0.43 | 0.13 | 0.05
$S_{\alpha}$ | 1.747 | 0.273 | 0.201
$Z_{\alpha}$ | 0.437 | 0.137 | 0.05
$E_{m}$ [MeV] | 15.96 | 18.08 | 20.97
Table 1:  Spectroscopic strengths and occupation probabilities of the valence
$p$-states of [12][]C at missing energy $E_{m}$, obtained from Eq. (12). The
experimental results reported in Ref. nikhef  are also listed, for comparison.

The total strength, $\sum_{\alpha}S_{\alpha}=2.22$, should be compared to
value, $S_{p}=2.8$ obtained from the LDA spectral function of Ref. LDA . This
$\sim$25% discrepancy is to be ascribed to additional contributions arising
from fragmentation of the shell-model states, clearly visible in the missing
energy spectra reported in Ref. nikhef  and not included in the determination
of the $S_{\alpha}$ from Eq. (12).

###  III.3 Improved spectral function model

The energy distribution assembled combining the fit to NIKHEF data, spanning
the region $E~{}\leq~{}E_{\rm match}=21.5$ MeV, and the LDA distribution of
Ref. LDA  at $E>E_{\rm match}$ is displayed by the thick solid line of Fig. 3.
Comparison with the histogram representing the full distribution of Ref. LDA ,
corresponding to Eq. (9), highlights the importance of taking into acount the
results of measurements performed with high missing energy resolution.

The integrated $s$-shell strength at $E\geq 21.5$ MeV discussed by Mougey et
al. saclay  has been also modeled, assuming a Maxwell-Boltzmann shape peaked
at $E=37$ MeV. This component, shown by the dot-dash line of Fig. 3,
contributes a spectroscopic strength $S_{s}=1.15$, which turns out to be
within $\sim 5$% of the corresponding result obtained from the data of Ref.
saclay .

Figure 3: The full line represents the energy distribution obtained combining
the NIKHEF data at $E\leq E_{\rm match}=21.5$ MeV and the results of Ref. LDA
at at $E>E_{\rm match}$, while the dot-dash line corresponds to the $s$-shell
contribution described in the text. For comparison, the full distribution of
Ref. LDA  is shown by the histogram.

A procedure analogue to that used for deriving the energy distribution has
been employed to obtain an improved model of the full spectral function.
Within this scheme, $P({\bf k},E>E_{\rm match})$ is the same as the LDA
spectral function of Ref. LDA , while in the region $E~{}\leq~{}E_{\rm match}$
we assume the factorised form

$\displaystyle P({\bf k},E<E_{\rm match})=n_{p}({\bf k})f(E)\ ,$ (13)

where $n_{p}({\bf k})=|\phi_{p}({\bf k})|^{2}$—with $\phi_{p}({\bf k})$ being
the $p$-state wave function in momentum space—and $f(E)$ is the distribution
represented by the solid line of Fig. 3. The energy and momentum dependence of
the resulting spectral function, normalised according to Eq. (3), is
illustrated in Fig. 4.

Figure 4: Energy and momentum dependence of the carbon spectral function
derived from the approach described in this work. Panels (a) and (b) have been
obtained with and without including the geometric factor $4\pi{\bf k}^{2}$,
respectively.

Integration of the spectral function over the energy range $\Delta E$ yields
the partial momentum distribution

$\displaystyle n_{\Delta E}({|\bf k|})=\int_{\Delta E}dE~{}P({\bf k},E)\ .$
(14)

The above definition implies that suitable choices of $\Delta E$ allow to pin
down the contributions of protons occupying different shell-model states, or
belonging to a correlated pair in the continuum. Figure 5 shows that the
partial distributions corresponding to $30<E<50$ MeV and $15<E<21.5$ MeV peak
at $|{\bf k}|=0$ and $\sim$100 MeV, as expected for $s$\- and $p$-shell
protons, respectively, and become vanishingly small at momenta larger than the
typical nuclear Fermi momentum $k_{F}\approx 250$ MeV. On the other hand,
integration over the full energy range, extending up to $E=300$ MeV, leads to
the appearance of a high momentum tail originating from strong short-range
correlations among the nucleons.

Figure 5: Partial proton momentum distributions of carbon, obtained from Eq.
(14) using the spectral function derived in this work and different
integration regions $\Delta E$. The spectroscopic factors resulting from
integration of the dashed and dot-dash lines, corresponding to $p$\- and
$s$-shell strength, respsctively, turn out to be $S_{p}=2.8$ and $S_{s}=1.2$.

In Fig. 6, the $p$-state momentum probability obtained multiplying by a factor
${\bf k}^{2}$ the corresponding distribution—shown by the dashed line of Fig.
5—is compared to the experimental result reported in Ref. saclay . Note that,
in order to perform a meaningful comparison, the distribution obtained from
our model has been modified to account for the effects of FSI, applying the
same correction employed in the analysis of Mougey et al. saclay . The
agreement between experimental data and the results of our approach turns out
to be remarkable.

Figure 6:  Momentum probability of the $p$-state of [12][]C, obtained from
integration of the spectral function described in this work in the range
$15<E<21.5$ MeV. The effects of FSI have been taken into account using the
same correction applied to the data by the authors of Ref. saclay .

##  IV Summary and Outlook

We have updated the carbon spectral function of Ref. LDA —obtained from a
semi-phenomenological model based on LDA—by including additional information
provided by the high-resolution [12][]C$(e,e^{\prime}p)$ experiment performed
at NIKHEF, Amsterdam, by Van der Steenhoven et al. nikhef .

The new spectral function features a modified energy dependence in the range
at $E\leq 21.5$ MeV, relevant to removal of a $p$-shell proton, and reduces to
the $P({\bf k},E)$ of Ref. LDA  at larger The primary improvement lies in the
capability to clearly identify the three previously unresolved states with
removal energies 15.96, 18.08, and 20.97 MeV, which turn out to account for
$\sim$55% of the $p$-state strength predicted by the nuclear shell model.

The strength in the region $E>21.5$ MeV is described in terms of an $s$-shell
contribution—extending at energies $21.5\lesssim E\lesssim 60$ MeV and modeled
by a Maxwell-Boltzmann distribution peaked at $E\sim 37$ MeV—and a smooth
background contribution, arising from fragmentation of shell-model states and
nucleon-nucleon correlations, extending to energies as high as $\sim 300$ MeV.

It should be pointed out that the modified energy dependence of the new
spectral function does not affect significantly the inclusive electron
scattering cross section, the calculation of which involves integrations of
$P({\bf k},E)$ over broad energy ranges.

Figure 7: Differential [12][]C$(e,e^{\prime})$ cross section at beam energy
$E_{e}=200$ MeV and electron scattering angle $\theta_{e}=60^{\circ}$, as a
function of electron energy loss $\omega$. Experimental data from Ref.
Barreau:1983ht .

To illustrate this point at quantitative level, we have compared the inclusive
electron-carbon cross section corresponding to beam energy $E_{e}=200$ MeV and
electron scattering angle $\theta_{e}=60^{\circ}$ obtained from the spectral
function derived in this work to the one computed by the authors of Ref.
Ankowski:2014yfa  using the spectral function of Ref. LDA . It appears that
the results of the two calculations are nearly indistinguishable from one
another over the whole range of electron energy loss $\omega$.

Note that the failure to reproduce the data at low $\omega$ signals the
breakdown of the factorisation approximation underlying the spectral function
formalism. In particular, the excitation of the giant dipole resonance at
$E_{X}=22.6$ MeV ref:nuclearDataSheets_12 —corresponding to $\omega=24.2$ MeV
due to nuclear recoil—is not accounted for in the calculations presented in
Fig. 7. In addition, the experimental data include contributions originating
from elastic scattering and transitions to discrete levels, leading to the
appearance of the jump of the cross section at the lowest value of $\omega$.

The potential of the improved $P({\bf k},E)$ can be fully exploited in studies
requiring a detailed description of the $E$-dependence, such as the analyses
of nuclear deexcitation processes. Notable example are the calculations of the
rate of $\gamma$-rays emitted by the residual nucleus in the aftermath of
neutral current neutrino interactions in the quasielastic channel gamma:16O ,
or the identification of the leading proton kinematics in Monte Carlo cascade
model simulations of FSI in neutrino-nucleus interactions Ershova:2023dbv .
The results reported in this article will allow to extend the study of Ref.
gamma:16O , relevant to measurements performed using water-Cherenkov
detectors, such as K2K Kameda:2006zn , to the analysis of experiments using
liquid scintillator detectors, such as KamLAND KamLAND:2022ptk  and JUNO
JUNO:2015zny . The need of accurate models of the carbon and oxygen spectral
functions for the description of deexcitation processes in Monte Carlo
simulations has been recently discussed in Ref. PhysRevD.109.036009 .

###### Acknowledgements.

A.M.A. and O.B. gratefully acknowledge partial support by grant NSF
PHY-1748958 to the Kavli Institute for Theoretical Physics (KITP), where this
research was initiated. The work of A.M.A. is partly supported by the National
Science Centre under grant UMO-2021/41/B/ST2/02778.

## References

  * (1) O. Benhar, “Exploring Nuclear Dynamics with $(e,e^{\prime}p)$ Reactions: From LNF to JLab”, Nuclear Physics News 26, 15-20 (2016). 10.1080/10619127.2016.1177334.
  * (2) S. Frullani and J. Mougey, “Single Particle Properties of Nuclei through $(e,e^{\prime}p)$ reactions”, Adv. Nucl. Phys. 14, 1-283 (1984).
  * (3) J. Mougey, M. Bernheim, A. Bussiere, A. Gillibert, X. H. Phan, M. Priou, D. Royer, I. Sick, and G. J. Wagner, “Quasifree ($e,e^{\prime}p$) Scattering on [12][]C, [28][]Si, [40][]Ca and [58][]Ni”, Nucl. Phys. A 262, 461-492 (1976). 10.1016/0375-9474(76)90510-8.
  * (4) O. Benhar, A. Fabrocini and S. Fantoni, “The Nucleon Spectral Function in Nuclear Matter”, Nucl. Phys. A 505, 267-299 (1989). 10.1016/0375-9474(89)90374-6.
  * (5) O. Benhar, A. Fabrocini, S. Fantoni and I. Sick, “Spectral function of finite nuclei and scattering of GeV electrons”, Nucl. Phys. A 579, 493-517 (1994). 10.1016/0375-9474(94)90920-2.
  * (6) O. Benhar, A. Fabrocini and S. Fantoni, “Nuclear-matter Green functions in correlated-basis theory”, Nucl. Phys. A 550, 201-222 (1992). 10.1016/0375-9474(92)90679-E.
  * (7) O. Benhar, N. Farina, H. Nakamura, M. Sakuda and R. Seki, “Electron- and neutrino-nucleus scattering in the impulse approximation regime”, Phys. Rev. D 72, 053005 (2005). 10.1103/PhysRevD.72.053005.
  * (8) O. Benhar and D. Meloni, “Total neutrino and antineutrino nuclear cross sections around 1 GeV”, Nucl. Phys. A 789 (2007), 379-402 10.1016/j.nuclphysa.2007.02.015
  * (9) O. Benhar, D. Day and I. Sick, “Inclusive quasi-elastic electron-nucleus scattering”, Rev. Mod. Phys. 80, 189-224 (2008). 10.1103/RevModPhys.80.189.
  * (10) A. M. Ankowski, O. Benhar, T. Mori, R. Yamaguchi and M. Sakuda, “Analysis of $\gamma$-ray production in neutral-current neutrino-oxygen interactions at energies above 200 MeV”, Phys. Rev. Lett. 108, 052505 (2012). 10.1103/PhysRevLett.108.052505.
  * (11) A. M. Ankowski, O. Benhar and M. Sakuda, “Improving the accuracy of neutrino energy reconstruction in charged-current quasielastic scattering off nuclear targets”, Phys. Rev. D 91, 033005 (2015). 10.1103/PhysRevD.91.033005.
  * (12) M. Bernheim, A. Bussiere, J. Mougey, D. Royer, D. Tarnowski, S. Turck-Chieze, S. Frullani, S. Boffi, C. Giusti and F. D. Pacati, et al. “The Influence of Bound State and Optical Potentials of 1p Momentum Distributions Obtained From [12][]C and [16][]O$(e,e^{\prime}p$) Reactions”, Nucl. Phys. A 375, 381-404 (1982). 10.1016/0375-9474(82)90020-3.
  * (13) M. Leuschner, J. R. Calarco, F. W. Hersman, E. Jans, G. J. Kramer, L. Lapikas, G. van der Steenhoven, P. K. A. de Witt Huberts, H. P. Blok and N. Kalantar-Nayestanaki, et al. “Quasielastic proton knockout from [16][]O”, Phys. Rev. C 49, 955-967 (1994). 10.1103/PhysRevC.49.955.
  * (14) G. Van der Steenhoven, H. P. Blok, E. Jans, M. De Jong, L. Lapikas, E. N. M. Quint and P. K. A. de Witt Huberts, “Knockout of 1p protons From [12][]C Induced by the ($e,e^{\prime}p$) Reaction”, Nucl. Phys. A 480, 547-572 (1988). 10.1016/0375-9474(88)90463-0.
  * (15) A.L Fetter and J.D. Walecka, Quantum Theory of Many-Particle Systems (McGraw-Hill, New York, NY, 1971).
  * (16) O. Benhar, A. Fabrocini and S. Fantoni, “Occupation probabilities and hole-state strengths in nuclear matter”, Phys. Rev. C 41, R24-R27 (1990). 10.1103/PhysRevC.41.R24.
  * (17) A. E. L. Dieperink, T. De Forest, I. Sick and R. A. Brandenburg, “Quasielastic electron scattering on [3][]He”, Phys. Lett. B 63 (1976), 261-264. 10.1016/0370-2693(76)90258-6.
  * (18) C. Ciofi degli Atti, E. Pace and G. Salme, “Realistic nucleon-nucleon interactions and the three-body electrodisintegration of [3][]He”, Phys. Rev. C 21 (1980), 805-815. 10.1103/PhysRevC.21.805.
  * (19) H. Meier-Hajduk, C. Hajduk, P. U. Sauer and W. Theis, “Quasielastic electron scattering on [3][]He”, Nucl. Phys. A 395 (1983), 332-348. 10.1016/0375-9474(83)90047-7.
  * (20) W. J. W. Geurts, K. Allaart, W. H. Dickhoff and H. Muther, “Spectroscopic factors for nucleon knock-out from [16][]O at small missing energy”, Phys. Rev. C 53 (1996), 2207-2212. 10.1103/PhysRevC.53.2207
  * (21) G. A. Rijsdijk, W. J. W. Geurts, K. Allaart and W. H. Dickhoff, “Hole spectral function and two-particle-one-hole response propagator”, Phys. Rev. C 53 (1996), 201-213. 10.1103/PhysRevC.53.201.
  * (22) A. Polls, M. Radici, S. Boffi, W. H. Dickhoff and H. Muther, “High momentum proton removal from [16][] O and the $(e,e^{\prime}p)$ cross section”, Phys. Rev. C 55 (1997), 810-819. 10.1103/PhysRevC.55.810.
  * (23) N. Rocco and C. Barbieri, “Inclusive electron-nucleus cross section within the Self Consistent Green’s Function approach”, Phys. Rev. C 98 (2018) no.2, 025501 10.1103/PhysRevC.98.025501.
  * (24) A. Ramos, A. Polls and W. H. Dickhoff, “Single-particle properties and short-range correlations in nuclear matter”, Nucl. Phys. A 503 (1989), 1-52. 10.1016/0375-9474(89)90252-2.
  * (25) C. Ciofi degli Atti and S. Liuti, “Realistic microscopic approach to deep inelastic scattering of electrons off few nucleon systems”, Phys. Rev. C 41 (1990), 1100-1114. 10.1103/PhysRevC.41.1100.
  * (26) H. Morita and T. Suzuki, “A realistic spectral function of [4][]He”, Prog. Theor. Phys. 86 (1991), 671-684. 10.1143/PTP.86.671.
  * (27) O. Benhar and V. R. Pandharipande, “Scattering of GeV electrons by light nuclei,” Phys. Rev. C 47 (1993), 2218-2227. 10.1103/PhysRevC.47.2218.
  * (28) C. Barbieri, N. Rocco and V. Somà, “Lepton Scattering from [40][]Ar and Ti in the quasielastic peak region,” Phys. Rev. C 100 (2019) no.6, 062501. 10.1103/PhysRevC.100.062501.
  * (29) C. Giusti and F. D. Pacati, “Electron distortion in quasifree $(e,e^{\prime}p)$ reactions”, Nucl. Phys. A 473 (1987), 717-735. 10.1016/0375-9474(87)90276-4.
  * (30) C. Giusti and F. D. Pacati, “Separation of structure functions and electron distortion in quasifree (e, e’p) reactions”, Nucl. Phys. A 485 (1988), 461-480. 10.1016/0375-9474(88)90548-9.
  * (31) C. Barbieri, D. Rohe, I. Sick and L. Lapikas, “Effect of kinematics on final state interactions in $(e,e^{\prime}p)$ reactions”, Phys. Lett. B 608 (2005), 47-52. 10.1016/j.physletb.2004.12.072.
  * (32) D. Rohe, C. S. Armstrong, R. Asaturyan, O. K. Baker, S. Bueltmann, C. Carasco, D. Day, R. Ent, H. C. Fenker and K. Garrow, et al. “Correlated strength in nuclear spectral function”, Phys. Rev. Lett. 93 (2004), 182501. 10.1103/PhysRevLett.93.182501.
  * (33) M. Alvioli, C. Ciofi degli Atti and H. Morita, “Universality of nucleon-nucleon short-range correlations: the factorization property of the nuclear wave function, the relative and center-of-mass momentum distributions, and the nuclear contacts”, Phys. Rev. C 94 (2016), 044309. 10.1103/PhysRevC.94.044309.
  * (34) P. A. Zyla et al. (Particle Data Group), “Review of Particle Physics”, Prog. Theor. Exp. Phys. 2020, 083C01 (2020). 10.1093/ptep/ptaa104.
  * (35) Nuclear Data Center, Japan Atomic Energy Agency (JAEA). https://wwwndc.jaea.go.jp/NuC/.
  * (36) R. Bevington and D.K. Robinson, Data reduction and error analysis for the physical sciences (McGraw-Hill, New York, NY, 2003).
  * (37) P. Barreau, M. Bernheim, J. Duclos, J. M. Finn, Z. Meziani, J. Morgenstern, J. Mougey, D. Royer, B. Saghai and D. Tarnowski, et al. “Deep Inelastic electron Scattering from Carbon,” Nucl. Phys. A 402, 515-540 (1983). 10.1016/0375-9474(83)90217-8.
  * (38) F. Ajzenberg-Selove, “Energy levels of light nuclei $A=11-12$,” Nuclear Physics A, 506, 1-158 (1990). 10.1016/0375-9474(90)90271-M.
  * (39) A. Ershova, K. Niewczas, S. Bolognesi, A. Letourneau, J. C. David, J. L. Rodríguez-Sánchez, J. T. Sobczyk, A. Blanchet, M. B. Avanzini and J. Chakrani, et al. “Role of deexcitation in the final-state interactions of protons in neutrino-nucleus interactions”, Phys. Rev. D 108, 112008 (2023). 10.1103/PhysRevD.108.112008.
  * (40) J. Kameda, “Observation of deexcitation gamma rays from nuclei in 1-kton detector in K2K experiment,” Nucl. Phys. B Proc. Suppl. 159, 44-49 (2006). 10.1016/j.nuclphysbps.2006.08.076.
  * (41) S. Abe et al. [KamLAND], “First measurement of the strange axial coupling constant using neutral-current quasielastic interactions of atmospheric neutrinos at KamLAND,” Phys. Rev. D 107, 072006 (2023). 10.1103/PhysRevD.107.072006.
  * (42) F. An et al. [JUNO], “Neutrino Physics with JUNO,” J. Phys. G 43, 030401 (2016). 10.1088/0954-3899/43/3/030401.
  * (43) S. Abe, “Nuclear deexcitation simulator for neutrino interactions and nucleon decays of ${}^{12}\mathrm{C}$ and ${}^{16}\mathrm{O}$ based on TALYS,” Phys. Rev. D 109, 036009 (2024). 10.1103/PhysRevD.109.036009.