††thanks: Present address: Graduate School of Science and Engineering, Saitama
University, Saitama 338-8570, Japan<EMAIL_ADDRESS>
# Negatively enhanced thermopower near a Van Hove singularity in electron-
doped Sr2RuO4

Rei Nishinakayama1<EMAIL_ADDRESS>Yoshiki J. Sato1 Takayoshi Yamanaka2
Yoshiteru Maeno3 Hiroshi Yaguchi1 Naoki Kikugawa4 Ryuji Okazaki1
<EMAIL_ADDRESS>1Department of Physics and Astronomy, Tokyo University of
Science, Noda 278-8510, Japan
2Institute for Materials Research, Tohoku University, Sendai 980-8577, Japan
3Toyota Riken-Kyoto University Research Center (TRiKUC), Kyoto 606-8501, Japan
4National Institute for Materials Science, Tsukuba 305-0003, Ibaraki, Japan

###### Abstract

The layered perovskite Sr2RuO4 serves as a model material of the two-
dimensional (2D) Fermi liquid but also exhibits various emergent phenomena
including the non-Fermi-liquid (NFL) behavior under external perturbations
such as uniaxial pressure and chemical substitutions. Here we present the
thermoelectric transport of electron-doped system Sr2-yLayRuO4, in which a
filling-induced Lifshitz transition occurs at the Van Hove singularity (VHS)
point of $y\approx 0.2$. We find that the sign of the low-temperature
thermopower becomes negative only near the VHS point, where the NFL behavior
has been observed in the earlier work. This observation is incompatible with
either a numerical calculation within a constant relaxation-time approximation
or a toy-model calculation for the 2D Lifshitz transition adopting an elastic
carrier scattering. As a promising origin of the observed negatively enhanced
thermopower, we propose a skewed NFL state, in which an inelastic scattering
with a considerable odd-frequency term plays a crucial role to negatively
enhance the thermopower.

##  I Introduction

Beyond the well-established Fermi-liquid (FL) picture, non-Fermi-liquid (NFL)
state has been intensively investigated as an essential concept for various
quantum phenomena [1, 2, 5, 6, 3, 4, 7]. Besides the well-known Tomonaga-
Luttinger liquid in one dimension, the prototypical NFL state appears in a
vicinity of the quantum critical point (QCP), in which a thermodynamic phase
transition into an ordered state is suppressed by tuning external parameters.
In the NFL state near QCP, the finite-temperature properties such as the
electronic specific heat and the electrical resistivity drastically deviate
from the FL behavior, as widely seen in correlated metals including
transition-metal oxides and heavy fermions.

A yet unsolved, fundamental issue is how the FL picture is modified in a
vicinity of the Lifshitz transition [8], an electronic topological transition
associated with the change in the topology of the Fermi surfaces. The Lifshitz
transition itself is ubiquitous; it is driven by various parameters such as
pressure [9, 10, 11], magnetic field [12, 13, 14, 15], band filling [16, 17,
18], and even temperature through the temperature-dependent chemical potential
[19, 20, 21, 22]. In spite of many observed examples, the kinetic properties
near the Lifshitz transition are complex and still controversial owing to a
peculiar energy-dependent relaxation time [23, 24], which should become more
complicated in correlated metals.

To tackle this problem, here we focus on the quasi-two-dimensional (q-2D) FL
material Sr2RuO4 [25, 26, 27, 28]. Indeed, its q-2D Fermi surfaces consisting
of hole-like $\alpha$ and electron-like $\beta$ and $\gamma$ sheets have been
accurately verified by the dHvA and the ARPES experiments [29, 30, 31, 32].
The calculated Fermi surfaces of Sr2RuO4 are shown in Fig. 1(a). The normal-
state nature in such a q-2D multi-band system is well understood within the FL
picture [26, 33, 34], whereas the pairing mechanism of the superconducting
state is still an unsolved issue [35, 36, 37, 38, 39]. Significantly, its
superconducting transition temperature $T_{\rm c}$ is enhanced to $T_{\rm
c}\approx 3.5$ K under compressional stress [40, 41, 42, 43], and near the
critical compression point in which $T_{\rm c}$ has a maximum value. Such
variation of $T_{\rm c}$ corresponds to the sharp peak in the density of
states (DOS) associated with a Van Hove singularity (VHS) point; such a
topology change in the $\gamma$ band is indeed observed by the ARPES [44].
Most importantly, the resistivity clearly deviates from the FL behavior near
the VHS point [45, 46], and also nontrivial electronic states such as entropic
anomaly [47] have been observed near the Lifshitz transition [48, 49, 50].
Thus, this layered material offers a suitable platform to investigate the NFL
nature near the Lifshitz transition.

In this paper, we report a thermopower study of the electron-doped system
Sr2-yLayRuO4 [51, 52, 53], in which a filling-induced Lifshitz transition
occurs near the critical concentration $y_{\rm c}\approx 0.2$ [31]. The
calculated Fermi surfaces for $y=y_{\rm c}$ and $y>y_{\rm c}$ are depicted in
Figs. 1(b) and 1(c), respectively; the topology of the $\gamma$ sheet changes
at $y=y_{c}$. In Sr2-yLayRuO4, similar electronic features including NFL
transport [51] and the effective mass enhancement [31] have been clearly
observed near $y=0.2$, offering a complementary approach toward such an
intriguing issue on the Lifshitz transition. We find that the low-temperature
thermopower depends on the La content $y$ and that the sign of the thermopower
becomes negative only near the VHS point. This is in contrast to the results
of the Hall effect measurements in which the Hall coefficient exhibits no
significant anomaly near $y_{\rm c}$ [52]. We also show that the present
experimental results cannot be explained either by numerical calculation
results within a constant relaxation-time approximation or by a simple model
for the 2D neck-disruption-type Lifshitz transition with an elastic carrier
scattering. Instead, we propose a skewed NFL state [54] as a promising
explanation for our results, in which an odd-frequency inelastic scattering is
considered. The skewed NFL state is indeed a unique state of matter, as it has
been studied as a nature of strange metal in cuprate superconductors and may
also have a relevance to the transport properties in twisted bilayer graphene
[54]. This skewed NFL state strengthens the electron-hole asymmetry owing to
the dominant odd-frequency term in the scattering rate, providing a crucial
role for negatively enhanced thermopower near the Lifshitz transition.

##  II Experimental

Single crystals of Sr2-yLayRuO4 were grown by a floating-zone method [51, 52,
53]. Typical dimension of the single crystals is $3\times 1\times 0.1$ mm3.
The in-plane thermopower was measured by a steady-state technique using a
manganin-constantan differential thermocouple in a closed-cycle refrigerator
[55, 56]. A typical temperature gradient of 0.5 K/mm, which is adjusted along
with the bath temperature, was applied along the in-plane direction using a
resistive heater and the distance between the thermocouple contacts is about 1
mm. The thermoelectric voltage from the wire leads was subtracted.

Figure 1:  The cross-sectional view of the calculated Fermi surfaces of
Sr2-yLayRuO4 at $k_{z}=0$ plane for (a) $y=0$, (b) $y=y_{\rm c}$, and (c)
$y>y_{\rm c}$, drawn by using FermiSurfer program [63].

##  III Results and discussion

Figure 2(a) shows the temperature dependence of the in-plane thermopower $S$
of Sr2-yLayRuO4. It is known that $T_{\rm c}$ decreases with La substitution
and is completely suppressed for $y$ greater than 0.04 [51]. For the parent
compound, the present results well agree with those of the thermopower in
previous reports [57, 58, 56, 59, 60]. The thermopower of Sr2RuO4 has also
been studied by the dynamical mean-field theory [61]. In the La-substituted
compounds, overall behavior of the thermopower is similar to that of the
parent crystal.

Figure 2:  (a) Temperature dependence of the thermopower $S$ of Sr2-yLayRuO4
$(0\leq y\leq 0.27)$ measured along the $ab$ plane direction. (b) Temperature
dependence of $dS/dT$. The arrows show the characteristic temperature $T^{*}$
below which the $dS/dT$ deviates from the linear temperature dependence of
$dS/dT$ at higher temperatures. The data are shifted vertically for clarity.
(c) Temperature dependence of $S/T$. (d) The La content $y$ dependence of
$S/T$ measured at 4 K (solid circles, left axis), $T^{*}$ (open squares, right
axis), and $T_{\rm M}$ (open circles, right axis). $T_{\rm M}$ is defined at
the peak temperature in the magnetic susceptibility [53]. The blue dotted line
is a guide to the eye to represent the $y$ dependence of these characteristic
temperatures. The vertical dashed line show the critical La content $y_{\rm
c}$.

The thermopower in Sr2RuO4 was analyzed in the differential form $dS/dT$ [Fig.
2(b)] to examine a characteristic temperature, and an anomaly was found near
$T^{*}\approx 25$ K, below which $dS/dT$ increases with decreasing temperature
[57]. Subsequently, through the Seebeck and the Nernst measurements, Xu et al.
have suggested that the coherence is developed below $T^{*}$ [58]. As
displayed in the right axis of Fig. 2(d), $T^{*}$ systematically decreases
with the La content $y$. It should be noted that the magnetic susceptibility
of Sr2RuO4 is Pauli paramagnetic, but the temperature dependence exhibits a
small peak structure at $T_{\rm M}\approx 30$ K [26], below which the FL
picture is well defined. The La content $y$ dependence of $T_{\rm M}$ taken
from Ref. [53] is also plotted in the right axis of Fig. 2(d). Notably, both
$T_{\rm M}$ and $T^{*}$ show similar $y$ dependence, indicating that the
characteristic temperature below which the coherence is formed in the
correlated carriers decreases with increasing $y$. This trend is consistent
with the NFL behavior near $y=0.2$ [51] and also signifies the inherent role
of the carrier scattering at the Lifshitz transition, as will be discussed
later.

In Fig. 2(c), we also plot $S/T$ of Sr2-yLayRuO4 as a function of $T$. The
low-temperature $S/T$ behavior notably depends on both temperature and the La
content $y$, suggesting the considerable change in the Fermi surfaces in the
present La content range as reported in earlier reports [51, 52, 53, 31]. It
should be noted that negative thermopower is found at low temperatures only
for the $y=0.2$ crystal, which is near the critical La content $y_{\rm c}$
[51, 52, 53, 31]. Figure 2(d) displays the La content $y$ dependence of $S/T$
obtained at 4 K for the left axis. With increasing $y$, $S/T$ slightly
increases and a singularity is clearly observed near $y\approx 0.2$. Note that
a positive value of $S/T$ is recovered at the higher content $y=0.27$. To
explain the observed La content dependence of $S/T$, we have examined the
chemical potential dependence of $S/T$ calculated within a constant
relaxation-time approximation [62] (Appendix), but the calculated data [Fig.
S2(b)] is positively enhanced near the VHS points. Obviously, this discrepancy
originates from the energy dependence of the relaxation time, which is ignored
in the constant relaxation-time approximation method.

Figure 3:  2D neck-disruption-type Lifshitz transitions and the physical
properties. (a) The Fermi surface shape and (b) the DOS are symmetric around
the critical energy $E_{\rm c}$ at which the Lifshitz transition occurs. (c)
The conductivity and (d) the thermopower as a function of the chemical
potential, which corresponds to the amount of electron doping by La
substitution, for several temperatures. HT and LT represent high and low
temperatures, respectively. The vertical dotted lines show $E_{\rm c}$ and the
horizontal dashed line in (d) represents the contributions from the regular
part of the Fermi surfaces, which is far from the VHS points.

Here we discuss the energy dependence of the relaxation time
$\tau(\varepsilon)$ near the Lifshitz transition. Similar to the case of
Sr2RuO4, the 2D neck-disruption-type Lifshitz transition for the cylindrical
Fermi surfaces [Fig. 3(a)] has been investigated [64, 65, 66, 67]. In this
case, the energy dependence of the DOS $D(\varepsilon)$ shows logarithmic
divergence near the critical energy $E_{\rm c}$ at which the Lifshitz
transition occurs [Fig. 3(b)]. Such logarithmic DOS behavior is also confirmed
by the numerical calculation for Sr2RuO4 [Fig. A2(a) in Appendix]. Then,
through an elastic impurity scattering, the scattering probability acquires a
correction of the energy dependence of $1/\tau(\varepsilon)\propto
D(\varepsilon)$ [23, 24], which significantly affects the energy dependence of
the conductivity function $\sigma(\varepsilon)\simeq D_{0}v_{0}^{2}\tau$. Note
that the DOS $D_{0}$ and the velocity $v_{0}$ in this conductivity function
exhibit weak energy dependence because these mainly come from the electrons in
the regular parts of the Fermi surfaces, which are far from the VHS points
[68]. Figure 3(c) shows the calculated electrical conductivity for this simple
model (Appendix). The horizontal axis is the chemical potential and
corresponds to the amount of electron doping by La substitution. As a
consequence, the electrical conductivity $\sigma$ decreases near $E_{\rm c}$
and such a modification leads to a NFL-like resistivity of
$\rho(T)=\rho_{0}+AT^{n}$ with $n<2$ for Sr2RuO4 [46], as observed near the
VHS point for both La-substituted [51] and uniaxially compressed [45] cases.

