# Single crystal growth and superconductivity in RbNi2Se2

Hui Liu  Xunwu Hu  Center for Neutron Science and Technology, Guangdong
Provincial Key Laboratory of Magnetoelectric Physics and Devices, School of
Physics, Sun Yat-Sen University, Guangzhou, 510275, China  Hanjie Guo  Neutron
Science Platform, Songshan Lake Materials Laboratory, Dongguan, Guangdong,
523808, China  Xiao-Kun Teng  Department of Physics and Astronomy, Rice
University, Houston, TX 77005, USA  Huanpeng Bu  Neutron Science Platform,
Songshan Lake Materials Laboratory, Dongguan, Guangdong, 523808, China  Zhihui
Luo  Lisi Li  Zengjia Liu  Mengwu Huo  Feixiang Liang  Hualei Sun  Bing Shen
Center for Neutron Science and Technology, Guangdong Provincial Key Laboratory
of Magnetoelectric Physics and Devices, School of Physics, Sun Yat-Sen
University, Guangzhou, 510275, China  Pengcheng Dai  Department of Physics and
Astronomy, Rice University, Houston, TX 77005, USA  Robert J. Birgeneau
Department of Physics, University of California, Berkeley, California 94720,
USA  Materials Science Division, Lawrence Berkeley National Laboratory,
Berkeley, California 94720, USA  Dao-Xin Yao  Center for Neutron Science and
Technology, Guangdong Provincial Key Laboratory of Magnetoelectric Physics and
Devices, School of Physics, Sun Yat-Sen University, Guangzhou, 510275, China
Ming Yi  Department of Physics and Astronomy, Rice University, Houston, TX
77005, USA  Meng Wang<EMAIL_ADDRESS>Center for Neutron Science
and Technology, Guangdong Provincial Key Laboratory of Magnetoelectric Physics
and Devices, School of Physics, Sun Yat-Sen University, Guangzhou, 510275,
China

###### Abstract

We report the synthesis and characterization of RbNi2Se2, an analog of the
iron chalcogenide superconductor RbxFe2Se2, via transport, angle resolved
photoemission spectroscopy, and density functional theory calculations. A
superconducting transition at $T_{c}$ = 1.20 K is identified. In normal state,
RbNi2Se2 shows paramagnetic and Fermi liquid behaviors. A large Sommerfeld
coefficient yields a heavy effective electron mass of $m^{*}\approx 6m_{e}$.
In the superconducting state, zero-field electronic specific-heat data
$C_{es}$ can be described by a two-gap BCS model, indicating that RbNi2Se2 is
a multi-gap superconductor. Our density functional theory calculations and
angle resolved photoemission spectroscopy measurements demonstrate that
RbNi2Se2 exhibits relatively weak correlations and multi-band characteristics,
consistent with the multi-gap superconductivity.

##  I Introduction

Since the discovery of copper oxide superconductors, researchers have
extensively searched for superconductivity in materials with transition
metalsBednorz et al. (1987); Takada et al. (2003). Significant progress has
been made in iron pnictide and chalcogenide compounds, where several
structural systems have been identified with the highest $T_{c}$ of 55 K
achieved in LaFeAsORen et al. (2008). Superconductivity was also observed in
the compounds consisting of Cr and Mn under pressure, such as CrAs, KCrAs, and
MnPWu et al. (2014); Mu et al. (2017); Cheng et al. (2015). Among all of them,
superconductivity has been found to be in the vicinity of an antiferromagnetic
(AF) order, suggesting that spin fluctuations may play an important role in
the mechanism of superconductivity. Nickel oxide materials have analogous
structures with copper oxide superconductors. Superconductivity with
$T_{c}=9-15$ K has also been observed in films of nickel-based compoundsLi et
al. (2019); Gu and Wen (2022).

$A$Ni2Se2 ($A=$ K, Cs, and Tl) crystalizes in the ThCr2Si2 structure and shows
metallic behavior and Pauli paramagnetism. At low temperatures,
superconductivity emerges with the superconducting (SC) transition
temperatures of $T_{c}$ $\approx$ 0.8 K for KNi2Se2 polycrystals, 2.7 K for
CsNi2Se2, and 3.7 K for TlNi2Se2 single crystalsNeilson et al. (2012); Chen,
Huimin and Yang, Jinhu and Cao, Chao and Li, Lin and Su, Qiping and Chen, Bin
and Wang, Hangdong and Mao, Qianhui and Xu, Binjie and Du, Jianhua and Fang,
Minghu (2016); Wang et al. (2013). While K0.95Ni1.86Se2 single crystals do not
show superconductivity down to 0.3 K, yielding that the superconductivity is
sensitive to the stoichiometry of the samplesLei et al. (2014). As a
comparison, AxFe2Se2 system exhibits $T_{c}$s ranging from 20 $-$ 30 K. With
different amount of iron vacancies, AxFe2-δSe2 exhibits a variety of AF orders
and iron vacancy ordersGuo et al. (2010); Fang et al. (2013); Dai (2015); Wang
et al. (2016). The replacement of Fe by Co suppresses the superconductivity
and induces a ferromagnetic (FM) order in RbCo2Se2Yang et al. (2013); Huang et
al. (2021). The $A$Ni2Se2 superconductors with a formal valence of Ni1.5+ have
been revealed to exhibit remarkable properties. In particular, they usually
exhibit a large Sommerfeld coefficient $\gamma$, suggesting a large density of
states and unconventional pairing at low temperatureWang et al. (2013); Chen,
Huimin and Yang, Jinhu and Cao, Chao and Li, Lin and Su, Qiping and Chen, Bin
and Wang, Hangdong and Mao, Qianhui and Xu, Binjie and Du, Jianhua and Fang,
Minghu (2016). One possibility that was proposed was that the large Sommerfeld
coefficient might be induced by local charge orderNeilson et al. (2012).
However, angle resolved photoemission spectroscopy (ARPES) measurements yield
weak electronic correlations in KNi2Se2 and the origin of the large Sommerfeld
coefficient may be driven by the large density states and the Van Hove
singularity in the vicinity of the Fermi energyFan et al. (2015).

In this work, we report the successful synthesis and characterization of
RbNi2Se2 single crystals. The crystal structure, electronic band structure,
and transport properties have been investigated. We find that RbNi2Se2 is a
Pauli paramagnetism and exhibits a SC transition at $T_{c}$ = 1.20 K. Normal-
state specific heat measurements suggest an effective electronic mass
enhancement with $m^{*}\approx 6m_{e}$. In SC state, a two-gap BCS model can
match well with the zero-field electronic specific heat, indicating that
RbNi2Se2 is a multi-gap superconductor. Comparison with the density functional
theory (DFT) calculations and ARPES measurements reveals that RbNi2Se2 is a
weakly correlated superconductor with multi bands crossing the Fermi level.

Figure 1: (a) XRD measurement on the $ab$ plane of RbNi2Se2 single crystal.
The inset shows the crystal structure. The pink, green, and gray balls
represent Rb, Se, and Ni ions. (b) A pattern of XRD measured on powder
samples. The inset is a photo of RbNi2Se2 single crystal.

##  II Experimental and calculation details

Single crystals of RbNi2Se2 were grown by the self-flux method. First, the
precursor NiSe was prepared by heating Ni powders and Se pellets at
$500^{\circ}$C. Then, NiSe powders and Rb were put into alumina crucibles
according to stoichiometry and sealed in an evacuated silica tube. The mixture
was kept at $150^{\circ}$C for 5 h, then heated to $760^{\circ}$C in 40 h and
kept for 5 h, after that heated to $1050^{\circ}$C in 40 h and held for 5 h.
Finally, the temperature was cooled down to $700^{\circ}$C at a rate of
$3.5^{\circ}$C/h. To prevent the reaction of Rb with water and air, all of the
processes were conducted in an argon-filled glove box. Shiny plate-like single
crystals with a typical size of 4 $\times$ 5 $\times$ 1 mm3 were grown as
shown in the inset of Fig. 1(b).

Table 1: Single-crystal of RbNi2Se2 refinement at 150 K. Formula weight | 223.14
---|---
Crystal system | Tetragonal
Space group | I4/mmm
Unit-cell parameters | a = b = 3.9272 (3) Å
| c = 13.8650(5) Å
| $\alpha=\beta=\gamma=90^{\circ}$
Atomic parameters |
Rb | 2b(0,0,1/2)
Ni | 4d(0,1/2,3/4)
Se | 4e(1/2,1/2,0.6502(1))
Density | 3.466 g/cm3
F(000) | 198
Radiation | Mo $\textit{K}\alpha$ ($\lambda$ = 0.7107Å)
2$\theta$ for data collection |  $10.794^{\circ}$ to $60.562^{\circ}$
Index ranges | -4 $\leq$ h $\leq$ 5, -5 $\leq$ k $\leq$ 5,
| -11 $\leq$ l $\leq$ 18
Reflections collected | 808
Independent reflections | 119
Data/restraints/parameters | 119/0/5
Goodness-of-fit on F2 | 1.204
Final R indexes [I $\geq$ 2$\sigma$(I)] |  $\textit{R}_{1}$ = 0.0378, $\textit{wR}_{2}$ =0.1001
Largest diff. peak/hole/e Å-3 | 3.06/-2.04

Single crystal x-ray diffraction (XRD) were conducted on a SuperNova (Rigaku)
x-ray diffractometer. The sample was blowed by N2 during the data collection
to avoid exposure to air. The elemental analysis was measured by using an
energy-dispersive x-ray spectroscopy (EDS) (EVO, Zeiss). Electrical transport,
magnetic and specific heat measurements were performed on a physical property
measurement system (PPMS, Quantum Design). The in-plane resistivity
$\rho_{ab}(T)$ was measured using the standard four-probe method on a
rectangular sheet crystal to keep current flowing in the ab-plane. The Vienna
Ab initio Simulation Package (VASP) was employed for the DFT
calculationsKresse and Hafner (1993). ARPES measurements were performed on a
helium-lamp based system with a DA30 electron analyzer. Single crystals were
cleaved _in-situ_ in ultra-high vacuum with a base pressure better than
$5\times 10^{-11}$ Torr at 30 K. Energy and angular resolutions were better
than 20 meV and $0.1^{\circ}$, respectively.

##  III Results and discussions

All peaks from single crystal XRD can be indexed with the ThCr2Si2-type
structure (space group: I4/mmm), which is illustrated in the inset of Fig.
1(a). The determined lattice parameters are $a=b=3.9272(3)$, and
$\textit{c}=13.8650(5)$ Å at 150 K with the volume of unit cell between that
of KNi2Se2 and CsNi2Se2. Details of the atom coordinates and other key
information are shown in Table 1. To show the quality of the samples, we
present a $\theta$$-$2$\theta$ scan of a single crystal along the
$(H=0,K=0,L)$ in Fig. 1(a) and an XRD pattern on a powder sample in Fig. 1(b),
where $(H,K,L)$ are Miller indices in reciprocal lattice units. No peaks from
impurity could be identified. The EDS results for several single crystals are
rather homogenous and the determined average atomic ratios are Rb:Ni:Se =
1.16(4):2.04(3):2.00(6) when the content of Se is normalized to be 2, close to
the stoichiometry of RbNi2Se2.

Figure 2: (a) Temperature dependence of in-plane resistivity $\rho_{ab}(T)$
with magnetic fields H = 0, 0.15, and 11 T along the c axis. The inset shows
the fitted result using $\rho_{ab}($T$)=\rho_{0}+A$T2 from 1.8 to 40 K. The
red line is a fitting curve. (b) Temperature dependence of ZFC and FC magnetic
susceptibility with applying various external magnetic fields along the $c$
axis. The inset shows a magnetization hysteresis loop at 3K. (c) Temperature
dependence of Cp(T ) from 2 to 300 K.

