A Gravitational Wave Background from Primordial Black Hole Lattices in Matter Dominated Era

TL;DR

This paper investigates gravitational wave production from primordial black hole lattices during matter domination using numerical simulations and analytical estimates, highlighting the role of black hole fraction in their evolution.

Contribution

It introduces a detailed simulation approach to study primordial black hole lattices and their gravitational wave emissions during matter-dominated era.

Findings

Black hole fraction significantly influences PBHL evolution

Primordial black holes move at relativistic speeds during expansion

Gravitational wave emission is confirmed through simulations and analysis

Abstract

We use the wide-used \textsf{Einstein Toolkit} to solve the Einstein constraints and then simulate the expansion of primordial black hole lattices (PBHLs) with different value of $f_{\mathrm{PBH}}$ and $m_{\mathrm{PBH}}$. We find that $f_{\mathrm{PBH}}$ plays an important role during the evolution of PBHLs. Since the motion of primordial black holes (PBHs) caused by the expansion of PBHLs occurs at speeds close to that of light, we expect the emission of gravitational waves (GWs) during the expansion of PBHLs. We use both analytical estimates and numerical simulations to cross check the production of GWs in expanding PBHLs.

Figures

Tables

Table 1. Table 1: The parameters setting of ten PBHLs.

Lattices	$f_{\mathrm{PBH}}$	$m_{\mathrm{PBH}}$	${\rho }_m$	${\rho }_t$	$L$	$\sigma$	$l=\frac{\mathrm{PBH}}{2}$
L0	$100\%$	$2$	$0$	$0.002$	$5$	$3.5$	$1$
L1	$75\%$	$1.5$	$0.0005$	$0.002$	$5$	$3.75$	$0.75$
L2	$50\%$	$1$	$0.001$	$0.002$	$5$	$4$	$0.5$
L3	$25\%$	$0.5$	$0.0015$	$0.002$	$5$	$4.25$	$0.25$
L4	$80\%$	$2$	$0.0005$	$0.0025$	$5$	$3.5$	$1$
L5	$66.7\%$	$2$	$0.001$	$0.003$	$5$	$3.5$	$1$
L6	$62.5\%$	$2$	$0.0015$	$0.0035$	$5$	$3.5$	$1$
L7	$100\%$	$1.5$	$0$	$0.0015$	$5$	$3.75$	$0.75$
L8	$100\%$	$1$	$0$	$0.001$	$5$	$4$	$0.5$
L9	$100\%$	$0.5$	$0$	$0.0005$	$5$	$4.25$	$0.25$

Table 2. Table 2: The solution of the Einstein constraints of ten PBHLs with spacing of the finest one of 5 5 5 refinement levels equal to 0.0625 0.0625 0.0625 .

Lattices	$f_{\mathrm{PBH}}$	$m_{\mathrm{PBH}}$	$K_c$	$D_{\mathrm{edge}}\left(\tau =0\right)$	$H_{\max }$	$H_{\min }$
L0	$100\%$	$2$	$-0.2480$	$14.16$	$0.0023$	$-0.0009$
L1	$75\%$	$1.5$	$-0.3413$	$12.03$	$0.0032$	$-0.0021$
L2	$50\%$	$1$	$-0.3864$	$10.76$	$0.0168$	$-0.0114$
L3	$25\%$	$0.5$	$-0.4000$	$10.07$	$0.2438$	$-0.2868$
L4	$80\%$	$2$	$-0.3711$	$12.63$	$0.0047$	$-0.0021$
L5	$66.7\%$	$2$	$-0.4869$	$11.21$	$0.0087$	$-0.0041$
L6	$62.5\%$	$2$	$-0.6200$	$9.84$	$0.0178$	$-0.0060$
L7	$100\%$	$1.5$	$-0.2338$	$13.12$	$0.0029$	$-0.0020$
L8	$100\%$	$1$	$-0.2120$	$12.07$	$0.0168$	$-0.0114$
L9	$100\%$	$0.5$	$-0.1698$	$11.01$	$0.2436$	$-0.2868$

Equations

$$R+K^2-K_{ij}{K}^{ij}$$

$$D_j{K}_i^j-D_{i}K$$

$$K_{ij}$$

$${\gamma }_{ij}$$

$$A_{ij}$$

