Secular evolution of compact binaries revolving around a spinning massive black hole

TL;DR

This paper investigates how the spin of a supermassive black hole influences the long-term orbital evolution of nearby compact binaries, revealing oscillating effects at the 1.5PN order that differ from non-spinning cases.

Contribution

It introduces a detailed analysis of the spin-induced effects on binary orbital evolution at the 1.5PN order, highlighting oscillatory behavior absent in non-spinning SMBH models.

Findings

Spin of SMBH causes oscillating secular evolution of binary eccentricity.

Spin effects appear at 1.5PN order through angular momentum coupling.

Long-term orbital dynamics are significantly affected by SMBH spin.

Abstract

The leading order effect of the spin of a distant supermassive black hole (SMBH) on the orbit of compact binary revolving around it appears at the 1.5 post Newtonian (PN) expansion through coupling with the inner and outer orbital angular momenta, and makes an oscillating characteristic contribution to the secular evolution of the compact binary eccentricity over a long time scale compared to the Kozai-Lidov dynamics caused by a non-spinning SMBH.

Equations

$${\mathbf{v}}_1=\frac{d{\mathbf{r}}_1}{dt}=\frac{m_2}{m}\mathbf{v}+\mathbf{V},{\mathbf{v}}_2=\frac{d{\mathbf{r}}_2}{dt}=-\frac{m_1}{m}\mathbf{v}+\mathbf{V},$$

$${\mathbf{a}}_{\beta }^{[1.5\text{PN,spin}]}=2a{m}_3{\mathbf{v}}_{\beta }\times \frac{({\mathbf{e}}_Z-3\frac{({\mathbf{e}}_Z\cdot {\mathbf{r}}_{\beta }){\mathbf{r}}_{\beta }}{r_{\beta }^2})}{r_{\beta }^3},$$

$$\mathbf{a}$$

$$\mathbf{A}$$

$${\mathbf{a}}^{[N]}$$

$${\mathbf{A}}^{[N]}$$

$${\mathbf{a}}^{[\text{N}]}$$

$${\mathbf{A}}^{[\text{N}]}$$

$${\mathbf{a}}^{[1.5\text{PN,spin}]}$$

$${\mathbf{A}}^{[1.5\text{PN,spin}]}$$

$$\ell$$

$$\mathbf{n}$$

$$L$$

$$\mathbf{N}$$

$$\frac{d\mathbf{h}}{dt}=\mathbf{\ell }\times \delta \mathbf{a},m\frac{d\mathbf{R}}{dt}=\delta \mathbf{a}\times \mathbf{h}+\mathbf{v}\times (\mathbf{\ell }\times \delta \mathbf{a}),$$

$$\begin{array}{rl}\frac{dp}{dt}=& 2\sqrt{\frac{p^3}{m}}\frac{\mathcal{S}}{1+e\cos \varphi },\\ \frac{de}{dt}=& \sqrt{\frac{p}{m}}(\sin \varphi \mathcal{R}+\frac{2\cos \varphi +e+e{\cos \varphi }^2}{1+e\cos \varphi }\mathcal{S}),\\ \frac{d\bar{\omega }}{dt}=& \frac{1}{e}\sqrt{\frac{p}{m}}(-\cos \varphi \mathcal{R}+\frac{2+e\cos \varphi }{1+e\cos \varphi }\sin \varphi \mathcal{S}),\\ \frac{d\iota }{dt}=& \sqrt{\frac{p}{m}}\frac{\cos (\omega +\varphi )}{1+e\cos \varphi }\mathcal{W},\\ \sin \iota \frac{d\Omega }{dt}=& \sqrt{\frac{p}{m}}\frac{\sin (\omega +\varphi )}{1+e\cos \varphi }\mathcal{W},\end{array}$$

$$⟨F⟩=\frac{1}{T_{\text{out}}}\frac{1}{T_{\text{in}}}{\int }_0^{T_{\text{out}}}{\int }_0^{T_{\text{in}}}Fdtd{t}^{\prime },$$

$$\begin{array}{rl}⟨F⟩=& \frac{1}{4{\pi }^2}(1-e^2{)}^{3/2}(1-E^2{)}^{3/2}\\ & {\int }_0^{2\pi }{\int }_0^{2\pi }\frac{F}{(1+e\cos \varphi {)}^2(1+E\cos \Phi {)}^2}d\varphi d\Phi ,\end{array}$$