Temperature dependence of the resistivity for RbNi2Se2 is shown in Fig. 2(a).
The electric current is applied in the ab-plane. The value of $\rho_{ab}$ is
about 92.1 $\mu\Omega\cdot$cm at 300 K and only about 1.5 $\mu\Omega\cdot$cm
at 1.8 K. The residual resistivity ratio (RRR) of 61.4 [$\rho_{ab}(300\
\text{K})/\rho_{ab}(1.8\ \text{K})$] suggests remarkable metallicity and high
quality of the single crystalsLei et al. (2014); Ehrlich and Rivier (1968);
Böhmer et al. (2016). The resistivity measured at 0, 0.15, and 11 T shows a
metallic behavior without any anomaly or strong magnetic field dependence. The
resistivity $\rho_{ab}$ below 40 K can be well described by the equation
$\rho_{ab}($T$)=\rho_{0}+A$T2 as shown in the inset of Fig. 2(a), where
$\rho_{0}$ = 1.211 $\mu\Omega\cdot$cm and A = 0.012 $\mu\Omega\cdot$cm/K2,
revealing a paramagnetic Fermi liquid behaviorAnalytis et al. (2014). The
magnetic susceptibility is nearly independent of temperature, yielding a Pauli
paramagnetic behaviour as shown in Fig. 2(b). A weak FM sign is revealed from
the hysteresis loop shown in the inset of Fig. 2(b), which could be ascribed
to a small amount of Ni impurityLi et al. (2020). The specific heat from 2 to
200 K shown in Fig. 2(c) also suggests that no phase transition occurs in this
temperature range.

To explore the possible superconductivity, we show measurements down to 50 mK
in Fig. 3. The magnetic susceptibility at low temperatures is shown in Fig.
3(a). A clear diamagnetic response appears below 1.20 K under zero field,
indicating a SC transition. The transition shifts to lower temperatures with
an increase of the $dc$ bias field, consistent with the Meissner effect of
superconductivity. With knowing $T_{c}$ $\approx$ 0.8 K for KNi2Se2Neilson et
al. (2012) and $T_{c}\approx$ 2.7 K for CsNi2Se2Chen, Huimin and Yang, Jinhu
and Cao, Chao and Li, Lin and Su, Qiping and Chen, Bin and Wang, Hangdong and
Mao, Qianhui and Xu, Binjie and Du, Jianhua and Fang, Minghu (2016), the
diamagnetic response at 1.20 K should correspond to the SC transition of
RbNi2Se2. We plot the onset SC transition temperatures at various magnetic
fields and fit the upper critical field $\mu_{0}H_{c2}(0)$ using the Ginzburg-
Landau theory with the formula
$\mu_{0}H_{c2}(T)=\mu_{0}H_{c2}(0)\times(1-t^{2})/(1+t^{2})$, where $t$ is the
reduced temperature $t=T/T_{c}$. As shown in Fig. 3(b), the resultant
$\mu_{0}H_{c2}(0)$ = 1.38 T is within the Pauli limit, $\mu_{0}H^{p}_{c2}(0)$
= 2.37 T, indicating a weak coupling behaviorBao et al. (2015).

Figure 3: (a) Diamagnetization of superconductivity under various magnetic
fields. (b) Upper critical field $H_{c2}$, as a function of temperature. The
red solid line shows a fit with the Ginzburg-Landau expression as discussed in
the text. (c) The low-temperature specific heat of RbNi2Se2, measured at
various fields near superconducting transition. The inset shows the Schotkky
anomaly at zero field and 5 T and a fitting of $\textit{C}_{p}$ in normal
state. (d) Reduced temperature $\textit{T}/\textit{T}_{c}$ dependence of
electronic specific heat divided by temperature and $\gamma_{N}$,
$\textit{C}_{es}/\gamma_{N}\textit{T}$, in the SC state at zero field, where
$\textit{C}_{es}=C_{p}-C_{lattice}-C_{N}$. The two solid lines show the
fitting curves of the one-gap BCS model and the two-gap model to
$\textit{C}_{es}/\gamma_{N}\textit{T}$, respectively. The dashed lines show
the contributions from two different gaps, $\Delta_{1}$ and $\Delta_{2}$,
respectively.

stop here

Figure 3(c) displays the specific heat measured at low temperatures. In normal
state, a fit to the specific heat $C_{p}$ from 1.5 to 3.8 K using $C_{p}/T$ =
$\gamma_{N}$+$\beta T^{2}$ results in the Sommerfeld coefficient $\gamma_{N}$
= 30.30 mJ$\cdot$mol-1K-2 and $\beta$ = 3.67 mJ$\cdot$mol-1K-4, as shown in
the inset. The Debye temperature $\Theta_{D}$ is estimated to be of 167 K from
the equation $\Theta_{D}$ = $(12\pi^{4}NR/5\beta)^{1/3}$, where $N=5$ is the
atomic number in each formula unit and $R$ is the ideal gas constant. The
$\gamma_{N}$ for RbNi2Se2 is comparable with that of KNi2Se2 ($\sim$44
mJ$\cdot$mol-1K-2)Neilson et al. (2012), CsNi2Se2 ($\sim$77
mJ$\cdot$mol-1K-2)Chen, Huimin and Yang, Jinhu and Cao, Chao and Li, Lin and
Su, Qiping and Chen, Bin and Wang, Hangdong and Mao, Qianhui and Xu, Binjie
and Du, Jianhua and Fang, Minghu (2016), and TlNi2Se2 ($\sim$40
mJ$\cdot$mol-1K-2)Wang et al. (2013).The effective mass of electrons $m^{*}$
can be estimated through Eq. 1DeLong et al. (1985):

${m^{*}}={\hbar^{2}k_{F}^{2}\gamma_{N}}/{\pi^{2}nk_{B}^{2}}$ (1)

where $k_{B}$ is the Boltzmann constant and the carrier density $n$ is
calculated by the number of electrons ($Z$) per cell volume ($V$). Using a
spherical Fermi surface approximation, the Fermi wave vector can be estimated
by $k_{F}$ = $(3\pi^{2}n)^{1/3}$. Assuming that Ni contributes 1.5 electrons
($Z$ = 6), we obtain $k_{F}$ = $9.4\times 10^{9}$ m-1. The estimated effective
mass of electrons $m^{*}/m_{e}$ = 6 for RbNi2Se2 is significantly enhanced
compared with the bare electron mass $m_{e}$. Combining the parameters from
fitting the electronic specific heat and the quadratic temperature dependent
regime of resistivity, the Kadowaki-Woods ratio, $A$/$\gamma_{N}^{2}$, is
calculated to be 0.94 $\times$ $10^{-5}\mu\Omega\cdot$cm(mol$\cdot$K2mJ)2,
close to $\sim$10${}^{-5}\mu\Omega\cdot$cm(mol$\cdot$K2mJ)2 of heavy fermion
systemsKadowaki and Woods (1986). This scaling relation yields RbNi2Se2 with
the heavy-fermion behavior.

In SC state, the specific heat data reveals a clear $\lambda$-anomaly under
zero field with a maximum at 0.94 K [Fig. 3(c)], suggesting bulk
superconductivity. Applying an external magnetic field, the SC transition
moves quickly to lower temperatures. As shown in the inset of Fig. 3(c), an
upturn appears on $C_{p}/T$ below $T^{2}<0.1$ K2. The small upturn could be
described by the Schottky anomaly for paramagnetic impurity spinsNeilson et
al. (2012), which can be well fitted by $C_{N}$ = $D(H)T^{-2}$ with $D(H=0)$ =
0.07 mJ$\cdot$K$\cdot$mol-1.

In order to get information about the SC gap, we extract the electronic
specific heat $C_{es}$ by subtracting the phonon contribution $C_{lattice}$
and Schottky anomaly $C_{N}$ from the total $C_{p}$,
$C_{es}=C_{p}-C_{lattice}-C_{N}$. The electronic specific heat data
$\textit{C}_{es}(\gamma_{N}\textit{T})^{-1}$ against $T/T_{c}$ are shown in
Fig. 3(d), where $\gamma_{N}$ has been revealed in the normal state. The
$T_{c}$ is determined to be 0.94 K using an equal-area entropy construction.
The heat jump at transition is 0.63, smaller than that of the theoretical
value of 1.43 in the BCS weak-coupling scenarioDaams and Carbotte (1981). We
employ both one-gap and two-gap BCS models $C_{es}=\sum
C_{i}exp(-\Delta_{i}/k_{B}T$) to fit the electronic specific heat, where
$\Delta_{i}$ is the size of the $i^{th}$ SC gap at 0 K. The one-gap model
reveals $\Delta_{0}/k_{B}T_{c}$ = 0.30 that is smaller than the value of 1.76
in the BCS theoryBardeen et al. (1957). In the two-gap model, the total
specific heat can be considered as the sum of electronic contributions from
two bands. Our fitting yields the sizes of two gaps of $\Delta_{1}/k_{B}T_{c}$
= 2.51 and $\Delta_{2}/k_{B}T_{c}$ = 0.25, respectively. The ratio of the
contributions from the two gaps is $\sim 1.5:1$ as presented by the dashed
lines in Fig. 3(d). The better fitting using the two-gap model indicates that
RbNi2Se2 may be a multi-gap superconductor.

Figure 4: (a) Projected band structure of Ni $3d$ in RbNi2Se2. Colors of the
circle represent the different nickel orbitals, and the weight of each orbital
is displayed by the size of the circles. The dashed lines are guide to eyes.
(b) Density of states (DOS) near the Fermi level of RbNi2Se2. The red, blue
and purple curves represent contributions from Ni, Rb and Se, respectively.
Inset: The coordinates of the high symmetry $\textbf{k}\textendash$path in
reciprocal space of the tetragonal unit cell. (c) Measured spectral images
along the high symmetry direction $\varGamma\sim M$. The black lines are band
structure calculated by the DFT. (d) The Fermi surface mapping at 35 K with
the Brillouin zone (BZ) marked by a red square.

The electron-phonon coupling strength is also calculated by employing the
inverted McMillan formulaMcMillan (1968):

${\lambda_{ep}}=\frac{1.04+\mu^{*}{\rm
ln}(\Theta_{D}/1.45T_{c})}{(1-0.62\mu^{*}){\rm
ln}(\Theta_{D}/1.45T_{c})-1.04}$ (2)

where $\mu^{*}$ represents the Coulomb repulsion pseudopotential, which we
adopt $\mu^{*}$ = 0.13 for this systemAmon et al. (2018). Generally, the
$\lambda_{ep}$ for strongly coupled superconductors are close to 1, and
$\lambda_{ep}\rightarrow$ 0.5 is viewed as weak coupled
superconductorsMcMillan (1968). The $\lambda_{ep}$ for RbNi2Se2 is 0.49,
suggesting that RbNi2Se2 is a weakly coupled superconductor, consistent with
the specific heat analysis.

To check the multi-band character, we conducted DFT calculations and ARPES
measurements on the electronic band structure. Figures 4(a) and 4(b) show the
calculated electronic bands of Ni ions and the integrated density of states of
Ni, Rb, and Se. The electronic states near the Fermi surface are governed by
the $3d$ orbitals of Ni ions. In Fig 4(c), the measured electronic structure
of RbNi2Se2 is compared with the band structure obtained from the DFT
calculations. The calculated curves (black lines) are scaled and overlaid onto
the photoemission intensity along the $\Gamma-M$ direction. The theoretical
and experimental data match qualitatively with a renormalization factor of
about 1.8. This factor is relatively small compared to other iron-based
superconductorsYi et al. (2017), indicating that the electronic correlations
are moderate. The moderate correlations suggest that the superconductivity is
most likely originated from the electron-phonon coupling, similar to the two-
gap superconductor MgB2Bouquet et al. (2001a, b). Figure 4(d) shows the
photoemission intensity map of the Fermi surface, overlaid with the DFT
calculated Fermi surfaces. The multiple Fermi pockets observed are consistent
with the multi-band behavior of this compound, reinforcing the condition for
the two-gap feature in the SC state. The band structure is reminiscent of the
closely-related compound RbCo2Se2Huang et al. (2021). However, the
$3d_{x^{2}+y^{2}}$ flat band that induces itinerant ferromagnetism in RbCo2Se2
is observed here to be well below the Fermi level, due to the electron doping
resulted from the replacement of Co by Ni.