In this model, however, as shown in Fig. 3(d), the thermopower should exhibit
positive and negative peaks with the same magnitudes below and above $E_{\rm
c}$, respectively [64, 65, 66], as indicated from the Mott formula [69]

$\displaystyle\frac{S}{T}\propto-\frac{1}{\sigma}\frac{\partial\sigma}{\partial\varepsilon}\sim-\frac{1}{\tau}\frac{\partial\tau}{\partial\varepsilon},$
(1)

where the energy dependence of $\tau$ is crucial as similar to the case of the
conductivity as mentioned before. It should be noted that such thermopower
behavior with positive and negative peaks is also obtained in the numerical
calculations for Sr2RuO4 in the elastic impurity scattering regime [46]. In
contrast, for Sr2-yLayRuO4, the low-temperature thermopower seems to be
enhanced only negatively near the critical content $y_{\rm c}\approx 0.2$
[Fig. 2(d)], which is difficult to explain with such a conventional neck
disruption case. Also, as indicated in Fig. 2(d), an inelastic electron-
electron scattering, not included in the model of Fig. 3, should be crucial
near $y_{\rm c}$.

Figure 4:  Scaling function $g(\omega/T)$ for an exponent $\nu=1$ and an
asymmetry parameter $\alpha=0$ and $1.5$ [54]. The inelastic scattering rate
$1/\tau_{\rm in}$ is given as $1/\tau_{\rm in}\propto(\pi
T)^{\nu}g(\omega/T)$.

To discuss the origin of the negatively enhanced thermopower near the VHS
point, we next consider a phenomenological model adopting a skewed NFL state
[54], in which an asymmetric inelastic scattering rate $1/\tau_{\rm
in}\propto(\pi T)^{\nu}g(\omega/T)$ characterized by a noneven scaling
function
$g(x)=|\Gamma(z)|^{2}\cosh(x/2)/[\cosh(\alpha/2)\Gamma\\{(1+\nu)/2\\}^{2}]$
becomes essential, where $\Gamma(z)$ is the $\Gamma$ function,
$z=(1+\nu)/2+i(x+\alpha)/2\pi$, $\omega=\varepsilon-\mu$ is a relative energy
from the chemical potential $\mu$, $\nu$ $(\leq 1)$ is an exponent, and
$\alpha$ is an parameter to induce the asymmetry in the scattering rate.
Figure 4 shows the examples of a scaling function $g(\omega/T)$ with $\nu=1$
for a symmetric (with an asymmetry parameter $\alpha=0$) and an asymmetric
($\alpha=1.5$) cases [54]. In the asymmetric case, owing to the odd-frequency
term of the scattering rate, the contribution of either electrons or holes
becomes stronger, and most importantly, the thermopower, a sensitive probe to
the electron-hole asymmetry, is enhanced either negatively or positively. We
infer that, along with the observation of the NFL resistivity near the VHS
point [51], the observed negative sign of the thermopower may be a hallmark of
such a skewed NFL state. Moreover, such unexpected sign change in the
thermopower has also been discussed in the NFL regime of cuprate
superconductors [70]. Note that the thermopower is also known as a sensitive
probe for the entropic properties such as magnetic fluctuations but the
present non-magnetic system may not be adapted to such situation. At this
stage, it is not easy to make a more quantitative calculation in this model,
and further theoretical study is necessary to quantitatively account for the
the sign change. Also, the frequency-dependent experiments such as the optical
conductivity measurement will be crucial to examine the skewed NFL state.

Lifshitz transitions in correlated matter thus bring an intriguing issue on
the NFL state. In this context, Sr2RuO4, in which a variety of electronic and
magnetic states emerges under external perturbations [71, 72, 73, 74, 75], is
of peculiar interest because it also exhibits the NFL behavior in Ti-
substituted system [76]. Moreover, recent experimental and theoretical studies
have revealed the appearance of enigmatic NFL states such as a Planckian metal
characterized by a linear temperature dependence of the resistivity [77, 78]
and the quantum critical phase in frustrated materials [79, 80], deepening
underlying physics of the NFL state of matter.

##  IV summary

In summary, we have measured the thermopower of the electron-doped system
Sr2-yLayRuO4 and observed an unusual sign change in the thermopower near the
Lifshitz transition at the VHS point $y_{\rm c}\approx 0.2$. We discuss the
thermopower in a skewed NFL state, in which an asymmetric frequency dependence
of the inelastic scattering rate is crucial, possibly leading to qualitative
explanation for the negative sign of the thermopower near $y_{\rm c}$. The
present results thus offer a fascinating playground to investigate a variety
of quantum phenomena of the correlated electrons near the Lifshitz transition.

## Acknowledgments

We appreciate Y. Fukumoto and R. Kurihara for discussion and R. Otsuki, H.
Shiina, and R. Taira for the assistance. This work was partly supported by
JSPS KAKENHI Grant No. 17H06136, No. 21H01033, No. 22H01166, No. 22H01168, and
No. 22K19093.

## Appendix

###  IV.1 First-principles calculations

In order to investigate the thermopower theoretically, we performed first-
principles calculations based on density functional theory (DFT) using Quantum
Espresso [81, 82, 83]. We used the projector-augmented-wave pseudopotentials
with the Perdew-Burke-Ernzerhof generalized-gradient-approximation (PBE-GGA)
exchange-correlation functional. The cut-off energies for plane waves and
charge densities were set to 70 and 560 Ry, respectively, and the $k$-point
mesh was set to $20\times 20\times 20$ uniform grid to ensure the convergence.
Using the obtained eigenvalues of the $n$-th band at $\bm{k}$ point
$E_{n,\bm{k}}$, the DOS $D(\varepsilon)=\sum_{n,\bm{k}}\delta(\varepsilon-
E_{n,\bm{k}})$ was obtained using the optimized tetrahedron method [84], where
$\delta$ is the delta function. We used on-site Coulomb energy $U=3.5$ eV and
exchange parameter $J=0.6$ eV for Ru ions [85], and performed full
relativistic calculations with spin-orbit coupling (${\rm DFT}+U+{\rm SOC}$).

Figure A1:  (a) Calculated band structure of Sr2RuO4 with ${\rm DFT}+U+{\rm
SOC}$ scheme (solid curves). The dashed curves represent the results of scalar
relativistic calculations and $U$ is not included. The Van Hove singularity is
at the M point, at which the $\gamma$ band is slightly above $E_{\rm F}$. (b)
High-symmetry points in the Brillouin zone.

Figure A1(a) shows the calculated electronic band structure near the Fermi
energy $E_{F}$ of the parent material, which well coincides with the results
in earlier studies [86, 87, 88, 89]. The depicted $k$ path is shown in Fig.
A1(b). The $\beta$ and $\gamma$ bands at the high-symmetry points $\Gamma$ and
Z split owing to the inclusion of SOC [88], and the $e_{g}$ bands are shifted
upward due to the on-site $U$. The calculated DOS is depicted in Fig. A2(a),
and the Van Hove singularity (VHS) point of the $\gamma$ band to show the cusp
anomaly is shifted slightly from $E\approx 50$ meV to $E\approx 30$ meV by
including $U+{\rm SOC}$. Such a trend is consistent with the results of ARPES
experiment [31, 32].

Figure A2:  (a) Calculated total and partial density of states with ${\rm
DFT}+U+{\rm SOC}$ scheme (solid filled curves). The dashed curves represent
the results of scalar relativistic calculations and $U$ is not included. (b)
Thermopower divided by temperature, $S/T$, calculated for several temperatures
within the relaxation time approximation. The horizontal axis shows the
chemical potential measured from the Fermi energy of the parent compound. The
black curves represent the total $S/T$ and the red, blue, and green curves
show the band-resolved data $S^{n}/T$ for the $\alpha$, $\beta$, and $\gamma$
bands, respectively. Varying $\mu$ across the VHS results in the enhancement
of the positive $S/T$.

###  IV.2 Calculated transport properties within a constant relaxation-time
approximation

To examine the experimentally observed singular behavior in $S/T$, we have
calculated the thermopower from the electronic band structure within the
constant relaxation time approximation. Note that the correlation effect of
$4d$ electrons is considerable in Sr2RuO4 and the calculation results with
local-density approximation are quantitatively different from the experimental
observations [61]. On the other hand, the thermopower behavior in correlated
metals is well described within a Fermi-liquid picture, where the correlation
effect is included in the carrier effective mass [69]. Also, in Sr2-yLayRuO4,
the electron-doping effect by La substitution is well explained within a rigid
band picture [51, 31]. We thus examine how the band structure affects the
thermopower.

The transport coefficients were calculated based on the linearized Boltzmann
equations under constant relaxation time approximation [62]. The transport
distribution function tensor $L_{ij}(\varepsilon)$ is calculated as

$\displaystyle
L_{ij}(\varepsilon)=\sum_{n}L_{ij}^{n}(\varepsilon)=\sum_{n}\sum_{\mathbb{k}}v_{i}v_{j}\tau\delta(\varepsilon-
E_{n,\mathbb{k}}),$ (2)

where $L_{ij}^{n}(\varepsilon)$ is the partial transport distribution function
tensor of the $n$-th band, $v_{i}$ is the $i$-th component of the band
velocity $\mathbb{v}=\frac{1}{\hbar}\nabla_{\mathbb{k}}E_{n,\mathbb{k}}$, and
$\tau$ is the relaxation time. We calculated the partial electrical
conductivity tensor of the $n$-th band of
$\sigma_{ij}^{n}(\mu)=e^{2}\int_{-\infty}^{\infty}d\varepsilon\left(-\frac{\partial
f_{0}}{\partial\varepsilon}\right)L_{ij}^{n}$, where $e$ is the elementary
charge and $f_{0}$ is the Fermi-Dirac distribution function for the chemical
potential $\mu$ and temperature $T$. The total electrical conductivity tensor
is $\sigma_{ij}(\mu)=\sum_{n}\sigma_{ij}^{n}(\mu)$. Similarly, the partial
Peltier conductivity tensor of $n$-th band $P_{ij}^{n}(\mu)=[\sigma
S]_{ij}^{n}(\mu)$ is

$\displaystyle
P_{ij}^{n}(\mu)=-\frac{e}{T}\int_{-\infty}^{\infty}d\varepsilon\left(-\frac{\partial
f_{0}}{\partial\varepsilon}\right)(\varepsilon-\mu)L_{ij}^{n},$ (3)

where $S_{ij}$ is the thermopower tensor. The total Peltier conductivity is
given as $P_{ij}(\mu)=\sum_{n}P_{ij}^{n}(\mu)$. Hereafter we consider the in-
plane component ($ij=aa$) only and the subscript will be omitted. The
thermopower of the $n$-th band is then obtained as $S^{n}=P^{n}/\sigma^{n}$,
and the total in-plane thermopower is given as
$S=P/\sigma=\sum_{n}P^{n}/\sum_{n}\sigma^{n}$ as is generally seen in multi-
band systems.

Figure A2(b) depicts the thermopower divided by temperature, $S/T$, calculated
for several temperatures within the constant relaxation time approximation.
The black curves represent the total $S/T$ and the red, blue, and green curves
show the band-resolved $S^{n}/T$ for $n=\alpha$, $\beta$, $\gamma$,
respectively. The horizontal axis is the chemical potential measured from the
Fermi energy of the parent compound and corresponds to the amount of electron
doping by La substitution. Note that the calculated values of $S/T$ are
significantly smaller than the experimental data because the electron
correlation effect is not accurately included in this scheme and should be
modified by using more precise methods such as the dynamical mean field theory
[61]. Nevertheless, characteristic features reflecting the Lifshitz transition
in Sr2RuO4 may be observed in the present calculations; the band-resolved
$S^{n}/T$ shows divergent behavior at low temperatures near $0.3$ eV ($\beta$
band) and 30 meV ($\gamma$ band) corresponding to the VHS points of DOS for
each band [Fig. A2(a)], while $S^{n}/T$ exhibits almost no temperature
dependence for the $\alpha$ band like a simple metal. In general, the VHS
points of DOS strongly affect the thermopower [90]. It is now important that
the experimental data of $S/T$ seems to negatively diverge near the critical
concentration [Fig. 2(d) in the main text], while the calculated data is
positively enhanced near the VHS points [Fig. A2(b)]. This discrepancy
obviously originates from the energy dependence of the relaxation time ignored
in the constant relaxation time approximation.

###  IV.3 Peculiar energy-dependent relaxation time and thermopower near the
Lifshitz transition