$$\begin{array}{rl}\frac{de}{d\tau }=& \frac{15\pi {\alpha }^{3}e\sqrt{1-e^2}{m}_3}{16A^3(1-E^2{)}^{3/2}(m_1+m_2)}({\sin }^2{\iota }_3(\cos 2\iota +3)\sin 2\omega \cos (2\Omega -2{\Omega }_3)+4{\sin }^2{\iota }_3\cos \iota \cos 2\omega \sin (2\Omega -2{\Omega }_3)\\ & -4\sin 2{\iota }_3\sin \iota \cos 2\omega \sin (\Omega -{\Omega }_3)-2\sin 2\iota \sin 2{\iota }_3\sin 2\omega \cos (\Omega -{\Omega }_3)+{\sin }^2\iota (3\cos 2{\iota }_3+1)\sin 2\omega ),\\ \frac{d\iota }{d\tau }=& \frac{3\pi {\alpha }^3{m}_3}{4A^3\sqrt{1-e^2}(1-E^2{)}^{3/2}(m_1+m_2)}(\sin \iota \sin {\iota }_3\cos (\Omega -{\Omega }_3)+\cos \iota \cos {\iota }_3)(\sin {\iota }_3\sin (\Omega -{\Omega }_3)(5e^2\cos 2\omega +3e^2+2)\\ & +5e^2\sin 2\omega (\sin {\iota }_3\cos \iota \cos (\Omega -{\Omega }_3)-\sin \iota \cos {\iota }_3)),\\ \frac{d\Omega }{d\tau }=& \frac{3\pi {\alpha }^3{m}_3}{4A^3\sqrt{1-e^2}(1-E^2{)}^{3/2}(m_1+m_2)}(\sin {\iota }_3\cos (\Omega -{\Omega }_3)+\cos {\iota }_3\cot \iota )(5e^2\sin {\iota }_3\sin 2\omega \sin (\Omega -{\Omega }_3)\\ & +(5e^2\cos 2\omega -3e^2-2)(\sin \iota \cos {\iota }_3-\sin {\iota }_3\cos \iota \cos (\Omega -{\Omega }_3))),\\ \frac{d\bar{\omega }}{d\tau }=& \frac{3\pi {\alpha }^3\sqrt{1-e^2}{m}_3}{8A^3(1-E^2{)}^{3/2}(m_1+m_2)}(10\sin \iota \sin 2{\iota }_3\sin 2\omega \sin (\Omega -{\Omega }_3)\\ & +{\sin }^2{\iota }_3\cos (2\Omega -2{\Omega }_3)(2{\sin }^2\iota (4-5{\cos }^2\omega )+20{\cos }^2\omega -10)-10{\sin }^2{\iota }_3\cos \iota \sin 2\omega \sin (2\Omega -2{\Omega }_3)\\ & +\sin 2\iota \sin 2{\iota }_3(3-5\cos 2\omega )\cos (\Omega -{\Omega }_3)+(3\cos 2{\iota }_3+1)({\sin }^2\iota (5{\cos }^2\omega -4)+1)),\\ \frac{dE}{d\tau }=& 0,\\ \frac{d{\iota }_3}{d\tau }=& -\frac{3\pi {\alpha }^{7/2}{m}_1{m}_2\sqrt{m_3}}{4A^{7/2}(1-E^2{)}^2(m_1+m_2{)}^{5/2}}(\sin {\iota }_3(\sin (2\Omega -2{\Omega }_3)({\sin }^2\iota (-5e^2{\cos }^2\omega +4e^2+1)+10e^2{\cos }^2\omega -5e^2)\\ & +5e^2\cos \iota \sin 2\omega \cos (2\Omega -2{\Omega }_3))+\cos {\iota }_3(\sin 2\iota \sin (\Omega -{\Omega }_3)(-5e^2{\cos }^2\omega +4e^2+1)-5e^2\sin \iota \sin 2\omega \cos (\Omega -{\Omega }_3))),\\ \frac{d{\Omega }_3}{d\tau }=& -\frac{3\pi {\alpha }^{7/2}{m}_1{m}_2\sqrt{m_3}\csc {\iota }_3}{8A^{7/2}(1-E^2{)}^2(m_1+m_2{)}^{5/2}}(\sin 2{\iota }_3(\frac{1}{2}{\sin }^2\iota (\cos (2\Omega -2{\Omega }_3)+3)(5e^2\cos 2\omega -3e^2-2)\\ & +5e^2\cos \iota \sin 2\omega \sin (2\Omega -2{\Omega }_3)-5e^2\cos 2\omega \cos (2\Omega -2{\Omega }_3)+3e^2+2)\\ & +\cos 2{\iota }_3(\sin 2\iota \cos (\Omega -{\Omega }_3)(5e^2\cos 2\omega -3e^2-2)-10e^2\sin \iota \sin 2\omega \sin (\Omega -{\Omega }_3))),\\ \frac{d{\bar{\omega }}_3}{d\tau }=& \frac{3\pi {\alpha }^{7/2}{m}_1{m}_2\sqrt{m_3}}{16A^{7/2}(1-E^2{)}^2(m_1+m_2{)}^{5/2}}(30e^2\sin \iota \sin 2{\iota }_3\sin 2\omega \sin (\Omega -{\Omega }_3)-30e^2{\sin }^2{\iota }_3\cos \iota \sin 2\omega \sin (2\Omega -2{\Omega }_3)\\ & +3{\sin }^2{\iota }_3\cos (2\Omega -2{\Omega }_3)({\sin }^2\iota (-5e^2\cos 2\omega +3e^2+2)+10e^2\cos 2\omega )\\ & +3\sin 2\iota \sin 2{\iota }_3\cos (\Omega -{\Omega }_3)(-5e^2\cos 2\omega +3e^2+2)+(2-3{\sin }^2{\iota }_3)({\sin }^2\iota (15e^2\cos 2\omega -9e^2-6)+6e^2+4))\end{array}$$