As $A$Ni2Se2 ($A=$ K, Cs, and Tl) superconductors, RbNi2Se2 also exhibits an
enhanced effective electron mass with a large Sommerfeld coefficient from
specific heat. However, ARPES data reveal that the electronic correlation is
not strong. To understand the reason of the large Sommerfeld coefficient,
$\gamma$ is estimated with $\lambda_{ep}$ and $N$(0) at the Fermi level
through the relationshipXiao et al. (2021):

${\gamma}=\frac{\pi^{2}k_{B}^{2}N(E_{F})(1+\lambda_{ep})}{3}$ (3)

From the DOS in Fig. 4(b), $N(E_{F})$ is estimated to be 11.68 states/eV per
formula, resulting in $\gamma$ = 40.98 mJ$\cdot$mol-1K-2 that is close to the
experimental value of 30.30 mJ$\cdot$mol-1K-2. Therefore, the result suggests
that the large $\gamma_{N}$ is related to the large DOS at the Fermi level as
proposed in KNi2Se2Fan et al. (2015), instead of the heavy fermion state. For
$A$Ni2Se2 ($A$ = K, Rb, and Cs), the increase of the atomic radiuses of alkali
metals works as applying a pressure to the Ni-Se layers. The DOS of Ni $3d$
orbitals at the Fermi surface, the Sommerfeld coefficient $\gamma_{N}$, and
the electronic correlations are all enhanced. In the framework of the BCS
theory, the SC transition temperature $T_{c}$ increases accordingly, as the
experimental observations in KNi2Se2, RbNi2Se2, and CsNi2Se2Neilson et al.
(2012); Chen, Huimin and Yang, Jinhu and Cao, Chao and Li, Lin and Su, Qiping
and Chen, Bin and Wang, Hangdong and Mao, Qianhui and Xu, Binjie and Du,
Jianhua and Fang, Minghu (2016); Wang et al. (2013).

##  IV Summary

In summary, we have successfully synthesized single crystals of RbNi2Se2 and
characterized the physical properties. RbNi2Se2 is found to be a weakly
coupled superconductor with $T_{c}=1.2$ K. In normal state, RbNi2Se2 exhibits
Fermi liquid behavior and Pauli paramagnetism. In the SC state, the zero-field
electronic specific heat can be well described with a two-gap BCS model,
indicating that RbNi2Se2 possesses multi-gap feature. DFT calculations and
ARPES measurements demonstrate that multi electronic bands of the $3d$
orbitals of Ni ions cross the Fermi level. Our analyses reveal that the large
Sommerfeld coefficient of RbNi2Se2 is originated from the large DOS at the
Fermi surface.

##  V acknowledgments

We thank Shiliang Li for fruitful discussions. Work at Sun Yat-Sen University
was supported by the National Natural Science Foundation of China (Grants No.
12174454, 11904414, 11904416, 11974432, U2130101), the Guangdong Basic and
Applied Basic Research Foundation (No. 2021B1515120015), National Key Research
and Development Program of China (No. 2019YFA0705702, 2018YFA0306001,
2017YFA0206203), and GBABRF-2019A1515011337. Work at SLAB was supported by NSF
of China with Grant No. 12004270, and GBABR-2019A1515110517. ARPES work at
Rice is supported by the Robert A. Welch Foundation Grant No. C-2024 (M. Y.).
P. D. is also supported by US Department of Energy, BES under Grant No. DE-
SC0012311. Work at Berkeley was funded by the U.S. Department of Energy,
Office of Science, Office of Basic Energy Sciences, Materials Sciences and
Engineering Division under Contract No. DE-AC02-05-CH11231 (Quantum Materials
program KC2202).

## References

  * Bednorz et al. (1987) J. G. Bednorz, M. Takashige, and K. A. Müller, Europhys. Lett. 3, 379 (1987), URL http://link.springer.com/10.1007/BF01303701.
  * Takada et al. (2003) K. Takada, H. Sakurai, E. Takayama-Muromachi, F. Izumi, R. A. Dilanian, and T. Sasaki, Nature 422, 53 (2003), ISSN 0028-0836, URL http://www.nature.com/articles/nature01450.
  * Ren et al. (2008) Z.-A. Ren, G.-C. Che, X.-L. Dong, J. Yang, W. Lu, W. Yi, X.-L. Shen, Z.-C. Li, L.-L. Sun, F. Zhou, et al., EPL (Europhysics Lett. 83, 17002 (2008), ISSN 0295-5075, URL https://iopscience.iop.org/article/10.1209/0295-5075/83/17002.
  * Wu et al. (2014) W. Wu, J. Cheng, K. Matsubayashi, P. Kong, F. Lin, C. Jin, N. Wang, Y. Uwatoko, and J. Luo, Nat. Commun. 5, 5508 (2014), ISSN 2041-1723, URL http://dx.doi.org/10.1038/ncomms6508http://www.nature.com/articles/ncomms6508.
  * Mu et al. (2017) Q.-G. Mu, B.-B. Ruan, B.-J. Pan, T. Liu, J. Yu, K. Zhao, G.-F. Chen, and Z.-A. Ren, Phys. Rev. B 96, 140504 (2017), ISSN 2469-9950, URL https://link.aps.org/doi/10.1103/PhysRevB.96.140504.
  * Cheng et al. (2015) J.-G. Cheng, K. Matsubayashi, W. Wu, J. P. Sun, F. K. Lin, J. L. Luo, and Y. Uwatoko, Phys. Rev. Lett. 114, 117001 (2015), ISSN 0031-9007, URL https://link.aps.org/doi/10.1103/PhysRevLett.114.117001.
  * Li et al. (2019) D. Li, K. Lee, B. Y. Wang, M. Osada, S. Crossley, H. R. Lee, Y. Cui, Y. Hikita, and H. Y. Hwang, Nature 572, 624 (2019), ISSN 0028-0836, URL http://dx.doi.org/10.1038/s41586-019-1496-5http://www.nature.com/articles/s41586-019-1496-5.
  * Gu and Wen (2022) Q. Gu and H.-H. Wen, Innov. 3, 100202 (2022), ISSN 26666758\.
  * Neilson et al. (2012) J. R. Neilson, A. Llobet, A. V. Stier, L. Wu, J. Wen, J. Tao, Y. Zhu, Z. B. Tesanovic, N. P. Armitage, and T. M. McQueen, Phys. Rev. B 86, 054512 (2012), ISSN 1098-0121, URL https://link.aps.org/doi/10.1103/PhysRevB.86.054512.
  * Chen, Huimin and Yang, Jinhu and Cao, Chao and Li, Lin and Su, Qiping and Chen, Bin and Wang, Hangdong and Mao, Qianhui and Xu, Binjie and Du, Jianhua and Fang, Minghu (2016) Chen, Huimin and Yang, Jinhu and Cao, Chao and Li, Lin and Su, Qiping and Chen, Bin and Wang, Hangdong and Mao, Qianhui and Xu, Binjie and Du, Jianhua and Fang, Minghu, Supercond. Sci. Technol. 29, 045008 (2016), URL https://iopscience.iop.org/article/10.1088/0953-2048/29/4/045008.
  * Wang et al. (2013) H. Wang, C. Dong, Q. Mao, R. Khan, X. Zhou, C. Li, B. Chen, J. Yang, Q. Su, and M. Fang, Phys. Rev. Lett. 111, 207001 (2013), ISSN 0031-9007, URL https://link.aps.org/doi/10.1103/PhysRevLett.111.207001.
  * Lei et al. (2014) H. Lei, M. Abeykoon, K. Wang, E. S. Bozin, H. Ryu, D. Graf, J. B. Warren, and C. Petrovic, J. Phys. Condens. Matter 26, 015701 (2014), ISSN 0953-8984, URL https://iopscience.iop.org/article/10.1088/0953-8984/26/1/015701.
  * Guo et al. (2010) J. Guo, S. Jin, G. Wang, S. Wang, K. Zhu, T. Zhou, M. He, and X. Chen, Phys. Rev. B 82, 180520 (2010), ISSN 1098-0121, URL https://link.aps.org/doi/10.1103/PhysRevB.82.180520.
  * Fang et al. (2013) M. Fang, H. Wang, C. Dong, and Q. Huang, J. Phys. Conf. Ser. 449, 012015 (2013), ISSN 1742-6596, URL https://iopscience.iop.org/article/10.1088/1742-6596/449/1/012015.
  * Dai (2015) P. Dai, Rev. Mod. Phys. 87, 855 (2015), ISSN 0034-6861, eprint 1503.02340, URL https://link.aps.org/doi/10.1103/RevModPhys.87.855.
  * Wang et al. (2016) M. Wang, M. Yi, W. Tian, E. Bourret-Courchesne, and R. J. Birgeneau, Phys. Rev. B 93, 075155 (2016), ISSN 2469-9950, eprint 1508.04482, URL https://link.aps.org/doi/10.1103/PhysRevB.93.075155.
  * Yang et al. (2013) J. Yang, B. Chen, H. Wang, Q. Mao, M. Imai, K. Yoshimura, and M. Fang, Phys. Rev. B 88, 064406 (2013), ISSN 1098-0121, URL https://link.aps.org/doi/10.1103/PhysRevB.88.064406.
  * Huang et al. (2021) J. Huang, Z. Wang, H. Pang, H. Wu, H. Cao, S.-K. Mo, A. Rustagi, A. F. Kemper, M. Wang, M. Yi, et al., Phys. Rev. B 103, 165105 (2021), ISSN 2469-9950, URL https://link.aps.org/doi/10.1103/PhysRevB.103.165105.
  * Fan et al. (2015) Q. Fan, X. P. Shen, M. Y. Li, D. W. Shen, W. Li, X. M. Xie, Q. Q. Ge, Z. R. Ye, S. Y. Tan, X. H. Niu, et al., Phys. Rev. B 91, 125113 (2015), ISSN 1098-0121, eprint arXiv:1411.6455v1, URL https://link.aps.org/doi/10.1103/PhysRevB.91.125113.
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993), ISSN 0163-1829, URL https://link.aps.org/doi/10.1103/PhysRevB.47.558.
  * Ehrlich and Rivier (1968) A. Ehrlich and D. Rivier, J. Phys. Chem. Solids 29, 1293 (1968), ISSN 00223697, URL https://linkinghub.elsevier.com/retrieve/pii/0022369768901819.
  * Böhmer et al. (2016) A. E. Böhmer, V. Taufour, W. E. Straszheim, T. Wolf, and P. C. Canfield, Phys. Rev. B 94, 024526 (2016), ISSN 2469-9950, eprint 1606.00500, URL https://link.aps.org/doi/10.1103/PhysRevB.94.024526.
  * Analytis et al. (2014) J. G. Analytis, H.-h. Kuo, R. D. McDonald, M. Wartenbe, P. M. C. Rourke, N. E. Hussey, and I. R. Fisher, Nat. Phys. 10, 194 (2014), ISSN 1745-2473, URL http://dx.doi.org/10.1038/nphys2869http://www.nature.com/articles/nphys2869.
  * Li et al. (2020) Q. Li, C. He, J. Si, X. Zhu, Y. Zhang, and H.-H. Wen, Communications Materials 1, 16 (2020), ISSN 2662-4443, eprint 1911.02420.
  * Bao et al. (2015) J.-K. Bao, J.-Y. Liu, C.-W. Ma, Z.-H. Meng, Z.-T. Tang, Y.-L. Sun, H.-F. Zhai, H. Jiang, H. Bai, C.-M. Feng, et al., Phys. Rev. X 5, 011013 (2015), ISSN 2160-3308, eprint 1412.0067, URL https://link.aps.org/doi/10.1103/PhysRevX.5.011013.
  * DeLong et al. (1985) L. E. DeLong, R. P. Guertin, S. Hasanain, and T. Fariss, Phys. Rev. B 31, 7059 (1985), ISSN 0163-1829, URL https://link.aps.org/doi/10.1103/PhysRevB.31.7059.
  * Kadowaki and Woods (1986) K. Kadowaki and S. Woods, Solid State Commun. 58, 507 (1986), ISSN 00381098, URL https://linkinghub.elsevier.com/retrieve/pii/0038109886907854.
  * Daams and Carbotte (1981) J. M. Daams and J. P. Carbotte, J. Low Temp. Phys. 43, 263 (1981), ISSN 0022-2291, URL http://link.springer.com/10.1007/BF00116155.
  * Bardeen et al. (1957) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957), ISSN 0031-899X, URL https://www.taylorfrancis.com/books/9780429964251%****␣RbNi2Se2_draft.bbl␣Line␣375␣****https://link.aps.org/doi/10.1103/PhysRev.108.1175.
  * McMillan (1968) W. L. McMillan, Phys. Rev. 167, 331 (1968), ISSN 0031-899X, URL https://link.aps.org/doi/10.1103/PhysRev.167.331.
  * Amon et al. (2018) A. Amon, E. Svanidze, R. Cardoso-Gil, M. N. Wilson, H. Rosner, M. Bobnar, W. Schnelle, J. W. Lynn, R. Gumeniuk, C. Hennig, et al., Phys. Rev. B 97, 014501 (2018), ISSN 2469-9950, URL https://link.aps.org/doi/10.1103/PhysRevB.97.014501.
  * Yi et al. (2017) M. Yi, Y. Zhang, Z. Shen, and D. Lu, npj Quantum Materials 2, 57 (2017).
  * Bouquet et al. (2001a) F. Bouquet, Y. Wang, R. A. Fisher, D. G. Hinks, J. D. Jorgensen, A. Junod, and N. E. Phillips, Europhys. Lett. 56, 856 (2001a), ISSN 0295-5075, eprint 0107196, URL https://iopscience.iop.org/article/10.1209/epl/i2001-00598-7.
  * Bouquet et al. (2001b) F. Bouquet, R. A. Fisher, N. E. Phillips, D. G. Hinks, and J. D. Jorgensen, Phys. Rev. Lett. 87, 047001 (2001b), ISSN 0031-9007, URL https://link.aps.org/doi/10.1103/PhysRevLett.87.047001.
  * Xiao et al. (2021) G. Xiao, Q. Zhu, Y. Cui, W. Yang, B. Li, S. Wu, G.-H. Cao, and Z. Ren, Sci. China Physics, Mech. Astron. 64, 107411 (2021), ISSN 1674-7348, URL https://link.springer.com/10.1007/s11433-021-1731-3.