$${\tilde{D}}_i{\tilde{D}}^i\Psi -\frac{1}{8}\tilde{R}\Psi +\frac{1}{8}{\hat{A}}_{ij}{\hat{A}}^{ij}{\Psi }^{-7}+2\pi E{\Psi }^5-\frac{1}{12}{K}^2{\Psi }^5$$

$${\tilde{D}}_j{\hat{A}}^{ij}-\frac{2}{3}{\Psi }^6{\tilde{D}}^{i}K$$

$${\hat{A}}^{ij}$$

$${\tilde{D}}_j{\tilde{A}}^{ij}$$

$${\tilde{\gamma }}^{ij}$$

$${\hat{A}}_{\text{TT}}^{ij}$$

$$E$$

$$p^i$$

$$K$$

$$△\psi +△\left(\frac{m_{\mathrm{PBH}}}{2r}W(r)\right)-\frac{1}{12}{K}^2{\Psi }^5+\frac{1}{8}{\hat{A}}_{ij}{\hat{A}}^{ij}{\Psi }^{-7}+2\pi {\rho }_m{\Psi }^5$$

$$△X^i+\frac{1}{3}{\partial }^i{\partial }_j{X}^j-\frac{2}{3}{\Psi }^6{\partial }^{i}K$$

$$\Psi =\psi +\frac{m_{\mathrm{PBH}}}{2r}(1-W(r))$$

$$2\pi m_{\mathrm{PBH}}-\frac{1}{12}{K}_c^2\int W^2{\Psi }^5+\frac{1}{8}\int {\hat{A}}_{ij}{\hat{A}}^{ij}{\Psi }^{-7}+2\pi \int {\rho }_m{\Psi }^5=0,$$

$$\psi (0)$$

$$X^i(0)$$

$$\tau (t,x)={\int }_0^t\alpha (t^{\prime },x)d{t}^{\prime },$$

$${\partial }_t\alpha -{\beta }^i{\partial }_i\alpha =-2K\alpha .$$

$$D(\tau )={\int }_{{\gamma }_{\tau }}[(-{\alpha }^2(\tau ,\ell )+{\beta }^2(\tau ,\ell ))({\partial }_{\ell }t{)}^2+2{\beta }_i(\tau ,\ell ){\partial }_{\ell }t{\partial }_{\ell }{x}^i+{\gamma }_{ij}(\tau ,\ell ){\partial }_{\ell }{x}^i{\partial }_{\ell }{x}^j{]}^{1/2}d\ell ,$$

$$E$$

$${\tilde{T}}^{ij}(\omega ,\omega \hat{\text{n}})=2\pi \delta (\omega -{\omega }_0)\frac{-i{m}_{\mathrm{PBH}}}{\omega }\left[\frac{\gamma (\tau )v^i(\tau )v^j(\tau )}{1-v(\tau )\cos \theta }-\frac{\gamma (\tau +\Delta \tau )v^i(\tau +\Delta \tau )v^j(\tau +\Delta \tau )}{1-v(\tau +\Delta \tau )\cos \theta }\right],$$

$${\Lambda }_{ij,kl}(\hat{\text{n}})={\delta }_{ik}{\delta }_{jk}-\frac{1}{2}{\delta }_{ij}{\delta }_{kl}-n_j{n}_l{\delta }_{ik}-n_i{n}_k{\delta }_{jl}+\frac{1}{2}{n}_k{n}_l{\delta }_{ij}+\frac{1}{2}{n}_i{n}_j{\delta }_{kl}+\frac{1}{2}{n}_i{n}_j{n}_k{n}_l,$$

Full text

I Introduction

Observations of gravitational waves (GWs) by LIGO/VIRGO Abbott:2016blz ; Abbott:2016nmj ; Abbott:2017vtc ; Abbott:2017oio ; TheLIGOScientific:2017qsa have proved there should be a gravitational wave background (GWB) produced by the black hole and neutron star binaries’ coalescence. Meanwhile, people also believe there is a primordial inflationary GWB from tensor perturbations generated by quantum fluctuations during inflation even though Planck didn’t measure it Ade:2015xua . Besides that, there are many other theories (or sources) which can produce GWBs, such as reheating after inflation GarciaBellido:2007af , thermal phase transitions from the decays of cosmic strings Maggiore:1999vm . Here we expect that there is a GWB from primordial black hole lattices (PBHLs) in matter dominated era.