$$\begin{array}{rl}\frac{de}{d\tau }=& \frac{dE}{d\tau }=\frac{dp}{d\tau }=\frac{dP}{d\tau }=\frac{d{\iota }_3}{d\tau }\supset 0,\\ \frac{d\iota }{d\tau }\supset & \frac{3\pi a{\alpha }^{3/2}{m}_3\sin 2{\iota }_3\sin (\Omega -{\Omega }_3)}{2A^3(1-E^2{)}^{3/2}\sqrt{m_1+m_2}},\\ \frac{d\Omega }{d\tau }\supset & -\frac{\pi a{\alpha }^{3/2}{m}_3}{2A^3(1-E^2{)}^{3/2}\sqrt{m_1+m_2}}\\ & \times (-3\sin 2{\iota }_3\cot \iota \cos (\Omega -{\Omega }_3)+3\cos 2{\iota }_3+1),\\ \frac{d\omega }{d\tau }\supset & -\frac{3\pi a{\alpha }^{3/2}{m}_3\sin 2{\iota }_3\csc \iota \cos (\Omega -{\Omega }_3)}{2A^3(1-E^2{)}^{3/2}\sqrt{m_1+m_2}},\\ \frac{d{\Omega }_3}{d\tau }\supset & \frac{4\pi a{\alpha }^{3/2}{m}_3}{A^3(1-E^2{)}^{3/2}\sqrt{m_1+m_2}},\\ \frac{d{\omega }_3}{d\tau }\supset & -\frac{12\pi a{\alpha }^{3/2}{m}_3\cos {\iota }_3}{A^3(1-E^2{)}^{3/2}\sqrt{m_1+m_2}},\end{array}$$

Full text

The first direct detection of gravitational wave (GW) by the Laser Interferometric Gravitational-wave Observatory (LIGO) Abbott:2016blz is a historical landmark. It provides a new way to explore the Universe, and opens the era of gravitational wave astronomy.

Up to now, several GW events produced by the inspiral and subsequent merger of two black holes (BHs) or two neutron stars were reported in Abbott:2016nmj ; Abbott:2017vtc ; Abbott:2017gyy ; Abbott:2017oio ; TheLIGOScientific:2017qsa ; TheLIGOScientific:2016pea ; LIGOScientific:2018mvr . The signal of GW150914 was strong enough to be apparent, without using any waveform model, in the filtered detector strain data. But the other GW signals were so weak that the matched filtering with waveform templates was essential for the detection. The templates adopted by LIGO in searching for GW signals assumed the circular obits for the compact binaries Abbott:2017vtc ; Abbott:2017gyy ; Abbott:2017oio ; TheLIGOScientific:2017qsa ; TheLIGOScientific:2016pea ; LIGOScientific:2018mvr . This assumption is reasonable for the isolated binary systems because the angular momentum that gravitational waves carry away causes the orbits to circularize faster than they shrink Peters:1963ux .