# The optimal weighting function for cosmic magnification measurement through
foreground galaxy-background galaxy (quasar) cross correlation

Xiaofeng Yang1,2, Pengjie Zhang1
$1$Key Laboratory for Research in Galaxies and Cosmology, Shanghai
Astronomical Observatory, Chinese Academy of Sciences,
Nandan Road 80, Shanghai, 200030, China
$2$Graduate University of Chinese Academy of Sciences, 19A, Yuquan Road,
Beijing, 100049, China<EMAIL_ADDRESS>
###### Abstract

Cosmic magnification has been detected through cross correlation between
foreground and background populations (galaxies or quasars). It has been shown
that weighing each background object by its $\alpha-1$ can significantly
improve the cosmic magnification measurement [16, 23]. Here, $\alpha$ is the
logarithmic slope of the luminosity function of background populations.
However, we find that this weighting function is optimal only for sparse
background populations in which intrinsic clustering is negligible with
respect to shot noise. We derive the optimal weighting function for general
case including scale independent and scale dependent weights. The optimal
weighting function improves the S/N (signal to noise ratio) by $\sim 20\%$ for
a BigBOSS-like survey and the improvement can reach a factor of $\sim 2$ for
surveys with much denser background populations.

###### keywords:

cosmology: weak lensing–cosmic magnification–theory

##  1 Introduction

Weak gravitational lensing directly probes the matter distribution of the
universe (e.g.[20]) and is emerging as one of the most powerful probes of dark
matter, dark energy [1] and the nature of gravity [12]. By far the most
sophisticated way to measure weak lensing is through cosmic shear, which is
the lensing induced coherent distortion in galaxy shape (Fu et al. 6 and
references therein). Coordinated projects on precision weak lensing
reconstruction through galaxy shapes have been carried out extensively (STEP,
Heymans et al. 8; STEP2, Massey et al. 15; GREAT08, Bridle et al. 3; GREAT10,
Kitching et al. 13).

Alternatively, one can reconstruct weak lensing through cosmic magnification,
namely the lensing induced coherent distortion in galaxy number density (e.g.
Gunn [7], Ménard & Bartelmann [16], Jain et al. [10] and references therein).
Neglecting high order corrections, the lensed galaxy (quasar) number
overdensity $\delta^{\rm L}_{\rm g}$ is related to the intrinsic overdensity
$\delta_{\rm g}$ by

$\delta^{\rm L}_{\rm g}=\delta_{\rm g}+2(\alpha-1)\kappa.$ (1)

Here, $\kappa$ is the lensing convergence and $\alpha(m,z)=2.5~{}{\rm dlog}\
n(m,z)/{\rm d}m-1$ is a function of the galaxy apparent magnitude $m$ and
redshift $z$. The number count of galaxy brighter than $m$ is
$N(m)=\int^{m}n(m)dm$. Throughout the paper we use the superscript “L” to
denote the lensed quantity.

Since cosmic magnification does not involve galaxy shape, weak lensing
reconstruction through cosmic magnification automatically avoids all problems
associated with galaxy shape. The key step in such reconstruction is to
eliminate $\delta_{\rm g}$, which is often orders of magnitude larger than the
lensing signal $\kappa$. Usually this is done by cross correlating foreground
population (galaxies) and background population (galaxies or quasars) with no
overlapping in redshift [23, 9, 24]. The resulting cross correlation is

$\displaystyle\langle\delta^{\rm L}_{\rm g,f}(\hat{\theta})\delta^{\rm L}_{\rm
g,b}(\hat{\theta^{\prime}})\rangle\approx 2(\alpha_{\rm
b}-1)\langle\delta_{\rm g,f}(\hat{\theta})\kappa_{\rm
b}(\hat{\theta^{\prime}})\rangle\ .$ (2)

Throughout the paper we use the subscript “b” for quantities of background
population and “f” for that of foreground galaxies. The above relation
neglects a term proportional to $\langle\kappa_{\rm f}\kappa_{\rm b}\rangle$,
which is typically much smaller than the $\langle\delta_{\rm g,f}\kappa_{\rm
b}\rangle$ term.111This term can be non-negligible or even dominant for
sufficiently high redshift foreground galaxy samples [27].

It is important to weigh the cross correlation measurement appropriately to
improve the S/N (signal to noise ratio). Since the signal scales as
$\alpha-1$, Ménard & Bartelmann [16] first suggested to maximize the S/N by
weighing each galaxy with its own $\alpha-1$. This weighting significantly
improves the measurement and a robust $8\sigma$ detection of the cosmic
magnification was achieved for the first time [23].

Nevertheless, we find that, the $\alpha-1$ weighting is optimal only in the
limit where the background galaxy (quasar) intrinsic clustering is negligible
with respect to the shot noise in background distribution. The statistical
errors (noises) are contributed by both shot noise and intrinsic clustering of
foreground and background galaxies. In this letter, we derive the exact
expression of the weighting function optimal for the cosmic magnification
measurement through cross correlation. The new weighting can further improve
the S/N by $\sim 20\%$ for galaxy-quasar cross correlation measurement in a
BigBOSS-like survey. We can also employ high redshift galaxies instead of
quasars as background sources which can have much larger number density and
smaller bias. Smaller shot noise results into better performance for the
derived optimal weighting. The improvement over the $\alpha-1$ weighting can
reach a factor of $\sim 2$ for surveys with background population density of
$\sim 2/$arcmin2.

Throughout the paper, we adopt the fiducial cosmology as a flat $\Lambda$CDM
universe with $\Omega_{\Lambda}=0.728$, $\Omega_{m}=1-\Omega_{\Lambda}$,
$\sigma_{8}=0.807$, $h=0.702$ and initial power index $n=0.961$, which are
consistent with WMAP seven years best-fit parameters [14].

##  2 THE OPTIMAL WEIGHTING FUNCTION

We are seeking for an optimal weighting function linearly operating on the
background galaxy (quasar) number overdensity in flux (magnitude) space. Let’s
denote the background galaxy number overdensity of the $i$-th magnitude bin as
$\delta^{(i)}_{\rm g,b}$ and the corresponding weighting function as $W_{i}$.

  * •

The simplest weighting function is scale independent, so the weighted
background galaxy overdensity is

$\sum_{i}W_{i}\delta^{(i)}_{\rm g,b}\ .$ (3)

  * •

We can further add scale dependence in $W_{i}$ to increase the S/N. The new
weighting function convolves the density field. For brevity, we express it in
Fourier space as $W_{i}(\ell)$. The Fourier transformation of the weighted
background overdensity is

$\sum_{i}\sum_{\ell}W_{i}(\ell)\tilde{\delta}^{(i)}_{\rm g,b}(\vec{\ell})\ .$
(4)

Here, $\tilde{\delta}_{\rm g,b}$ is the Fourier component of the overdensity
$\delta_{\rm g,b}$. The weighting function $W_{i}(\ell)$ is real and only
depends on the amplitude of the wavevector $\ell\equiv|\vec{\ell}|$. It
guarantees the weighted overdensity to be real.

The S/N of the background-foreground galaxy cross-correlation depends on the
weighting function, so we use the subscript “W” to denote the S/N after
weighting. The overall S/N can be conveniently derived in the Fourier space,

$\displaystyle\left(\frac{S}{N}\right)_{W}^{2}=$
$\displaystyle\sum_{\ell}\left[\frac{\langle C^{\rm
CM-g}_{\ell}\rangle_{W}}{\langle\Delta C^{\rm
CM-g}_{\ell}\rangle_{W}}\right]^{2}$ $\displaystyle=$
$\displaystyle\sum_{\ell}\frac{(2\ell+1)\Delta\ell f_{\rm sky}\langle C^{\rm
CM-g}_{\ell}\rangle_{W}^{2}}{\langle C^{\rm
CM-g}_{\ell}\rangle_{W}^{2}+(\langle C_{\rm g,b}\rangle_{W}+\langle C_{\rm
s,b}\rangle_{W})(C_{\rm g,f}+C_{\rm s,f})}$ $\displaystyle=$
$\displaystyle\sum_{\ell}\frac{(2\ell+1)\Delta\ell f_{\rm sky}}{1+(C_{\rm
g,f}+C_{\rm s,f})\frac{\langle b_{\rm g,b}W\rangle^{2}C_{\rm m,b}+\langle
W^{2}\rangle C_{\rm s,b}}{4\langle W(\alpha_{b}-1)\rangle^{2}C_{\kappa_{b}\rm
g}^{2}}}\ .$ (5)

Here, $C^{\rm CM-g}$ is the cosmic magnification-galaxy cross correlation
power spectrum and $\Delta C^{\rm CM-g}$ is the associated statistical error.
$\langle\cdots\rangle_{W}$ denotes the weighted average of the corresponding
quantity. We then have

$\langle C^{\rm CM-g}\rangle_{W}=\langle 2(\alpha_{\rm b}-1)W\rangle
C_{\kappa_{\rm b}\rm g}\ .$ (6)

Here, $C_{\kappa_{\rm b}\rm g}$ is the cross correlation power spectrum
between background lensing convergence and foreground galaxy overdensity.
$\langle uv\rangle$ is the averaged product of $uv$,

$\displaystyle\langle uv\rangle\equiv\frac{\sum_{i}u(m_{i})v(m_{i})N_{{\rm
b},i}}{\sum_{i}N_{{\rm b},i}}\ .$ (7)

Here, $N_{{\rm b},i}$ is the number of background galaxies (quasars) in the
given magnitude bin $m_{i}-\Delta m_{i}/2<m<m_{i}+\Delta m_{i}/2$.