To what extend can the configuration of primordial black holes (PBHs) just after the epoch of matter-radiation equality be considered as a PBHL? As we known, a pair of PBHs would decouple from the expansion of the Universe and form a gravitationally bound system if the average energy density of PBHs over the volume $R^{3}$ is larger than the total background cosmic energy density $\rho_{t}$, that is, if $m_{\mathrm{PBH}}R^{-3}>\rho_{t}$. Since $\rho_{t}\approx m_{\mathrm{PBH}}\bar{R}^{-3}f_{\mathrm{PBH}}^{-1}$, we know that most of PBHs are gravitationally free on average in matter dominated era. On the other hand, due to the cosmological principle that the Universe is homogeneous and isotropic on large scales, we can say that the PBHL is the most reasonable configuration of PBHs with a monochromatic mass distribution on our hands.

In fact, there is a cosmological model, black hole lattices (BHLs) Lindquist ; Yoo:2012jz ; Yoo:2013yea ; Bentivegna:2013jta , to investigate the so-called “backreaction” that the local inhomogeneities can affect the global expansion of the Universe. For review of BHLs, see Bentivegna:2018koh . It’s worth pointing out that the PBHLs doesn’t serve as a cosmological model like BHLs does but is an early epoch of the Universe. In practice, however, we do apply the same technologies developed by Yoo:2012jz ; Yoo:2013yea ; Bentivegna:2013jta when they evolve the BHLs in numerical relativity to the PHBLs.

Here we will turn to the wide-used Einstein Toolkit Loffler:2011ay to solve the Einstein constraints and simulate the spacetime and hydrodynamical evolution. More precisely: the thorn CT_MultiLevel (and CT_Analytic) Bentivegna:2013xna gets the initial data of PBHLs; the thorn McLachlan Brown:2008sb ; Reisswig:2010cd ; code evolves spacetime using the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism Baumgarte:1998te ; Shibata:1995we ; Alcubierre:2000xu ; the thorn GRHydro evolves the hydrodynamical system Moesta:2013dna ; Baiotti:2004wn ; Hawke:2005zw .

This paper is organized as follows. In Subsec. II.1, we give the Einstein constraints of PBHLs with dust. In Subsec. II.2, we solve the Einstein constraints of PBHLs with different value of $f_{\mathrm{PBH}}$ and $m_{\mathrm{PBH}}$. In Subsec. III.1, we show the expansion of PBHLs with different value of $f_{\mathrm{PBH}}$ and $m_{\mathrm{PBH}}$. In Subsec. III.2, we use both analytical estimates and numerical simulations to derive the production of GWB in PBHLs. At last, a brief summary and discussion are included in Sec. IV.

In this paper, we adopt the following conventions: Greek indices run in {0, 1, 2, 3}, Latin indices run in {1, 2, 3} and repeated indices implies summation and we are in a geometric unit system with $G=c=M_{\odot}=1$.

II Initial Data of Primordial Black Hole Lattices

To mock an infinite 3-dimensional lattice in our simulations, we impose the periodic boundary conditions on a cubic with $x^{i}$ in $[-5,5]$. A PBH is located at the center of cubic and surrounded by flat-distributed dust with energy density $\rho_{m}$. It’s worth noting that our simulation is not done in Schwarzschild coordinates but in its isotropic coordinates. For simplicity, we just add the dust outside the PBH horizon $r_{h}=0.5m_{\mathrm{PBH}}$ so that the inner boundary conditions of PBHLs can be set as that of a single puncture Brandt:1997tf . And there are no outer boundary conditions for PBHLs. Although the set up of initial data is usually a non-trivial task, here we can solve PBHLs’ coupled system of one non-linear and three linear elliptic partial differential equations thanks to Yoo:2012jz ; Bentivegna:2013jta ; Bentivegna:2013xna .

II.1 Einstein Constraints of Primordial Black Hole Lattices

The initial data of PBHLs must satisfy their Hamiltonian constraint and momentum constraint

[TABLE]

where the 3-metric $\gamma_{ij}$ is the intrinsic metric, $K_{ij}$ is the extrinsic curvature, $R$ is 3-Ricci scalar, $D_{i}$ is the covariant derivative associated with $\gamma_{ij}$, $E$ is the matter energy density and $p_{i}$ is the matter momentum density as measured by the Eulerian observer. After conformal decomposition of $\gamma_{ij}$ and $K_{ij}$ with the conformal factor $\Psi$