The real environment where binary BHs live in are often complicated and they are mostly living in multiple systems. Hierarchical triple system is the simplest and stable one, where binary BHs might obtain large eccentricity perturbed by the third body through Kozai-Lidov mechanism Kozai:1962zz ; Lidov:1962zs . The formation channel of binary BHs roughly fall into two categories: they are from the remanent of evolvement of isolated stellar binaries in “field” chanel Voss:2003ep , or they are formed by “dynamically” channels through three body encounters in dense star clusters Sigurdsson:1993 . Triple system of stars are believed to be common in universe Raghavan:2010ds ; Tokovinin:2014ds ; Fuhrmann:2017df ; Tokovinin:2006jm ; Pribulla:2006gk , it is possible that triple BHs form in the end of the evolution of triple star systems. On the other hand, the dynamical formation requires a very high stellar density which are believed to exist within the cores of globular clusters (GCs), and the BHs being more massive than the average stars sink to the center of GC until the the majority of the BHs reside in the cluster core Spitzer:1969df . This “mass segregation” process guaranteed that three body encounters in the core can frequently occur Ivanova:2005mi , producing binary BHs at high rates Rodriguez:2015oxa .

Kozai-Lidov mechanism has an important secular effect in hierarchical triple systems, and plays an important role in dynamical evolution of triples. GWs emitted by highly eccentric orbits of compact binaries excited from the Kozai-Lidov mechanism might be detectable by LIGO and VIRGO Wen:2002km ; OLeary:2005vqo ; Mandel:2007hi , pulsar timing arrays Finn:2010ph ; AmaroSeoane:2009af ; Kocsis:2011ch , and also future space-based GW observatories such as LISA Robson:2018svj . Even though Kozai-Lidov mechanism has been extended to more general cases Katz:2011hn ; Naoz:2012bx ; Naoz:2010xi ; Will:2017vjc in the past decades, the effect of spin is absent in the literature. Actually the black holes in our Universe universally have spins. In particular, the observations mainly based on X-ray Reflection Spectroscopy Reynolds:2013qqa ; Reynolds:2013rva indicate that most of SMBHs especially with masses larger than $10^{7}M_{\odot}$ have large spins (e.g. the dimensionless spin parameter larger than 0.9).

In this letter, we will explore the secular evolution of a compact BH binary revolving around a spinning SMBH. At the leading order the effect of spin of SMBH can be interpreted by the gravitomagnetic force Nichols:2011pu which induces couplings between spin of SMBH and the inner and outer orbital angular momenta, and leads to an oscillating characteristic correction to the secular evolution of the compact binary eccentricity over a long time scale compared to the Kozai-Lidov dynamics caused by a non-spinning SMBH.

Let’s start with a hierarchical triple system illustrated in Fig. 1, where two binary BHs with mass $m_{1}$ and $m_{2}$ revolve around a spinning SMBH with mass $m_{3}(\gg m_{1},m_{2})$ and spin $a$. Here the spin of SMBH is taken to align the $Z$ direction in the fundamental reference frame whose origin is located at the position of SMBH. ${\bf{S}}=m_{3}a{\bf{e}}_{Z}$ denotes the spin vector of SMBH, and ${\bf{J}}_{\text{in}}$ and ${\bf{J}}_{\text{out}}$ are the orbital angular momentum of the inner and outer orbits respectively. The green lines are the lines of ascending nodes, and the blue lines are the lines of pericenters. The coordinates of these two small BHs are denoted by $\mathbf{r}_{\beta}$, $\mathbf{r}_{\beta\gamma}\equiv\mathbf{r}_{\beta}-\mathbf{r}_{\gamma}$, and $\mathbf{n}_{\beta\gamma}={\mathbf{r}_{\beta\gamma}/r_{\beta\gamma}}$ with $r_{\beta\gamma}\equiv|\mathbf{r}_{\beta\gamma}|$, where $\beta,\gamma=(1,2)$. The coordinate $\bf{L}$ of the mass center of the inner binary satisfies $m_{1}\mathbf{r}_{1}+m_{2}\mathbf{r}_{2}=m\bf{L}$, where $m=m_{1}+m_{2}$. Introducing $\bm{\ell}=\mathbf{r}_{12}$, we have $\mathbf{r}_{1}={m_{2}\over m}\bm{\ell}+\mathbf{L}$, $\mathbf{r}_{2}=-{m_{1}\over m}\bm{\ell}+\mathbf{L}$, and then the velocities of these two small black holes read