The S/N scales with $f^{1/2}_{\rm sky}$ and $f_{\rm sky}$ is the fractional
sky coverage. $C_{\rm s}$ is the shot noise power spectrum and the weighted
one is $\langle C_{\rm s,b}\rangle_{W}=\langle W^{2}\rangle C_{\rm s,b}$.
$C_{\rm g}$ is the galaxy power spectrum. We adopt a simple bias model for the
foreground and background galaxies. We then have $C_{{\rm g},i}=b_{{\rm
g},i}^{2}C_{\rm m}$ where $b_{{\rm g},i}$ is the bias of the $i$-th magnitude
bin and $C_{\rm m}$ is the corresponding matter angular power spectrum. The
weighted background galaxy power spectrum is

$\langle C_{\rm g,b}\rangle_{W}=\langle b_{\rm g,b}W\rangle^{2}C_{\rm m,b}\ .$
(8)

###  2.1 The scale independent optimal weighting function

The optimal weighing function $W$ can be obtained by varying the S/N (Eq. 2)
with respect to $W$ and maximizing it. The derivation is lengthy but
straightforward, so we leave details in the appendix and only present the
final result here.

The optimal weighting function is of the form222The derived scale independent
weighting function implicitly assumes no scale dependence in the galaxy bias.
In reality, the galaxy bias is scale dependent and the application of Eq. 9 is
limited. The exact optimal weighting function applicable to scale dependent
bias is given by Eq. 10.

$\displaystyle W=(\alpha_{\rm b}-1)+\varepsilon b_{\rm g,b}\ \ .$ (9)

where the scale independent constant $\varepsilon$ is determined by Eq. 12. It
is a fixed number for the given redshift bin of the given survey. In the limit
that shot noise of background galaxies overwhelms their intrinsic clustering
($C_{\rm s,b}\gg C_{\rm g,b}$), $\varepsilon\rightarrow 0$. In this case, the
weighting function $\alpha-1$ proposed by Ménard & Bartelmann [16] becomes
optimal.

###  2.2 The scale dependent optimal weighting function

The weighting function $W$ (Eq. 9) is optimal under the condition that $W$ is
scale independent. If we relax this requirement and allow for scale dependence
in $W$, we are able to maximize the S/N of the cross power spectrum
measurement at each $\ell$ bin. Clearly, this further improves the overall
S/N.

In this case, $W$ of different $\ell$ bins are independent to each other. This
significantly simplifies the derivation and we are now able to obtain an
analytical expression for $W$,

$\displaystyle W(\ell)=(\alpha_{\rm b}-1)+\left[-\frac{\langle(\alpha_{\rm
b}-1)b_{\rm g,b}\rangle C_{\rm m,b}(\ell)/C_{\rm s,b}}{1+\langle b_{\rm
g,b}^{2}\rangle C_{\rm m,b}(\ell)/C_{\rm s,b}}\right]b_{\rm g,b}(\ell).$ (10)

This form is similar to Eq. (9), except that it is now scale dependent. Here
again, in the limit $C_{\rm s,b}\gg C_{\rm m,b}\sim C_{\rm g,b}$,
$W\rightarrow\alpha-1$ and we recover the result of Ménard & Bartelmann [16].
This is indeed the case for SDSS background quasar sample.

###  2.3 The applicability

Are the derived weighting functions (Eq. 9 & 10) directly applicable to real
surveys? From Eq. 9 & 10, it seems that we need to figure out $b_{\rm g,b}$
first. Since $b^{2}_{\rm g,b}\equiv C_{\rm g,b}/C_{\rm m,b}$ and $C_{\rm m,b}$
is not directly given by the observation, cosmological priors or external
measurements (e.g. weak lensing) are required to obtain the absolute value of
$b_{\rm g,b}$. Hence it seems that the applicability of Eq. 9 & 10 is limited
by cosmological uncertainty.

However, this is not the case. Eq. 10 shows that, it is the combination
$b_{\rm g,b}^{2}C_{\rm m,b}$ determines $W$. Since $C_{\rm g,b}\equiv b_{\rm
g,b}^{2}C_{\rm m,b}$ and $\alpha_{b}$ are directly observable, Eq. 10 is
determined completely by observations. Closer look shows that Eq. 9 is also
determined by the combination $b_{\rm g,b}^{2}C_{\rm m,b}$, so the
corresponding weighting is determined completely by observations, too. Hence
the derived optimal weighting functions are indeed directly applicable to real
surveys.

For ongoing and planned surveys such as CFHTLS, COSMOS, DES, BigBOSS, LSST,
SKA, etc., the number density of background populations (galaxies) can be high
and the intrinsic clustering can be non-negligible or even dominant comparing
to shot noise. In next section we will quantify the improvement of the optimal
weighting functions for a BigBOSS-like survey and briefly discuss implications
to surveys with even denser background populations.

##  3 The improvement

Table 1: Improving the cosmic magnification measurement by the optimal weighting function. The target survey is BigBOSS. The terms on the left side of the arrows are the estimated S/N using the weighting $\alpha-1$. The ones on the right side are the S/N using the optimal weighting function, where the ones on the left side of $``/"$ are what expected using the scale independent weighting (Eq. 9), and the ones on the right are what expected using the scale dependent weighting (Eq. 10). The improvement depends on the bias dependence on galaxy luminosity. To demonstrate such dependence, we adopt a toy model $b_{Q}\propto F^{\beta}$ and investigate different values of the parameter $\beta$. In general, the optimal weighting function improves the S/N by $10\%$-$20\%$ for BigBOSS, whose background quasar density is $\sim 0.02/$arcmin2. The improvement can reach a factor of $\sim 2$ for surveys with background (galaxy) populations reaching surface density of $\sim 2/$arcmin2.  Flux dependence of quasar bias model | $\beta=0$ | 0.1 | 0.2 | 0.3
---|---|---|---|---
Detection significance of $\rm LRG\times quasar$ | 111.2$\rightarrow$129.0/136.8 | 109.8$\rightarrow$126.3/133.9 | 108.3$\rightarrow$123.7/131.1 | 106.5$\rightarrow$119.3/126.7
Detection significance of $\rm ELG\times quasar$ | 94.3$\rightarrow$106.4/110.3 | 93.2$\rightarrow$104.5/108.6 | 92.0$\rightarrow$102.2/106.7 | 90.5$\rightarrow$99.1/104.1

BigBOSS333http://bigboss.lbl.gov/ is a planned spectroscopic redshift survey
of $24000\ {\rm deg}^{2}$ (BigBOSS-North plus South). Cosmic magnification can
be measured by BigBOSS through LRG (luminous red galaxy)-quasar and ELG
(emission line galaxy)-quasar cross correlations. In principle, it can also be
measured through LRG-ELG cross correlation. But the interpretation of the
measured cross correlation signal would be complicated by the intricate
selection function of ELGs [28]. In the current paper, we only consider the
LRG-quasar and ELG-quasar cross correlations.

There are some uncertainties in the BigBOSS galaxy (quasar) redshift
evolution, flux distribution and intrinsic clustering. To proceed, we will
take a number of simplifications. So the absolute S/N of cross correlation
measurement that we calculate is by no means accurate. But our calculation
should be sufficiently robust to demonstrate the relative improvement of the
exact optimal weighting function over the previous one.

The LRG and ELG luminosity functions are calculated based on the BigBOSS white
paper [21]. The comoving number density of LRG and ELG is $3.4\times
10^{-4}(h/\rm Mpc)^{3}$, then we have $1.1\times 10^{7}$ LRGs in the redshift
range of $z=0.2-1$ and $3.3\times 10^{7}$ ELGs in the redshift range of
$z=0.7-1.95$. Clustering of LRGs evolves slowly, so we adopt LRG bias as
$b_{\rm g,f}(z)=1.7/D(z)$ [19]. Here $D(z)$ is the linear density growth
factor and is normalized such that $D(z=0)=1$. Existing knowledge on
clustering of ELGs is rather limited. So we simply follow Padmanabhan et al.
[19] and approximate the ELG bias as $b_{\rm g,f}=0.8/D(z)$.

For the luminosity function of background quasars, we adopt the LDDE
(Luminosity dependent density evolution) model with best fit parameters from
Croom et al. [5]. The magnitude limit is $g=23$, then we have $2.1\times
10^{6}$ quasars in the redshift range of $z=2-3.5$. We choose a redshift gap
($z_{\rm b,min}-z_{\rm f,max}=0.05$) such that the intrinsic cross correlation
between foreground and background populations can be safely neglected. We
adopt a bias model for quasar clustering, with $b_{\rm
Q}(z)=0.53+0.289(1+z)^{2}$ from the analysis of $3\times 10^{5}$quasars [18].

The S/N depends on many issues and can vary from $90$-$140$ (Table 1). The
S/Ns of LRG-quasar and ELG-quasar correlations are comparable because of a
consequence of several competing factors including the lensing efficiency,
galaxy surface density and clustering. Nevertheless, a robust conclusion is
that BigBOSS can measure cosmic magnification through galaxy-quasar cross
correlation measurement with high precision. Given such high S/N and accurate
redshift available in BigBOSS, it is feasible to directly measure the angular
diameter distance from such measurement by the method of Jain & Taylor [11],
Zhang et al. [25], Bernstein [2].

Unambiguous improvement in the cross correlation measurement by using our
optimal weighting (Eq. 9 & 10) is confirmed, as shown in Table 1. The S/N is
improved by $\sim 15\%$ by using the scale independent optimal weighting (Eq.
9) and by $\sim 20\%$ by using the scale dependent one (Eq. 10).

We further investigate the dependence of the above improvement on the flux
dependence of quasar bias. We adopt a toy model with $b_{\rm Q}(z,F)=b_{\rm
Q}(z)(F/F^{*})^{\beta}$. Here, $F^{*}$ is the flux corresponding to that of
the $M^{*}$ in the quasar luminosity function. $\beta$ is an adjustable
parameter and we will try $\beta=0,0.1,0.2,0.3$, then the corresponding
parameters in scale independent weighting are $\varepsilon=0.049,0.048,0.047$
and $0.045$. Larger value of $\beta$ ($\geq 0.4$) leads to too large quasar
bias ($\geq 10$) and hence will not be investigated here. Table 1 shows
consistent improvement by our optimal weighting functions. Hence, despite
uncertainties in quasar modeling, we expect our optimal weighting function to
be useful to improve cosmic magnification measurement in the BigBOSS survey.

Nevertheless the improvement is only moderate. The major reason is that, even
for BigBOSS, the quasar sample is still sparse, with a surface number density
$\sim 0.02/$arcmin2. Hence shot noise dominates over the intrinsic clustering.
Indeed, we find that typically $C_{\rm m,Q}/C_{\rm s,Q}\la 0.1$. For imaging
surveys like CFHTLS, COSMOS, DES and LSST, we can also use high redshift
galaxies as background galaxies to correlate with low redshift foreground
galaxies. For these surveys, high redshift galaxy population (e.g. with
$z>1$-$2$) can reach surface number density $\sim 0.2$-$2/$arcmin2 or even
higher. So the shot noise in these surveys can be suppressed by a factor of
$10$-$100$ or more. The overall improvement of our optimal weighting would be
larger.

To demonstrate these further improvements, we hypothetically increase the
surface density of BigBOSS quasars by a factor of 10 and 100 respectively, but
keep $\beta=0$ and all other parameters unchanged. Shot noise will be
decreased by a factor of 10 and 100 correspondingly. The scale independent
weighting parameter $\varepsilon$ can reach 0.12 and 0.22 respectively. For
the first case, the S/N is improved by $\sim 38\%$ for the scale independent
optimal weighting and by $\sim 51\%$ for the scale dependent one. For the
second case, the improvement is $\sim 72\%$ for the scale independent optimal
weighting and is $\sim 94\%$ for the scale dependent one. It is now clear that
for measuring cosmic magnification through cross correlation between
foreground and background galaxies of many existing and planned surveys, one
should adopt the optimal weighting function derived in this letter.

##  4 Summary

We have derived the optimal weighting functions for cosmic magnification
measurement through cross correlation between foreground and background
populations, for scale independent and scale dependent weights respectively.
Our weighting functions outperform the commonly used weighting function
$\alpha-1$ by $\sim 20\%$ for a BigBOSS-like survey and by larger factors for
surveys with denser background populations. Hence we recommend to use our
optimal weighting function for cosmic magnification measurement in BigBOSS,
CFHTLS, COSMOS, DES, Euclid, LSST, SKA, WFIRST, etc.