[TABLE]

we can get the conformal Einstein constraints

[TABLE]

where $\tilde{D}_{i}$ is the covariant derivative associated with the conformal metric $\tilde{\gamma}_{ij}$. Adopting the conformal transverse traceless method

[TABLE]

and choosing the free data as following

[TABLE]

we can get the conformal transverse traceless form of Einstein constraints for PBHLs

[TABLE]

where we have expanded the conformal factor as

[TABLE]

according to the inner boundary conditions of a single puncture.

II.2 Solving the Einstein Constraints

We use a multigrid approach Bentivegna:2013xna to solving the Einstein constraints (II.1). But we impose the integrability condition as Yoo:2012jz

[TABLE]

at the end of each relaxation step to determine the negative parameter $K_{c}$. Finally we perform the inner resetting

[TABLE]

at the end of every relaxation step.

Tab. 1 shows the parameters setting of ten PBHLs and Tab. 2 shows the solution of the Einstein constraints of ten PBHLs when all of $5$ refinement levels cover the whole domain with spacing $1$, $0.5$, $0.25$, $0.125$, and $0.0625$ respectively. Comparing L0 to L7, L8 and L9, we can see that a smaller PBH mass $m_{\mathrm{PBH}}$ gives a shorter initial proper cubic edge $D_{\mathrm{edge}}(\tau=0)$ and leads a lower initial expansion rate $-K_{c}$. Comparing L0 to L4, L5 and L6, we can see that a higher matter energy density $\rho_{m}$ has a higher initial expansion rate $-K_{c}$ but produces a shorter initial proper cubic edge $D_{\mathrm{edge}}(\tau=0)$. Comparing L0 to L1, L2 and L3, we can see that a larger PBH mass $m_{\mathrm{PBH}}$ and a smaller matter energy density $\rho_{m}$ still produce a larger initial proper cubic edge $D_{\mathrm{edge}}(\tau=0)$ but lead a lower initial expansion rate $-K_{c}$. That is to say, the initial expansion rate $-K_{c}$ is more sensitive to the matter energy density $\rho_{m}$ than the PBH mass $m_{\mathrm{PBH}}$. From the Tab. 2, we can also see that a smaller PBH mass $m_{\mathrm{PBH}}$ violates the Hamiltonian constraint more severely.

In Fig. 1, the $\Psi=\psi>1$ at the boundary means that each PBH is gravitationally coupled to its neighbors; the boundary $\psi$ of L0 is larger than that of L9 just because the PBH with larger mass will be gravitationally coupled to its neighbors more tightly; the boundary $\psi$ of L0 is larger than that of L6, which means the matter surrounding the PBH suppresses the gravitational correlation among PBHs; the boundary $\psi$ of L0 is larger than that of L3, which is consistent with above two cases.

In Fig. 2, we give the solutions of $X_{1}$ of PBHLs.

In Fig. 3, we shows the initial Hamiltonian constraint violation of PBHLs.

The last plot in Fig. 1, Fig. 2 and Fig. 3 is a convergence test for solving the initial data and shows the $\psi$, $X_{1}$ and $H$ of L0 with spacing of the finest one of $5$ refinement levels equal to $0.0625$ (L$0_{f}$), $0.078125$ (L$0_{m}$) and $0.1041667$ (L$0_{c}$) respectively.

III Time Evolution of Primordial Black Hole Lattices

Given the initial data of PBHLs, we can simulate the time evolution of PBHLs. Due to the negative trace of the extrinsic curvature $K_{ij}$ around the boundary, we can expect an expansion of PBHLs during the time evolution. Meanwhile, since the motion of PBHs caused by the expansion of PBHLs occurs at speeds close to that of light, we also expect the emission of GWs during the time evolution.

III.1 Expansion of Primordial Black Hole Lattices

In order to circumvent the central puncture, here we will study the expansion of PBHLs through the rescaling of the proper length of lattice cells’ edge with proper time. The proper time at anywhere is given by integrating the corresponding lapse function $\alpha(t,x)$ from the beginning of simulation

[TABLE]

where the evolution of $\alpha$ is given by

[TABLE]

And the proper length measured in a cctk_delta_time (or a timestep) is then given by