Figure A3:  Two types of the Lifshitz transitions for 3D case and the physical
properties: (a-d) void formation and (e-h) neck disruption. For each panel,
the dashed lines represent the contributions from the regular part of the
Fermi surfaces. The vertical dotted lines show the critical energy $E_{c}$ at
which the Lifshitz transition occurs. HT and LT represent high and low
temperatures, respectively.

Here we briefly review the three-dimensional (3D) case to see the significance
of the scattering process. In the 3D case, the Lifshitz transition is
categorized into two types of topological change in the Fermi surface called
void formation and neck disruption, which are schematically shown in Figs.
A3(a) and A3(e), respectively [8]. In the void formation case, for example,
the DOS behaves as $D(\varepsilon)-D_{0}\sim|\varepsilon-E_{c}|^{1/2}$, where
$D_{0}$ is the DOS from the regular part of the Fermi surface (large void) and
$E_{c}$ is the critical energy above which a new void appears [Fig. A3(b)]. At
first glance, such a small void with almost zero carrier velocity seems to
give no contribution to the transport properties. Through the scattering
process, however, electrons in the regular part exchange the momenta with
those in the singular part (small void) and get virtually into the singular
part [68, 23, 24]. As a result, according to the golden rule, the scattering
probability acquires a correction of the energy dependence of
$1/\tau(\varepsilon)\propto D(\varepsilon)$, and then the electrical
conductivity $\sigma\propto\tau$ decreases above $E_{c}$ [Fig. A3(c)]. Note
that the energy dependence of $\tau$ is essential here as an approximate form
$\sigma(\mu)\simeq
e^{2}\int_{-\infty}^{\infty}d\varepsilon\left(-\frac{\partial
f_{0}}{\partial\varepsilon}\right)Dv^{2}\tau\simeq
e^{2}D_{0}v_{0}^{2}\tau(\mu)$, because the DOS $D_{0}$ and the velocity
$v_{0}$ in the conductivity mainly come from the electrons in the large void
[68].

According to the Mott formula, the thermopower $S$ is given as

$\displaystyle\frac{S}{T}\propto-\frac{1}{\sigma}\frac{\partial\sigma}{\partial\varepsilon}\sim-\frac{1}{\tau}\frac{\partial\tau}{\partial\varepsilon},$
(4)

where the energy dependence of $\tau$ is crucial as similar to the case of the
conductivity, and thus it shows a sharp peak structure at low temperatures
[Fig. A3(d)]. The similar situation occurs in the case of the neck disruption
[Figs. A3(e-h)] and the thermopower is also enhanced positively near the
critical point. Note that the singularities in $\sigma$ and $S$ are smeared
with increasing temperature. These thermoelectric singularities in 3D case
have been experimentally observed in the Li-Mg alloy [91].

###  IV.4 Calculations of the transport properties for the Lifshitz
transitions

Here we show the calculation details for the transport coefficients near the
Lifshitz transition by using the energy-dependent scattering time. The
electrical conductivity $\sigma$ and the Peltier conductivity $P=\sigma S$ are
given as

$\displaystyle\begin{bmatrix}\sigma\\\ P\\\ \end{bmatrix}$
$\displaystyle=\begin{bmatrix}e^{2}\\\ -\frac{e}{T}\\\
\end{bmatrix}\int_{-\infty}^{\infty}d\varepsilon\left(-\frac{\partial
f_{0}}{\partial\varepsilon}\right)\begin{bmatrix}1\\\ \varepsilon-\mu\\\
\end{bmatrix}L$ (5)
$\displaystyle=\frac{e}{4k_{B}T^{2}}\int_{-\infty}^{\infty}\frac{d\omega}{\cosh^{2}(\beta\omega/2)}\begin{bmatrix}eT\\\
-\omega\\\ \end{bmatrix}L,$ (6)

where $\omega=\varepsilon-\mu$ is the relative energy measured from the
chemical potential. Using the energy-dependent scattering time, the transport
function near the Lifshitz transition is approximately given as

$\displaystyle L=D_{0}v_{0}^{2}\tau(\varepsilon),$ (7)

and the transport coefficients are given as

$\displaystyle\begin{bmatrix}\sigma\\\ P\\\ \end{bmatrix}$
$\displaystyle=\frac{eD_{0}v_{0}^{2}}{4k_{B}T^{2}}\int_{-\infty}^{\infty}\frac{d\omega}{\cosh^{2}(\beta\omega/2)}\begin{bmatrix}eT\\\
-\omega\\\ \end{bmatrix}\tau(\omega),$ (8)

where the scattering time is model-dependent as described below.

For the 3D void formation case in Figs. A3(a-d), the density of states is
expressed as

$\displaystyle D(\varepsilon)$ $\displaystyle\sim D_{0}+a|\varepsilon-
E_{c}|^{1/2}\theta(\varepsilon-E_{c}),$ (9)

where $a$ $(>0)$ is a constant and $\theta$ is the Heaviside step function, as
is shown in Fig. A3(b). The scattering time is then given as

$\displaystyle\tau(\varepsilon)$ $\displaystyle\sim D(\varepsilon)^{-1}$ (10)
$\displaystyle\sim D_{0}-a|\varepsilon-E_{c}|^{1/2}\theta(\varepsilon-E_{c})$
(11) $\displaystyle=D_{0}-a|\omega+Z|^{1/2}\theta(\omega+Z),$ (12)

where $Z=\mu-E_{c}$ is the chemical potential measured from the critical
energy. The transport coefficients are now calculated and the thermopower $S$
is give as $S=P/\sigma$. Note that although the regular part $D_{0}$ also
depends on the energy [23, 24], the energy dependence of the singular part is
much significant. Indeed, the calculation results shown in Fig. A3 are similar
to the earlier works.

For the 3D neck disruption case [Figs. A3(e-h)], the scattering time is given
as

$\displaystyle\tau(\varepsilon)$ $\displaystyle\sim
D_{0}+a|E_{c}-\varepsilon|^{1/2}\theta(E_{c}-\varepsilon)$ (13)
$\displaystyle=D_{0}+a|-\omega-Z|^{1/2}\theta(-\omega-Z),$ (14)

which is similar to the case of void formation.

For the symmetric 2D neck disruption case [Figs. 3(a-d) in the main text], the
density of states near the Lifshitz transition is given as a logarithmic form
of

$\displaystyle D(\varepsilon)$ $\displaystyle\sim\ln\frac{t}{|\varepsilon-
E_{c}|},$ (15)

where $t$ $(>0)$ is a constant as is shown in Fig. 3(b) in the main text. The
scattering time is given as

$\displaystyle\tau(\varepsilon)$
$\displaystyle\sim\left(\ln\frac{t}{|\varepsilon-
E_{c}|}\right)^{-1}=\left(\ln\frac{t}{|\omega+Z|}\right)^{-1}.$ (16)