##  5 Acknowledgment

This work is supported in part by the one-hundred talents program of the
Chinese Academy of Sciences, the national science foundation of China (grant
No. 10821302, 10973027 & 11025316), the CAS/SAFEA International Partnership
Program for Creative Research Teams and the 973 program grant No.
2007CB815401.

## References

  * Albrecht et al. [2006] Andreas Albrecht, et al. The Dark Energy Task Force Final Report. arXiv:astro-ph/0609591
  * Bernstein [2006] Bernstein, G. 2006, ApJ, 637, 598
  * Bridle et al. [2009] Bridle, S., et al. 2009, Annals of Applied Statistics, Vol. 3, No. 1, 6
  * Bridle et al. [2009] Bridle, S., et al. arXiv:astro-ph/0908.0945
  * Croom et al. [2009] Croom,S.M., et al. arXiv:astro-ph/0907.2727
  * Fu et al. [2008] Fu, L., et al. 2008, AAP, 479, 9
  * Gunn [1967] Gunn, J. E. 1967, ApJ, 150, 737
  * Heymans et al. [2006] Heymans,C., et al. 2007, MNRAS, 368,1323
  * Hildebrandt et al. [2009] Hildebrandt,H., et al. arXiv:astro-ph/0906.1580v2
  * Jain et al. [2003] Jain, B., Scranton, R., & Sheth, R. K. 2003, MNRAS, 345, 62
  * Jain & Taylor [2003] Jain, B., & Taylor, A. 2003, PRL, 91, 141302
  * Jain & Zhang [2008] Jain, B.,& Zhang, P. 2008, PRD, 78,063503
  * Kitching et al. [2010] Kitching, T., 2010, arXiv:astro-ph/1009.0779
  * Komatsu et al. [2011] Komatsu, E., et al. 2011, Astrophys.J.Suppl.192,18
  * Massey et al. [2007] Massey,R., 2007, MNRAS, 376,13
  * Ménard & Bartelmann [2002] Ménard, B. & Bartelmann, M., 2002, AAP, 386, 784
  * Menard et al. [2003] Ménard, B., Hamana, T., Bartelmann, M., & Yoshida, N. 2003, AAP, 403, 817
  * Myers et al. [2007] Myers,A.D.,et al. 2007, ApJ, 658, 85
  * Padmanabhan et al. [2006] Padmanabhan,N., et al. 2006, MNRAS,359,237
  * Refregier [2003] Refregier,A. 2003, Annual Review of Astronomy & Astrophysics, 41, 645
  * Schlegel et al. [2009] Schlegel,D., et al. arXiv:astro-ph/0904.0468
  * Schrabback et al. [2009] Schrabback,T., et al. arXiv:astro-ph/0911.0053
  * Scranton et al. [2005] Scranton,R.,et al. 2005, ApJ, 633, 589
  * Wang et al. [2011] Wang,L., et al. arXiv:astro-ph/1101.4796
  * Zhang et al. [2005] Zhang, J., Hui, L., & Stebbins, A. 2005, ApJ, 635, 806
  * Zhang & Pen [2005] Zhang, P.,& Pen, U.-L. 2005, PRL, 95, 241302
  * Zhang & Pen [2006] Zhang, P.,& Pen, U.-L. 2006, MNRAS, 367, 169
  * Zhu et al. [2009] Zhu,G., et al. 2009, ApJ, 701, 86

##  Appendix A Deriving the optimal weighting function

In Section 2, we give the optimal weighting function without derivation. Here,
we present a brief derivation for the scale independent weighting function.
The derivation of the scale dependent weighting function is similar, but
simpler.

Maximizing the S/N requires the variation of $(S/N)^{2}$ with respect to $W$
to be zero. From this condition, we obtain

$\displaystyle W$
$\displaystyle=\frac{\sum_{\ell}(2\ell+1)\Delta\ell\eta(\nu\langle b_{\rm
g,b}W\rangle^{2}+\langle W^{2}\rangle)\big{/}(1+F)^{2}}{\langle W(\alpha_{\rm
b}-1)\rangle\sum_{\ell}(2\ell+1)\Delta\ell\eta\big{/}(1+F)^{2}}(\alpha_{\rm
b}-1)$ (11) $\displaystyle-\frac{\langle W(\alpha_{\rm b}-1)\rangle\langle
b_{\rm
g,b}W\rangle\sum_{\ell}(2\ell+1)\Delta\ell\eta\nu\big{/}(1+F)^{2}}{\langle
W(\alpha_{\rm
b}-1)\rangle\sum_{\ell}(2\ell+1)\Delta\ell\eta\big{/}(1+F)^{2}}b_{\rm g,b}.$

Here, for brevity, we have denoted $\eta=C_{\rm s,b}(C_{\rm g,f}+C_{\rm
s,f})/C_{\kappa_{b}g}^{2}$, $\nu=C_{\rm m,b}/C_{\rm s,b}$ and
$F_{(W)}=\eta[\nu\langle b_{\rm g,b}W\rangle^{2}+\langle
W^{2}\rangle]/[4\langle W(\alpha_{\rm b}-1)\rangle^{2}]$.

Noticing that the coefficients of $\alpha_{\rm b}-1$ and $b_{\rm g,b}$ only
involve the weighted average of $W$ and taking the advantage that the optimal
$W$ remains optimal after a constant (flux independent) scaling, we are able
to seek for the solution of the form $W=(\alpha_{\rm b}-1)+\varepsilon b_{\rm
g,b}$, with $\varepsilon$ a constant to be determined. Plugging it into the
above equation, we obtain the equation of $\varepsilon$,

$\displaystyle\varepsilon=-\frac{\langle(\alpha_{\rm b}-1)b_{\rm
g,b}\rangle}{\sum_{\ell}\frac{(2\ell+1)\Delta\ell\eta}{(1+F)^{2}}\bigg{/}\sum_{\ell}\frac{(2\ell+1)\Delta\ell\eta\nu}{(1+F)^{2}}+\langle
b_{\rm g,b}^{2}\rangle}\ .$ (12)

Noticing that $F$ depends on $\varepsilon$. The above equation can be solved
numerically to obtain the solution of $\varepsilon$.

In the case of scale dependent weighting, each $W(\ell)$ are independent to
each other. Through similar derivation, one can show that

$\displaystyle W(\ell)=(\alpha_{\rm b}-1)+\left[-\frac{\nu\langle(\alpha_{\rm
b}-1)b_{\rm g,b}\rangle}{1+\nu\langle b_{\rm g,b}^{2}\rangle}\right]b_{\rm
g,b}\ .$ (13)


# A Novel Approach for Earthquake Early Warning System Design using Deep
Learning Techniques

Tonumoy Mukherjee
Advanced Technology Developement Center
Indian Institute of Technology Kharagpur
Kharagpur 721302, India
<EMAIL_ADDRESS>
&Chandrani Singh
Geology & Geophysics Department
Indian Institute of Technology Kharagpur
Kharagpur 721302, India
<EMAIL_ADDRESS>
&Prabir Kumar Biswas
Electronics & Electrical Communication Department
Indian Institute of Technology Kharagpur
Kharagpur 721302, India
<EMAIL_ADDRESS>

###### Abstract

Earthquake signals are non-stationary in nature and thus in real-time, it is
difficult to identify and classify events based on classical approaches like
peak ground displacement, peak ground velocity. Even the popular algorithm of
STA/LTA requires extensive research to determine basic thresholding parameters
so as to trigger an alarm. Also, many times due to human error or other
unavoidable natural factors such as thunder strikes or landslides, the
algorithm may end up raising a false alarm. This work focuses on detecting
earthquakes by converting seismograph recorded data into corresponding audio
signals for better perception and then uses popular Speech Recognition
techniques of Filter bank coefficients and Mel Frequency Cepstral Coefficients
(MFCC) to extract the features. These features were then used to train a
Convolutional Neural Network(CNN) and a Long Short Term Memory (LSTM) network.
The proposed method can overcome the above-mentioned problems and help in
detecting earthquakes automatically from the waveforms without much human
intervention. For the 1000Hz audio data set the CNN model showed a testing
accuracy of 91.1% for 0.2-second sample window length while the LSTM model
showed 93.99% for the same. A total of 610 sounds consisting of 310 earthquake
sounds and 300 non-earthquake sounds were used to train the models. While
testing, the total time required for generating the alarm was 1.68 seconds
which included individual times for data collection, processing, and
prediction. Taking into consideration the processing and prediction delays,
the total time is thus considered to be approximately 2 seconds. This shows
the effectiveness of the proposed method for EEW applications. Since the input
of the method is only the waveform, it is suitable for real-time processing,
thus, the models can very well be used also as an onsite earthquake early
warning system requiring a minimum amount of preparation time and workload.

##  1 Introduction

Earthquakes have been an integral part of the planet earth since time
immemorial. Any movement in the tectonic plates of the earth’s crust releases
massive amounts of energy which passes through the earth’s surface as seismic
waves and results in mild to tremendous shaking. This whole process although
seems to be very long but happens within a few minutes. As a result, the time
to respond to these occurrences have increased many folds. Various scales of
1-10 have been adopted as the magnitude indicator of earthquakes. For
addressing a wider range of earthquake sizes, the moment magnitude scale,
abbreviated $M_{W}$, is preferred and is applicable globally [1]. Magnitudes
of earthquakes are exponential. To put it simply, for each whole number that
we go up on a magnitude scale, the amplitude of the ground motion goes up by a
factor of 10 when recorded by a seismograph [2]. Thus, by using this scale as
a reference, it is realized that the level of ground shaking caused by a
magnitude 2 earthquake would be ten times more than a magnitude 1 earthquake
(and 32 times as much energy would be released). To put that into context, if
a magnitude 1 earthquake releases as much energy as blowing up 6 ounces of
TNT, a magnitude 8 earthquake would release as much energy as detonating 6
million tons of TNT. Major earthquakes can cause significant damage to life as
well as property. Preventing an earthquake from occurring is an illogical
thing to do. We need to focus on how to mitigate the devastation caused by
these events. For that, we need some robust and reliable prediction methods.
There were a few successful predictions. The Haicheng earthquake(1975) is a
perfect example of a successful prediction. Precursors for this event included
a foreshock sequence, peculiarity in animal behavior, and anomalies like
geodetic deformation, groundwater level differences[3]. But these
abnormalities were not present in other major earthquakes. Hence there is no
generic precursor for earthquakes. This is where modern approaches of machine
learning and deep learning comes into play. Although researches in seismology
using machine learning and deep learning are quite limited, with the huge
amount of data that are accessible to researchers, good researches like
foreshock identification in real-time using deep learning [4], Earthquake
detection and location using convolutional neural network [5] and Machine
Learning Seismic Wave Discrimination [6] have been done. Thus having an idea
about the occurrence of an earthquake as early as possible will result in
developing a certain amount of alertness to respond to such events with ease.
Unlike various methods of detecting an earthquake such as STA/LTA (Short-Term
Average/Long-Term Average) algorithm as shown in Figure: 1 [7], another
possible way of detecting and classifying earthquakes using their sounds, have
been explored in this research. Sound is an effect produced by anything
physical including earthquakes [8] and is one of the most common effects
reported during or immediately after the felt tremors caused by them[9]. An
earthquake of magnitude 5, will have a sound or vibration different from
another earthquake whose magnitude is 4 [10]. Also, similar magnitude
earthquakes should have similar vibrations or sounds that are
indistinguishable by us but can be picked up very efficiently by neural
networks.
In this paper, two types of machine learning models were used to train a
system with data consisting of earthquake and non-earthquake sounds.
Convolutional Neural Network (CNN) and Long Short Term Memory (LSTM) models
were used for the purpose. The former model showed a testing accuracy of
91.102% for a 0.2-second audio data sample and the latter showed a testing
accuracy of 93.999% for the same.

Figure 1: Figure showing a 6.8 magnitude earthquake acceleration signal. The
red line represents the first arrival of P-wave detected by STA/LTA algorithm
marked by the letter ’P’ and the S-wave arrival marked by the letter ’S’.

