# Weyl phonons in chiral crystals

Tiantian Zhang<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical Physics,
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China  Zhiheng Huang  This author contribute equally to this work. Beijing
National Laboratory for Condensed Matter Physics, Institute of Physics,
Chinese Academy of Sciences, Beijing 100190, China  University of Chinese
Academy of Sciences, Beijing 100049, China  Zitian Pan  Beijing National
Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy
of Sciences, Beijing 100190, China  University of Chinese Academy of Sciences,
Beijing 100049, China  Luojun Du<EMAIL_ADDRESS>Beijing National
Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy
of Sciences, Beijing 100190, China  Guangyu Zhang  Beijing National Laboratory
for Condensed Matter Physics, Institute of Physics, Chinese Academy of
Sciences, Beijing 100190, China  Songshan Lake Materials Laboratory, Dongguan
523808, China  Shuichi Murakami  Department of Physics, Tokyo Institute of
Technology, Okayama, Meguro-ku, Tokyo 152-8551, Japan

###### Abstract

### Abstract

Chirality is an indispensable concept that pervades fundamental science and
nature, manifesting itself in diverse forms, $e.g.$, quasiparticles and
crystal structures. Of particular interest are Weyl phonons carrying specific
Chern numbers and chiral phonons doing circular motions. Up to now, they have
been studied independently and the interpretations of chirality seem to be
different in these two concepts, impeding our understanding. Here, we
demonstrate that they are entangled in chiral crystals. Employing a typical
chiral crystal of elementary tellurium (Te) as a case study, we expound on the
intrinsic relationship between Chern number of Weyl phonons and pseudo-angular
momentum (PAM, $l_{ph}$) of chiral phonons. We propose Raman scattering as a
new technique to demonstrate the existence of Weyl phonons in Te, by detecting
the chirality-induced energy splitting between the two constituent chiral
phonon branches for Weyl phonons. Moreover, we also observe the obstructed
phonon surface states for the first time.

Weyl phonon, chiral phonon, chiral crystal, Raman scattering, Chern number,
pseudo-angular momentum

The breaking of continuous space-time symmetries in condensed matter solids
has led to the emergence of novel quasiparticle excitations beyond the
conventional Weyl-Dirac-Majorana classification. Among them, Weyl
quasiparticles with definite chirality have garnered significant attention.
For example, unconventional Weyl fermions with various band degeneracies and
topological charges have been uncovered Kruthoff et al. (2017); Bradlyn et al.
(2016); Huang et al. (2016); Wieder et al. (2016); Zhu et al. (2016); Bouhon
and Black-Schaffer (2017); Tang et al. (2017); Chang et al. (2017); Rao et al.
(2019); Takane et al. (2019); Schröter et al. (2019); Lv et al. (2019);
Schröter et al. (2020); Miao et al. (2018); Stenull et al. (2016); Lu et al.
(2015); Chang et al. (2017), strengthening our understanding of elementary
particles. In addition to Weyl fermions, recent advances have revealed
topological elementary excitations for bosonic systems, such as Weyl phonons
in the THz energy range ($\sim$10 meV) Zhang et al. (2018); Li et al. (2018);
Xia et al. (2019); Zhang et al. (2020); Liu et al. (2019, 2020); Ding et al.
(2022). These excitations have drawn significant attention both theoretically
and experimentally due to their potential to trigger novel quantum physical
phenomena and applications, such as phonon Hall effects Zhang et al. (2010);
Qin et al. (2012); Tang and Cao (2021); Im et al. (2022) and topological
thermal devices Liu et al. (2017). However, unlike Weyl fermions, which can be
detected through various approaches such as angle-resolved photoemission
spectroscopy Lv et al. (2015); Xu et al. (2015) and electronic transport
measurements Hosur et al. (2012); Hosur and Qi (2013); Zhang et al. (2016);
Liu et al. (2019), the experimental probe of Weyl phonons is challenging and
has only been achieved by inelastic x-ray spectroscopy (IXS) on their bulk
wave functions in conjunction with first-principle calculations Miao et al.
(2018); Li et al. (2021).

Apart from Weyl phonons, recent studies also uncover the concept of “chiral
phonon”, characterized by negative circular polarization and pseudoangular
momentum (PAM) Zhang and Niu (2015); Zhang and Murakami (2022); Chen et al.
(2021); Streib (2021); Tsunetsugu and Kusunose (2023). Chiral phonons can
interact strongly with other circularly polarized excitations (e.g., photons)
and thus be easily detected experimentally Malard et al. (2009); Garanin and
Chudnovsky (2015); Nakane and Kohno (2018); Long et al. (2018); Tatsumi et al.
(2018); Tatsumi and Saito (2018); Drapcho et al. (2017); Chen et al. (2021,
2019); Li et al. (2019, 2019); Zhu et al. (2018); Ishito et al. (2023, 2023).
Up to now, Weyl and chiral phonons are studied independently, and the concept
of “chirality” seems to be different between these two concepts. The former is
related to its (pseudo-)spin and momentum coupling patterns and is defined in
reciprocal space, as shown in Figs. 1(a-b), while the latter corresponds to
phonon modes with atoms vibrating in a circular way in real space, as shown in
Figs. 1(c-d). If chiral phonons can entangle with Weyl phonons, they may
provide the possibility to demonstrate the existence of Weyl phonons.

Figure 1: Spin-momentum locking configurations for Weyl points with (a) C = +1
and (b) C = $-1$. (c-d) Illustration for phonon modes with $l_{ph}^{s}$ = + 1
and $l_{ph}^{s}$ = $-1$, respectively. The insert figure shows the top view of
Te crystal with a right-hand structure.

In this study, we demonstrate that Weyl and chiral phonons are entangled in Te
and further experimentally validate our theory by utilizing helicity-resolved
Raman spectra. The evidence for the existence of Weyl phonons is revealed by
directly tracking the chirality-induced energy splitting between the two
constituent chiral phonon branches. Additionally, obstructed surface phonon
states have been uncovered by Raman spectra for the first time. All the
results are in perfect agreement with first-principles calculations,
establishing Raman spectroscopy as a fast, non-destructive, and high spatial
resolution probe for exploring topological phonons in chiral crystal
structures, as well as the obstructed surface modes.

The Hamiltonian in the vicinity of conventional Weyl points can be succinctly
characterized by $H(\mathbf{k})\propto\mathbf{S}\cdot\mathbf{k}$, where
$\mathbf{S}$ denotes the rotation generator for spin-$\frac{1}{2}$ particles.
Two bands of Weyl points display spin components of
$\mathbf{S}\cdot\hat{k}=\pm\frac{\hbar}{2}$, leading to an inward or outward
hedgehog configuration on each equal-energy surface. These distinct spin-
momentum-locking properties determine the topological invariant referred to as
the Chern number, which assumes the values of C = $\pm 1$ and serves to
differentiate between those two cases, as shown in Figs. 1 (a) and (b),
respectively. In addition, the sign of the Chern number corresponds to the
chirality of the Weyl points. The conventional Weyl point with C = 1 may
transform into an unconventional one when additional crystalline symmetries
exist, such as a double Weyl point with C = 2 under threefold (screw) rotation
symmetry Liu and Vanderbilt (2014); Fang et al. (2015); Zhang et al. (2020).

Figure 2: (a) Side (left panel) and top (right panel) views of elemental
crystal tellurium with right-hand structure. (b-c) Brillouin zone and phonon
band structure of tellurium. The twofold degenerate band marked by the light
grey box at $\Gamma$ is a double Weyl phonon with C = +2. (d) Dynamical
structural factor simulation on the double Weyl phonon with C = +2 by
following the $k$-path of $(4.0,5.0,-0.5)$-$(4.0,5.0,0.0)$-$(4.0,5.0,0.5)$.
(e1) Illustration for the bulk dispersion (solid cone) and surface states
(yellow sheets) for the double Weyl phonon at $\Gamma$ on the (010) surface.
(e2) Surface arcs for the double Weyl phonon in tellurium.

Tellurium crystals have two enantiomers, classified as the right-hand
structure with space group $P3_{1}21$ or the left-hand structure with space
group $P3_{2}21$. Figures 2 (a-c) depict the crystal structure, Brillouin zone
(BZ), and phonon spectra for the right-hand tellurium crystal with space
group, where the three Raman active branches are denoted by $E^{1}$, $A_{1}$,
and $E^{2}$ with increasing frequency at $\Gamma$. According to first-
principles calculations, the twofold phonon modes with $E^{1}$ and $E^{2}$
irreducible representations are both double Weyl phonons with $|$C$|$ = 2,
protected by threefold screw rotation symmetry ($C_{3}^{\prime}$) and time-
reversal symmetry ($\mathcal{T}$) Fang et al. (2012); Tsirkin et al. (2017);
Yang and Nagaosa (2014); Jian and Yao (2015), and they are beyond the
traditional classification. The low-energy effective
$\mathbf{k}\cdot\mathbf{p}$ model for $E^{1}$ can be written as:

$H(\mathbf{k})=\begin{pmatrix}Ak_{z}&B(k_{x}-ik_{y})^{2}\\\
B(k_{x}+ik_{y})^{2}&-Ak_{z}\end{pmatrix}$ (1)

where A and B are real constants. The Hamiltonian gives rise to a Weyl point
with C = $+2$.

Figure 2 (d) illustrates the phonon dynamical structure factor (DSF) of Te
simulated by first-principle calculations for the $E^{1}$ mode along the
$k$-path of A$(4.0,5.0,-0.5)$-$\Gamma(4.0,5.0,0.0)$-A$(4.0,5.0,0.5)$. The
nontrivial DSF intensity distribution arising from the topological bulk wave
functions provides the possibility to verify the double Weyl phonon using IXS
Miao et al. (2018); Zhang et al. (2019); Li et al. (2021); Jin et al. (2022).
In addition to the bulk topology property, Weyl phonons can also show
topological surface states due to the bulk-surface correspondence (BSC).
Figure 2 (e1) shows a double Weyl phonon with C = +2 (solid cone) associated
with two helical surface states (yellow sheets), rotating in a right-hand way
with increasing frequency and connecting to the bulk Weyl bands, sharing the
same chirality with the Weyl point. Figure 2 (e2) displays the iso-energetic
contour for the $E^{1}$ Weyl phonon along the (010) direction, where the
double Weyl phonon is projected onto the center of the surface BZ, associated
with two surface arcs. It is noteworthy that $E^{1}$ Weyl modes carry opposite
Chern numbers between two enantiomers of chiral crystal tellurium, and the
Chern number alters from +2 in the right-hand crystal to $-2$ in the left-hand
crystal, as shown in Figs. 3 (a4) and (b4).

Figure 3: Entanglement between the chirality of crystal structure, Weyl
phonon, and PAM. (a1) and (b1) are tellurium with right-hand and left-hand
chiral crystal structures, respectively. (a2) and (b2) are the atomic
vibration for phonon modes with the same energy in different chiral crystal
structures, showing opposite chiral motions in real space, and their PAMs
shown in (a3) and (b3) are defined in terms of the threefold screw symmetry
with respect to $k_{z}$ direction. (a4) and (b4) are double Weyl phonons
carrying opposite Chern numbers for (a1) and (b1), respectively.