## References

  * [1] J. A. Hertz, Quantum Critical Phenomena, Phys. Rev. B 14, 1165 (1976).
  * [2] A. J. Millis, Effect of a Nonzero Temperature on Quantum Critical Points in Itinerant Fermion Systems, Phys. Rev. B 48, 7183 (1993).
  * [3] G. R. Stewart, Non-Fermi-liquid behavior in $d$\- and $f$-electron metals, Rev. Mod. Phys. 73, 797 (2001).
  * [4] H. v. Löhneysen, A. Rosch, M. Vojta, and Peter Wölfle, Fermi-liquid instabilities at magnetic quantum phase transitions, Rev. Mod. Phys. 79, 1015 (2007).
  * [5] S. Sachdev, Quantum Magnetism and Criticality, Nat. Phys. 4, 173 (2008).
  * [6] P. Gegenwart, Q. Si, and F. Steglich, Quantum Criticality in Heavy-Fermion Metals, Nat. Phys. 4, 186 (2008).
  * [7] D. Chowdhury, A. Georges, O. Parcollet, and S. Sachdev, Sachdev-Ye-Kitaev models and beyond: Window into non-Fermi liquids, Rev. Mod. Phys. 94, 035004 (2022).
  * [8] I. Lifshitz, Anomalies of electron characteristics of a metal in the high pressure region, Sov. Phys. JETP 11, 1130-1135 (1960).
  * [9] T. Nishimura, H. Sakai, H. Mori, K. Akiba, H. Usui, M. Ochi, K. Kuroki, A. Miyake, M. Tokunaga, Y. Uwatoko, K. Katayama, H. Murakawa, and N. Hanasaki, Large Enhancement of Thermoelectric Efficiency Due to a Pressure-Induced Lifshitz Transition in SnSe, Phys. Rev. Lett. 122, 226601 (2019).
  * [10] M. Krottenmüller, M. Vöst, N. Unglert, J. Ebad-Allah, G. Eickerling, D. Volkmer, J. Hu, Y. L. Zhu, Z. Q. Mao, W. Scherer, and C. A. Kuntscher, Indications for Lifshitz transitions in the nodal-line semimetal ZrSiTe induced by interlayer interaction, Phys. Rev. B 101, 081108(R) (2020).
  * [11] S. Sen and G.-Y. Guo, Pressure induced Lifshitz transition in ThFeAsN, Phys. Rev. Materials 4, 104802 (2020).
  * [12] H. Pfau, R. Daou, S. Friedemann, S. Karbassi, S. Ghannadzadeh, R. Küchler, S. Hamann, A. Steppke, D. Sun, M. König, A. P. Mackenzie, K. Kliemt, C. Krellner, and M. Brando, Cascade of Magnetic-Field-Induced Lifshitz Transitions in the Ferromagnetic Kondo Lattice Material YbNi4P2, Phys. Rev. Lett. 119, 126402 (2017).
  * [13] D. Schulze Grachtrup, N. Steinki, S. Süllow, Z. Cakir, G. Zwicknagl, Y. Krupko, I. Sheikin, M. Jaime, and J. A. Mydosh, Magnetic phase diagram and electronic structure of UPt2Si2 at high magnetic fields: A possible field-induced Lifshitz transition, Phys. Rev. B 95, 134422 (2017).
  * [14] A. Pourret, S. G. Sharapov, T. D. Matsuda, G. Knebel, G. Zwicknagl, and A. A. Varlamov, Transport Spectroscopy of the Field Induced Cascade of Lifshitz Transitions in YbRh2Si2, J. Phys. Soc. Jpn. 88, 104702 (2019).
  * [15] L. Wu, S. Chi, H. Zuo, G. Xu, L. Zhao, Y. Luo, and Z. Zhu, Field-induced Lifshitz transition in the magnetic Weyl semimetal candidate PrAlSi, npj Quantum Materials 8, 4 (2023).
  * [16] Y. Okamoto, A. Nishio, and Z. Hiroi, Discontinuous Lifshitz transition achieved by band-filling control in NaxCoO2, Phys. Rev. B 81, 121102(R) (2010).
  * [17] J. Kwon, M. Kim, D. Song, Y. Yoshida, J. Denlinger, W. Kyung, and C. Kim, Lifshitz-Transition-Driven Metal-Insulator Transition in Moderately Spin-Orbit-Coupled Sr2-xLaxRhO4, Phys. Rev. Lett. 123, 106401 (2019).
  * [18] N. Ito, M. Ishii, and R. Okazaki, Enhanced Seebeck coefficient by a filling-induced Lifshitz transition in KxRhO2, Phys. Rev. B 99, 041112(R) (2019).
  * [19] Y. Wu et al., Temperature-Induced Lifshitz Transition in WTe2, Phys. Rev. Lett. 115, 166602 (2015).
  * [20] Y. Zhang et al., Electronic evidence of temperature-induced Lifshitz transition and topological nature in ZrTe5, Nat. Commun. 8, 15512 (2017).
  * [21] F. C. Chen, Y. Fei, S. J. Li, Q. Wang, X. Luo, J. Yan, W. J. Lu, P. Tong, W. H. Song, X. B. Zhu, L. Zhang, H. B. Zhou, F. W. Zheng, P. Zhang, A. L. Lichtenstein, M. I. Katsnelson, Y. Yin, N. Hao, and Y. P. Sun, Temperature-Induced Lifshitz Transition and Possible Excitonic Instability in ZrSiSe, Phys. Rev. Lett. 124, 236601 (2020).
  * [22] Y. Tian, Y. Zhu, R. Li, Z. Mao, and J. H. Ross, Jr., NMR determination of Van Hove singularity and Lifshitz transitions in the nodal-line semimetal ZrSiTe, Phys. Rev. B 104, L041105 (2021).
  * [23] A. A. Varlamov, V. S. Egorov, and A.V. Pantsulaya, Kinetic properties of metals near electronic topological transitions (2 1/2-order transitions), Adv. Phys. 38, 469 (1989).
  * [24] A. A. Varlamov, Y. M. Galperin, S. G. Sharapov,, and Y. Yerin, Concise guide for electronic topological transitions, Low Temp. Phys. 47, 672 (2021).
  * [25] Y. Maeno, H. Hashimoto, K. Yoshida, S. Nishizaki, T. Fujita, J. G. Bednorz, and F. Lichtenberg, Superconductivity in a layered perovskite without copper, Nature 372, 532 (1994).
  * [26] A. P. Mackenzie and Y. Maeno, The superconductivity of Sr2RuO4 and the physics of spin-triplet pairing, Rev. Mod. Phys. 75, 657 (2003).
  * [27] A. P. Mackenzie, T. Scaffidi, C. W. Hicks, and Y. Maeno, Even odder after twenty-three years: The superconducting order parameter puzzle of Sr2RuO4, npj Quant. Mater. 2, 40 (2017).
  * [28] S. A. Kivelson, A. C. Yuan, B. Ramshaw, and R. Thomale, A proposal for reconciling diverse experiments on the superconducting state in Sr2RuO4, npj Quant. Mater. 5, 43 (2020).
  * [29] A. P. Mackenzie, S. R. Julian, A. J. Diver, G. J. McMullan, M. P. Ray, G. G. Lonzarich, Y. Maeno, S. Nishizaki, and T. Fujita, Quantum Oscillations in the Layered Perovskite Superconductor Sr2RuO4, Phys. Rev. Lett. 76, 3786 (1996).
  * [30] A. Damascelli, D. H. Lu, K. M. Shen, N. P. Armitage, F. Ronning, D. L. Feng, C. Kim, Z.-X. Shen, T. Kimura, Y. Tokura, Z. Q. Mao, and Y. Maeno, Fermi Surface, Surface States, and Surface Reconstruction in Sr2RuO4, Phys. Rev. Lett. 85, 5194 (2000).
  * [31] K. M. Shen, N. Kikugawa, C. Bergemann, L. Balicas, F. Baumberger, W. Meevasana, N. J. C. Ingle, Y. Maeno, Z.-X. Shen, and A. P. Mackenzie, Evolution of the Fermi Surface and Quasiparticle Renormalization through a van Hove Singularity in Sr2-yLayRuO4, Phys. Rev. Lett. 99, 187001 (2007).
  * [32] A. Tamai, M. Zingl, E. Rozbicki, E. Cappelli, S. Riccò, A. de la Torre, S. McKeown Walker, F. Y. Bruno, P. D. C. King, W. Meevasana, M. Shi, M. Radović, N. C. Plumb, A. S. Gibbs, A. P. Mackenzie, C. Berthod, H. U. R. Strand, M. Kim, A. Georges, and F. Baumberger, High-Resolution Photoemission on Sr2RuO4 Reveals Correlation-Enhanced Effective Spin-Orbit Coupling and Dominantly Local Self-Energies, Phys. Rev. X 9, 021048 (2019).
  * [33] D. Stricker, J. Mravlje, C. Berthod, R. Fittipaldi, A. Vecchione, A. Georges, and D. van der Marel, Optical Response of Sr2RuO4 Reveals Universal Fermi-Liquid Scaling and Quasiparticles Beyond Landau Theory, Phys. Rev. Lett. 113, 087404 (2014).
  * [34] C. Bergemann, A. P. Mackenzie, S. R. Julian, D. Forsythe, and E. Ohmichi, Quasi-two-dimensional Fermi liquid properties of the unconventional superconductor Sr2RuO4, Adv. Phys. 52, 639 (2003).
  * [35] A. Pustogow, Y. Luo, A. Chronister, Y.-S. Su, D. A. Sokolov, F. Jerzembeck, A. P. Mackenzie, C. W. Hicks, N. Kikugawa, S. Raghu, E. D. Bauer, and S. E. Brown Constraints on the superconducting order parameter in Sr2RuO4 from oxygen-17 nuclear magnetic resonance, Nature 574, 72 (2019).
  * [36] K. Ishida, M. Manago, K. Kinjo, and Y. Maeno, Reduction of the 17O Knight Shift in the Superconducting State and the Heat-up Effect by NMR Pulses on Sr2RuO4, J. Phys. Soc. Jpn. 89, 034712 (2020).
  * [37] S. Benhabib, C. Lupien, I. Paul, L. Berges, M. Dion, M. Nardone, A. Zitouni, Z. Q. Mao, Y. Maeno, A. Georges, L. Taillefer, and C. Proust Ultrasound evidence for a two-component superconducting order parameter in Sr2RuO4, Nat. Phys. 17, 194 (2021).
  * [38] S. Ghosh, A. Shekhter, F. Jerzembeck, N. Kikugawa, D. A. Sokolov, M. Brando, A. P. Mackenzie, C, W. Hicks, and B. J. Ramshaw, Thermodynamic evidence for a two-component superconducting order parameter in Sr2RuO4, Nat. Phys. 17, 199 (2021).
  * [39] V. Grinenko, S. Ghosh, R. Sarkar, J.-C. Orain, A. Nikitin, M. Elender, D. Das, Z. Guguchia, F. Brückner, M. E. Barber, J. Park, N. Kikugawa, D. A. Sokolov, J. S. Bobowski, T. Miyoshi, Y. Maeno, A. P. Mackenzie, H. Luetkens, C. W. Hicks, and H.-H. Klauss, Split superconducting and time-reversal symmetry-breaking transitions in Sr2RuO4 under stress, Nat. Phys. 17, 748 (2021).
  * [40] C. W. Hicks, D. O. Brodsky, E. A. Yelland, A. S. Gibbs, J. A. N. Bruin, M. E. Barber, S. D. Edkins, K. Nishimura, S. Yonezawa, Y. Maeno, and A. P. Mackenzie, Strong Increase of $T_{c}$ of Sr2RuO4 Under Both Tensile and Compressive Strain, Science 244, 283 (2014).
  * [41] H. Taniguchi, K. Nishimura, S. K. Goh, S. Yonezawa, and Y. Maeno, Higher-$T_{c}$ Superconducting Phase in Sr2RuO4 Induced by In-Plane Uniaxial Pressure, J. Phys. Soc. Jpn. 84, 014707 (2015).
  * [42] A. Steppke, L. Zhao, M. E. Barber, T. Scaffidi, F. Jerzembeck, H. Rosner, A. S. Gibbs, Y. Maeno, S. H. Simon, A. P. Mackenzie, and C. W. Hicks, Strong peak in $T_{c}$ of Sr2RuO4 under uniaxial pressure, Science 355, eaaf9398 (2017).
  * [43] Y.-S. Li, N. Kikugawa, D. A. Sokolov, F. Jerzembeck, A. S. Gibbs, Y. Maeno, C. W. Hicks, J. Schmalian, M. Nicklas, and A. P. Mackenzie, High-sensitivity heat-capacity measurements on Sr2RuO4 under uniaxial pressure, Proc. Natl. Acad. Sci. U.S.A. 118, e2020492118 (2021).
  * [44] V. Sunko, E. A. Morales, I Marković, M. E. Barber, D. Milosavljević, F. Mazzola, D. A. Sokolov, N. Kikugawa, C. Cacho, P. Dudin, H. Rosner, C. W. Hicks, P. D. C. King, and A. P. Mackenzie, Direct observation of a uniaxial stress-driven Lifshitz transition in Sr2RuO4, npj Quantum Mater. 4, 46 (2019).
  * [45] M. E. Barber, A. S. Gibbs, Y. Maeno, A. P. Mackenzie, and C. W. Hicks, Resistivity in the Vicinity of a van Hove Singularity: Sr2RuO4 under Uniaxial Pressure, Phys. Rev. Lett. 120, 076602 (2018).
  * [46] F. Herman, J. Buhmann, M. H. Fischer, and M. Sigrist, Deviation from Fermi-liquid transport behavior in the vicinity of a Van Hove singularity, Phys. Rev. B 99, 184107 (2019).
  * [47] Y.-S. Li, M. Garst, J. Schmalian, S. Ghosh, N. Kikugawa, D. A. Sokolov, C. W. Hicks, F. Jerzembeck, M. S. Ikeda, Z. Hu, B. J. Ramshaw, A. W. Rost, M. Nicklas, and A. P. Mackenzie, Elastocaloric determination of the phase diagram of Sr2RuO4, Nature 607, 276 (2022).
  * [48] Yongkang Luo, A. Pustogow, P. Guzman, A. P. Dioguardi, S. M. Thomas, F. Ronning, N. Kikugawa, D. A. Sokolov, F. Jerzembeck, A. P. Mackenzie, C. W. Hicks, E. D. Bauer, I. I. Mazin, and S. E. Brown, Normal State 17O NMR Studies of Sr2RuO4 under Uniaxial Stress, Phys. Rev. X 9, 021044 (2019).
  * [49] E. A. Morales, G.-R. Siemann, A. Zivanovic, P. A. E. Murgatroyd, I. Markovic, B. Edwards, C. A. Hooley, D. A. Sokolov, N. Kikugawa, C. Cacho, M. D. Watson, T. K. Kim, C. W. Hicks, A. P. Mackenzie, and P. D. C. King, Hierarchy of Lifshitz Transitions in the Surface Electronic Structure of Sr2RuO4 under Uniaxial Compression, Phys. Rev. Lett. 130, 096401 (2023).
  * [50] B. Kim, S. Khmelevskyi, C. Franchini, and I. I. Mazin, Suppressed Fluctuations as the Origin of the Static Magnetic Order in Strained Sr2RuO4, Phys. Rev. Lett. 130, 026702 (2023).
  * [51] N. Kikugawa, A. P. Mackenzie, C. Bergemann, R. A. Borzi, S. A. Grigera, and Y. Maeno, Rigid-band shift of the Fermi level in the strongly correlated metal: Sr2-yLayRuO4, Phys. Rev. B 70, 060508(R) (2004).
  * [52] N. Kikugawa, A. P. Mackenzie, C. Bergemann, and Y. Maeno, Low-temperature Hall effect in substituted Sr2RuO4, Phys. Rev. B 70, 174501 (2004).
  * [53] N. Kikugawa, C. Bergemann, A. P. Mackenzie, and Y. Maeno, Band-selective modification of the magnetic fluctuations in Sr2RuO4: A study of substitution effects, Phys. Rev. B 70, 134520 (2004).
  * [54] A. Georges and J. Mravlje, Skewed non-Fermi liquids and the Seebeck effect, Phys. Rev. Research 3, 043132 (2021).
  * [55] N. Kouda, K. Eguchi, R. Okazaki, and M. Tamura, Anomalous scaling law for thermoelectric transport of two-dimension-confined electrons in an organic molecular system, Phys. Rev. Research 4, 043050 (2022).
  * [56] T. Yamanaka R. Okazaki, and H. Yaguchi, Enhanced Seebeck coefficient through magnetic fluctuations in ${\mathrm{Sr}}_{2}{\mathrm{Ru}}_{1-x}{M}_{x}{\mathrm{O}}_{4}$ ($M=$ Co, Mn), Phys. Rev. B 105, 184507 (2022).
  * [57] H. Yoshino, K. Murata, N. Shirakawa, Y. Nishihara, Y. Maeno, and T. Fujita, Thermopower of a Layered Perovskite Superconductor, Sr2RuO4, J. Phys. Soc. Jpn. 65, 1548 (1996).
  * [58] X. F. Xu, Z. A. Xu, T. J. Liu, D. Fobes, Z. Q. Mao, J. L. Luo, and Y. Liu, Band-Dependent Normal-State Coherence in Sr2RuO4: Evidence from Nernst Effect and Thermopower Measurements, Phys. Rev. Lett. 101, 057002 (2008).
  * [59] R. Daou, S. Hébert, G. Grissonnanche, E. Hassinger, L. Taillefer, H. Taniguchi, Y. Maeno, A. S. Gibbs, and A. P. Mackenzie, Anisotropic Seebeck coefficient of Sr2RuO4 in the incoherent regime, Phys. Rev. B 108, L121106 (2023).
  * [60] R. Otsuki, Y. J. Sato, R. Okazaki, T. Komine, R. Kurihara, and H. Yaguchi, Carrier filtering effect for enhanced thermopower in a body-centered tetragonal ruthenate, Phys. Rev. Mater. 7, 125401 (2023).
  * [61] J. Mravlje and A. Georges, Thermopower and Entropy: Lessons from Sr2RuO4. Phys. Rev. Lett. 117, 036401 (2016).
  * [62] G. K.H. Madsen, BoltzTraP. A code for calculating band-structure dependent quantities, Comput. Phys. Commun. 175, 67 (2006).
  * [63] M. Kawamura, FermiSurfer: Fermi-surface viewer providing multiple representation schemes, Comput. Phys. Commun. 239, 197 (2019).
  * [64] Y. M. Blanter, A.V. Pantsulaya, and A. A. Varlamov, Thermoelectric power and topological transitions in quasi-two-dimensional electronic systems, Phys. Rev. B 45, 6267 (1992).
  * [65] Y. Yamaji, T. Misawa, and M. Imada, Quantum and Topological Criticalities of Lifshitz Transition in Two-Dimensional Correlated Electron Systems, J. Phys. Soc. Jpn. 75, 094719 (2006).
  * [66] J. Lin, Lifshitz transition in two-dimensional spin-density wave models, Phys. Rev. B 82, 195110 (2010).
  * [67] J. M. Buhmann and M. Sigrist, Thermoelectric effect of correlated metals: Band-structure effects and the breakdown of Mott’s formula, Phys. Rev. B 88, 115128 (2013).
  * [68] A. A. Abrikosov, Fundamentals of the Theory of Metals (Dover, 1988), pages 111-114.
  * [69] K. Behnia, D. Jaccard, and J. Flouquet, On the thermoelectricity of correlated electrons in the zero-temperature limit, J. Phys.: Condens. Matter 16, 5187 (2004).
  * [70] A. Gourgout, G. Grissonnanche, F. Laliberté, A. Ataei, L. Chen, S. Verret, J.-S. Zhou, J. Mravlje, A. Georges, N. Doiron-Leyraud, and L. Taillefer, Seebeck Coefficient in a Cuprate Superconductor: Particle-Hole Asymmetry in the Strange Metal Phase and Fermi Surface Transformation in the Pseudogap Phase, Phys. Rev. X 12, 011037 (2022).
  * [71] S. Nakatsuji and Y. Maeno, Quasi-Two-Dimensional Mott Transition System Sr2-xCaxRuO4, Phys. Rev. Lett. 84, 2666 (2000).
  * [72] M. Minakata and Y. Maeno, Magnetic ordering in Sr2RuO4 induced by nonmagnetic impurities, Phys. Rev. B 63, 180504(R) (2001).
  * [73] J. P. Carlo, T. Goko, I. M. Gat-Malureanu, P. L. Russo, A. T. Savici, A. A. Aczel, G. J. MacDougall, J. A. Rodriguez, T. J. Williams, G. M. Luke, C. R. Wiebe, Y. Yoshida, S. Nakatsuji, Y. Maeno, T. Taniguchi, and Y. J. Uemura, New magnetic phase diagram of (Sr,Ca)2RuO4, Nat. Mater. 11, 323 (2012).
  * [74] J. E. Ortmann, J. Y. Liu, J. Hu, M. Zhu, J. Peng, M. Matsuda, X. Ke, and Z. Q. Mao, Competition Between Antiferromagnetism and Ferromagnetism in Sr2RuO4 Probed by Mn and Co Doping, Sci. Rep. 3, 2950 (2013).
  * [75] B. Zinkl and M. Sigrist, Impurity-induced magnetic ordering in Sr2RuO4, Phys. Rev. Research 3, 023067 (2021).
  * [76] N. Kikugawa and Y. Maeno, Non-Fermi-Liquid Behavior in Sr2RuO4 with Nonmagnetic Impurities, Phys. Rev. Lett. 89, 117001 (2002).
  * [77] A. Legros, S. Benhabib, W. Tabis, F. Laliberté, M. Dion, M. Lizaire, B. Vignolle, D. Vignolles, H. Raffy, Z. Z. Li, P. Auban-Senzier, N. Doiron-Leyraud, P. Fournier, D. Colson, L. Taillefer, and C. Proust, Universal $T$-linear resistivity and Planckian dissipation in overdoped cuprates, Nat. Phys. 15, 142 (2019).
  * [78] C. M. Varma, Colloquium: Linear in temperature resistivity and associated mysteries including high temperature superconductivity, Rev. Mod. Phys. 92, 031001 (2020).
  * [79] H. Zhao, J. Zhang, M. Lyu, S. Bachus, Y. Tokiwa, P. Gegenwart, S. Zhang, J. Cheng, Y.-f. Yang, G. Chen, Y. Isikawa, Q. Si, F. Steglich, and P. Sun, Quantum-critical phase from frustrated magnetism in a strongly correlated metal, Nat. Phys. 15, 1261 (2019).
  * [80] S. Paschen and Q. Si, Quantum phases driven by strong correlations, Nat. Rev. Phys. 3, 9 (2021).
  * [81] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, m. Cococcioni, I. Dabo, A. Dal Corso, S. Fabris, G. Fratesi, S. de Gironcoli, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter 21, 395502 (2009).
  * [82] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio Jr., A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter 29, 465901 (2017).
  * [83] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, Quantum ESPRESSO toward the exascale, J. Chem. Phys. 152, 154105 (2020).
  * [84] M. Kawamura, Y. Gohda, and S. Tsuneyuki, Improved tetrahedron method for the Brillouin-zone integration applicable to response functions, Phys. Rev. B 89, 094515 (2014).
  * [85] H.-L. Huang and H.-t. Jeng, Orbital ordering and magnetism in layered perovskite Ruthenate Sr2RuO4, Sci. Rep. 10, 7089 (2020).
  * [86] T. Oguchi, Electronic band structure of the superconductor Sr2RuO4, Phys. Rev. B 51, 1385 (1995).
  * [87] D. J. Singh, Relationship of Sr2RuO4 to the superconducting layered cuprates, Phys. Rev. B 52, 1358 (1995).
  * [88] M.W. Haverkort, I. S. Elfimov, L. H. Tjeng, G. A. Sawatzky, and A. Damascelli, Strong Spin-Orbit Coupling Effects on the Fermi Surface of Sr2RuO4 and Sr2RhO4, Phys. Rev. Lett. 101, 026406 (2008).
  * [89] M. Behrmann, C. Piefke, and F. Lechermann, Multiorbital physics in Fermi liquids prone to magnetic order, Phys. Rev. B 86, 045130 (2012).
  * [90] Y. J. Sato, A. Nakamura, R. Nishinakayama, R. Okazaki, H. Harima, and D. Aoki, Fermi surface topology and electronic transport properties of chiral crystal NbGe2 with strong electron-phonon interaction, Phys. Rev. B 108, 235115 (2023).
  * [91] V. S. Egorov and A. N. Fedorov, Thermopower of lithium-magnesium alloys at the 2 1/2-order transition, Zh. Eksp. Teor. Fiz. 85, 1647 (1983).