##  2 Deep Learning and Seismology

With the rapid increase in the seismic data quantity, major challenges are
faced by modern seismology in the fields of data analyzing and processing
techniques. Most of the many popular techniques that are used in major data
centers date back to the time when the amount of seismic data was small and
the computational power was limited.
Today with the advancements in the fields of machine learning and deep
learning, scientists and researchers can very easily extract useful
information from voluminous data as they provide a large collection of tools
to work with. Once trained with sufficient data, deep learning models just
like humans, can acquire their knowledge by extracting features from raw data
[11] to recognize natural objects and make expert-level decisions in various
disciplines. Besides, the high computational costs for training such networks
are balanced by their low-cost online operation [5]. Advantages like these,
make deep learning suitable for applications in real-time seismology and
earthquake early warning (EEW).
This paper would look into two commonly used and widely known Deep learning
models, CNN and LSTM respectively and the accuracy achieved by them in
detecting and classifying earthquakes from their sounds. The architecture of
both the networks would be looked at in detail in the following sections.

###  2.1 Convolutional Neural Network (CNN) Architecture

The proposed architecture for the CNN model has 4 convolutional layers, 1
maxpool layer, and 3 dense layers. 2(b) shows the first layer of the 4
convolutional layers consisting of 16 filters built with 3x3 convolution with
’Relu’ as the activation function and 1x1 stride. All the parameters for the
second, third, and fourth layers remain the same except for the number of
filters in each layer multiplies two times with the number of filters in the
previous layer. Or in other words, 16 filters in the first layer, 32 filters
in the second layer, 64 in the third and 128 in the final layer (shown in
2(a)).

(a) Block Diagram of the CNN Model used

(b) 16 Filters of the First Convolutional Layer arranged in a 4X4 Matrix

Figure 2: CNN Block Diagram with First Convolutional Layer Filters

The idea behind increasing the number of filters in each layer is to be more
specific about the features as the data starts to convolve through each layer.
A kernel of 2x2 has been used for the maxpool layer. The three dense layers
after the maxpool layer consist of 128, 64, and 2 neurons so as to pool down
the features for the final 2 class classification. The first two dense layers
use ’Relu’ as their activation function whereas the last dense layer uses
’Softmax’ as its activation function as we use categorical cross-entropy for
multi-class classification purposes and ’adam’ as the optimizer.

###  2.2 Long Short Term Memory (LSTM) Architecture

The proposed architecture for the LSTM model (shown in Figure: 3) has 2 LSTM
layers consisting of 128 neurons each. 4 time distributed fully connected
layers of 64, 32, 16, and 8 neurons respectively are added after the 2 LSTM
layers with ’relu’ as their activation function. Lastly, a dense layer
consisting of 2 neurons is added for the final 2 class classification with
’softmax’ as its activation function and ’adam’ as the optimizer.

Figure 3: Block Diagram of the LSTM Model used

##  3 Data

Indian Earthquake acceleration data of the past 10 years and magnitudes
approximately between 2-8$M_{W}$ was collected from PESMOS (Program for
Excellence in Strong Motion Studies, IIT Roorkee) [12] and EGMB (Eastern Ghats
Mobile Belt) Data of earthquakes with magnitudes ranging between 2-5$M_{W}$
Collected at Geology and Geophysics Lab IIT Kharagpur, were also used.

Table 1: Training Data for Deep Learning Models (CNN & LSTM)

Data Source | Data Type & Format | Feature Vector Dimensions | No. of Earthquake Data | No. of Non-Earthquake Data
---|---|---|---|---
PESMOS (IIT Roorkee) | Audio, .wav | 9 x 13 | 212 | 0
EGMB (IIT Kharagpur) | Audio, .wav | 9 x 13 | 68 | 280

Table 2: Testing Data for Deep Learning Models (CNN & LSTM)

Data Source | Data Type & Format | Feature Vector Dimensions | No. of Earthquake Data | No. of Non-Earthquake Data
---|---|---|---|---
PESMOS (IIT Roorkee) | Audio, .wav | 9 x 13 | 11 | 0
EGMB (IIT Kharagpur) | Audio, .wav | 9 x 13 | 21 | 18

Table 1 and Table 2 shows the Train-Test Data split for both the Deep Learning
models. The collected dataset was replicated twice. The first copy of the
dataset was converted into corresponding audio signals of .wav format by
keeping the original sensor sampling rate 200Hz and the second copy of the
same dataset was converted into corresponding audio signals of the same format
whose sampling rate was increased to 1000Hz using upsampling techniques. The
sampling rate of the upsampled audio signals was estimated using trial and
error method so as to hear and clearly notice the change of the signal with
respect to time. A total of 610 sounds consisting of 310 earthquake sounds and
300 non-earthquake sounds were used to train the models. Both the audio signal
datasets of 200Hz and 1000Hz were fed as inputs to both the above models.

##  4 Methodology

###  4.1 Data Preparation

Before training the models, proper cleaning of the audio data was done and
essential features from the data were then extracted using popular speech
processing methodologies of Filter Banks and Mel-Frequency Cepstral
Coefficients (MFCCs).

###  4.2 Pre-Emphasis

The very first step of data processing included filtering the data using a
Pre-Emphasis filter. This was mainly done to amplify the high frequencies.
Apart from amplification, the filter helped in balancing the frequency
spectrum since higher frequencies usually have smaller magnitudes compared to
lower frequencies. The filter also was able to improve the Signal-to-Noise
Ratio (SNR). The first order filter is represented by the equation :

$y(t)=x(t)-\alpha x(t-1)~~$ (1)

was used to apply the pre-emphasis filter over the range of audio signal data.
From the typical values of 0.95 and 0.97, the former was used for the filter
coefficient $\alpha$.

###  4.3 Framing and FFT

After pre-emphasis filtering, the data was divided into short-time frames to
avoid the loss of frequency contours of the signal over time for performing
Fourier transform. A good approximation of the frequency contours of the
signal was achieved by concatenating adjacent frames and applying the Fourier
transform over those short-time frames. Popular settings of 25ms for the
frame-size and 10ms stride were used for framing the data. A Hamming window
function (shown in Eq.2) was applied after the signal was sliced into frames.

$w[n]=0.54-0.46\cos\frac{2\pi n}{N-1},\text{for}\ 0\leq n\leq N-1$ (2)

where N represents the number of frames in which the signal was divided.
After dividing the signal into frames, an N-point FFT was performed on each
frame to calculate the frequency spectrum, which also happens to be the Short-
Time Fourier-Transform (STFT), where N is typically 256 or 512 (256 in this
case).

###  4.4 Filter Banks and Mel-Frequency Cepstral Coefficients (MFCCs)

The rationale behind using filter banks was to separate the input signal into
its multiple components such that each of those components carries a single
frequency sub-band of the original signal. Triangular filters, typically 26,
were applied on a Mel-scale to the power spectrum of the short-time frames to
extract the frequency bands as shown in Figure: 4.

The formula for converting from frequency to Mel scale is given by:

$M(f)=1125\times ln(1+\frac{f}{700})$ (3)

To go back from Mels to frequency, the formula used is given by:

$M^{-1}(m)=700\times(e^{m/1125}-1)$ (4)

Figure 4: Features extracted from the signal for training the models after
applying the filter bank to the power spectrums of the short-time frames of
the signal. Figure 5: Mel-Frequency Cepstral Coefficients (MFCCs) as extracted
from the short frames in which the signal was divided.

MFCC is a biologically inspired and by far the most successful and most used
feature in the area of speech processing [13]. The algorithm was used in
volcano classification also [14]. For speech signals, the mean and the
variance changes continuously with time, and thus it makes the signal non-
stationary [15]. Similarly like speech, earthquake signals are also non-
stationary [16] as each of them has different arrival of P, S, and surface
waves. Therefore, normal signal processing techniques like Fourier transform,
cannot be directly applied to it. But, if the signal is observed in a very
small duration window (say 25ms ), the frequency content in that small
duration appears to be more or less stationary. This opened up the possibility
of short-time processing of the earthquake sound signals. The small duration
window is called a frame, discussed in section 4.3. For processing the whole
sound segment, the window was moved from the beginning to the end of the
segment consistently with equal steps, called shift or stride. Based on the
frame-size and frame-stride, it gave us M frames. Now, for each of the frames
as depicted in Figure: 5, MFCC coefficients were computed.

Figure 6: Block diagram of the entire System

Moreover, the filter bank energies computed were highly correlated since all
the filterbanks were overlapping. This becomes a problem for most of the
machine learning algorithms. To reduce autocorrelation between the filterbank
coefficients and get a compressed representation of the filter banks, a
Discrete Cosine Transform (DCT) was applied to the filterbank energies. This
also allowed the use of diagonal covariance matrices to model the features for
training. Also, 13 of the 26 DCT coefficients were kept and the rest were
discarded due to the fact that fast changes in the filterbank energies are
represented by higher DCT coefficients. These fast changes resulted in
degrading the model performances. Thus a small improvement was observed by
dropping them.
The overall system representation could be better understood from Figure: 6.

##  5 Results and Discussion

The CNN and the LSTM models performed almost similarly for the 200Hz audio
data set, but significant improvements in the train-test accuracy percentages
are observed for the 1000Hz data set. For 1000Hz audio data set the CNN model
showed a testing accuracy of 91.102% for 0.2-second sample window length while
the LSTM model showed 93.999% for the same (shown in Figure: 7). This
observation can be backed by the fact that LSTMs performs better for
sequential or time-series data classifications [17].

Figure 7: Testing and Training Accuracy comparisons of CNN and LSTM Models
over several sampling window time frames for the audio signal data sets of
1000Hz & 200Hz sampling frequencies

The Kappa statistics(values), generally used for comparing an Observed
Accuracy with an Expected Accuracy (random chance), was used for validating
the model accuracies for both the data sets (shown in Figure: 8).
For both data classes, activations by random 5 out of 16 filters of the first
layer of the CNN Model along with their inputs is represented by Figure: 9 and
Figure: 10 respectively.

Figure 8: Testing Accuracy and Cohen Kappa Score comparisons of the CNN and
the LSTM Models over several sampling window time frames for the audio signal
data sets of 1000Hz and 200Hz sampling frequencies
Figure 9: 1st Convolutional layer activations for Earthquake Data with
corresponding Filters
Figure 10: 1st Convolutional layer activations for Non- Earthquake Data with
the same Filters as used for Earthquake data

The time required for generating the alarm by the CNN model includes 1.28
seconds to gather enough data to start the MFCC computations, 8ms for
processing and a prediction time of 0.2 seconds. The summation gives a result
of 1.68 seconds. Taking into consideration the processing and prediction
delays, the total time for the CNN model is thus considered to be
$\simeq{2}$seconds.

LSTMs being computationally a bit more expensive, the processing and the
prediction times were 10ms and 0.5 seconds respectively, giving a total of
2.19 seconds. Taking into consideration the processing and prediction delays,
the total time for the LSTM model is thus considered to be
$\simeq{2.5}$seconds.

Table 3: Overall Comparision between standard STA/LTA algorithm and the
proposed algorithms

Algorithm | Data Source | Data Sampling Frequency | Time to Alarm | Accuracy | Prerequisites
---|---|---|---|---|---
STA/LTA |  | PESMOS (IIT Roorkee)
---
&
EGMB (IIT Kharagpur)
200 Hz | 3 Seconds after P-wave arrival |  | 95.43%
---
(Heavily dependant on Prerequisites)
| 1\. Proper user-defined Thresholds
---
2\. Different Thresholds for different Regions
3\. Threshold kept high for Strong Motion Events
(more earthquakes missed, lesser false alarms)
4\. Threshold kept low for Weak Motion Events
(fewer earthquakes missed, more false alarms)
| MFCC with CNN Model
---
(Proposed Method)
| PESMOS (IIT Roorkee)
---
&
EGMB (IIT Kharagpur)
1000 Hz | 2 Seconds after P-wave arrival | 91.1% | None
| MFCC with LSTM Model
---
(Proposed Method)
| PESMOS (IIT Roorkee)
---
&
EGMB (IIT Kharagpur)
1000 Hz | 2.5 Seconds after P-wave arrival | 93.99% | None