[TABLE]

where $\mathbf{V}={d\bm{L}/dt}$, $\mathbf{v}={d\bm{\ell}/dt}$.

The interactions induced by spin in binary or multiple systems have been well investigated in post Newtonian (PN) expansion in the weak field and low velocity limit in Thorne:1984mz ; Barker:1975ae ; DEath:1975wqz ; Barker:1979 ; Thorne:1980ru ; Poisson:1997ha . At the leading order the spin of SMBH makes a contribution to the equation of motion of the compact binary BHs in 1.5 PN which can be explained as the gravitomagnetic force Nichols:2011pu from the space-time of spinning SMBH. The accelerations of $m_{1}$ and $m_{2}$ BHs due to the gravitomagnetic force caused by the spin of SMBH take the form

[TABLE]

where $a\equiv J_{3}/m_{3}$ and $J_{3}$ is the angular momentum of black hole $m_{3}$. In order to figure out the leading order contribution to the dynamics of inner binary from spin of SMBH, we decompose the equations of motion of inner and outer orbits up to 1.5 PN as follows

[TABLE]

where $\mathbf{a}\equiv{d^{2}\bm{\ell}/dt^{2}}$, $\mathbf{A}\equiv{d^{2}\bm{L}/dt^{2}}$. Here $\mathbf{a}^{[\text{N}]}$ and $\mathbf{A}^{[\text{N}]}$ are the Newtonian accelerations for the inner and outer orbits, namely

[TABLE]

Up to the quadrupole order, we have

[TABLE]

where $\bf{n}={\bm{\ell}/\ell}$, ${\bf{N}}={{\bf{L}}/L}$, and $N_{n}=\bf{N}\cdot\bf{n}$. From Eq. (2), the leading contributions to the accelerations of the relative motion and center of mass of inner binary BHs from the spin of SMBH are given by

[TABLE]

where the higher order corrections of ${\cal O}({\ell/L})$ are ignored.

The inner and outer orbits are dominated by the Keplerian orbits, and the acceleration in (3) and (4) are dominated by $a^{[\text{N}]}$ and $A^{[\text{N}]}$ respectively. The inner orbital elements $p,e,\omega,\Omega,\iota$ are defined as

[TABLE]

and the semi-major axes of inner orbit is $\alpha=p(1-e^{2})$. Similarly, for the outer orbit, we have

[TABLE]

and the semi-major axes of outer orbit is $A=P(1-E^{2})$.

According to Eqs. (3), (4), (7) and (8), we introduce the perturbing accelerations $\delta\mathbf{a}=\mathbf{a}+{m\over\ell^{2}}\mathbf{n}$ and $\delta\mathbf{A}=\mathbf{A}+{m_{3}\over L^{2}}\mathbf{N}$. The dynamical evolution of the inner orbital elements are govern by the following equations of motion

[TABLE]

where $\mathbf{h}\equiv\bm{\ell}\times\mathbf{v}=\sqrt{mp}\hat{\mathbf{h}}$, $\mathbf{R}$ is the Runge-Lenz vector defined by $\mathbf{R}\equiv{\mathbf{v}\times\mathbf{h}/m}-\mathbf{n}=e(\cos{\phi}\mathbf{n}-\sin{\phi}\bm{\lambda})$. Then the evolution of the orbital elements can be obtained by solving the above two equations, namely

[TABLE]