For 3D chiral systems with (screw) rotation symmetry, where both mirror and
inversion symmetries are absent, the lattice chirality is closely related to
the circular polarization and propagation direction of phonons Chen et al.
(2022). Thus, in addition to describing topological properties, chirality can
be also used to catalog chiral phonons with left/right-hand circular vibrating
motions, carrying nonzero phonon circular polarization and a noninteger PAM
$l_{ph}$. Since chiral phonons that we study in this paper are located in the
vicinity of $\Gamma$, PAM $l_{ph}$ can be simplified as the intracell part
$l_{ph}^{s}$, which describes the phase difference of rotational motions
between two neighboring atoms under the (screw) rotation symmetry Zhang and
Murakami (2022). As the example shown in Fig. 1 (c), Te1${}^{{}^{\prime}}$ is
rotated from Te1 by $C_{3}^{{}^{\prime}}$ symmetry in the right-hand crystal,
having a $\frac{2}{3}\pi$ phase advancing to Te2, and therefore this case
corresponds to $l_{ph}^{s}$ = +1. Likewise, Fig. 1 (d) corresponds to
$l_{ph}^{s}$ = $-1$.

Figures 3 (a2) and (b2) show the two-dimensional projections of atomic motions
for the lower branches of $E^{1}$ modes along $\Gamma$-A path (i.e., $k_{z}$
direction) in right- and left-hand crystals, respectively. The lower branch of
$E^{1}$ mode along $\Gamma$-A path in the right-hand (left-hand) crystal does
counter-clockwise (clockwise) circular motions, which indicates the opposite
$l_{ph}^{s}$ between right- and left-hand crystals. According to the three-
fold rotation symmetry operation, $l_{ph}^{s}$ for the lower branches of
$E^{1}$ modes along $\Gamma$-A path is +1 ($-1$) in right-hand (left-hand)
crystal, as shown in Fig. 3 (a3) [Fig. 3 (b3)]. By comparing Fig. 3 (a4) and
Fig. 3 (a3) [Fig. 3 (b4) and Fig. 3 (b3)], we conclude that Weyl phonons in Te
are entangled with chiral phonons, as well as their chiral crystal structures,
i.e., the Weyl phonon bands with C = +2 (C = $-2$) carries $l_{ph}^{s}$ = +1
($-1$) for the lower phonon mode of $E^{1}$ along $\Gamma$-A path with the
right(left)-hand crystal structure.

Figure 4: Weyl phonons detected by helicity-resolved Raman scattering. (a)
Backscattering geometry of helicity-resolved Raman spectroscopy on Te, where
the incident/scattered photon propagates along the c axis. (b1)/(c1) Schematic
diagram of $E^{1}$ phonon dispersion along the $\Gamma-$A path, which is
parallel to the c axis. The red (green) branch corresponds to chiral phonon
mode with $l_{ph}^{s}$ = +1 ($-1$), which can be selectively detected under
incident/scattered polarization configuration of $\sigma^{+}$/$\sigma^{-}$
($\sigma^{-}$/$\sigma^{+}$) by backscattering Raman spectroscopy. $\sigma^{+}$
($\sigma^{-}$) denotes left (right) circularly polarized light and propagates
along the screw rotation axis. (b2)-(c2) Helicity-resolved Raman spectra in
the $\sigma^{+}$/$\sigma^{-}$ (green) and $\sigma^{-}$/$\sigma^{+}$ (red)
polarization conﬁgurations at 10 K for 633 nm excitation. (b3)-(c3) Zoom-in
Raman spectra of $E^{1}$ mode, which clearly show the chirality-induced band
splitting for the Weyl phonon. (b1)-(b3) are the results for Te with a right-
hand crystal structure, while (c1)-(c3) are the results for Te with a left-
hand crystal structure.

As discussed above, Weyl phonons in chiral crystal Te are entangled with two
chiral phonon branches holding opposite $l^{s}_{ph}$. This entanglement
enables the detection of Weyl phonons through the fast, non-destructive, and
high spatial resolution probe of helicity-resolved Raman scattering along the
screw rotation axis (i.e., c axis), as illustrated in Fig. 4 (a). Typically,
three conservation rules on the PAM, the pseudo-momentum, and the energy
should be satisfied in the Raman scattering: (i) $l_{i}\ -$ $l_{s}$ =
$l_{ph}^{s}$ $modulo$ 3, where $l_{i}$ ($l_{s}$) is psedoangular momentum for
the incident (scattered) photons, and modulo three can be understood by
considering the threefold screw rotation symmetry of Te and the Umklapp
process Ishito et al. (2023); Du et al. (2019); (ii) The momentum $N$ for
phonons detected by backscattering Raman spectroscopy should arise from the
wavenumber of photons Pine and Dresselhaus (1971); (iii) Raman shift equal to
the energy of the observed phonon modes at momentum $N$ shown in Figs. 4 (b1)
and (c1). Here, with the wavelength $\lambda=633$ nm for the incident light,
the corresponding wavenumber of momentum $N$ is
$k=\frac{4\pi|n+i\eta|}{\lambda}\simeq 0.8\times 10^{6}$ cm-1, which is about
$\frac{1}{120}$ of the Brillouin zone along the c axis direction. $n=2.3$ and
$\eta=3.5$ are the real and imaginary parts of the complex refractive index of
Te, respectively Pine and Dresselhaus (1971). Since the nonzero momentum leads
to the energy splitting of the Weyl phonon carrying opposite chirality, Weyl
phonons can be diagnosed by Raman spectroscopy.

As depicted in Figure 3 (a3), $l^{s}_{ph}$ is equal to $\pm 1$ for the lower
and upper branches of $E^{1}$ in Te with a right-hand crystal structure,
respectively. In light of conservation rules, the upper (lower) branch can be
selectivity detected under incident/scattered polarization configurations of
$\sigma^{+}$/$\sigma^{-}$ ($\sigma^{-}$/$\sigma^{+}$), as shown in Fig. 4
(b1). Figure 4 (b2) is the polarized Raman spectra of right-hand Te, observed
by using cross-circularly polarized detection with different circularly
polarized light at 10 K, recorded with 633 nm excitation. The spectra exhibit
two $E$ phonon modes ($E^{1}$: $\sim 99$ cm-1 and $E^{2}$: $\sim 146$ cm-1)
and one $A_{1}$ mode ($\sim 127$ cm-1), which are in good agreement with
first-principles calculations and previous results Torrie (1970); Pine and
Dresselhaus (1971). Notably, as a nonchiral phonon mode, $A_{1}$ peaks appear
almost identical in both configurations, while $E^{1}$ and $E^{2}$ phonon
modes show significant energy splitting. Close-ups of Raman spectra for
$E^{1}$ phonon mode (hollow circles) are shown in Fig. 4 (b3), fitted by the
Voigt function (solid lines). The chirality-induced band splitting extracted
from fitting is $\sim 1.15$ cm-1 for $E^{1}$ phonon mode, demonstrating the
existence of two chiral phonons carrying opposite $l_{ph}^{s}$ and showing the
existence of Weyl phonons in Te.

Figure 4 (c2) is the polarized Raman spectra of left-hand Te, which are in
close resemblance to those of right-hand Te but exhibit opposing energy splits
between $\sigma^{+}$/$\sigma^{-}$ and $\sigma^{-}$/$\sigma^{+}$
configurations. These results give information of $l_{ph}^{s}$ for the two
chiral phonons and, consequently, the Chern number of the Weyl phonon in Te.
Furthermore, both $l_{ph}^{s}$ and the Chern number of each band will change
their sign if the chirality of the crystal is altered, showing good agreement
with Figure 3.

Figure 5: Demonstration for obstructed surface phonon modes. (a) First-
principles calculations on right-hand Te with a slab configuration along [010]
direction. Black solid lines at $\sim$100 cm-1 are the obstructed surface
phonon modes. (b) Backscattering geometry of helicity-resolved Raman
spectroscopy on Te, where the incident and scattered photons are perpendicular
to the c axis. (c) Co- (black) and cross- (red) circularly polarized Raman
spectra under $\sigma^{-}$ excitation at 10 K. Inset: Fitting plot for the
measured and calculated phonon energies, which shows a perfect match.

The symmetry data (irreducible representation) for the lowest six phonon bands
of Te is equal to that of an atomic insulator, so topological surface states
are not expected in the frequency gap around 100 $cm^{-1}$. However, the
symmetry data for the lowest six phonon bands is induced from the atomic
positions where no real atoms reside, which makes Te belong to the recently
proposed “obstructed atomic insulator”, which can be diagnosed by the real
space topological invariant Po et al. (2018); Cano et al. (2018); Khalaf et
al. (2021); Zhang et al. (2023). As a consequence, “metallic” surface states
separate from the bulk bands will appear when the crystal is cleaved at these
virtual sites on the surface.

Figure 5 (a) illustrates the surface phonon spectra obtained from first-
principles calculations along the [010] direction, wherein the obstructed
surface states marked by the black solid lines originate from the mismatch
between the atomic position and the vibrating mass center. Since they spread
over the whole surface Brillouin zone and are in the gap of bulk phonon
spectra, they may therefore be observed by Raman scattering.

Figure 5 (c) presents the co- (black) and cross- (red) circularly polarized
Raman spectra under $\sigma{{}^{-}}$ excitation at 10 K, with the incident and
scattered photons perpendicular to the c axis, as depicted in Figure 5 (b).
Remarkably, a new phonon mode $S$ at $\sim 115$ cm-1 is detected under the
cross-circularly polarized configuration ($\sigma^{-}$/$\sigma^{+}$) in
addition to three bulk phonon modes. The energies of the four experimentally
observed modes are in excellent agreement with the first-principles calculated
ones, as shown in the inset of Figure 5 (c), which confirms the $S$ mode being
the obstructed surface phonon state. In the same setup, when the surface is
destroyed by a high-power laser, the $S$ mode will be efficiently quenched,
which further proves that $S$ mode is a surface mode (See details in the
Supplementary Material).

We highlight that, by combining Raman scattering and first-principle
calculations, the existence of obstructed surface phonon modes is demonstrated
for the first time. Note that surface phonon modes associated with Weyl
phonons also spread widely in the surface Brillouin. However, they have the
degeneracy in energy at $\Gamma$, which can not be distinguished by the Raman
scattering.

The interplay between Weyl phonons and chiral phonons has received limited
attention due to their independent investigations in different systems. In
this study, we propose that the Weyl phonons carrying nonzero Chern numbers,
including $\pm 1$, $\pm 2$, $et\ al.$, are intricately entangled with chiral
phonons doing circular motions, with specific propagation directions and PAM,
in systems with $C_{3/4/6}$ (screw) rotation symmetries. We present a novel
experimental approach to probing and characterizing the chirality of Weyl
phonons in chiral crystals that possess (screw) rotation symmetry, i.e.,
circularly polarized Raman scattering. Our approach is testified in chiral
crystal Te, combined with first-principles calculations. Furthermore, we
report the observation of obstructed phonon surface states for the first time
by using Raman scattering and first-principles calculations, providing new
insights into the phonon physics of chiral materials.

##  I Acknowledge

This work is supported by the Strategic Priority Research Program of Chinese
Academy of Sciences (CAS) under the grant No. XDB30000000, the National
Natural Science Foundation of China (Grant No. 12047503), and Japan Society
for the Promotion of Science (JSPS), KAKENHI Grant No. JP21K13865, No.
JP22K18687 and JP22H00108, MEXT X-NICS under Grant No. JPJ011438, the National
Key R&D Program of China under grant No. 2021YFA1202900, the National Science
Foundation of China (Grant Nos. 12274447, 12074412, 61888102, and 11834017).

##  II References