It was also observed that by increasing the sampling rate of the audio
signals, the training and the testing accuracies of both CNN and LSTM model
increased. This is one of the very first applications of Machine learning in
the field of earthquake detection using MFCC’s and Filterbank Coefficients
which are generally used in the field of Speech Recognition. An earthquake is
basically understood well by three types of waves namely P-wave, S-wave, and
surface wave. Interaction of these waves with the surrounding medium gives an
earthquake its intensity. Any wave is a vibration. Any vibration has some
sound associated with it. It might be inaudible to the human ears, but the
sound remains. In real-time, it is difficult to identify and classify events
based on classical approaches like peak ground displacement, peak ground
velocity, or even the widely recognized algorithm of STA/LTA as they require
extensive research to determine basic thresholding parameters so as to trigger
an alarm. Many times due to human error or other unavoidable natural factors
such as thunder strikes or landslides, the conventional algorithms may end up
raising a false alarm (shown in Figure: 11).

Figure 11: Figure Showing 3 subplots where the first plot shows incorrect
P-wave arrival detection by STA/LTA, marked by the red line, due to poor
threshold parameters, leading to late detection of P-wave and in-turn leading
to late alarm generation. The second plot shows the spectrogram with almost
appropriate detection of the P-wave first arrival marked by the green line and
the third plot shows the same spectrogram with a -50dB threshold applied, so
as to get a clear view of the first P-wave arrival.

Table 3 shows a detailed comparision between the standard STA/LTA algorithm
and the proposed Deep Learning models. The main disadvantage of STA/LTA
algorithm is that it has to be tweaked differently for different types of
event detections. The user defined threshold which acts as the trigger, varies
from region to region due to its direct dependence on the geographical and
topological features of a particular region. Also, the threshold is kept high
for strong motion events and low for weak motion events. This results in
serious ambiguity as both the features can’t be used together. If the
threshold is high, more earthquakes are missed but lesser false alarms are
generated and if the threshold is low, less earthquakes are missed but more
false alarms generated. The proposed method can overcome these problems as it
can extract essential features of a raw sound or vibrational data without
manually hand engineering them and doesn’t require any knowledge or expertise
in the relevant field. It is also invariant to small variations in occurrence
in time or position and can understand representations of data with multiple
levels of abstraction. Since the input of the method is only the waveform, it
is suitable for real-time processing, thus, the models can very well be used
also as an onsite earthquake early warning system requiring a minimum amount
of preparation time and workload. Until now, the use of earthquake early
warning systems for earthquake disasters are mainly limited by false alarm
generation and delay in detection. By using the suggested approach, these
problems can be overcome, leading to automatic, fast, and accurate detection
of earthquake seismic signals.

###  5.1 Seismic Sensor Prototype Hardware Design and Experimental Setup for
Model Testing

The main reason for developing the sensor hardware prototype was to test the
validity and robustness of the proposed models. Also, it was very crucial to
understand and evaluate the complexities and solve the challenges at a system
level as it is targeted to deploy the system in earthquake prone areas. It was
also essential for comparing the system results with already present
commercial solutions.
The hardware is developed with the help of commercial equipments to mimic the
vibrations from an actual earthquake. Figure: 12 depicts a schematic of the
prototype sensor hardware used. Figure: 13(a) shows the actual PCB with boxes
representing the utility of sections. On the hard left we have the power
management block to supply a steady source of power to the circuit
(3.3V,80mA). On the top right we have the accelerometer sensor acting as the
principle of detection (seismometer) having an output dynamic range of 0.1V to
2.5V and a sensitivity of 300mV/g. In between the Power Management and the
accelerometer sensor blocks we have the microcontroller and wifi communication
blocks to receive and transmit the sensor data wirelessly. The Wifi has a
bandwidth of 100mHz to 10Hz, RF carrier frequency of 2.4GHz with a data
transmission rate upto 256Kbps. The processor is of a 16 bit RISC architecture
with an operating frequency of 16MHz.

Figure 12: Schematic of an accelerometer with front-end amplifiers (box in
black), and a microcontroller with wireless transceiver (box in red).

(a) Photo of the prototype circuit board for implementing a wireless seismic
sensor node.

(b) A 3D printed packaging scheme for the wireless seismic sensor prototype.

(c) Specifications of the Seismic sensor Prototype.

Figure 13: Prototype Sensor PCB with its 3D printed packaging and
Specifications

Figure: 13(b) and Figure: 13(c) gives us the image of the 3D printed packaging
solution used to house the sensor board and the exact specification details of
the components used to design the sensor PCB respectively. After the Deep
learning models were validated and finalized, an interesting experimental
setup was made so as to test the effectiveness and robustness of the models.
The idea was to mimic the vibrations from an actual earthquake signal, sense
it with the help of the custom designed sensor PCB, transmit the data
wirelessly to a remote computer where the data gets plotted, processed and
stored in real time. Also, as soon as the Deep Learning models detects the
required data, it activates and predicts the label so as to classify whether
or not the recorded data corresponds to an earthquake or a non earthquake
sound. The experimental setup can be visualized from Figure: 14.

Figure 14: Test Setup for the entire System.

The test setup and the overall system included the following steps:

  * •

An Arbitary Waveform Generator (AWG). Real earthquake data was stored in it so
that it could generate the output from that data.

  * •

A 4 ohm 40 Watt Speaker. The output of the AWG was connected to the speaker so
as to make the speaker vibrate according to the output waveform and thus mimic
an earthquake.

  * •

The Prototype circuit board consisting of the accelerometer sensor to sense
the speaker vibrations, microcontroller and an ESP8266 wifi communication
block (transmitter built within the PCB and receiver is coonected to the
remote PC via a USB to UART cable).

  * •

MATLAB code running in the remote PC to receive, process, plot and store the
data in real time.

  * •

Python code running simultaneously along with the MATLAB code in the remote
PC, monitoring the stored data directory to activate the deep learning model
for prediction as soon as it detects any data in that directory.

The test run of the system was a success as it was able to produce the desired
results by accurately predicting the class of the vibration signals from the
speaker received by the wireless ESP8266 Wi-Fi module.

##  6 Conclusion

In this paper, a new way of automatic classification of earthquake signals is
presented based on CNN and LSTM by using only MFCC features extracted from the
waveform. The performance of this algorithm has been tested by its application
to regional and local earthquake events selected from PESMOS (Program for
Excellence in Strong Motion Studies, IIT Roorkee) and EGMB (Eastern Ghats
Mobile Belt, Geology and Geophysics Lab IIT Kharagpur) data sets. Using
optimal parameters, for 1000 Hz audio data set the CNN model showed a testing
accuracy of 91.102% for a 0.2-second sample window while the LSTM model
obtained an accuracy of 93.999% for the same. This brings down the standard
alarm generation time to approximately 2 seconds after P-wave arrival. The
results outperform the conventional algorithm of STA/LTA in terms of
classification accuracy with respect to classification speed and constraint
requirements. While the models were mainly aimed at classifying earthquake
events as quickly as possible, they can also be easily tweaked to give early
estimates of magnitudes, ground velocities, shaking intensity, and many more
useful parameters.
Also, the uniquely innovative experimental setup helped in creating real time
earthquake simulations in a cheap, precise and effective way.
MFCC for earthquake detection is like putting our ears to hear and listen the
sounds inside the earth which we otherwise could not hear and thus treating
the earth sounds as earth’s speech signals. The most interesting and effective
part of this model is that it can be trained on various classes of sounds
other than earthquake sounds alone, to classify every signal that the sensor
detects, using their sound signatures. This can be of enormous help in
military applications also. If trained with human movement sounds, the model
could be deployed in the border and other high-security areas so as to provide
us instant information regarding trespassing and other such unlawful
activities, releasing the burden to some extent from the security officials
and the soldiers. This could be a solution not just for earthquake detection
alone but for many other such applications also.

## Acknowledgments

This work is funded by the Ministry of Human Resource Development(MHRD),
Government of India. The authors are thankful to the Department of Earthquake
Engineering, IIT Roorkee, and the Department of Geology and Geophysics, IIT
Kharagpur, for proving Program for Excellence in Strong Motion Studies(PESMOS)
and Easter Ghats Mobile Belt(EGMB) earthquake datasets respectively for the
research. This work would not have been possible without the research
facilities provided by IIT Kharagpur. Finally, authors are very thankful to
all the members of Image Processing and Computer Vision Laboratory, Department
of Electronics and Electrical Communication Engineering and Computational
Laboratory, Department of Geology and Geophysics, Indian Institute of
Technology Kharagpur, for their kind help and support.

## References

  * [1] Thomas C Hanks and Hiroo Kanamori.  A moment magnitude scale.  Journal of Geophysical Research: Solid Earth, 84(B5):2348–2350, 1979.
  * [2] Hiroo Kanamori.  Quantification of earthquakes.  Nature, 271(5644):411–414, 1978.
  * [3] Anshu Jin and Keiiti Aki.  Temporal change in coda q before the tangshan earthquake of 1976 and the haicheng earthquake of 1975.  Journal of Geophysical Research: Solid Earth, 91(B1):665–673, 1986\.
  * [4] K Vikraman.  A deep neural network to identify foreshocks in real time.  arXiv preprint arXiv:1611.08655, 2016.
  * [5] Thibaut Perol, Michaël Gharbi, and Marine Denolle.  Convolutional neural network for earthquake detection and location.  Science Advances, 4(2):e1700578, 2018.
  * [6] Zefeng Li, Men-Andrin Meier, Egill Hauksson, Zhongwen Zhan, and Jennifer Andrews.  Machine learning seismic wave discrimination: Application to earthquake early warning.  Geophysical Research Letters, 45(10):4773–4779, 2018.
  * [7] Wu Yih-Min, Hiroo Kanamori, Richard M Allen, and Egill Hauksson.  Determination of earthquake early warning parameters, $\tau$ c and pd, for southern california.  Geophysical Journal International, 170(2):711–717, 2007.
  * [8] PL Bragato, M Sugan, P Augliera, M Massa, A Vuan, and A Saraò.  Moho reflection effects in the po plain (northern italy) observed from instrumental and intensity data.  Bulletin of the Seismological Society of America, 101(5):2142–2152, 2011.
  * [9] Patrizia Tosi, Valerio De Rubeis, Andrea Tertulliani, and Calvino Gasparini.  Spatial patterns of earthquake sounds and seismic source geometry.  Geophysical research letters, 27(17):2749–2752, 2000.
  * [10] Patrizia Tosi, Paola Sbarra, and Valerio De Rubeis.  Earthquake sound perception.  Geophysical Research Letters, 39(24), 2012.
  * [11] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep learning.  MIT press, 2016.
  * [12] Himanshu Mittal, Ashok Kumar, and Rebecca Ramhmachhuani.  Indian national strong motion instrumentation network and site characterization of its stations.  International Journal of Geosciences, 3(06):1151, 2012.
  * [13] Md Afzal Hossan, Sheeraz Memon, and Mark A Gregory.  A novel approach for mfcc feature extraction.  In 2010 4th International Conference on Signal Processing and Communication Systems, pages 1–5. IEEE, 2010.
  * [14] I. Ã?lvarez, G. CortÃ©s, A. de la Torre, C. BenÃtez, L. GarcÃa, P. Lesage, R. ArÃ¡mbula, and M. GonzÃ¡lez.  Improving feature extraction in the automatic classification of seismic events. application to colima and arenal volcanoes.  In 2009 IEEE International Geoscience and Remote Sensing Symposium, volume 4, pages IV–526–IV–529, 2009.
  * [15] Mark G Frei and Ivan Osorio.  Intrinsic time-scale decomposition: time–frequency–energy analysis and real-time filtering of non-stationary signals.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 463(2078):321–342, 2006.
  * [16] JK Hammond and PR White.  The analysis of non-stationary signals using time-frequency methods.  Journal of Sound and vibration, 190(3):419–447, 1996.
  * [17] Sepp Hochreiter and Jrgen Schmidhuber.  Long short-term memory.  Neural Computation, 9(8):1735–1780, 1997.