* (42) A. Joyal and R. Street “Braided Tensor Categories,” Advances in Mathematics 102, no.1, 20-78 (1993).
  * (43) J. Kaidi, Z. Komargodski, K. Ohmori, S. Seifnashri and S. H. Shao, “Higher central charges and topological boundaries in 2+1-dimensional TQFTs,” hep-th/2107.13091.
  * (44) A. Kapustin and N. Saulina, “Topological boundary conditions in abelian Chern-Simons theory,” Nucl. Phys. B 845, 393-435 (2011), hep-th/1008.0654.
  * (45) A. Kapustin and N. Saulina, “Surface operators in 3d Topological Field Theory and 2d Rational Conformal Field Theory,” hep-th/1012.0911.
  * (46) A. Kawauchi, S. Kojima, “Algebraic classification of linking pairings on 3-manifolds,” Math. Ann. 253, 29–42 (1980).
  * (47) M. V. Kneser, P. Dieter, “Quadratische Formen und Verschlingungsinvarianten von Knoten,” Mathematische Zeitschrift 58, no.1, 376-384 (1953).
  * (48) L. Kong, “Anyon condensation and tensor categories,” Nucl. Phys. B 886, 436-482 (2014), cond-mat.str-el/1307.8244.
  * (49) L. Kong and Z. H. Zhang, “An invitation to topological orders and category theory,” cond-mat.str-el/2205.05565.
  * (50) T. Kohno, T. Takata, “Level-rank duality of Witten’s 3-manifold invariants,” Progress in algebraic combinatorics. 24, 243-26 (1996).
  * (51) Y. Lee and Y. Tachikawa, “A study of time reversal symmetry of abelian anyons,” JHEP 07, 090 (2018), hep-th/1805.02738.
  * (52) M. Levin and A. Stern, “Classification and analysis of two-dimensional abelian fractional topological insulators.,” hysical Review B 86, 115-131 (2012), cond-mat.str-el/1205.1244.
  * (53) M. Mignard and P. Schauenburg, “Modular categories are not determined by their modular data,” Letters in Mathematical Physics 111, no.3, 1-9 (2017), math.QA/1708.02796.
  * (54) R. Miranda, “Nondegenerate symmetric bilinear forms on finite abelian 2-groups,” Transactions of the American Mathematical Society 284, 535-542 (1984).
  * (55) R. Miranda, “Embeddings of Integral Quadratic Forms.”
  * (56) G. W. Moore and N. Seiberg, “Classical and Quantum Conformal Field Theory,” Commun. Math. Phys. 123, 177 (1989).
  * (57) G. W. Moore and N. Seiberg, “LECTURES ON RCFT,” In _Strings ’89_ ,Proceedings of the Trieste Spring School on Superstrings, 3-14 April 1989, M. Green, et. al. Eds. World Scientific, 1990; Also in Lee, H.C. (eds) Physics, Geometry and Topology. NATO ASI Series, vol 238. Springer, Boston, MA.
  * (58) G. W. Moore and N. Seiberg, “Naturality in Conformal Field Theory,” Nucl. Phys. B 313, 16-40 (1989).
  * (59) M. Mueger, “From Subfactors to Categories and Topology II. The quantum double of tensor categories and subfactors,” J. Pure Appl. Alg. 180, 159-219 (2003), math.QA/0601012.
  * (60) S. G. Naculich, H. A. Riggs, and H. J. Schnitzer “Group-level duality in WZW models and Chern-Simons theory,” Physics Letters B 246, no. 3, 417-422 (1990).
  * (61) T. Neupert, L. Santoz, S. Ryu, C. Chamon, and C. Mudry “Fractional topological liquids with time-reversal symmetry and their lattice realization,” Physical Review B 84, no.16, (2011), cond-mat.str-el/1106.3989.
  * (62) T. Neupert, L. Santoz, S. Ryu, C. Chamon, and C. Mudry “Time-reversal symmetric hierarchy of fractional incompressible liquids,” Physical Review B 84, no.16, (2011), cond-mat.str-el/1108.2440.
  * (63) T. Neupert, L. Santoz, S. Ryu, C. Chamon, and C. Mudry “Wire deconstructionism of two-dimensional topological phases,” Physical Review B 90, no.20, (2014), cond-mat.str-el/1403.0953.
  * (64) S. H. Ng, P. Schauenburg, “Frobenius–Schur indicators and exponents of spherical categories,” Advances in Mathematics 70, no.1, 34–71 (2007), math.QA/0601012.
  * (65) S. H. Ng, A. Schopieray, Y. Wang, “Higher Gauss sums of modular categories,” Selecta Mathematica 25, no.4, 329-452 (2019), math.QA/1812.11234.
  * (66) S. H. Ng, E. Rowell, Y. Wang, Q. Zhang, “Higher central charges and Witt groups,” Advances in Mathematics 404, (2022), math.QA/2002.03570.
  * (67) V. V. Nikulin, “Integral symmetric bilinear forms and some of their applications,” Mathematics of the USSR-Izvestiya 14, (1980).
  * (68) V. Ostrik, M. Sun, “Level-rank duality via tensor categories,” Commun. Math. Phys. 326, 49–61 (2014), math-ph/1208.5131.
  * (69) E. Rowell, Y. Ruan, and Y. Wang “The Witt classes of $SO(2r)_{2r}$,” math.QA/2107.06746.
  * (70) E. Rowell, Z. Wang, “Mathematics of Topological Quantum Computing,” Bulletin of the American Mathematical Society 55, no. 2, 183–238 (2017), math.QA/1705.06206.
  * (71) S. Sanford, “Fusion Categories over Non-Algebraically Closed Fields”, PhD thesis.
  * (72) N. Seiberg and E. Witten, “Gapped Boundary Phases of Topological Insulators via Weak Coupling,” PTEP 2016, no.12, (2016), cond-mat.str-el/1602.04251.
  * (73) H. Seifert, “Verschlingungsinvarianten,” Sitzungsber. Preuß. Akad. Wiss. 70, no.26-29, 811-828 (1933).
  * (74) A. Sertöz, “Which singular K3 surfaces cover an Enriques surface,” Proceedings of the American Mathematical Society. 133, no.1, 43-50 (2005), math.AG/0205282.
  * (75) S. D. Stirling, “Abelian Chern-Simons theory with toral gauge group, modular tensor categories, and group categories,” hep-th/0807.2857.
  * (76) L. R. Taylor, “Gauss Sums in Algebra and Topology,” math.AT/2208.06319.
  * (77) C. T. C. Wall, “Quadratic forms on finite groups, and related topics,” Topology 2, 281-298 (1963).
  * (78) C. Wang and T. Senthil “Boson topological insulators: A window into highly entangled quantum phases,” Physical Review B, 87, no.23, (2013), cond-mat.str-el/1302.6234.
  * (79) L. Wang and Z. Wang, “In and around Abelian anyon models,” J. Phys. A 53 no.50, (2020), math.QA/2004.12048.
  * (80) E. P. Wigner, “Group Theory and its Application to the Quantum Mechanics of Atomic Spectra,” Academic Press Inc., New York, (1959).
  * (81) E. Witten, “Quantum field theory and the Jones polynomial,” Communications in Mathematical Physics 121, no.3, 351-399 (1989).
  * (82) E. Witten, “SL(2,Z) action on three-dimensional conformal field theories with Abelian symmetry,” hep-th/0307041.
  * (83) D. Yetter, “Framed tangles and a theorem of Deligne on braided deformations of Tannakian categories,” Deformation Theory and Quantum Groups with Applications to Mathematical Physics, 325–349 (1992)