## References

  * Kruthoff et al. (2017) Kruthoff, J.; de Boer, J.; van Wezel, J.; Kane, C. L.; Slager, R.-J. _Phys. Rev. X_ 2017, _7_ , 041069.
  * Bradlyn et al. (2016) Bradlyn, B.; Cano, J.; Wang, Z.; Vergniory, M. G.; Felser, C.; Cava, R. J.; Bernevig, B. A. _Science_ 2016, _353_.
  * Huang et al. (2016) Huang, S.-M. et al. _Proceedings of the National Academy of Sciences_ 2016, _113_ , 1180–1185.
  * Wieder et al. (2016) Wieder, B. J.; Kim, Y.; Rappe, A. M.; Kane, C. L. _Phys. Rev. Lett._ 2016, _116_ , 186402.
  * Zhu et al. (2016) Zhu, Z.; Winkler, G. W.; Wu, Q.; Li, J.; Soluyanov, A. A. _Phys. Rev. X_ 2016, _6_ , 031003.
  * Bouhon and Black-Schaffer (2017) Bouhon, A.; Black-Schaffer, A. M. _Phys. Rev. B_ 2017, _95_ , 241101\.
  * Tang et al. (2017) Tang, P.; Zhou, Q.; Zhang, S.-C. _Phys. Rev. Lett._ 2017, _119_ , 206402.
  * Chang et al. (2017) others,, et al. _Physical Review Letters_ 2017, _119_ , 206401.
  * Rao et al. (2019) others,, et al. _Nature_ 2019, _567_ , 496–499.
  * Takane et al. (2019) others,, et al. _Physical Review Letters_ 2019, _122_ , 076402.
  * Schröter et al. (2019) others,, et al. _Nature Physics_ 2019, _15_ , 759–765.
  * Lv et al. (2019) others,, et al. _Physical Review B_ 2019, _99_ , 241104\.
  * Schröter et al. (2020) Schröter, N. B. M. et al. _Science_ 2020, _369_ , 179–183.
  * Miao et al. (2018) Miao, H.; Zhang, T.; Wang, L.; Meyers, D.; Said, A.; Wang, Y.; Shi, Y.; Weng, H.; Fang, Z.; Dean, M. _Physical Review Letters_ 2018, _121_ , 035302.
  * Stenull et al. (2016) Stenull, O.; Kane, C.; Lubensky, T. _Physical Review Letters_ 2016, _117_ , 068001.
  * Lu et al. (2015) Lu, L.; Wang, Z.; Ye, D.; Ran, L.; Fu, L.; Joannopoulos, J. D.; Soljačić, M. _Science_ 2015, _349_ , 622–624.
  * Chang et al. (2017) Chang, M.-L.; Xiao, M.; Chen, W.-J.; Chan, C. T. _Physical Review B_ 2017, _95_ , 125136.
  * Zhang et al. (2018) Zhang, T.; Song, Z.; Alexandradinata, A.; Weng, H.; Fang, C.; Lu, L.; Fang, Z. _Phys. Rev. Lett._ 2018, _120_ , 016401.
  * Li et al. (2018) Li, J.; Xie, Q.; Ullah, S.; Li, R.; Ma, H.; Li, D.; Li, Y.; Chen, X.-Q. _Physical Review B_ 2018, _97_ , 054305.
  * Xia et al. (2019) Xia, B.; Wang, R.; Chen, Z.; Zhao, Y.; Xu, H. _Physical Review Letters_ 2019, _123_ , 065501.
  * Zhang et al. (2020) Zhang, T.; Takahashi, R.; Fang, C.; Murakami, S. _Physical Review B_ 2020, _102_ , 125148.
  * Liu et al. (2019) Liu, J.; Hou, W.; Wang, E.; Zhang, S.; Sun, J.-T.; Meng, S. _Physical Review B_ 2019, _100_ , 081204.
  * Liu et al. (2020) Liu, Q.-B.; Qian, Y.; Fu, H.-H.; Wang, Z. _npj Computational Materials_ 2020, _6_ , 95.
  * Ding et al. (2022) Ding, G.; Zhou, F.; Zhang, Z.; Yu, Z.-M.; Wang, X. _Physical Review B_ 2022, _105_ , 134303.
  * Zhang et al. (2010) Zhang, L.; Ren, J.; Wang, J.-S.; Li, B. _Physical review letters_ 2010, _105_ , 225901.
  * Qin et al. (2012) Qin, T.; Zhou, J.; Shi, J. _Physical Review B_ 2012, _86_ , 104305\.
  * Tang and Cao (2021) Tang, D.-S.; Cao, B.-Y. _Journal of Applied Physics_ 2021, _129_ , 085102.
  * Im et al. (2022) Im, J.; Kim, C. H.; Jin, H. _Nano Letters_ 2022, _22_ , 8281–8286.
  * Liu et al. (2017) Liu, Y.; Xu, Y.; Zhang, S.-C.; Duan, W. _Physical Review B_ 2017, _96_ , 064106.
  * Lv et al. (2015) others,, et al. _Physical Review X_ 2015, _5_ , 031013\.
  * Xu et al. (2015) others,, et al. _Science_ 2015, _349_ , 613–617.
  * Hosur et al. (2012) Hosur, P.; Parameswaran, S.; Vishwanath, A. _Physical Review Letters_ 2012, _108_ , 046602.
  * Hosur and Qi (2013) Hosur, P.; Qi, X. _Comptes Rendus Physique_ 2013, _14_ , 857–870.
  * Zhang et al. (2016) others,, et al. _Nature communications_ 2016, _7_ , 1–9.
  * Liu et al. (2019) others,, et al. _Science_ 2019, _365_ , 1282–1285.
  * Li et al. (2021) Li, H.; Zhang, T.; Said, A.; Fu, Y.; Fabbris, G.; Mazzone, D. G.; Zhang, J.; Lapano, J.; Lee, H. N.; Lei, H. C.; Dean, M. P. M.; Murakami, S.; Miao, H. _Phys. Rev. B_ 2021, _103_ , 184301.
  * Zhang and Niu (2015) Zhang, L.; Niu, Q. _Physical Review Letters_ 2015, _115_ , 115502.
  * Zhang and Murakami (2022) Zhang, T.; Murakami, S. _Phys. Rev. Res._ 2022, _4_ , L012024.
  * Chen et al. (2021) Chen, H.; Wu, W.; Zhu, J.; Yang, S. A.; Zhang, L. _Nano Letters_ 2021, _21_ , 3060–3065.
  * Streib (2021) Streib, S. _Phys. Rev. B_ 2021, _103_ , L100409.
  * Tsunetsugu and Kusunose (2023) Tsunetsugu, H.; Kusunose, H. _Journal of the Physical Society of Japan_ 2023, _92_ , 023601.
  * Malard et al. (2009) Malard, L.; Pimenta, M. A.; Dresselhaus, G.; Dresselhaus, M. _Physics reports_ 2009, _473_ , 51–87.
  * Garanin and Chudnovsky (2015) Garanin, D. A.; Chudnovsky, E. M. _Phys. Rev. B_ 2015, _92_ , 024421\.
  * Nakane and Kohno (2018) Nakane, J. J.; Kohno, H. _Phys. Rev. B_ 2018, _97_ , 174403\.
  * Long et al. (2018) Long, Y.; Ren, J.; Chen, H. _Proceedings of the National Academy of Sciences_ 2018, _115_ , 9951–9955.
  * Tatsumi et al. (2018) Tatsumi, Y.; Kaneko, T.; Saito, R. _Phys. Rev. B_ 2018, _97_ , 195444.
  * Tatsumi and Saito (2018) Tatsumi, Y.; Saito, R. _Physical Review B_ 2018, _97_ , 115407\.
  * Drapcho et al. (2017) Drapcho, S. G.; Kim, J.; Hong, X.; Jin, C.; Shi, S.; Tongay, S.; Wu, J.; Wang, F. _Physical Review B_ 2017, _95_ , 165417.
  * Chen et al. (2021) Chen, C.; Chen, X.; Deng, B.; Watanabe, K.; Taniguchi, T.; Huang, S.; Xia, F. _Physical Review B_ 2021, _103_ , 035405.
  * Chen et al. (2019) Chen, X.; Lu, X.; Dubey, S.; Yao, Q.; Liu, S.; Wang, X.; Xiong, Q.; Zhang, L.; Srivastava, A. _Nature Physics_ 2019, _15_ , 221–227.
  * Li et al. (2019) others,, et al. _ACS Nano_ 2019, _13_ , 14107–14113.
  * Li et al. (2019) others,, et al. _Nature Communications_ 2019, _10_ , 1–7.
  * Zhu et al. (2018) Zhu, H.; Yi, J.; Li, M.-Y.; Xiao, J.; Zhang, L.; Yang, C.-W.; Kaindl, R. A.; Li, L.-J.; Wang, Y.; Zhang, X. _Science_ 2018, _359_ , 579–582.
  * Ishito et al. (2023) Ishito, K.; Mao, H.; Kousaka, Y.; Togawa, Y.; Iwasaki, S.; Zhang, T.; Murakami, S.; Kishine, J.-i.; Satoh, T. _Nature Physics_ 2023, _19_ , 35–39.
  * Ishito et al. (2023) Ishito, K.; Mao, H.; Kobayashi, K.; Kousaka, Y.; Togawa, Y.; Kusunose, H.; Kishine, J.-i.; Satoh, T. _Chirality_ 2023, _35_ , 338–345.
  * Liu and Vanderbilt (2014) Liu, J.; Vanderbilt, D. _Phys. Rev. B_ 2014, _90_ , 155316\.
  * Fang et al. (2015) Fang, C.; Chen, Y.; Kee, H.-Y.; Fu, L. _Physical Review B_ 2015, _92_ , 081201.
  * Fang et al. (2012) Fang, C.; Gilbert, M. J.; Dai, X.; Bernevig, B. A. _Phys. Rev. Lett._ 2012, _108_ , 266802.
  * Tsirkin et al. (2017) Tsirkin, S. S.; Souza, I.; Vanderbilt, D. _Phys. Rev. B_ 2017, _96_ , 045102.
  * Yang and Nagaosa (2014) Yang, B.-J.; Nagaosa, N. _Nature Communications_ 2014, _5_ , 4898.
  * Jian and Yao (2015) Jian, S.-K.; Yao, H. _Physical Review B_ 2015, _92_ , 045121\.
  * Zhang et al. (2019) others,, et al. _Physical Review Letters_ 2019, _123_ , 245302.
  * Jin et al. (2022) others,, et al. _Physical Review B_ 2022, _106_ , 224304\.
  * Chen et al. (2022) Chen, H.; Wu, W.; Zhu, J.; Yang, Z.; Gong, W.; Gao, W.; Yang, S. A.; Zhang, L. _Nano Letters_ 2022, _22_ , 1688–1693.
  * Du et al. (2019) others,, et al. _Advanced Functional Materials_ 2019, _29_ , 1904734.
  * Pine and Dresselhaus (1971) Pine, A. S.; Dresselhaus, G. _Phys. Rev. B_ 1971, _4_ , 356–371.
  * Torrie (1970) Torrie, B. H. _Solid State Communications_ 1970, _8_ , 1899–1901.
  * Po et al. (2018) Po, H. C.; Watanabe, H.; Vishwanath, A. _Physical review letters_ 2018, _121_ , 126402.
  * Cano et al. (2018) Cano, J.; Bradlyn, B.; Wang, Z.; Elcoro, L.; Vergniory, M. G.; Felser, C.; Aroyo, M. I.; Bernevig, B. A. _Physical Review B_ 2018, _97_ , 035139.
  * Khalaf et al. (2021) Khalaf, E.; Benalcazar, W. A.; Hughes, T. L.; Queiroz, R. _Physical Review Research_ 2021, _3_ , 013239.
  * Zhang et al. (2023) Zhang, R.; Deng, J.; Sun, Y.; Fang, Z.; Guo, Z.; Wang, Z. _Phys. Rev. Res._ 2023, _5_ , 023142.


Impact of big bang nucleosynthesis on the $H_{0}$ tension