[TABLE]

where $\beta_{i}(\tau,x)$ is the shift vector obeying a Gamma driver and $\gamma_{\tau}$ is a constant-$\tau$ edge parameterized by $\ell$.

Fig. 4 shows the expansion of four PBHLs with the smallest initial Hamiltonian constraint violation in each group, where a prime represents a derivative with respect to the proper time $\tau$. We can see that $L0$ and $L1$ evolve differently even though they share a same $\rho_{t}$. That is to say $f_{\mathrm{PBH}}$ will play an important role during the evolution of PBHLs. Comparing $L4$ to $L0$ and $L1$, we find that the PBHL with a larger $\rho_{t}$ no matter due to an extra $\rho_{m}$ or a larger $m_{\mathrm{PBH}}$ expands faster. There is no intersection point between the evolutions of $D_{\mathrm{edge}}(\tau)^{\prime}$, which means a PBHL with a larger $-K_{c}$ will keep expanding faster for ever. Finally, we also find that the motion of PBHs caused by the expansion of PBHLs does occur at speeds close to that of light.

The right plot in Fig. 4 shows the expansion of $L0$ and Fig. 5 shows the $L_{2}$ norms of the Hamiltonian and momentum constraint for $L0$ with spacing of the finest one of $5$ refinement levels equal to $0.0625$, $0.078125$ and $0.1041667$ respectively. They serve as a convergence test for evolving PBHLs.

III.2 Gravitational Waves in Primordial Black Hole Lattices

We have shown that the motion of PBHs caused by the expansion of PBHLs occurs at speeds close to that of light. Here we will use both analytical estimates and numerical simulations to cross check the production of GWs in expanding PBHLs.

III.2.1 Theoretical estimation

The total energy radiated by one PBH in an expanding PBHL between $\tau$ and $\tau+\Delta\tau$ can be considered as the gravitational radiation from an accelerated mass estimated by Maggiore:1900zz

[TABLE]

where the energy-momentum tensor is

[TABLE]

the Lambda tensor $\Lambda_{ij,kl}(\hat{\textbf{n}})$ is

[TABLE]

the scalar product of the direction of gravitational radiation $\hat{\textbf{n}}$ and the velocity of PBH v is $n_{i}v^{i}=v\cos\theta$, $\omega_{0}$ is the frequency of the gravitational radiation at $\tau$, $T=2\pi\delta(0)\approx\Delta\tau$ and $\gamma=(1-v^{2})^{-1/2}$. The term of $[v(\tau)-v(\tau+\Delta\tau)]^{2}$ predicts that, in our PBHLs, the total energy radiated when $\tau\lesssim 14$ is larger than that when $\tau\gtrsim 14$ because $|D_{\mathrm{edge}}(\tau\lesssim 14)^{\prime\prime}|$ is much larger than $|D_{\mathrm{edge}}(\tau\gtrsim 14)^{\prime\prime}|$ as show in Fig. 4. The term of $[v(\tau)+v(\tau+\Delta\tau)]^{2}$ predicts that when $\tau\lesssim 14$ the total energy radiated in L4 should be larger than that in L0 and L7 because $D_{\mathrm{edge}}(\tau)^{\prime}$ of L4 is larger than that of L0 and L7. Similarly, according to the initial distribution of $K=K_{c}W(r)$, the total energy radiated in PBHLs should not be uniform and the locations far away from the center are full of more radiation energy.

III.2.2 Simulation results

Fig. 6 shows the waveforms of GWs, as estimated by the Newman-Penrose scalar $\Psi_{4}(\tau)=h^{\prime\prime}_{+}-ih^{\prime\prime}_{\times}$, produced at several distances $r$ in expanding PBHLs with different $\rho_{t}$. For $\tau\lesssim 14$, there is an obvious production of GWs and we also find an obvious deceleration in Fig. 4; for $\tau\gtrsim 14$, the amplitude of waveforms decreases and there is a plateau for $D_{\mathrm{edge}}(\tau)^{\prime}$ in Fig. 4. That is to say, there are some gravitational potential energy among PBHs converted to the GWs radiation and Re$[\Psi_{4}^{2,0}]\propto|D_{\mathrm{edge}}(\tau)^{\prime\prime}|$. The amplitude of waveforms increases with $\rho_{t}$ increasing, which also means Re$[\Psi_{4}^{2,0}]\propto D_{\mathrm{edge}}(\tau)^{\prime}$ since the $D_{\mathrm{edge}}(\tau)^{\prime\prime}$ is not sensitive to $\rho_{t}$ as shown in Fig. 4. The amplitude of GWs increases with $r$ increasing, which means the expansion of PBHLs is not uniform and the initial distribution of $K=K_{c}W(r)$ keeps a faster expansion at positions far away from PBHs. All of above features are consistent with theoretical predictions.