where $\mathcal{R}=\mathbf{n}\cdot\delta\mathbf{a}$, $\mathcal{S}=\bm{\lambda}\cdot\delta\mathbf{a}$, $\mathcal{W}=\hat{\mathbf{h}}\cdot\delta\mathbf{a}$, $\bar{\omega}$ is defined by $\dot{\omega}=\dot{\bar{\omega}}-\dot{\Omega}\cos{\iota}$. Similarly, the outer orbital elements evolution can be expressed by replacing all the elements of inner binary by the the elements of outer binary, like $e\to E,p\to P,m\to m_{3},\phi\to\Phi$, $\dot{\omega}_{3}=\dot{\bar{{\omega}}}_{3}-\dot{\Omega}_{3}\cos{\iota_{3}}$, etc., and $\mathcal{R}_{3}=\mathbf{N}\cdot\delta\mathbf{A}$, $\mathcal{S}_{3}=\mathbf{\Lambda}\cdot\delta\mathbf{A}$, $\mathcal{W}_{3}=\hat{\mathbf{H}}\cdot\delta\mathbf{A}$.

The secular evolution of orbital elements are double averaged results through integrating over one orbital period of both inner and outer orbits as follows

[TABLE]

where $T_{\text{in}}$ and $T_{\text{out}}$ are the periods of inner and outer orbits. For convenience, the time integration can be replaced by phase angle integration according to the relation of $dt=\sqrt{p^{3}/m}(1+e\cos{\phi})^{-2}d\phi$ and $dt^{\prime}=\sqrt{P^{3}/m_{3}}(1+E\cos{\Phi})^{-2}d\Phi$, and thus the average in Eq. (17) becomes

[TABLE]

For simplicity, we also convert the time derivation ${d/dt}$ to a dimensionless one ${d/d\tau}$ by rescaling time compared to the inner orbital period with $\tau\equiv t/T_{\text{{in}}}={t\over 2\pi}\sqrt{m\over{\alpha}^{3}}$.

The Newtonian quadrupole perturbing accelerations in Eqs. (7) and (8) result in the well-known Kozai-Lidov effect as follows

[TABLE]

Here both the inner orbital angular momentum $\mathbf{J}_{\text{in}}={m_{1}m_{2}\over m}\sqrt{pm}\hat{\mathbf{h}}$ and outer orbital angular momentum $\mathbf{J}_{\text{out}}={m}\sqrt{Pm_{3}}\hat{\mathbf{H}}$ can be in arbitrary directions, and the above results come back to the familiar Kozai-Lidov formular when the total orbital angular momentum $\mathbf{J}=\mathbf{J}_{\text{in}}+\mathbf{J}_{\text{out}}$ is along the $Z$ direction. According to Eqs. (9) and (10), the spin of SMBH $\mathbf{S}$ couples with $\mathbf{J}_{\text{in}}$ and $\mathbf{J}_{\text{out}}$ in 1.5 PN at the leading order. The perturbing forces in Eqs. (9) and (10) from the leading order effect of spin contribute to the double averaged results of secular evolution of orbital elements as follows

[TABLE]

where “$\supset$” indicates the contribution from the spin only. It apparently seems that the spin of SMBH does not lead to the secular evolution of the eccentricity of the compact binary. However, Eqs. (20) couples to the those in Eqs. (19), and the spin of SMBH influence the eccentricity of the compact binary through the effects on the other orbital elements, such as $\iota$, $\omega$, $\Omega-\Omega_{3}$ and etc., which are closely related to the eccentricity $e$ of the inner orbit in Eqs. (19). For a numerical example, see Fig. 2.

Fig. 2 indicates that the spin of SMBH leads to a correction to the evolution of inner orbital eccentricity in both the amplitude and phase. Furthermore, such a correction oscillates over a long time scale.

In this letter we sketch out the leading order contribution to the secular evolution of the compact binary revolving around a spinning SMBH. The leading order interactions induced by the spin of SMBH appear at the 1.5 PN through coupling between spin and the inner and outer orbital angular momenta. The relative angle between the inner orbital longitude of ascending node and the outer orbital longitude of ascending node becomes dynamical. Through affecting the evolutions of some orbital elements, the spin of SMBH finally influence the secular evolution of the orbit eccentricity of compact binary. Compared to the Kozai-Lidov mechanism caused by a non-spinning SMBH, the spin of SMBH makes an oscillating characteristic contribution to the secular evolution of the compact binary eccentricity over a long time scale.

Finally, the features on the secular evolution of the orbit elements of compact binary from the spin of SMBH leaves some fingerprints in the GW waveform which might be used to measure the spin of SMBH in the future.

Acknowledgments. This work is supported by grants from NSFC (grant No. 11690021, 11575271, 11747601), the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDB23000000), Top-Notch Young Talents Program of China, and Key Research Program of Frontier Sciences of CAS.