Tomo Takahashi$\,{}^{1}$ and Yo Toda$\,{}^{2}$

1Department of Physics, Saga University, Saga 840-8502, Japan
2Department of Physics, Hokkaido University, Sapporo 060-0810, Japan

We investigate the impact of big bang nucleosynthesis (BBN) on the Hubble
tension, focusing on how the treatment of the reaction rate and observational
data affect the evaluation of the tension. We show that the significance of
the tension can vary by $0.8\sigma$ in some early dark energy model, depending
on the treatment of the reaction rate and observational data. This indicates
that how we include the BBN data in the analysis can give a significant impact
on the Hubble tension, and we need to carefully consider the assumptions of
the analysis to evaluate the significance of the tension when the BBN data is
used.

EPHOU-23-011

##  1 Introduction

The so-called $\Lambda$CDM model has now been recognized as the standard
paradigm of cosmology. Indeed almost all observations can be understood in the
framework of the $\Lambda$CDM model, however, there have been several tensions
reported, among which the Hubble tension ($H_{0}$ tension) is one of the most
pronounced one where the values of the Hubble constant $H_{0}$ are
inconsistent at almost 5$\sigma$ level between the ones obtained from direct
measurements such as Cepheid calibrated distance ladder $H_{0}=73.30\pm
1.04~{}{\rm km/s/Mpc}$ [1]#1#1#1 Here we quote the value based on the analysis
including high-redshift supernovae.  and from indirect observations such as
cosmic microwave background (CMB) from Planck $H_{0}=67.66\pm 0.42~{}{\rm
km/s/Mpc}$ [2]. Even if different data are used, local direct observations
consistently infer a relatively higher value of $H_{0}$ compared to the one
obtained from Planck, on the other hand, indirect observations, even without
Planck data, gives $H_{0}$ consistent with that obtained from Planck#2#2#2
Even without CMB data, the combination of data from baryon acoustic
oscillation, type Ia supernovea and big bang nucleosynthesis also gives a
relatively low value of $H_{0}$ in the framework of $\Lambda$CDM and other
extended models [3, 4, 5].  (see, e.g., [6, 7] for reviews of the $H_{0}$
tension).

A lot of works have been devoted to resolving the $H_{0}$ tension by extending
the $\Lambda$CDM framework (for various models proposed to solve the tension,
see, e.g., reviews [6, 8]). In those models, the value of $H_{0}$ can become
higher than $\Lambda$CDM case when they are fitted to CMB in combination with
some other data such as baryon acoustic oscillation (BAO), type Ia supernova
(SNeIa), and so on. It would be important to notice that there exists the
correlation between $H_{0}$ and other cosmological parameters in the fit to
CMB#3#3#3 One can easily recognize such a correlation from 2D constraints in
the plane of $H_{0}$ and some other parameters (e.g., [2]). See also, e.g.,
[9, 10, 11] to see how $H_{0}$ is correlated with other cosmological
parameters in the structure of acoustic peaks in CMB power spectrum. .
Although such a correlation is somewhat involved in models proposed to resolve
the $H_{0}$ tension, the increase of the value of $H_{0}$ certainly affects
other aspects of cosmology. Indeed the baryon energy density generally tends
to become larger than that for the case with the $\Lambda$CDM in models
proposed to resolve the $H_{0}$ tension when fitted to CMB in combination with
other observations [12, 13, 14, 15]#4#4#4 It has also been argued that
constraints on the neutrino mass [16], the spectral index for primordial power
spectrum $n_{s}$ [17, 18, 19, 20] are also modified in models proposed to
resolve the tension such as time-varying electron mass [21] and early dark
energy models [22]. .

The baryon density can also be well determined by big bang nucleosynthesis
(BBN) and hence we can also include BBN data in the fitting to derive the
cosmological parameters. However, as we discuss in this paper, we need to be
cautious regarding what observational data we use and what nuclear reaction
rates are adopted, depending on which the derived values of the cosmological
parameters can change. Therefore when the BBN data is combined with CMB and
some other data, it can affect the value of $H_{0}$ in models to resolve the
$H_{0}$ tension through the correlations among the cosmological parameters in
the fitting. Since, among the light elements, the baryon density is mainly
determined by the primordial deuterium abundance, which we denote $D_{p}$ in
the following, we mainly investigate the impact of the treatment of the
primordial deuterium on the $H_{0}$ tension#5#5#5 The primordial helium
abundance $Y_{p}$ is also precisely measured [23, 24, 25, 26], however $Y_{p}$
is not so severely constrained by CMB data (see, e.g., [2, 27, 28]), and hence
the inclusion of $Y_{p}$ data would not affect $H_{0}$ much. Therefore we
focus on the treatment of the deuterium abundance in this paper although the
$Y_{p}$ data of [23] is included in our analysis. . To this end, we consider
several different treatments to include the primordial deuterium abundance in
the analysis: For observational data, we use two different data for $D_{p}$
from Cooke et al. [29] and the weighted mean of the recent 11 measurements
compiled by Particle Data Group [30]. For the theoretical calculations of
$D_{p}$, we adopt the nuclear reaction rate for $d(p,\gamma)^{3}{\rm He}$
obtained from the extrapolation of experiments [31, 32] and the theoretically
evaluated one via an ab-initio approach [33]. As we argue in this paper, the
significance of the $H_{0}$ tension can vary up to 0.8$\sigma$ depending on
the treatment of how we include $D_{p}$ in the analysis and models assumed,
which shows that BBN has a significant impact on the $H_{0}$ tension.

The organization of this paper is as follows. In the next section, we
summarize our analysis method and models assumed in the analysis. Then, in
Section 3, we show our results and argue that the treatment of BBN in the
analysis can give a significant impact on the $H_{0}$ tension depending on
what observational data and the reaction rate are adopted. We give the summary
of this paper in the final section.

##  2 Models and analysis method

###  2.1 Models

To investigate the impact of the treatment of the primordial deuterium
abundance on the $H_{0}$ tension, we consider two model frameworks: a flat
$\Lambda$CDM and early dark energy (EDE) models [22]. In the framework of a
flat $\Lambda$CDM model, we vary 6 standard cosmological parameters: baryon
density $\Omega_{b}h^{2}$, cold dark matter (CDM) density $\Omega_{c}h^{2}$,
the acoustic scale $\theta$, the reionization optical depth $\tau$, the
amplitude of primordial power spectrum $A_{s}$ and its spectral index $n_{s}$.
Since we assume a flat Universe, the Hubble constant $H_{0}$ is derived once
these parameters are given. The amplitude $A_{s}$ is defined at the reference
scale $k_{{\rm ref}}=0.05~{}{\rm Mpc}^{-1}$.

EDE model [34, 22, 35, 36] has been extensively studied as a possible
candidate to resolve the $H_{0}$ tension (for recent review on early dark
energy, see [37, 38]). There are several variations of EDE model discussed in
the literature depending on the potential of a scalar field and/or the actual
implementation of the model, here we adopt the one based on an axionlike field
[34]#6#6#6 The cosmic birefringence suggested by the recent analyses of the
Planck data [39, 40] also motivates to introduce such an axionlike field
$\phi$ [41, 42].  whose potential is assumed as

$V(\phi)=\Lambda^{4}\left(1-\cos(\phi/f)\right)^{n}\,,$ ( 2.1)

where $\Lambda$ sets the energy scale of the potential and $f$ is a parameter
which could be related to a breaking scale of the model. To describe the
energy fraction of the EDE field at the redshift when $\phi$ starts to
oscillate, denoted as $z_{c}$, one usually introduces the following quantity:

$f_{\rm de}(z_{c})=\frac{\rho_{\rm de}(z_{c})}{\rho_{\rm total}(z_{c})}\,,$ (
2.2)

where $\rho_{\rm de}$ and $\rho_{\rm total}$ are energy densities of the EDE
field and the total component, respectively. When $n=2$, the energy density of
EDE scales as $\rho_{\rm de}\propto a^{-4}$. When $n=\infty$, $\rho_{\rm de}$
decreases as $\rho_{\rm de}\propto a^{-6}$. We use the camb [43] where
axionlike EDE is already implemented#7#7#7 The fluid approximation as in Eq.
(16) in the paper [34] is used for an axionlike field, which is implemented in
the background and perturbation equations.  and perform the MCMC analysis
sampling $f_{\rm de}(z_{c})\in[0.00001,0.15]$, $z_{c}\in[1000,50000]$, and
$\Theta_{i}\equiv\phi_{\rm ini}/f\in[0.01,3.14]$ (the initial value of $\phi$)
in addition to the 6 standard parameters.

###  2.2 Analysis method

To evaluate the impact of the treatment of BBN on the Hubble tension, we
perform Markov Chain Monte Carlo analysis by using a modified version of
CosmoMC [44], which accommodates the EDE model discussed in [34], by using the
data from Planck (TT,TE,EE+lowE) [45] and BAO [46, 47, 48]. In some cases, we
assume a prior on the Hubble constant as $H_{0}=73.30\pm 1.04$ [1]
(hereinafter, we refer this prior as R21, after the author and the year of the
paper). In addition to these, we include the light element data of the
primordial helium abundance $Y_{P}=0.2449\pm 0.0040$ [23] and two different
data for the primordial deuterium abundance from Cooke et al. $D_{p}=(2.527\pm
0.030)\times 10^{-5}$ [29] and the weighted mean of the recent 11 measurements
compiled by Particle Data Group $D_{p}=(2.547\pm 0.025)\times 10^{-5}$ [30].

For the theoretical calculations of light elements, we consider three
different nuclear rate for $d(p,\gamma)^{3}{\rm He}$, which is important in
computing the deuterium abundance. The one is from Adelberger et al. [31]
which is evaluated from the extrapolation of experimental results. Another
reaction rate is from Marcucci et al. [33] which is an ab-initio theoretical
calculation of the reaction rate, which shows a larger $d(p,\gamma)^{3}{\rm
He}$ rate, more deuterium burning, and less abundance of deuterium $D_{p}$.
The other is taken from Pisanti et al. [32] which is the extrapolation of the
most recent experimental results and their value is somewhat in between the
two rates mentioned above. The value of the theoretical uncertainties are
given by $\sigma_{Y_{p},\rm th}=0.0003$ [49] for the all analysis,
$\sigma_{D_{p},\rm th}=0.06\times 10^{-5}$ [31] for the analysis using
Adelberger et al., $\sigma_{D_{p},\rm th}=0.03\times 10^{-5}$ [33] for the
analysis using Marcucci et al., and $\sigma_{D_{p},\rm th}=0.06\times 10^{-5}$
[50] for the analysis using Pisanti et al. In total, we perform 7 different
analyses depending on the assumption for the data and the theoretical
calculation of $D_{p}$, which is summarized in Table 1. “offset” in Table 1
indicates that we adopt the theoretical bias offset $-0.091\times 10^{-5}$ and
corrected theoretical error $\sigma_{D_{p},\rm th}=0.089\times 10^{-5}$. These
values are introduced in the standard CosmoMC analysis so that the BBN
constraints correspond to the prior $\Omega_{b}h^{2}=0.0222\pm 0.005$#8#8#8
See Sec. 7.6.1 of the Planck paper [2] for details. . We note that Case II or
III is commonly adopted one in the literature when the BBN data is included in
the CMB analysis.