# Muon radiography experiments on the subway overburden structure detection

Xin Mao<EMAIL_ADDRESS>Beijing Institute of Spacecraft Environment
Engineering, 104 YouYi Road, Beijing 100094, China  Zhiwei Li
<EMAIL_ADDRESS>Beijing Institute of Spacecraft Environment Engineering, 104
YouYi Road, Beijing 100094, China  Shuning Dong  Beijing Institute of
Spacecraft Environment Engineering, 104 YouYi Road, Beijing 100094, China
Jingtai Li  Beijing Institute of Spacecraft Environment Engineering, 104 YouYi
Road, Beijing 100094, China  Jianming Zhang  School of Nuclear Science &
Engineering, North China Electric Power University, Beijing 102206, China  Jie
Pang  School of Nuclear Science & Engineering, North China Electric Power
University, Beijing 102206, China  Yaping Chen  School of Nuclear Science &
Engineering, North China Electric Power University, Beijing 102206, China  Bin
Liu  College of Nuclear Science and Technology, Beijing Normal University, 19
Xinjiekou outer St., Beijing 13 100875, China  Xiaoping Ouyang  College of
Nuclear Science and Technology, Beijing Normal University, 19 Xinjiekou outer
St., Beijing 13 100875, China  Ran Han<EMAIL_ADDRESS>Beijing Institute
of Spacecraft Environment Engineering, 104 YouYi Road, Beijing 100094, China

###### Abstract

Muon radiography is an innovative and non-destructive technique for internal
density structure imaging, based on measuring the attenuation of cosmic-ray
muons after they penetrate the target. Due to the strong penetration ability
of muons, the detection range of muon radiography can reach the order of
hundreds of meters or even kilometers. Using a portable muon detector composed
of plastic scintillators and silicon photomultipliers, we performed a short-
duration(1h) flux scanning experiment of the overburden above the platform and
tunnel of the Xiaoying West Road subway station under construction. With the
observation direction facing up, the detector is placed on the north side of
the track and moved eastward from the platform section inside the station to
the tunnel section. The scanning length is 264m and a total of 21 locations
are observed. By comparing the observed and predicted values of the muon
survival ratio at different locations, the experiment accurately detects the
anomalous density loss caused by the lobby on the first floor of the platform
section. Furthermore, unknown anomalies caused by random placed light brick
piles and side passage entrance above the observation locations are detected
and confirmed later. This experiment verifies the feasibility of using natural
muons to quickly detect abnormal structures of the upper layer of tunnel, and
shows that muon radiography has broad application prospects in tunnel safety
and other similar aspects.

##  1 Introduction

By the end of 2021, China’s total subway mileage has ranked first in the world
with 7253.73 $km$, and the construction of municipal subway tunnels is also in
a period of rapid development[1], therefore, the safety of tunnel construction
and operation cannot be ignored. The compactness and stability of the rock and
soil above the tunnel are directly related to the safety of the subway
construction and operation. Quicksand, cavities or water bladders formed in
the overburden due to reasons such as leakage of urban underground pipelines
can put subway tunnels at risk of water inrush and collapse[2]. Between 2003
and 2016, there were 49 subway construction collapse accidents in China,
including 28 collapse accidents caused caused by pipeline leakage[3].
Reasonable internal inspection of the upper layer of the tunnel and timely
detection of the abnormal density structure can help avoid risks in the
construction and operation of the subway in advance.

Traditional geophysical methods, including electromagnetism, seismic and
gravity have their own advantages and limitations. The electromagnetic method
can realize the detection of geological bodies by obtaining the conductivity
structure of the underground medium, which has the advantages of large
detection depth and high efficiency, but is insensitive to the density
structure and is easily affected by large metal structures or mineralized
zones[4]. The seismic method has a large detection range and can obtain high-
resolution wave velocity structures, but the reliability of the density
structure obtained according to the wave velocity-density empirical
relationship is affected to a certain extent, and it needs to rely on natural
or artificial seismic sources, which will have a certain impact on the target
structure[5, 6]. The gravitational method has problems such as low
spatiotemporal resolution, insufficient vertical resolution, and strong multi-
solution[7]. In contrast, Muon radiography, also called muography in brief,
has the characteristics of high precision, large detection depth, non-
destructive by using of natural sources, and can obtain the internal density
structure of the detection target, which is more suitable for urban subway
tunnels, ancient buildings and other large targets that are not suitable for
intrusion and destruction.

Muons are secondary particles produced by the interaction of high-energy
cosmic rays from the universe, mainly protons, with atomic nuclei in the
atmosphere, with an average lifetime of 2 microseconds and traveling at the
speed of light. The muons reaching sea level have an energy of more than 1
$GeV/c$ and a flux of about $170Hz/m^{2}$[8]. On the contrary, the number of
other particles such as protons, mesons and electrons reaching the surface can
be ignored. Muons will undergo a process of energy loss and deflection of the
incident direction during the penetration of the target, due to inelastic
collisions with electrons and elastic scattering from nuclei[9]. High-energy
muons can penetrate hundreds or even thousands of meters of rock
formations[10] while low-energy muons will be stopped. The density structure
inside the target can be indirectly obtained by measuring the attenuation of
the muon flux before and after passing through the target object, which is the
basic principle of muography.

Research on muography began in 1950 and was first applied in 1955 to measure
overburden thickness above underground laboratories[11]. With the improvement
of detector level and electronic technology, muography has been successfully
applied to geology and archaeology[12, 13, 14, 15, 16]. Muography is more
suitable for detecting internal chambers in large buildings or underground
structures. In 2017, an unknown cave is detected by performing a muography
experiment 40 meters underground in Mt. Echia, Naples[17]. In the same year,
an unknown space hidden inside the Pyramid of Khufu is discovered by the same
technique[16]. In 2020, a muon flux observation experiment of the 770m-long
Alfreton tunnel in the United Kingdom was conducted and an unknown overburden
void was identified, which was subsequently officially confirmed[18].

In order to verify the feasibility of muography in detecting structural
anomalies of the overburden of the tunnel and explore the potential of the
muography in engineering applications, we conducted a muon observation
experiment of the overburden of the tunnel at the Xiaoying West Road subway
station under construction. We analyzed the abnormalities in the overburden
layer, and derived the average density distribution and average thickness
distribution. In this paper, after an introduction on the principles of
determining anomalies, the detection system and experimental process will be
presented, and finally the experimental results will be analyzed and
summarized.

##  2 Methods

###  2.1 Determination method of anomalies

The degree of attenuation of muons through the substance is defined as the
ratio of the flux reaching the detector after penetrating the overburden
(underground case) to the amount of flux reaching the detector before
penetrating the overburden (open-sky case), here referred to as
survival/transmission ratio, expressed by ratio, abbreviated as R. The
principle for judging the existence of anomalies is that the observation value
of ratio is deviated from the prediction value of ratio without abnormal
structure, also the deviation is outside the statistical error range of the
data.

When there is no anomalous structure inside the object, the number $N$ of
muons from a specific direction $\theta_{0}$ detected by the detector in
$\Delta T$ time can be calculated by the following formula:

$N(\theta_{0})=\tau(\theta_{0})\cdot\Delta{T}\cdot\varepsilon\cdot\int^{\infty}_{E_{min}(\rho)}\phi(\theta_{0},E)dE$
(1)

Among them, $\tau(\theta_{0})$ is the geometric acceptance of the muon
detector, representing the influence of the detector geometry on the detection
performance, the unit is $cm^{2}\cdot{sr}$, depending on the number of
detector pixels $N_{x}\cdot{N_{y}}$, pixel size $d$, and the separation
distance $D$ of the top and bottom layers of the detector. $\tau(\theta_{0})$
can be obtained by integrating the detection area with solid angular
resolution[13, 14]. The geometric acceptance in different directions is
different and is the largest in the vertical direction. Futhermore, the
greater the deviation from the vertical direction, the smaller the acceptance.
$\Delta{T}$ is the detection time, in second. The $\varepsilon$ is the
detection efficiency, which can be obtained by scaling the efficiency of each
scintillator unit[14]. $\phi(\theta_{0},E)$ is the differential energy
spectrum of the incident muon in
$cm^{-2}\cdot{sr^{-1}}\cdot{s^{-1}}\cdot{GeV^{-1}}$, and when predicting the
muon flux, it is usually approximated by using empirical formulas obtained by
fitting the muon data measured at sea level, such as Gaisser’s formula
[27][28], or directly using some software packages that simulate the
interaction of primary cosmic rays with the atmosphere, such as Sybill[19],
CRY[20], FLUKA[21], etc. $E_{min}$ is the minimum energy required for the muon
to penetrate a substance of a certain thickness, depending on the density
length $\rho$ of the penetrating object. $\rho$ is defined as the integral of
the density of matter penetrated by the muon along the path, or can be
approximated by multiplying the average density along the path by the length
of the path. The relationship between the $E_{min}$ and $\rho$ in the energy
range of $10MeV\sim 100TeV$ in different common substances can be obtained by
looking up the summarized tables[22], or by simulating with Monte Carlo
simulation tools such as Geant4.

In the absence of an anomalous structure inside the object, within the
$\theta$ cone angle, the predicted value of the muon survival ratio can be
expressed by the following formula.

$R(\theta)=\frac{\int_{\theta}N(\theta_{0})d\theta_{0}}{\int_{\theta}N_{opensky}(\theta_{0})d\theta_{0}}=\frac{\int_{E_{min}(\rho)}\phi(\theta_{0},E)dE}{\int_{E_{0}}\phi(\theta_{0},E)dE}$
(2)

Among them, $E_{0}$ is the minimum energy threshold of the muon that the
detector can detect in the open sky case. It can be seen from Equation (2)
that using the survival ratio to judge the anomaly of the density structure
can eliminate the influence of the detector effect on imaging.

###  2.2 The relationship between the ratio and density length

According to the definition of density length, if its relationship to the
ratio is known, the average value of either quantity in density and length can
be derived by assuming in advance that the other quantity is constant. To get
the relationship, the CRY software package is used to generate cosmic ray
muons, and Geant4 is used to simulate the interaction between muons and
matter, and finally the relationship between ratio and penetration density
length is obtained. CRY is short for Cosmic-ray shower generator, and the CRY
software package can generate the distribution of cosmic ray particle beams at
three altitudes, that is sea level, 2100m and 11300m, as a simulated input to
the detector response. Geant4 is a Monte Carlo simulation tool commonly used
to simulate the transportation and interaction of particles in matter, and is
widely used in high-energy physics, nuclear physics, medicine, aerospace and
other fields. The main process of this method is to establish the CAD model of
the target, import it into Geant4, use the muon source generated by CRY,
establish the detector in Geant4, select the built-in physical process in
Geant4 for simulation, and finally obtain the muon flux rate received by the
detector. Based on this method, the penetration of muons through different
thicknesses of concrete of which the average density is assumed to be
$2.3g/cm^{3}$ is simulated, and the ratios of muons received in different
angle ranges are obtained as shown in Figure 1. The solid black line in the
figure represents the relationship between the ratio of the muons in the range
of 50° cone angle and the density length of the penetrated object, and the
remaining dashed lines represent the relative deviation of the ratio value of
muons in the range of other cone angles to the result of 50°, which is
represented by the following formula.

$\Delta(\theta)=\frac{|R(\theta)-R(50^{\circ})|}{R(50^{\circ})}$ (3)