The right plots in Fig. 6 show the waveforms of GWs produced at several distances $r$ in $L0$. The oscillations due to numerical error at $\tau\approx 35$ in $L0_{c}$ disappear in $L0_{f}$, which guarantees the other oscillations in Fig. 6 are the waveforms of GWs produced in expanding PBHLs.

IV Summary and discussion

In this paper, we first use the wide-used Einstein Toolkit to solve the Einstein constraints of PBHLs with different value of $m_{\mathrm{PBH}}$ and $f_{\mathrm{PBH}}$ which is determined by the surrounding flat-distributed dust with energy density $\rho_{m}$. From the solutions, we find that a smaller PBH mass $m_{\mathrm{PBH}}$ gives a shorter initial proper cubic edge $D_{\mathrm{edge}}(\tau=0)$ and leads a lower initial expansion rate $-K_{c}$; a higher matter energy density $\rho_{m}$ has a higher initial expansion rate $-K_{c}$ but produces a shorter initial proper cubic edge $D_{\mathrm{edge}}(\tau=0)$; a larger PBH mass $m_{\mathrm{PBH}}$ and a smaller matter energy density $\rho_{m}$ still produce a larger initial proper cubic edge $D_{\mathrm{edge}}(\tau=0)$ but lead a lower initial expansion rate $-K_{c}$; a smaller PBH mass $m_{\mathrm{PBH}}$ violates the Hamiltonian constraint more severely; the PBH with larger mass will be gravitationally coupled to its neighbors more tightly; the matter surrounding the PBH suppresses the gravitational correlation among PBHs.

And then we simulate the expansion of PBHLs with the smallest initial Hamiltonian constraint violation in each group. From the evolutions, we find that $f_{\mathrm{PBH}}$ plays an important role during the evolution of PBHLs; the PBHL with a larger $\rho_{t}$ no matter due to an extra $\rho_{m}$ or a larger $m_{\mathrm{PBH}}$ expands faster; a PBHL with a larger $-K_{c}$ will keep expanding faster for ever; the motion of PBHs caused by the expansion of PBHLs does occur at speeds close to that of light. Meanwhile, we use both analytical estimates and numerical simulations to cross check the production of GWs in expanding PBHLs and find that Re$[\Psi_{4}^{2,0}]\propto D_{\mathrm{edge}}(\tau)^{\prime}~{}D_{\mathrm{edge}}(\tau)^{\prime\prime}$.

Comparing with the black hole and neutron star binaries’ coalescence, there is an abnormal feature that the amplitude of GWs doesn’t decrease with radius but increase with radius in expanding PBHLs. Our explanation is that the final waveform of GWs at any location is the superposition of eight waveforms of GWs produced by eight adjacent PBHs in expanding PBHLs and every waveform of GWs obeys the theoretical estimation (13).

There is one caveat: the total energy density $\rho_{t}$ of the real Universe in matter dominated era is much smaller than our exaggerated initial values $\sim 10^{-3}$. As we known, in SI units, this value corresponds to $6.1727\cdot 10^{17}\,\mathrm{kg/m^{3}}$ hence $H^{2}=3.45\cdot 10^{47}\mathrm{km^{2}/s^{2}/Mpc^{2}}$. Although we can suppress $\rho_{t}$ through enlarging the volume of our cubic until a reasonable $\rho_{t}$ compared to the real Universe in dominated matter era, we can’t afford the accompanying high computational cost. Therefore, we wouldn’t give the forecast whether GWs detectors can detect such signals presented in our paper or not in future.

Acknowledgments We would like to thank Xiao Guo and You-Jun Lu for their helpful discussions and advices on this paper. This work is partly supported by the National Natural Science Foundation of China under grant No. 11690024, the Strategic Priority Program of the Chinese Academy of Sciences (Grant No. XDB 23040100).