| Reaction rate for $d(p,\gamma)^{3}{\rm He}$ |  $D_{p}$ data
---|---|---
Case I | Adelberger et al. | no-data
Case II | Adelberger et al. | Cooke et al. + offset
Case III | Adelberger et al. | Cooke et al.
Case IV | Marcucci et al. | Cooke et al.
Case V | Marcucci et al. | recent 11 weighted mean
Case VI | Pisanti et al. | Cooke et al.
Case VII | Pisanti et al. | recent 11 weighted mean
Table 1:  Combinations of the data of $D_{p}$ and theoretical treatment of the
reaction rate for $d(p,\gamma)^{3}{\rm He}$ used in our analysis.

##  3 Results

Now we present the results of our analysis. In Tables 2, 3 and 4, the derived
values of $H_{0}$ and $\Omega_{b}h^{2}$, as well as the Gaussian tension are
summarized for the framework models of $\Lambda$CDM, EDE ($n=2$) and EDE
($n=\infty$), respectively. The Gaussian tension is evaluated between the
direct measurement of $H_{0}$ from Cepheid calibrated supernova distance
ladder [1] $H_{0}=73.30\pm 1.04~{}{\rm km/s/Mpc}$ and the one obtained in our
analysis using CMB+BAO (with/without BBN and the $H_{0}$ prior). For the cases
of EDE, we also give the value of $f_{{\rm de}}(z_{c})$. We have performed 10
analyses including Case I - VII in Table 1 and the cases with and without the
prior on $H_{0}$ (R21) for some cases.

For the case of $\Lambda$CDM, depending on the treatment of the observational
data and theoretical calculation for $D_{p}$, the significance of the tension
varies from $4.97\sigma$ to $5.38\sigma$ when R21 is not included. As already
mentioned, the setup of Case II or III is commonly adopted when one includes
the BBN data in combination with CMB. It is also worth noting that in Case III
the Hubble tension is slightly eased with respect to Case I (no-BBN data
analysis) and Case II (standard BBN analysis including the offset). In any
case, just by changing the assumption on the observational data for $D_{p}$
and the reaction rate, the significance is increased by $0.41\sigma$. The
difference mainly comes from the treatment of the reaction rate. In Figure 1,
we show the posterior distributions and two dimensional constraints for
$\Omega_{b}h^{2},H_{0}$ and $D_{p}$, which shows that the deuterium abundance
$D_{p}$ decrease in the order of Case III (Adelberger et al.), Case VI
(Pisanti et al.), and Case IV (Marcucci et al.) for a given baryon fraction
$\Omega_{b}h^{2}$ (all these cases use Cooke et al. [29] for the observational
data of $D_{p}$). In order to be consistent with $D_{p}$ observations, it is
important to notice that the decrease of $D_{p}$ leads to a lower baryon
fraction, a larger sound horizon $r_{*}$, and a lower Hubble constant $H_{0}$,
which makes the tension worse. The observational data also affects slightly
although two data from Cooke et al. [29] and the averaged mean value from 11
recent measurements of $D_{p}$ [30] are almost consistent. For example, by
comparing Case IV and V in which Cooke et al. and the recent 11 weighted mean
value are adopted for the $D_{p}$ data with the same reaction rate, the
tension in Case V is slightly larger than that in Case IV.

| $\Omega_{b}h^{2}$ | $H_{0}$ | Tension
---|---|---|---
Case I | $0.02242\pm 0.00014$ | $67.67\pm 0.43$ | 5.00$\sigma$
Case I +R21 | $0.02259\pm 0.00013$ | $68.49\pm 0.40$ | 4.32$\sigma$
Case II | $0.02240\pm 0.00013$ | $67.63\pm 0.43$ | 5.04$\sigma$
Case III | $0.02246\pm 0.00013$ | $67.73\pm 0.42$ | 4.97$\sigma$
Case III+R21 | $0.02260\pm 0.00012$ | $68.51\pm 0.39$ | 4.31$\sigma$
Case IV | $0.02230\pm 0.00012$ | $67.42\pm 0.40$ | 5.28$\sigma$
Case V | $0.02226\pm 0.00011$ | $67.33\pm 0.39$ | 5.38$\sigma$
Case V +R21 | $0.02238\pm 0.00011$ | $68.09\pm 0.38$ | 4.71$\sigma$
Case VI | $0.02242\pm 0.00013$ | $67.66\pm 0.41$ | 5.05$\sigma$
Case VII | $0.02240\pm 0.00013$ | $67.62\pm 0.42$ | 5.06$\sigma$
Table 2:  Derived values of $\Omega_{b}h^{2}$, $H_{0}$ and the Gaussian
tension for $\Lambda$CDM model.  Figure 1:  Posterior distributions and two
dimensional constraints for $\Omega_{b}h^{2},H_{0}$ and $D_{p}$ in the
$\Lambda$CDM model.

The impact of the treatment of BBN is more significant when we study the issue
in the EDE model which has been proposed to relax the $H_{0}$ tension. In
Tables 3 and 4, derived values of $\Omega_{b}h^{2},f_{{\rm de}}(z_{c})$,
$H_{0}$ and the Gaussian tension are listed for the cases of $n=2$ and
$n=\infty$. In Figs. 2 and 3, the posterior distribution and two dimensional
constraints for $\Omega_{b}h^{2},f_{\rm DE}(z_{c})$ and $H_{0}$ are depicted
for the cases of $n=2$ and $n=\infty$, respectively. We first note that, in
the EDE model we adopted here, the significance of the tension is $3.6\sigma$
for Case I where we do not include the BBN data, which can be compared to
$5.0\sigma$ in the counterpart of the $\Lambda$CDM. However, when the BBN is
included in the analysis for the case of $n=2$, it increases to $3.83\sigma$
for the standard treatment of Case II and even worse to $4.41\sigma$ for Case
V. It should be noticed that the significance of the tension varies from
$3.61\sigma$ to $4.41\sigma$ ($0.8\sigma$ difference) just by changing the
assumption on the treatment of $D_{p}$. The tendency of which reaction rate
makes the tension larger is the same as that in the $\Lambda$CDM model,
however, the change of the tension is more pronounced in the EDE model. For
the case of $n=\infty$, the Gaussian tension is increased from $4.58\sigma$ in
Case II to $5.23\sigma$ in Case V, which is $0.65\sigma$ increase. Although
the change is less significant compared to the case of $n=2$, BBN still have
an impact on the $H_{0}$ tension.

| $\Omega_{b}h^{2}$ | $f_{{\rm de}}(z_{c})$ | $H_{0}$ | Tension
---|---|---|---|---
Case I | $0.02252\pm 0.00016$ | $0.0106_{-0.0088}^{+0.0043}$ | $68.64_{-0.86}^{+0.65}$ | $3.61\sigma$
Case I +R21 | $0.02277\pm 0.00016$ | $0.0277_{-0.0092}^{+0.0071}$ | $70.64\pm 0.73$ | $2.09\sigma$
Case II | $0.02245\pm 0.00015$ | $0.0098_{-0.0082}^{+0.0039}$ | $68.49_{-0.79}^{+0.60}$ | $3.83\sigma$
Case III | $0.02254\pm 0.00014$ | $0.0111_{-0.010}^{+0.0036}$ | $68.68_{-0.86}^{+0.61}$ | $3.61\sigma$
Case III+R21 | $0.02276_{-0.00014}^{+0.00016}$ | $0.0261_{-0.0081}^{+0.0073}$ | $70.54\pm 0.70$ | $2.20\sigma$
Case IV | $0.02235\pm 0.00012$ | $0.0085_{-0.0078}^{+0.0028}$ | $68.13_{-0.68}^{+0.57}$ | $4.26\sigma$
Case V | $0.02229\pm 0.00012$ | $0.0081_{-0.0079}^{+0.0021}$ | $67.96_{-0.72}^{+0.50}$ | $4.41\sigma$
Case V +R21 | $0.02238\pm 0.00015$ | $0.033_{-0.015}^{+0.012}$ | $70.22_{-0.88}^{+0.76}$ | $2.32\sigma$
Case VI | $0.02250\pm 0.00014$ | $0.0099_{-0.0089}^{+0.0034}$ | $68.53_{-0.81}^{+0.58}$ | $3.80\sigma$
Case VII | $0.02247\pm 0.00014$ | $0.0096_{-0.0086}^{+0.0032}$ | $68.46_{-0.80}^{+0.58}$ | $3.86\sigma$
Table 3:  Derived values of $\Omega_{b}h^{2},f_{{\rm de}}(z_{c})$, $H_{0}$ and the Gaussian tension for the EDE model with $n=2$.  | $\Omega_{b}h^{2}$ | $f_{{\rm de}}(z_{c})$ | $H_{0}$ | Tension
---|---|---|---|---
Case I | $0.02249_{-0.00017}^{+0.00014}$ | $0.0174_{-0.017}^{+0.0037}$ | $67.93_{-0.62}^{+0.42}$ | $4.60\sigma$
Case I +R21 | $0.02302\pm 0.00021$ | $0.079_{-0.019}^{+0.021}$ | $70.85\pm 0.77$ | $1.89\sigma$
Case II | $0.02245_{-0.00015}^{+0.00013}$ | $0.0165_{-0.016}^{+0.0034}$ | $67.82_{-0.54}^{+0.42}$ | $4.78\sigma$
Case III | $0.02251_{-0.00016}^{+0.00013}$ | $0.0177_{-0.018}^{+0.0035}$ | $67.99_{-0.62}^{+0.38}$ | $4.58\sigma$
Case III+R21 | $0.02293\pm 0.00019$ | $0.075_{-0.018}^{+0.022}$ | $70.65\pm 0.76$ | $2.06\sigma$
Case IV | $0.02232\pm 0.00012$ | $0.0131_{-0.013}^{+0.00226}$ | $67.52\pm 0.42$ | $5.15\sigma$
Case V | $0.02227\pm 0.00011$ | $0.0135_{-0.014}^{+0.0026}$ | $67.44_{-0.45}^{+0.38}$ | $5.23\sigma$
Case V +R21 | $0.02243\pm 0.00015$ | $0.044_{-0.020}^{+0.027}$ | $69.31_{-0.80}^{+0.92}$ | $2.95\sigma$
Case VI | $0.02246_{-0.00015}^{+0.00013}$ | $0.0174_{-0.017}^{+0.0038}$ | $67.87_{-0.57}^{+0.40}$ | $4.72\sigma$
Case VII | $0.02243\pm 0.00014$ | $0.0149_{-0.015}^{+0.0029}$ | $67.76_{-0.48}^{+0.41}$ | $4.90\sigma$
Table 4:  Derived values of $\Omega_{b}h^{2},f_{{\rm de}}(z_{c})$, $H_{0}$ and
the Gaussian tension for the DE model with $n=\infty$.  Figure 2:  Posterior
distributions and two dimensional constraints for $\Omega_{b}h^{2},f_{\rm
DE}(z_{c})$ and $H_{0}$ for the EDE model with $n=2$. Figure 3:  Posterior
distributions and two dimensional constraints for $\Omega_{b}h^{2},f_{\rm
DE}(z_{c})$ and $H_{0}$ for the EDE model with $n=\infty$.

##  4 Conclusion and discussion

We have investigated the impact of the BBN on the Hubble tension, focusing on
to what extent the evaluation of the tension is affected by how we include the
BBN in the analysis. In models proposed to resolve the $H_{0}$ tension, not
only the derived value of $H_{0}$ tends to be increased in the fitting to
cosmological data such as CMB, the determination of other cosmological
parameters is also affected, among which the baryon density is significantly
affected. Indeed the baryon density can also be well determined by the BBN
data, in particular the primordial deuterium abundance. We have investigated
this issue in the $\Lambda$CDM model, and the EDE as an example of models
which has been proposed to resolve the $H_{0}$ tension.