Figure 1: The relationship between ratio and density length. (solid black
line: ratio of muons within $50^{\circ}$ cone angle; Blue, pink, yellow,
green, red dotted line: $5^{\circ}$, $10^{\circ}$, $20^{\circ}$, $30^{\circ}$,
$40^{\circ}$ relative deviation from the case of
$50^{\circ},{\Delta(\theta)=\frac{|R(\theta)-R(50^{\circ})|}{R(50^{\circ})}}$.

##  3 Experiment and results

###  3.1 Detection system

The detector, based on plastic scintillators and silicon photomultiplier
devices (SIPM), is an upgraded version of the previous generation
[20][19][21]. The main improvements are higher position resolution, higher
detection efficiency, smaller size and less power consumption, and more
portable and convenient for moving and observing in complex scenes. As shown
in Figure 2, the detector consists of four layers of detection unit plates
arranged equidistant vertically to form $Z$ coordinates, and the plate spacing
is $20cm$. Each layer consists of two $2500cm^{2}$ planes, each consisting of
25 plastic scintillator arrangements, each coupled with a single-ended SIPM.
The scintillators of the two planes are orthogonally arranged together to form
$X$ and $Y$ coordinates. The geometric size of the detector determines that
its angular aperture of about $79^{\circ}$, and the maximum acceptance is
$2.78cm^{2}sr$. The bottom of the detector is mounted with pulleys for
convenient movement. In addition, the detector consumes $20W$ and is equipped
with a $2000Wh$ mobile lithium battery, thus can continuously operate in
places without mains power.

Figure 2: Schematic diagram of the structure of the detector.

The most basic unit of the detector is the strip plastic scintillator with a
rectangular cross-section. When the muon hits the scintillator, there will be
ionization, the atomic deexcitation will produce a light signal, and the light
signal will be detected by SIPM when it propagates along the scintillator to
the side. The scintillator density used is $1.02g/cm^{3}$, the light output
relative to anthracene crystals is $50\%\sim 60\%$, the emission wavelength is
between $395nm\sim 425nm$, the attenuation time is $2.4ns$, the attenuation
length is $2.1m$, and the measurement range is between $50keV\sim{3MeV}$. Each
scintillator is $2cm\times 1cm\times 50cm$ in size, resulting a high average
angular resolution of $33.3mrad$, and is coated with a silver reflective film
(ESR), which can increase the reflectivity of photons inside and shield
foreign light sources to improve the efficiency of photon collection. At the
scintillator end face, the J series SIPM sensor produced by onsemi is used to
detect the optical signal, the size is $6mm\times 6mm$, and it has a high
quantum efficiency (PDE) for weak light, the maximum is $50\%$, and the
corresponding wavelength is $420nm$, which corresponds to the maximum emission
wavelength of the scintillator.

The front-end and back-end electronics (FEE) is used to process and upload the
electrical signal, converted by SIPM from the detected optical signal, in
16-base form via TCP/IP protocol. Finally the data is analyzed and stored by
the host computer software. The FEE consists of 8 amplification screening
modules and 1 data acquisition board. Each amplification screening module can
perform year-over-year amplification and over-threshold screening of 25
signals, and if the threshold of $32mV$ is exceeded after triggering signal
amplification, it is considered "fire", and the over-threshold screening
circuit then outputs a $3.3V$ TTL level signal to the detector data
acquisition board. The acquisition board is based on the Xilinx ZYNQ7035
series FPGA (Field Programable Gate Array), which can simultaneously read a
total of 200 level signals from four layers of detection units. In addition,
the board integrates multi-function peripherals such as temperature control
voltage module, FPGA start-stop control module, and data transmission
interface module, which can provide SiPM bias and temperature compensation for
the detector at the same time. The sampling frequency of the acquisition card
is $100MHZ$. When the optical signal is transmitted from the hitting point
position to the SIPM detection, there is a delay in time. In order to prevent
the signal acquisition is incomplete, the signal collected in 20 clock cycles
is taken by the same logical bit "or" into a "real signal". When the real
signals of X and Y "fire" at the same time, the acquisition board will package
the detection unit number, temperature, signal into a 16-base data of 19 byte
and upload it to the host computer through the TCP/IP protocol. Finally, the
host computer selects, displays and saves the data of the muon events. Among
them, the case can be considered as muon case only if the hitting time of the
four detection units are in the same window ($200ns$) and the hitting points
of the four detection units are in the same straight line. The detection
framework of the system is shown Figure 3.

Figure 3: Block diagram of the detection system. Figure 4: 1:1 CAD model of
the Xiaoying West Road subway station and tunnel for Geant4 simulation.

###  3.2 Experiment

The experimental site is located at Xiaoying West Road subway station in
Beijing. The station has two reverse lines, the left line is on the north side
and heads from west to east, the right line is on the south side and heads
from east to west. This experiment is located on the left line. The subway
station and tunnel of the left line are drawn at the actual 1:1 scale CAD
model shown in Figure 4. The subway track extends from west to east through
the subway station to the tunnel and is represented in pink. The platform
inside the station is $165m$ long (platform section), and the tunnel outside
the station is about $1.17km$ (tunnel section) east to the next station.
Inside the station is a two-story hall, the top of the hall is an arch with a
radius of $7.6m$ (dark blue in the picture), the vault is $0.8m$ thick, the
thickness of the overburden between the vault and the ground is $8.7m$,
between the two floors of the hall is a $0.4m$ thick concrete floor, and the
hall is arranged with 26 marble columns with a diameter of $0.8m$ along the
east-west platform. On the north edge of the hall is a $0.8m$ thick concrete
wall (brown in the picture). Outside of the hall, the shield of the tunnel
part has an inner diameter of $2.7m$ and a lining thickness of $0.3m$, which
is represented by a purple hollow cylinder in the figure, and the thickness of
the overburden between the top of the lining and the surface is $17m$. A 1:1
size detector is placed in the model, and the red cone in the figure is the
detection range of the detector. It can be seen from the figure that the
observation range of the two observation points has a large combination,
indicating that the entire line can be scanned at this observation distance
interval. In addition, since the subway line is located directly below
Xiaoying West Road, there is no influence of pavement buildings above the
observation sites, which is conducive to simulation verification.

In order to eliminate the error of the detector itself and measure the
stability of the detector, it is necessary to conduct a ground/open-sky
experiment at first, as shown by equation (2). Therefore, the detector was
placed skyward on the ground and measured for 114 hours before entering the
tunnel. After this, the underground observation started. The detector panel is
placed horizontally, facing the upper overlay. In order not to affect the
engineering vehicles moving on the track, the detector was placed on the north
side of the track against the wall, and the detector center is about $0.6m$
from the track. Started at the position of the marble column No. 11 in the
middle of the platform (the first red cone point in Figure 4), the detector
was moved from west to east for point-by-point scanning, with each point
observed for 1 hour. A total of 21 points were observed, the cumulative time
was 210 hours and the scanning distance was $264m$. Among them, each point is
observed every $9.56m$ in the platform section, a total of 11 points are
observed, and the scanning distance is $90.72m$. It is worth noting that
observation point 11 is located directly below the interface between the
eastern edge of the hall and the tunnel section. A total of 10 points were
observed in the tunnel section, with the first 7 points spaced $12m$ apart and
the last 3 points spaced $30m$ apart. Figure 5 shows the observation scenes of
the detector in the platform section and tunnel section.

(a)

(b)

Figure 5: Detection scenes in the platform section (a) and tunnel section (b).

##  4 Results

###  4.1 Open-sky muon flux measurement

The open-sky flux measurement result before entering the subway tunnel is
shown in Figure 6. The histogram shows the statistical distribution and the
Gaussian fitting result of muon rates obtained from 114-hour measurement. It
can be seen that the observed muon rate under the open-sky case is $15670\pm
264h^{-1}$, and the data stability of the detector is $1.69\%$.

Figure 6: Open-sky histogram distribution and statistical fitting result of
the muon rates by 114-hour measurement.

###  4.2 Underground muon flux measurement

Muon event rate measurements at 21 locations within the platform section and
tunnel section are shown in Figure 7, of which the error bars represent the
statistical errors and the subway tunnel model profile is drawn at the bottom
of the diagram according to the coordinate relationship. It can be seen that
the average muon rates are $3795.4\pm 62.68h^{-1}$ in the platform section,
$3027\pm 55h^{-1}$ at the intersection face and $2137.3\pm 105.96h^{-1}$ in
the tunnel section.

Figure 7: Underground measurement results of muon rates and corresponding
statistical errors as a function of distance. Figure 8: The observed and
predicted ratios vary with distance (a); Assuming that the structure size is
fixed, the average density of the 10-degree-cone overburden varies with
distance (b).Assuming an average density of $2.3g/cm^{3}$, the average
thickness of the 10-degree-cone overburden varies with distance (c).

###  4.3 Ratio, density and thickness

The observed and predicted values of the ratios at 21 locations within the
platform section and tunnel section are distributed as a function of distance,
as shown in Figure 8(a). The data points in different colors are measured
ratios within different observation angle ranges, and the corresponding color
ramp range is the statistical error of the data. The dotted lines
corresponding to the colors in the figure are the predicted ratios obtained
from the Geant4 simulation using the CAD model in Figure 4. It can be seen
from the figure that there is a big difference in the ratio points within the
platform section and the tunnel section in the figure, point 11 is located at
the intersection interface and its ratio is in the middle position. In
addition, the ratio value of the tunnel section has a tendency to increase.
Secondly, by comparing with the predicted values, points 7 and 8 deviate
greatly from the simulated data, indicating that there are anomalies in the
corresponding observation range, resulting in a low observed number at
location 7 and a large number at location 8.

According to the relationship between ratio and density length in Figure 1,
since the structural dimensions of the subway station and tunnel are known,
assuming that the thickness of the platform and tunnel part of the overburden
remains constant, using the observation results in the near vertical direction
of 10 degrees, the average density distribution of the overburden above the
detection location is inverted as a function of distance, as shown in Figure
8(b).

Similarly, assuming that the average density for the overburden above the
detector is fixed at $2.3g/cm^{3}$, the average value of the thickness of the
overburden can be calculated according to the density length. The average
thickness distribution of the 10 degree cone angle overburden above the
detector is obtained as shown in Figure 8(c). It can be seen that the
thickness of overburden within the platform section is relatively stable, with
an average thickness of $11.11m$. The average thickness of the overburden
within the tunnel section is $18.5m$, and there is a trend of gradual
thinning.

##  5 Discussion

###  5.1 Anomaly Analysis

As can be seen from Figure 8(a), the change in the muon ratio detected by the
detector is highly sensitive to anomalies. Due to the characteristics of large
geometric acceptance of detector in the vertical direction but small in the
inclined direction, the muons with large inclination angle would have fewer
event statistics and large statistical errors. In order to reduce the
influence of the muon with a large angle tilt, the muon events in the range of
10 degrees, 20 degrees, and all observation angles were selected and analyzed
separately, with muons in the 10-degree cone angle range was approximately
considered to be vertical. According to the observation results in Figure
8(a), the major anomaly is the difference between the platform section and the
tunnel, where the ratio values at observation locations in the two sections in
the figure are quite different. Secondary anomalies within the platform
section and tunnel section can be analyzed separately. In the platform
section, under the full observation range of the detector, the observed values
of points 7 and 8 deviate from the predicted values and are outside the error
range. Within the tunnel section, the observed values deviate from the
predicted values and the amount of deviation gradually increase with distance.
Here, we mainly discuss these three part of anomalies, other abnormalities
caused by tiny structures in the subway station overburden or other reasons
will not be discussed. Such as point 4, the ratio observation value within the
full observation angle range is consistent with the predicted value. However,
when the angle range decreases, there is an anomaly. It is guessed that there
is a density loss caused by structure such as cable channels directly above
this location. Due to its small size, it will not be discussed in the article.

For the major anomaly, the difference in the ratios between the platform
section and the tunnel is very significant, nearly doubling. The interface
ratio value is somewhere in between. Since the experimental observations and
simulated predictions in the Figure 8(a) are basically consistent within the
statistical error range, it can be seen from the structure of the subway model
that the anomaly is caused by the density loss of the lobby, around $8m$ in
height, on the first floor of the station, which is known from the
architectural drawings.

For the secondary anomalies in the platform section, according to Figure 8(a),
there are two anomalies by the full-range observation, location 7 and 8. As
the observation angle range decreases, the anomaly at location 7 has always
existed, while the anomaly at location 8 has gradually disappeared, indicating
that the anomaly at location 7 may be located directly above the detector,
while the anomaly at location 8 is above the side of the detector. In the
lobby on the first floor within the platform section, it was observed that
there was a pile of light bricks above location 7, and a passage mouth with a
depth of about $2m$ on the wall side above above location 8, as shown in
Figure 9(a). Consulting the engineering drawings of the subway station shows
that the passage entrance is the entrance of the subway station passage B,
which is currently in the stage of unexcavated. The lightweight brick pile
makes the muon flux at location 7 small, while the passage mouth makes the
muon flux at location 8 large, resulting in abnormalities. In order to confirm
the guess of the anomaly at locations 7 and 8, the light brick pile and the
passage mouth were added to the original subway tunnel model in Figure 4 for
simulation, as shown in Figure 9(b). The light body brick pile was $1.6m$
high, $14m$ long, $2.8m$ wide, and the density was set to $1.7g/cm^{3}$. The
entrance of passage B is $7m$ wide and $3.6m$ high, and is topped by a round
arch with a radius of $6.5m$. Through the Geant4 simulation, the prediction
value of ratio in the full observation angle range is shown by the dotted line
in figure 10. It can be seen that the predicted values of the anomalies at
locations 7 and 8 are consistent with the actual observed values.

(a)

(b)

Figure 9: Light body brick pile and passage mouth above locations 7 and 8 (a).
Model diagram of the light body brick pile and passage mouth (b). Figure 10:
Simulation results for the anomalies at locations 7 and 8 in the platform
section.

For secondary anomalies in the tunnel section, it can be seen from Figure 8(a)
that from the position of $174m$, the observed ratio within the tunnel section
gradually becomes larger and deviates from the predicted value. It is
speculated that the overburden above the tunnel section is gradually thinning.
Due to the lack of actual data on the thickness of the overburden in the
tunnel section, further observation of the muons at longer time and denser
points are required to confirm this hypothesis.

###  5.2 Error source analysis

For the observation values, the detector itself has good stability and the
cosmic ray muon flux fluctuates little [23] at the same altitude, under the
same pressure and temperature conditions. Moreover, for the muon radiography
of large scale target, the effect of multiple coulomb scattering is negligible
[24]. Therefore, the system error of detector is neglected [25] and only
statistical error is considered in the paper. The statistical error of the
muon rate is estimated using $\sigma=1/\sqrt{N}$, and the ratio error is
calculated using the error transfer formula of division, as shown in the error
band of Figure 8(a). For the inverted average density (Figure 8(b)) and
average thickness (Figure 8(c)) of the overburden, the main error comes from
the inaccuracy of the structural CAD model as well as the assumed density
value. Since the subway station in the actual situation is much more complex
than the model, there are many structures in the overburden that are
relatively small compared to the light brick pile and the passage mouth, such
as ventilation ducts, cable sandwiches, etc., which will cause the thickness
and density of the inversion to deviate from the actual value more or less. In
addition, since the relationship between ratio and density length is used to
derive the average density and average thickness (Figure 1), the error caused
by this conversion relationship, whether estimated using empirical formulas or
simulated by software packages, cannot be ignored. Therefore, the anomaly
analysis is mainly carried out by comparing the difference between the
observed value and the predicted value of the ratio, and the average thickness
and average density distribution of the conversion are only used as reference.

##  6 Conclusion

Muography is an innovative way to study the internal structure of large
buildings non-intrusively, non-destructively and conveniently. Nowadays,
muography has many successful research cases in volcanoes, ancient building
archaeology, nuclear safety and other aspects at both home and abroad. The
main principle of muography is to use the muon detector to observe the change
of muon rate before and after muons penetrate the object, and reflect the
density or thickness inside the target according to the relationship between
the survival ratio of the muon rate and density length of the matter. For
subway tunnels, structural anomalies such as cavities in the overburden are
safety factors that have to be considered. Compared with traditional detection
methods, muography has no damage and deep detection thickness, which is more
suitable for the detection of the overburden of urban subway tunnels. In this
paper, a muon detector composed of four layers of unit plates sized
$50cm\times 50cm$ based on plastic scintillator and SIPM readout is developed,
with a pixel size of $2cm$. The detector has a broad angular aperture of about
$79^{\circ}$, a high average angular resolution of $33.3mrad$, and a maximum
acceptance of $2.78cm^{2}sr$. The detector is equipped with a $2000Wh$ mobile
power supply and pulleys at the bottom, which is convenient for moving and
continuous acquisition in the tunnel. Based on this detector, this paper
conducts a muon rate scanning experiment for the overburden of Xiaoying West
Road railway station under construction in Beijing, and analyzes the main
structural anomalies.

The first measurements were made on the ground before the underground
experiment. 114 hours of data have been obtained in the case of open sky and
the data stability was $1.69\%$. The underground experiment can be divided
into two parts, the platform section where there is a two-story hall inside
the station and the tunnel section outside the station. The detector scanned
the overburden along the track eastward from the mid-position of the platform
section, observing at each location for 1 hour. A total of 21 points were
scanned, with a total length of about $264m$.

By comparing the observed and predicted values of the ratios, the detected
anomalies are divided into major anomaly, secondary anomalies in the platform
section, and secondary anomalies in the tunnel section. The major anomaly was
caused by the lobby on the first floor of the subway station in the platform
section, resulting in the ratio values within the platform section being
nearly double those within the tunnel section. As for the secondary anomalies
in the platform section, the main ones are at locations 7 and 8 where the
ratios deviate from the predicted values within the full observation angle
range. The ratio value at location 7 is low and does not change with the
narrowing of the observation angle range. The ratio value at location 8 is
large, and the difference gradually decreases as the observation angle range
decreases. Through observation, the speculative anomalies were found in the
lobby on the first floor, which were the pile of light bricks directly above
location 7 and the passage mouth above the side of location 8. In order to
confirm this conjecture, two structures, the light brick pile and the passage
mouth, were added to the tunnel model according to their actual sizes. The
Geant4 simulation results were consistent with the observed values within the
error range. For the tunnel section, the ratio observation deviates from the
predicted value from the location at $174m$ and gradually becomes larger, it
is guessed that the overburden of the tunnel in this section gradually becomes
thinner. If the anomaly is to be further analyzed, further long-term, denser-
point observation is required.

The systematic error is independent of the statistic and is relatively small
and negligible. For the muon rate and ratio observations, only the impact of
the statistical error is considered, and the estimation is made using
$\sigma=1/\sqrt{N}$. For the inverted average density and average thickness of
the overburden, the main error comes from the inaccuracy of the model size and
the assumed density. Secondly, the error from the conversion relationship
between ratio and density length cannot be ignored either.

In general, through this muon observation experiment, not only the density
loss caused by the lobby located on the first floor of the platform section
was accurately detected, but also secondary anomalies such as light brick
piles and unexcavated passage mouth were detected. It shows that the muon
event rate is very sensitive to the change of the structure of the overburden
layer, and muography may be a very promising technology for tunnel safety
detection.

##  7 Acknowledge

This work is supported by National Natural Science Foundation of China (Grant
No. 42174076, 41974064 and U1865206), Xi’an Research Institute of China Coal
Technology & Engineering Foundation (Grant No. 2021XAYJC01). We also express
our gratitude to YanHeng Li for his great recommendations of experiment site.

## References

  * [1] Editorial Department of China Journal of Highway and Transport.  Review on china’s traffic tunnel engineering research: 2022.  China Journal of Highway and Transport, 35(4):1–40, 2022.
  * [2] Chengping Zhang, Yuejing Yue, and Yi Cai.  Influence of pipeline leakage range on ground deformation and failure duration shallow tunneling.  Chinese Journal of Rock Mechanics and Engineering, 34(02):392–400, 2015.
  * [3] Gong Wang.  Study on Accidents Statistical Analysis and Risk Assessment Methods of Urban Subway Tunnel.  PhD thesis, Beijing Jiaotong University, Beijing, 2018.  https://cdmd.cnki.com.cn/Article/CDMD-10004-1018144591.htm.
  * [4] Thomas Weiland.  Time domain electromagnetic field computation with finite difference methods.  International Journal of Numerical Modelling: Electronic Networks, 9(4):295–319, 1996.
  * [5] Nikolai M. Shapiro, Michel Campillo, Laurent Stehly, and Michael H. Ritzwoller.  High-resolution surface-wave tomography from ambient seismic noise.  Science, 307(5715):1615–1618, 2005.
  * [6] Chenglong Duan, Changhong Yan, Baotian Xu, and Yinkang Zhou.  Crosshole seismic ct data field experiments and interpretation for karst caves in deep foundations.  Engineering Geology, page S0013795217300121, 2017.
  * [7] Anne Barnoud, Valérie Cayol, Valentin Niess, Cristina Cârloganu, Peter Lelièvre, Philippe Labazuy, and Eve LeMénédeu.  Bayesian joint muographic and gravimetric inversion applied to volcanoes.  Geophysical Journal International, (3):2179–2194, 2019.
  * [8] M. Tanabashi, K. Nakamura, K. Hikasa, and K. Hagiwara.  Review of particle physics.  Journal of Physics G: Nuclear and Particle Physics, (10):88, 2016.
  * [9] R. K. Adair and H. Kasha.  Cosmic-ray muons.  Muon Physics, pages 323–385, 1977.
  * [10] O.C. Allkofer, G. Bella, W.D. Dau, H. Jokisch, G. Klemke, Y. Oren, and R. Uhr.  Cosmic ray muon spectra at sea-level up to 10 tev.  Nuclear Physics B, 1985.
  * [11] E.P. George.  Cosmic rays measure overburden of tunnel.  Commonwealth Engineer, pages 455–457, 01 1955.
  * [12] H. K. M. Tanaka, T. Nakano, S. Takahashi, J. Yoshida, M. Takeo, J. Oikawa, T. Ohminato, Y. Aoki, E. Koyama, and H. and Tsuji.  Radiographic imaging below a volcanic crater floor with cosmic-ray muons.  American Journal of Science, 308(7):843–850, 2008.
  * [13] Ya Ping Cheng, Ran Han, Zhi Wei Li, Jing Tai Li, Xin Mao, Wen Qiang Dou, Xin Zhuo Feng, Xiao Ping Ou-Yang, Bin Liao, and Fang and Liu.  Imaging internal density structure of the laoheishan volcanic cone with cosmic ray muon radiography.  Nuclear Science and Techniques, 33(7):1–10, 2022.
  * [14] R. Han, Q. Yu, Z. Li, J. Li, Y. Cheng, B. Liao, L. Jiang, S. Ni, Z. Yi, T. Liu, and Z. Wang.  Cosmic muon flux measurement and tunnel overburden structure imaging.  Journal of Instrumentation, 15(06):P06019, jun 2020.
  * [15] R. Nishiyama, A. Ariga, T. Ariga, A. Lechmann, and Fritz Schlunegger.  Bedrock sculpting under an active alpine glacier revealed from cosmic- ray muon radiography.  Scientific Reports, 9(6970), 2019.
  * [16] Kunihiro Morishima, Mitsuaki Kuno, Akira Nishio, Nobuko Kitagawa, Yuta Manabe, Masaki Moto, Fumihiko Takasaki, Hirofumi Fujii, Kotaro Satoh, and Hideyo and Kodama.  Discovery of a big void in khufu’s pyramid by observation of cosmic-ray muons.  Nature, 2017\.
  * [17] G Saracino, L Amato, F Ambrosino, G Antonucci, L Bonechi, L Cimmino, L Consiglio, R. D. ’ Alessandro, E. De Luzio, and G and Minin.  Imaging of underground cavities with cosmic-ray muons from observations at mt. echia (naples).  Scientific Reports, 7(1):1181, 2017.
  * [18] L. F. Thompson, J. P. Stowell, S. J. Fargher, C. A. Steer, K. L. Loughney, E. M. O’Sullivan, J. G. Gluyas, S. W. Blaney, and R. J. Pidcock.  Muon tomography for railway tunnel imaging.  Phys. Rev. Res., 2:023017, Apr 2020.
  * [19] Eun-Joo Ahn, Ralph Engel, Thomas K. Gaisser, Paolo Lipari, and Todor Stanev.  Cosmic ray interaction event generator sibyll 2.1.  Phys. Rev. D, 80:094003, Nov 2009.
  * [20] Christian Hagmann, David Lange, and David M. Wright.  Cosmic-ray shower generator (cry) for monte carlo transport codes.  2007 IEEE Nuclear Science Symposium Conference Record, 2:1143–1146, 2007.
  * [21] G. Battistoni, A. Ferrari, S. Muraro, and P.R. Sala.  Atmospheric muon simulation using the fluka mc model.  Nuclear Physics B - Proceedings Supplements, 168:286–288, 2007.  Proceedings of the Neutrino Oscillation Workshop.
  * [22] Donald E. Groom, Nikolai V. Mokhov, and Sergei I. Striganov.  Muon Stopping Power and Range Tables 10 MeV-100 TeV.  Atomic Data and Nuclear Data Tables, 78(2):183–356, July 2001.
  * [23] N. Lesparre, D. Gibert, J. Marteau, Y. Déclais, D. Carbone, and E. Galichet.  Geophysical muon imaging: feasibility and limits.  Geophysical Journal International, 183(3):1348–1361, 2010.
  * [24] Jianming Zhang, Zhiwei Li, Fang Liu, Jintai Li, Xin Mao, Yaping Cheng, Jie Pang, Xinzhuo Feng, Sidao Ni, Xiaoping Ouyang, and Ran Han.  Influence of multiple coulomb scattering on accuracy of muon transmission imaging of small-scale matter.  Acta Physica Sinica, 72(2):021401, may 2023.
  * [25] J. Li, Z. Li, R. Han, Y. Cheng, X. Mao, L. Yu, X. Feng, B. Liu, L. Jiang, and X. Ouyang.  Investigation of structures in tunnel overburdens by means of muon radiography.  Journal of Instrumentation, 17(05):P05029, may 2022.