When one includes the BBN data in addition to CMB and BAO in the analysis,
there are several ways to take account of it as to which observational data is
used and which nuclear reaction rate, especially for $d(p,\gamma)^{3}{\rm He}$
which is particularly relevant to the deuterium abundance, is adopted. We have
made the analysis adoping two different observational data for $D_{p}$ and
three separate treatments for the reaction rate of $d(p,\gamma)^{3}{\rm He}$.
We found that the significance of the $H_{0}$ tension between the direct
measurement of Cepheid calibrated supernova distance ladder from [1] and
indirect ones from Planck+BAO in combination with BBN data varies depending on
which observational data and the nuclear reaction rate are used in the
analysis.

In $\Lambda$CDM model, the tension is $5.04\sigma$ for Case II where the value
of $D_{p}$ from Cooke et al. [29] and the reaction rate of Adelberger et al.
[31] for $d(p,\gamma)^{3}{\rm He}$ are assumed, which is usually adopted in
the literature. However, it increases to $5.38\sigma$ for Case V where the
weighted average of 11 recent measurements of $D_{p}$ [30] and the reaction
rate from Marcucci et al. [33] are assumed.

The impact is more significant in the EDE model, which has been extensively
studied as a possible solution to the $H_{0}$ tension. For the case of the EDE
model with $n=2$, the significance is increased from $3.83\sigma$ (Case II) to
$4.41\sigma$ (Case V), which indicates that the assumption in the treatment of
BBN would significantly affect the tension. Our analysis shows that we need to
carefully consider which data and the reaction rate are adopted in the BBN
analysis to investigate the Hubble tension.

Since the Hubble tension is now one of the most serious problem in cosmology
today, one needs to carefully quantify the tension particularly when we
analyze models to resolve the $H_{0}$ tension. We have various accurate
cosmological observations such as CMB, BAO and BBN and the combinations of
these data would constrain cosmological models even further. However, we may
have some options as to which data and how the theoretical calculation are
performed, which can appreciably affect the significance of the $H_{0}$
tension. We have explicitly demonstrated this issue for the BBN, which gives
some cautions to future research in the $H_{0}$ tension.

## Acknowledgements

This work was supported in part by JSPS KAKENHI Grant Number 19K03874 (TT),
MEXT KAKENHI Grant Number 23H04515 (TT), and JST SPRING, Grant Number
JPMJSP2119 (YT).

## References

  * [1] A. G. Riess et al., A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km s-1 Mpc-1 Uncertainty from the Hubble Space Telescope and the SH0ES Team, Astrophys. J. Lett. 934 (2022), no. 1 L7, [arXiv:2112.04510].
  * [2] Planck Collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6, [arXiv:1807.06209]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [3] N. Schöneberg, J. Lesgourgues, and D. C. Hooper, The BAO+BBN take on the Hubble tension, JCAP 10 (2019) 029, [arXiv:1907.11594].
  * [4] F. Okamatsu, T. Sekiguchi, and T. Takahashi, $H_{0}$ tension without CMB data: Beyond the $\Lambda$CDM, Phys. Rev. D 104 (2021), no. 2 023523, [arXiv:2105.12312].
  * [5] N. Schöneberg, L. Verde, H. Gil-Marín, and S. Brieden, BAO+BBN revisited – Growing the Hubble tension with a 0.7km/s/Mpc constraint, arXiv:2209.14330.
  * [6] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, In the realm of the Hubble tension—a review of solutions, Class. Quant. Grav. 38 (2021), no. 15 153001, [arXiv:2103.01183].
  * [7] L. Perivolaropoulos and F. Skara, Challenges for $\Lambda$CDM: An update, New Astron. Rev. 95 (2022) 101659, [arXiv:2105.05208].
  * [8] N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, The H0 Olympics: A fair ranking of proposed models, Phys. Rept. 984 (2022) 1–55, [arXiv:2107.10291].
  * [9] K. Ichikawa, M. Fukugita, and M. Kawasaki, Constraining neutrino masses by CMB experiments alone, Phys. Rev. D 71 (2005) 043001, [astro-ph/0409768].
  * [10] K. Ichikawa, T. Sekiguchi, and T. Takahashi, Probing the Effective Number of Neutrino Species with Cosmic Microwave Background, Phys. Rev. D 78 (2008) 083526, [arXiv:0803.0889].
  * [11] T. Sekiguchi, K. Ichikawa, T. Takahashi, and L. Greenhill, Neutrino mass from cosmology: Impact of high-accuracy measurement of the Hubble constant, JCAP 03 (2010) 015, [arXiv:0911.0976].
  * [12] T. Takahashi and S. Yamashita, Big bang nucleosynthesis and early dark energy in light of the EMPRESS Yp results and the H0 tension, Phys. Rev. D 107 (2023), no. 10 103520, [arXiv:2211.04087].
  * [13] O. Seto and Y. Toda, Comparing early dark energy and extra radiation solutions to the Hubble tension with BBN, Phys. Rev. D 103 (2021), no. 12 123501, [arXiv:2101.03740].
  * [14] O. Seto and Y. Toda, Big bang nucleosynthesis constraints on varying electron mass solution to the Hubble tension, Phys. Rev. D 107 (2023), no. 8 083512, [arXiv:2206.13209].
  * [15] O. Seto and Y. Toda, Hubble tension in lepton asymmetric cosmology with an extra radiation, Phys. Rev. D 104 (2021), no. 6 063019, [arXiv:2104.04381].
  * [16] T. Sekiguchi and T. Takahashi, Cosmological bound on neutrino masses in the light of $H_{0}$ tension, Phys. Rev. D 103 (2021), no. 8 083516, [arXiv:2011.14481].
  * [17] G. Ye, B. Hu, and Y.-S. Piao, Implication of the Hubble tension for the primordial Universe in light of recent cosmological data, Phys. Rev. D 104 (2021), no. 6 063510, [arXiv:2103.09729].
  * [18] F. Takahashi and W. Yin, Cosmological implications of $n_{s}\approx 1$ in light of the Hubble tension, Phys. Lett. B 830 (2022) 137143, [arXiv:2112.06710].
  * [19] J.-Q. Jiang and Y.-S. Piao, Toward early dark energy and ns=1 with Planck, ACT, and SPT observations, Phys. Rev. D 105 (2022), no. 10 103514, [arXiv:2202.13379].
  * [20] G. Ye, J.-Q. Jiang, and Y.-S. Piao, Towards hybrid inflation with $n_{s}=1$ in light of Hubble tension and primordial gravitational waves, arXiv:2205.02478.
  * [21] T. Sekiguchi and T. Takahashi, Early recombination as a solution to the $H_{0}$ tension, Phys. Rev. D 103 (2021), no. 8 083507, [arXiv:2007.03381].
  * [22] V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Early Dark Energy Can Resolve The Hubble Tension, Phys. Rev. Lett. 122 (2019), no. 22 221301, [arXiv:1811.04083].
  * [23] E. Aver, K. A. Olive, and E. D. Skillman, The effects of He I $\lambda$10830 on helium abundance determinations, JCAP 07 (2015) 011, [arXiv:1503.08146].
  * [24] T. Hsyu, R. J. Cooke, J. X. Prochaska, and M. Bolte, The PHLEK Survey: A New Determination of the Primordial Helium Abundance, Astrophys. J. 896 (2020), no. 1 77, [arXiv:2005.12290].
  * [25] O. A. Kurichin, P. A. Kislitsyn, V. V. Klimenko, S. A. Balashev, and A. V. Ivanchik, A new determination of the primordial helium abundance using the analyses of H II region spectra from SDSS, Mon. Not. Roy. Astron. Soc. 502 (2021), no. 2 3045–3056, [arXiv:2101.09127].
  * [26] A. Matsumoto et al., EMPRESS. VIII. A New Determination of Primordial He Abundance with Extremely Metal-poor Galaxies: A Suggestion of the Lepton Asymmetry and Implications for the Hubble Tension, Astrophys. J. 941 (2022), no. 2 167, [arXiv:2203.09617].
  * [27] K. Ichikawa and T. Takahashi, Revisiting the constraint on the helium abundance from cmb, Phys. Rev. D 73 (2006) 063528, [astro-ph/0601099].
  * [28] K. Ichikawa, T. Sekiguchi, and T. Takahashi, Primordial Helium Abundance from CMB: a constraint from recent observations and a forecast, Phys. Rev. D 78 (2008) 043509, [arXiv:0712.4327].
  * [29] R. J. Cooke, M. Pettini, and C. C. Steidel, One Percent Determination of the Primordial Deuterium Abundance, Astrophys. J. 855 (2018), no. 2 102, [arXiv:1710.11129].
  * [30] Particle Data Group Collaboration, R. L. Workman et al., Review of Particle Physics, PTEP 2022 (2022) 083C01.
  * [31] E. G. Adelberger et al., Solar fusion cross sections II: the pp chain and CNO cycles, Rev. Mod. Phys. 83 (2011) 195, [arXiv:1004.2318].
  * [32] O. Pisanti, G. Mangano, G. Miele, and P. Mazzella, Primordial Deuterium after LUNA: concordances and error budget, JCAP 04 (2021) 020, [arXiv:2011.11537].
  * [33] L. E. Marcucci, G. Mangano, A. Kievsky, and M. Viviani, Implication of the proton-deuteron radiative capture for Big Bang Nucleosynthesis, Phys. Rev. Lett. 116 (2016), no. 10 102501, [arXiv:1510.07877]. [Erratum: Phys.Rev.Lett. 117, 049901 (2016)].
  * [34] V. Poulin, T. L. Smith, D. Grin, T. Karwal, and M. Kamionkowski, Cosmological implications of ultralight axionlike fields, Phys. Rev. D 98 (2018), no. 8 083525, [arXiv:1806.10608].
  * [35] M. Braglia, W. T. Emond, F. Finelli, A. E. Gumrukcuoglu, and K. Koyama, Unified framework for early dark energy from $\alpha$-attractors, Phys. Rev. D 102 (2020), no. 8 083513, [arXiv:2005.14053].
  * [36] P. Agrawal, F.-Y. Cyr-Racine, D. Pinner, and L. Randall, Rock ’n’ Roll Solutions to the Hubble Tension, arXiv:1904.01016.
  * [37] M. Kamionkowski and A. G. Riess, The Hubble Tension and Early Dark Energy, arXiv:2211.04492.
  * [38] V. Poulin, T. L. Smith, and T. Karwal, The Ups and Downs of Early Dark Energy solutions to the Hubble tension: a review of models, hints and constraints circa 2023, arXiv:2302.09032.
  * [39] Y. Minami and E. Komatsu, New Extraction of the Cosmic Birefringence from the Planck 2018 Polarization Data, Phys. Rev. Lett. 125 (2020), no. 22 221301, [arXiv:2011.11254].
  * [40] P. Diego-Palazuelos et al., Cosmic Birefringence from the Planck Data Release 4, Phys. Rev. Lett. 128 (2022), no. 9 091302, [arXiv:2201.07682].
  * [41] L. M. Capparelli, R. R. Caldwell, and A. Melchiorri, Cosmic birefringence test of the Hubble tension, Phys. Rev. D 101 (2020), no. 12 123529, [arXiv:1909.04621].
  * [42] K. Murai, F. Naokawa, T. Namikawa, and E. Komatsu, Isotropic cosmic birefringence from early dark energy, Phys. Rev. D 107 (2023), no. 4 L041302, [arXiv:2209.07804].
  * [43] http://camb.info.
  * [44] A. Lewis and S. Bridle, Cosmological parameters from CMB and other data: A Monte Carlo approach, Phys. Rev. D 66 (2002) 103511, [astro-ph/0205436].
  * [45] Planck Collaboration, N. Aghanim et al., Planck 2018 results. V. CMB power spectra and likelihoods, Astron. Astrophys. 641 (2020) A5, [arXiv:1907.12875].
  * [46] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, The 6dF Galaxy Survey: Baryon Acoustic Oscillations and the Local Hubble Constant, Mon. Not. Roy. Astron. Soc. 416 (2011) 3017–3032, [arXiv:1106.3366].
  * [47] A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, The clustering of the SDSS DR7 main Galaxy sample – I. A 4 per cent distance measure at $z=0.15$, Mon. Not. Roy. Astron. Soc. 449 (2015), no. 1 835–847, [arXiv:1409.3242].
  * [48] BOSS Collaboration, S. Alam et al., The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample, Mon. Not. Roy. Astron. Soc. 470 (2017), no. 3 2617–2652, [arXiv:1607.03155].
  * [49] R. Consiglio, P. F. de Salas, G. Mangano, G. Miele, S. Pastor, and O. Pisanti, PArthENoPE reloaded, Comput. Phys. Commun. 233 (2018) 237–242, [arXiv:1712.04378].
  * [50] S. Gariazzo, P. F. de Salas, O. Pisanti, and R. Consiglio, PArthENoPE revolutions, Comput. Phys. Commun. 271 (2022) 108205, [arXiv:2103.05027].


# Modeling meteorite craters by impacting melted tin on sand

H. Y. Huang1, P. R. Tsai1, C. Y. Lu1, H. Hau1, Y. L. Chen2, Z. T. Ling3, Y. R.
Wu4 and Tzay-Ming Hong1∗ 1Department of Physics, National Tsing Hua
University, Hsinchu, Taiwan 30013, Republic of China  2Department of
Mechanical Engineering, National Taiwan University, Taipei, Taiwan 10617,
Republic of China  3Department of New Media Art, Taipei National University of
the Arts, Taipei, Taiwan 110024, Republic of China  4National Hsinchu Senior
High School, Hsinchu, Taiwan 30013, Republic of China

###### Abstract

To simulate the heated exterior of a meteorite, we impact a granular bed by
melted tin. The morphology of tin remnant and crater is found to be sensitive
to the temperature and solidification of tin. By employing deep learning and
convolutional neural network, we can quantify and map the complex impact
patterns onto network systems based on feature maps and Grad-CAM results. This
gives us unprecedented details on how the projectile deforms and interacts
with the granules, which information can be used to trace the development of
different remnant shapes. Furthermore, full dynamics of granular system is
revealed by the use of Particle Image Velocimetry. Kinetic energy, temperature
and diameter of the projectile are used to build phase diagrams for the
morphology of both crater and tin remnant. In addition to successfully
reproducing key features of simple and complex craters, we are able to detect
a possible artifact when compiling crater data from field studies. The depth
of craters from high-energy impacts in our work is found to be independent of
their width. However, when mixing data from different energy, temperature and
diameter of projectile, a bogus power-law relationship appears between them.
Like other controlled laboratory researches, our conclusions have the
potential to benefit the study of paint in industry and asteroid sampling
missions on the surface of celestial bodies.

††preprint: APS/123-QED

##  I Introduction

Meteorite impact events have appeared in many disaster movies. However,
Chelyabinsk meteor impact event and 2019 OK asteroid sweeping past the Earth
remind us that meteorite impacts remain a menace to life on earth. Since the
days of Galileo, geologists and astronomers have accumulated a lot of
knowledge on the meteorite craters via telescope and field studies. This
information can help us not only prepare for and mitigate potential risks
posed by future impacts, but also understand the history and geology of
habitable exoplanets. Besides, impact craters provide valuable insights into
the formation and location of recoverable mineral resources both on Earth and
beyond.

Starting roughly 30 years ago, scientists began to smash various projectiles
on a granular bed to simulate the impact event [3, 4, 5]. This endows precious
information, esp. on how the impact energy $E$, proportional to drop height
$h$, affects the width $w$ and depth $d$ of crater. Roughly proportional to
$d$, $w$ scales as $E^{\alpha}$ with $\alpha=0.25$ [6] or 0.17 [7] depending
on whether a hard (e.g., steel ball) or deformable (water droplet) bullet is
used. Meanwhile, it was found that some specific features, e.g., the central
peak in complex craters, are better reproduced by adopting a granular
projectile [8]. As opposed to field studies, knowledge accumulated from
controlled laboratory research like these also benefits the study of soil
erosion in agriculture [9, 10], paint in industry [11], and asteroid sampling
missions on the surface of celestial bodies [12, 13]. It is thus important to
continue lab experiments of impacting granular beds.

The important role of heat transfer and solidification for melted tin
impacting on a solid surface has been discussed [14] to gain insight on
thermal spraying and additive manufacturing, and extreme ultraviolet
lithography in chip production. In this Letter we adopt melted tin droplets as
the bullet to better simulate the high-temperature environment typical of
meteorite impact events. Not limited to $h$ in Fig. 1, the tin remnant also
changes its shape with temperature $T$. This information will be summarised in
diagrams and shown to be correlated to the morphology of crater and tin
remnant. As in previous studies [15], we also investigate $w(E)$ and $d(E)$ to
see whether the solidification of melted tin will create any new feature.
Particle image velocimetry (PIV) and image processing are employed to analyze
the flow field of granules, while image processing enables us to project the
tin shape on a 2-D plane and witness its deformation in real time.
Additionally, we use deep learning - convolutional neural network (CNN) and
reverse backpropagation Grad-CAM to detect features of tin during the impact
process and discuss their network properties.

Figure 1:  (a$\sim$d) Snapshots of melted tin of $R=6.4$ mm at 380 ℃ impacting
sand for $h$ = 100, 250, 400, and 700 mm. Scale bars: 6.0 mm. False colors are
used as an aid in visualizing features. The shape of tin remnant, highlighted
by dark solid lines, is categorized into sphere, disk, shotgun, and snowflake,
as labelled by red circle, black rhombus, green square, and blue triangle.
Figure 2: (a) Schematic of experimental setup. (b) Phase diagram with $h$ and
$T$ as axes shows the morphology of impact crater created by melted tin of $R$
= 6.4 mm can mimic either letter U (in blue) or V (brown). Green area denotes
an intermediate phase, and labels for the remnant shape follow those in Fig.
1. (c) Similar phase diagram with $h$ and $R$ at 380 ℃. (d) Schematic profile
of U- and V-shape craters. (e) Shape ratio $r_{s}$ vs $h$ at 300 and 420 °C,
represented by red square and black rhombus.

##  II Experimental Setup

As shown in Fig. 2(a), our setup is typical of impact experiments [16], except
that it comprises of an elevated station with two soldering irons controlled
at $T$ = 300, 340, 380, 420, and 460 °C, respectively. The spacing and
relative orientation between irons can be adjusted to produce tin droplets of
mass 0.30, 0.32, 0.34, 0.36, and 0.38 g or equivalently 6.4, 6.5, 6.6, 6.8,
and 6.9 mm in diameter $R$. For the height $h$ of release from 0.05 to 1 m,
the decrease in temperature during free fall is negligible. To eliminate the
effect of granular compaction, we follow the practice of Walsh by reloading
the sand and scraping off extra sands by a ruler to flatten its surface before
each trial [3]. The impacts are characterized by the inverse Frode number
${\rm Fr}^{-1}=gD_{b}/2v^{2}$, the ratio of gravitational to dynamic pressures
[3]. Here $g$ denotes the gravitational acceleration, $D_{b}$ the projectile
diameter, and $v$ the impact velocity. In our system, ${{10}^{-3}\leq{\rm
Fr}^{-1}}\leq{10}^{-2}$ overlaps with ${{10}^{-6}\leq{\rm
Fr}^{-1}}\leq{10}^{-2}$ for meteorite craters, which justifies the analogy to
astronomical impact events [5].

##  III Experimental Results

Figure 3: (a) The top-view area $A$ of tin remnant versus time $t$ for sphere,
sphere-II, disk, shotgun, and snowflake shapes is denoted by red circle,
orange triangle, black rhombus, green square, and blue triangle, respectively.
(b) The total kinetic energy $K$ of sand for different tin shapes is plotted
as a function of $t$. (c) PIV results for snapshots show the evolution of
melted tin with $R=6.4$ mm, 380 ℃ and $h$ = 50 mm where velocity colorbars for
sand on the left and tin on the right. Image of tin remnant after texture
segmentation [22] is inserted on the upper right corner of each photo.

The morphology of melted tin in Fig. 1 is recorded by a high-speed camera of
8146 fps. Irrespective of $h$, the process can be roughly separated into four
stages based on the configuration of tin droplet: (1) circumference increases
with time, but remains approximately circular, (2) irregular shape emerges,
(3) its boundary retracts, and (4) static configuration is reached. As the
impact energy increases, the final shape of tin remnant can roughly be
categorized into four shapes: sphere, disk, shotgun, and snowflake, as
illustrated in Fig. 1.

The sequence of appearance for these shapes can be understood physically.
Since the distortion of tin melt raises its surface energy, we expect the
aspect ratio of disk to increase with impact energy. As the thickness of disk
reaches the threshold of Rayleigh-Plateau (RP) instability [17], clumping
ensues and produces many small and smooth tin balls in the stage-4 photo of
Fig. 1(c), which mimics the shotgun pattern. As more kinetic energy is fueled,
fast spreading of melted tin hastens the conduction rate of its heat to the
sand bed. As a result, the transition from liquid to solid phase precedes and
prohibits the RP instability. This creates a splash-like [18] irregular
configuration in one piece, which we term snowflake in Fig. 1(d).

###  III.1 Phase Diagrams for Morphology

####  III.1.1 Tin Remnants

Figure 2(b, c) show the phase diagrams for the morphology of both the crater
and tin remnant as a function of parameters ($T$, $h$) and ($R$, $h$). At low
impact energies or for $h$ = 50$\sim$250 mm, a higher $T$ renders the melted
tin more susceptible to flattening upon impact. At slightly higher energies or
$h$ = 400$\sim$500 mm, $h_{c}$ roughly marks the emergence of shotgun pattern
at high temperatures, again because the fluid behavior is less restricted
which allows the RP instability to enter the picture. Within the same scenario
that high temperature promotes fluidity, the transition from shotgun to
snowflake around $h_{s}$ is analogous to the effect of increasing $h$ at high
energies. We do not understand why two shapes can be observed in some
intervals of $R$ in Fig. 2(c), as marked by crosses, daggers, and stars.

Figure 4: (a) Grad-CAM results for different shapes of tin remnant. Network
density $D$ vs $L$ are shown in (b$\sim$e) for stage 1$\sim$4 respectively
where the symbols are defined in (a).

The trend of phase boundary in Fig. 2(c) for remnant shape at low energies
mimics that in Fig. 2(b). It is because the dissipation rate per unit mass is
proportional to the surface-area-to-volume ratio which for a sphere drops as
$R$ increases. A larger melted tin thus gets to retain its malleability longer
and better. This effect of preferring disk over sphere at low impact energy is
similar to raising $T$. Ostensibly, the trend of $h_{c}$ and $h_{s}$ in Fig.
2(b) are both reversed in Fig. 2(c). Their physical explanations are in fact
similar and consistent. The reduction in thickness as we decrease $R$ favors
the occurrence of RP instability at medium $h$, while making it easier to lose
heat and solidify before RP occurs at high $h$ \- both effects are similar to
those by increasing $T$.

For complex meteorite craters, due to extremely large mass, more energy will
transform into heat during the impact process [19]. As a result, the material
on the surface of celestial bodies will change to liquid or gas phases. It was
observed that the edge of these large craters is so steep that it often
collapses under its own weight and concentrates at the bottom of craters [20].
This is one of several mechanisms of the central peak in the majority of
complex craters [21]. Note that the spherical shape reappears at higher-$h$ in
Fig. 2(b, c), which is associated with a steeper impact crater in stage 2 and
will be denoted by sphere-II. The appearance of this shape is always
accompanied by a quasi-bump forms in the middle of crater as the sand on the
wall and rim of the transient cavity collapses and collides at the center in
stage 3. This mimics the central-peak structure in complex craters.

####  III.1.2 Crater

We first try to reproduce the general observation that the complex craters
tend to be more plate-like with $d/w=0.05\sim 0.1$, as opposed to bowl-shaped
with $d/w=0.14\sim 0.2$ for simple ones [21]. With laser displacement sensors
and a homemade translation stage, we can obtain detailed information on the
profile of our impact crater. In Fig. 2(b, c), two critical heights, $h_{c}$
and $h_{s}$, are defined to manifest the shape change of tin remnant.
Physically we expect the transition to leave its mark on the morphology of
crater and, thus, introduce the shape ratio $r_{s}\equiv w_{r}/w$ to quantify
the profile where $w_{r}$ is the distance between inflection points as shown
in Fig. 2(d). According to Fig. 2(e), it turns out that $h<h_{c}$ indeed
corresponds to a smaller $r_{s}$, i.e., resembling V-shape or a bowl. On the
other hand, U-shape or plate-like craters are more likely to occur at
$h>h_{s}$. This is consistent with the distinctive shapes for simple and
complex craters. We speculate that this dichotomy in shape has to do with the
viscoelastic property of sand and soil [4] \- a high-energy collision triggers
more elastic response that is more effective at halting the forward motion of
tin and forcing it to spread out horizontally.

Figure 5: (a) Crater width $w$ vs $h$ for steel balls (in black circle) and
melted tin with 300, 340, 380, 420 and 460 ℃ (red star, orange triangle,
yellow square, green rhombus and cyan triangle). (b) Crater depth $d$ vs $h$
with the same symbols as (a). (c) $w$ vs $h$ for steel balls (black circle)
and melted tin with 6.36, 6.5, 6.63, 6.76 and 6.88 mm (red star, orange
triangle, yellow square, green rhombus and cyan triangle). (d) $d$ vs $h$ with
the same symbols as (c). The insets in (a$\sim$d) show the same data in full-
log plots. (e) $d$ vs $w$ for V- (red square) and U-shape craters (black
rhombus).

###  III.2 Deep-learning-augmented Analysis

####  III.2.1 Expanded Area Projection and PIV

In order to confirm our speculations in Sec. III.A.1 regarding the different
shapes of tin remnant, we employ image processing techniques and a high-speed
camera to measure the top-view area $A$ of tin at successive times in Fig.
3(a). Note that the exhibition of a peak value for all lines signals the time
beyond which the melted tin starts to retract. Two things worth mentioning.
First, the maximum $A$ value of sphere-II is smaller than that of sphere. This
is against our intuition because the former is dropped from a higher $h$. It
turns out that sphere-II causes a steeper cavity in stage 2 which hinders the
spread of tin. Second, there exists a plateau immediately following the peak
for shotgun and snowflake. This is ascribed to the faster solidification due
to their lesser thickness when spread out to a larger area. As for the sudden
drop following the plateau for shotgun, it is due to the breakup of melted tin
into several small spheres.

Using PIV analysis and image processing, we can measure the speed of moving
granules and calculate their total kinetic energy $K$, as shown in Fig. 3(b).
Physically $K$ has to compete with the surface energy of melted tin for the
energy converted from the gravitational potential during the impact. The
reason why shotgun allocates lesser $K$ for sand than other shapes from a
lower height is that the surface energy already gets the lion’s share. Aside
from the leftmost image that illustrates the initial impact, the successive
images in Fig. 3(c) correspond to the four stages.

####  III.2.2 Network of Tin Remnants

Creating the equivalent network for the tin remnant allows us to distinguish
the highly complex pattern of tin remnants quantitatively and more
objectively. This is made possible by employing CNN and Grad-CAM [23]. The
structure and learning curve of CNN can be found in Ref.[22]. We adopt the
features, used by Grad-CAM in distinguishing remnant shapes, to create
distance-based networks [22]. The Grad-CAM results are shown in Fig. 4(a).

Based on the network density $D$ and the average distance $L$ between nodes,
we are able to not only pinpoint when the morphology starts to develop
different shapes, but also quantify their relative proximity. According to
Fig. 4(b$\sim$e), we learn that the shape branches off in stage 3 or 4. In
contrast to Fig. 2(b,c) using $h$, $T$ and $R$ as parameters, the relative
position of shapes in Fig. 4(e) is based on the network properties. The
proximity of shapes, say, disk and sphere-II, disk and shotgun, and shotgun
and snowflake, revealed by both figures is consistent. What is surprising is
that snowflake and sphere-II that are separated by a third shape in Fig.
2(b,c) turn out to share a very similar network, whiles the opposite happens
to sphere and disk, i.e., bordering shapes in Fig. 2(b,c) turn out to drift
apart in Fig. 4(e).

###  III.3 Diameter and Depth of Crater

Figure 5(a$\sim$d) reveals how the width $w$ and depth $d$ of craters are
determined by the impact energy or equivalently the height $h$. Consistent
with previous studies, $w(h)$ is found to obey a power-law relation, i.e.,
$w\propto h^{\alpha}$ [24]. However, in contrast to $\alpha$=0.17 for water
droplet [25] and 0.25 for steel ball [6], both values occur in Fig. 5(a, c),
depending on whether $h$ is smaller or greater than $h_{c}$ as defined in Fig.
2(b, c). This change of behavior was reminiscent of the observation for
granular ball [8], except that the latter exhibits a discontinuous increase
and the $\alpha$ on both sides are close to 0.25.

The $d\propto h^{\beta}$ is monotonically increasing with $\beta=1/3$ [26] or
0.25 [6] for steel ball and $\beta=0.17$ [25] for water droplet, but exhibits
a discontinuous drop for granular ball [8]. In Fig. 5(b, d) for melted tin,
$d(h)$ is not monotonic - increase at $h<h_{c}$, and decrease continuously
when $h_{c}<h<h_{s}$ before finally saturate at $h>h_{s}$. Note that the data
points in Fig. 5(b, d) for melted tin at $h<h_{c}$ are higher than that for
steel ball with the same weight. At first glance, this is counter-intuitive
because melted tin ought to have less energy to dig the crater while expanding
its area at the same time. The answer to this puzzle is that the crater is
largely formed in stage 2 before melted tin has fully retracted to become a
sphere or disk. It was checked [22] that while the bottom of this V-shape tin
sheet or, equivalently, $d$ for crater is deeper than steer ball, their
centers of mass are in fact at roughly the same depth in stage 2. This
explanation is consistent with the fact that the discrepancy of $d$ for melted
tin and steer ball increases with temperature.

By compiling the data in Fig. 5(a$\sim$d), we are able to obtain $d\propto
w^{0.61}$ for $h>h_{s}$ in Fig. 5(e) which is in line with the consensus that
there exists a power law and the exponent is greater than 0, but less than 1
for complex meteorite craters [27]. However, this consistency is alarming
because our data points for $h>h_{s}$ in fact consist of several constant
functions of $w$. In other words, the non-zero exponent in Fig. 5(e) is
spurious because it is an artifact of compiling data from different $T$ and
$R$, which inevitably is the case for meteorite scholars.

##  IV Conclusion and Discussions

In contrast to solid and liquid projectiles used in previous studies of impact
experiments, we adopt melted tin to better simulate the high temperature
involved in real meteorite events. Modifying tin temperature enables us to
discern V- and U-shape craters, analogous to the distinction of simple and
complex craters in real events. Similar to previous research, we also found
the crater width to follow a power-law relation with the impact energy, except
that low and high-energy impacts turn out to exhibit different exponents. The
fact that the correlation of depth and width is lost as the depth saturates at
high energies allows us to uncover a possible misinterpretation when combining
data from different parameters, as is inevitable in field studies.
Specifically, we were misled to conclude that depth vs width follows a power
law when data from different tin diameter and temperature are plotted on the
same graph, which in fact consists of several constant functions of $d(w)$.

To aid the analysis and visualization of the complex interactions between
melted tin and granules, we employed a couple of technical methods, including
deep learning and image processing. Mapping the complex patterns generated by
tin in its phase transition to a network allows us to quantitatively compare
and elucidate the solidification process. We expect our conclusions are not
unique to the use of melted tin, but general to other substances that
experience a similar phase transition during the impact process, such as
thermoplastic materials. It is our hope that the advantages of recruiting such
a material will encourage more applications. For instance, the fluid-
solidification mechanisms in porous granules, which are highly applicable to
additive manufacturing or metal production processes [28] and extreme
ultraviolet (EUV) lithography in chip production [29], where molten tin is
used to generate EUV radiation. In a nutshell, our work is not limited to
academic interests, but of potential industrial applications.

## References

  * [1]
  * [2] [∗]<EMAIL_ADDRESS>  * [3] D. J. Roddy, Impact and Explosion Cratering, edited by D. J. Roddy, R. O. Pepin and R. B. Merril (Pergamon Press, New York, 1978) p. 185.
  * [4] H. J. Melosh, Impact Cratering: A Geologic Process (Oxford University Press, New York, 1989).
  * [5] K. Holsapple et al., Annu. Rev. Earth Planet. Sci. 21, 333 (1993).
  * [6] A. M. Walsh et al., Phys. Rev. Lett. 91, 104301 (2003).
  * [7] G. Delon et al., Phys. Rev. E. 84, 046320 (2011).
  * [8] F. Pacheco-Va´zquez et al., Phys. Rev. Lett. 107, 218001 (2011).
  * [9] D. Pimentel et al., Science 267, 1117 (1995).
  * [10] D. J. Furbish et al., J. Geophys. Res. 112, F01001. (2007).
  * [11] L. Zhang et al., J, Thermal Spray Technology 30, 181 (2021).
  * [12] Y. Chen et al., Scientia Sinica Physica, Mechanica & Astronomica 41, 1104 (2011).
  * [13] S. R. Schwartz et al., Granular Matter 14, 363 (2012).
  * [14] J. de Ruiter et al., Nature Phys. 14, 35 (2018).
  * [15] Hiroaki Katsuragi, Phys. Rev. Lett. 104, 218001 (2010).
  * [16] E. Nefzaoui et al. Experimental Thermal and Fluid Science 41, 43 (2012).
  * [17] T. R. Powers et al., Phys. Fluids 10, 1052 (1998).
  * [18] M. M. Driscoll et al., Phys. Rev. Lett. 107, 154502 (2011).
  * [19] B. M. French, LPI Contribution No. 954, Lunar and Planetary Institute (Houston, 1998), pp. 120.
  * [20] R. S. Dietz, Science 105, 42 (1947).
  * [21] R. S. Dietz, Popular Astronomy 54, 465 (1946).
  * [22] Supplemental Information
  * [23] R. R. Selvaraju, IEEE International Conference on Computer Vision (ICCV) (Venice, Italy, 2017) pp. 618.
  * [24] Yu Matsuda et al., Phys. Rev. E 99, 032906 (2019).
  * [25] Runchen Zhao et al., PNAS 112, 342 (2014).
  * [26] J. S. Uehara et al., Phys. Rev. Lett. 90, 194301 (2003).
  * [27] J. R. Stuart et al., Meteorit. Planetary Sci. 53, 583 (2018).
  * [28] R. O. David et al., PNAS 116, 8184 (2019).
  * [29] P. Kant et al., PNAS 117, 2788 (2020).