Orbit design and thruster requirement for various constant-arm space mission concepts for gravitational-wave observation

TL;DR

This paper evaluates orbit design and thruster needs for various gravitational-wave space missions, highlighting fuel requirements and proposing solutions for maintaining constant arm lengths in different orbital configurations.

Contribution

It extends previous work by applying orbit and thruster analysis to B-DECIGO and DECIGO missions, assessing fuel needs and proposing alternatives for proof mass actuation.

Findings

Heliocentric orbits require manageable fuel for constant arm maintenance.

Geocentric orbits pose significant fuel mass challenges.

Alternative proof mass solutions are necessary for adequate dynamical range.

Abstract

In previous papers, we have addressed the issues of orbit design and thruster requirement for the constant arm versions of AMIGO (Astrodynamical Middle-frequency Interferometric Gravitational-wave Observatory) mission concept and for the constant arm GW (gravitational wave) mission concept of AIGSO (Atom Interferometric Gravitational-wave Space Observatory). In this paper, we apply similar methods to the orbit design and thruster requirement for the constant arm GW missions B-DECIGO and DECIGO, and estimate the yearly propellant requirements at the specific impulse Isp = 300 sec and Isp = 1000 sec. For the geocentric orbit options of B-DECIGO which we have explored, the fuel mass requirement is a concern. For the heliocentric orbit options of B-DECIGO and DECIGO, the fuel requirement to keep the arm equal and constant should be easily satisfied. Furthermore, we explore the thruster and propellant requirements for constant arm versions of LISA and TAIJI missions and find the fuel mass requirement is not a show stopper either. The proof mass actuation noise is a concern. To have enough dynamical range, an alternate proof mass is required. Detailed laboratory study is warranted

Equations

$$\begin{array}{rl}{X}_k& =R(\cos {\psi }_k+e)\cos ϵ\\ Y_k& =R\sqrt{1-e^2}\sin {\psi }_k(k=1,2,3)\\ Z_k& =R(\cos {\psi }_k+ϵ)\sin ϵ\end{array}$$

$${\psi }_k+e\sin {\psi }_k=\Omega (t-t_0)-(k-1)\frac{2\pi }{3},$$

$$\begin{array}{rl}{x}_k& =X_k\cos \left[\frac{2\pi }{3}(k-1)+{\phi }_0\right]-Y_k\sin \left[\frac{2\pi }{3}(k-1)+{\phi }_0\right],\\ y_k& =X_k\cos \left[\frac{2\pi }{3}(k-1)+{\phi }_0\right]+Y_k\sin \left[\frac{2\pi }{3}(k-1)+{\phi }_0\right],\\ z_k& =Z_k,\end{array}$$

$$\begin{array}{rl}{\mathbf{r}}_{S/Ck}& =[x_k,y_k,z_k](k=1,2,3)\end{array}$$

$$\begin{array}{rl}{\mathbf{n}}_{23}(t)& =\frac{{\mathbf{r}}_{S/C3}-{\mathbf{r}}_{S/C2}}{|{\mathbf{r}}_{S/C3}-{\mathbf{r}}_{S/C2}|},\\ {\mathbf{n}}_{21}(t)& =\frac{{\mathbf{r}}_{S/C1}-{\mathbf{r}}_{S/C2}}{|{\mathbf{r}}_{S/C1}-{\mathbf{r}}_{S/C2}|},\\ {\mathbf{n}}_z(t)& =\frac{{\mathbf{n}}_{23}\times {\mathbf{n}}_{21}}{|{\mathbf{n}}_{23}\times {\mathbf{n}}_{21}|},\end{array}$$

$$\begin{array}{rl}{\mathbf{n}}_1& ={\mathbf{n}}_{23},\\ {\mathbf{n}}_2& ={\mathbf{n}}_z\times {\mathbf{n}}_1.\end{array}$$

$$\begin{array}{rl}{\mathbf{r}}_{\mathrm{traj},S/C1}& ={\mathbf{r}}_{S/C1},\\ {\mathbf{r}}_{\mathrm{traj},S/C2}& ={\mathbf{r}}_{S/C1}-b{\mathbf{n}}_1-(a+c){\mathbf{n}}_2,\\ {\mathbf{r}}_{\mathrm{traj},S/C3}& ={\mathbf{r}}_{S/C1}+b{\mathbf{n}}_1-(a+c){\mathbf{n}}_2,\end{array}$$

$$\begin{array}{rl}{\mathbf{r}}_{\mathrm{traj},S/C1b}& ={\mathbf{r}}_{S/C1}-2c{\mathbf{n}}_2,\\ {\mathbf{r}}_{\mathrm{traj},S/C2b}& ={\mathbf{r}}_{S/C1}-b{\mathbf{n}}_1+(a-c){\mathbf{n}}_2,\\ {\mathbf{r}}_{\mathrm{traj},S/C3b}& ={\mathbf{r}}_{S/C1}+b{\mathbf{n}}_1+(a-c){\mathbf{n}}_2.\end{array}$$

$${\mathbf{a}}_{\mathrm{traj}}={\ddot{\mathbf{r}}}_{\mathrm{traj}}.$$

$${\mathbf{a}}_{\mathrm{eph}}({\mathbf{r}}_{\mathrm{traj}},{\dot{\mathbf{r}}}_{\mathrm{traj}})={\mathbf{a}}_{\mathrm{Newton}}+{\mathbf{a}}_{1PN}+{\mathbf{a}}_{\mathrm{fig}}+{\mathbf{a}}_{\mathrm{asteroid}},$$

$${\mathbf{a}}_{\mathrm{thruster}}={\mathbf{a}}_{\mathrm{traj}}-{\mathbf{a}}_{\mathrm{eph}}.$$

$$F_{\mathrm{thruster}}(t)=g_0{I}_{\mathrm{sp}}\dot{m}(t),$$

$$\Delta m=M\left[1-\exp \left(-\frac{1}{g_0{I}_{\mathrm{sp}}}{\int }_0^T{a}_{\mathrm{thruster}}dt\right)\right].$$

$$\Delta m\approx \frac{1}{g_0{I}_{\mathrm{sp}}}\int F_{\mathrm{thruster}}dt\simeq \frac{M}{g_0{I}_{\mathrm{sp}}}{\int }_0^T{a}_{\mathrm{thruster}}dt,$$

$${\overline{\Delta m}}_{\mathrm{yr}}=\Delta m\times \frac{1\mathrm{yr}}{T}$$

Full text

{history}

\ccode

PACS numbers: 04.80.Nn, 04.80.-y, 95.10.Eg, 95.55.Ym

1 Introduction

O1 and O2 observations of advanced LIGO and advanced Virgo have detected 10 GW (gravitational wave) events from stellar-size binary black hole mergers together with a binary neutron stars coalescence[1, 2, 3]. The frequencies of these GW events are in the high-frequency band. Efforts of observations have been made in all other frequency bands – from below Hubble frequency to over terahertz frequency also[4].

With the successful demonstration of drag-free technology by LISA Pathfinder[5, 6], the technology for space GW detection is becoming mature in low-frequency and middle-frequency GW bands. In the low-frequency band, LISA is under development[7]. In China, two low-frequency GW mission proposals – TAIJI[8] and TianQin[9] are under actively study. TAIJI is a mission concept in heliocentric orbit[10]. TianQin is a mission concept in geocentric orbit[11]. Other low or middle frequency space detection methods under conceptual study are: AIGSO (Atom Interferometric Gravitational-wave Space Observatory)[12, 13], AMIGO (Astrodynamical Middle-frequency Interferometric GW Observatory)[14, 15], ASTROD-GW[16, 17, 18, 19, 20, 21, 22, 23], BBO[24], B-DECIGO[25, 26], DECIGO[26, 27], Super-ASTROD[28], other AI (atom interferometry) proposals[29, 30, 31], and optical clock tracking proposals[32, 33, 34, 35]. In the middle frequency band, there are a few ground-based proposals — MIGA[36, 37], SOGRO[38, 39], TOBA[40, 41], and ZAIGA[42].

To have significant sensitivity in the frequency band 0.1–10 Hz and yet to be a first-generation candidate for space GW missions, we have proposed a middle-frequency GW mission AMIGO[14, 15]. The mission concept is to use three drag-free spacecraft to form a triangular formation with nominal arm length 10,000 km, the first-generation TDI (time delay interferometry), laser power 2-10 W and telescope diameter 300-360 mm[14]. The targeting sensitivity in the middle frequency band is $3\times 10^{-21}$ Hz*-1*. Four options of orbits have been studied: (i) Earth-like heliocentric orbits (3-20 degrees behind the Earth); (ii) 600,000 km high orbit formation around the Earth; (iii) 100,000 km-250,000 high orbit formation around the Earth; (iv) near Earth-Moon L4 (or L5) halo orbit formation. All four options have LISA-like formations, that is the triangular formation is nearly 60∘ inclined to the orbit plane. In 2017, we proposed this mission concept with the aim of either as a stand-alone mission or as a pathfinder with only two spacecraft. The heliocentric formation is easy to obtain; along with it, the deployment method is also obtained. When we looked deeper into orbit simulation for the geocentric orbits, we found that it was difficult to find a hoped triangular formation[15]. Since the arm length is small, we tried to use thruster to keep the arms at constant equal length, and found that the fuel requirement make the geocentric AMIGO orbits technologically not feasible currently, but perfectly feasible for the heliocentric orbits. Therefore, for the heliocentric mission orbit choice of AMIGO, there can be both geodesic implementation and constant arm implementation of orbit configuration. In the AMIGO constant arm implementation for heliocentric orbit, the acceleration to maintain the formation can be designed to be less than 15 nm/s2 and the thruster requirement can be smaller than 15 $\mu$N[15].

AIGSO is a mission concept using atom interferometry to detect the GWs mainly in the middle frequency band (0.1–10 Hz)[12]. AIGSO proposes to have three spacecraft in linear formation with 10 km baseline. The three spacecraft maintain 5 km + 5 km constant arm-length formation. In a previous paper, we have addressed the issue of orbit design and thruster requirement for the constant arm AIGSO mission concept. The acceleration to maintain the formation can be designed to be less than 30 pm/s2 and the required amplitude of thruster force will be smaller than 30 nN[13].

B-DECIGO of arm length 100 km and DECIGO of arm length 1000 km are constant arm GW missions due to its Fabry-Perot implementation. Since arm length of B-DECIGO is shorter than AMIGO by 2 orders of magnitude, we consider three shrunk orbit configurations of AMIGO as the possible orbits of B-DECIGO. By using these downscaled orbits, we estimate the propellant requirement to maintain the constant arm-length for B-DECIGO mission concept.

Our paper is organized as follows. In Section 2, we briefly describe the algorithm to attain the selected orbits. In Section 3, we expatiate the method to obtain the trajectories of the S/Cs to maintain the constant arm interferometry, and to estimate the acceleration requirement for the thrusters. In Section 4.1, we obtain the thruster acceleration requirement for the B-DEICGO at different orbital configurations. In Section 4.2, we find that the acceleration to maintain the formation can be designed to be less than 0.2 nm/s2 and the thruster requirement be less than 0.2 nN for DECIGO. As far as these requirements are concerned, DECIGO is perfect feasible. In Section 4.3, we explore a new scheme of formation control for LISA and TAIJI. In this scheme of control, a constant-arm equilateral triangle is formed by one geodesic spacecraft followed by two non-geodesic spacecraft instead of three geodesic motion spacecraft forming a nearly equilateral triangle. In Section 5, we obtain the propellant requirement for various missions at specific impulse 300 sec and 1000 sec, and compile the results in Table 5. We present feasibility discussions and deliberate on the implementation requirements in Section 6.

2 Mission Orbit Selection

2.1 Orbit Selection

In this subsection, we briefly state the initial orbit selection algorithm used in our previous works \refcitewang&ni2017,dnw2013,wangni2013CQG, which is specifically discussed in \refciteDhurandhar+etal+2005. For the B-DECIGO configurations, as we will describe in the following section, their orbits are downscaled from the selected AMIGO configurations. The initial orbit configuration of the DECIGO used in this work are generated from the Eqs. (1)-(4) below without any optimization. In subsection 4.3, The initial LISA orbit configuration is picked from Ref. \refcitewang&ni2017, while the orbit of TAIJI is recalculated and re-optimized according to the method of Ref. \refcitewang&ni2017 for the configuration ahead of Earth by $20^{\circ}$.

For a LISA-like orbit configuration which has a $60^{\circ}$ inclination with respect to the ecliptic plane, the orbit of each S/C has an eccentricity $e$ and inclination $\iota$. The first order parameter is $\alpha\ [=\iota=l/(2R)]$, where $l$ is the nominal arm length and $R$ is the orbital radius. A set of S/C initial conditions in the heliocentric elliptical coordinate is defined in \refciteDhurandhar+etal+2005,

[TABLE]

where $\epsilon\simeq 3.34\times 10^{-9}\times l/$km; orbital eccentricity $e\simeq 1.93\times 10^{-9}\times l$/km; $R=1$ AU; and $\psi_{k}$ is the eccentric anomaly which is related to the mean anomaly $\Omega(t-t_{0})$ by equation

[TABLE]

where $\Omega$ is 2$\pi$/(one sidereal year). The $x_{k},y_{k},z_{k}(k=1,2,3)$ is defined as

[TABLE]

where $\varphi_{0}=\psi_{E}-\theta$ and $\psi_{E}$ is the position angle of Earth with respect to the X-axis at initial time in the ecliptic plane. The initial positions of the S/Cs in the heliocentric coordinate are

[TABLE]

2.2 Vector Defination in the Instantaneous Plane

Follow the procedures in Section 2.1, we can obtain the geodesic orbits of the three S/Cs. To identify the instantaneous plane formed by the three S/Cs, the direction $\mathbf{n}_{z}$ of the plane can be defined from the instantaneous positions of the S/Cs by

[TABLE]

where $\mathbf{r}_{\mathrm{S/C}i}$ is the instantaneous positions of the S/C$i$ ($i$=1,2,3) at $t$.

From the vectors obtained in Eq. (5), we define the two orthogonal unit vectors in a instantaneous plane as shown in Fig 1,

[TABLE]

3 Trajectory Choices and Thruster Calculation Algorithm

3.1 Trajectory Choices

From the two vectors in Eq. (6), we calculate the trajectories of the S/Cs which keep the constant equal-arm configuration as follows: one S/C follows its geodesic orbit as selected in Section 2.1, and other two S/Cs follow this S/C to form an equilateral triangle with desired constant arm length by using the thruster compensation.

Without loss of generality, here we let the S/C2 and S/C3 follow the geodesic motion of S/C1. The trajectories of the three S/Cs are calculated by

[TABLE]

where $a={l}/{(2\sqrt{3})},b=l/2,c=l/\sqrt{3}$, the $l$ is the proposed arm length. $\mathbf{r}_{{\mathrm{S/C}}i}$ is the geodesic orbit of S/C$i$. Therefore, the three trajectories, $\mathbf{r}_{\mathrm{traj,S/C}i}\ (i=1,2,3)$ will form a constant arm equilateral triangle.

For the DECIGO configuration, there is cluster b (30-degree-b/60-degree-b) as shown in Fig. 3. We assume the three spacecraft of cluster b will follow the geodesic S/C1 in cluster a. The trajectories of the three S/Cs are calculated by

[TABLE]

3.2 Thruster Acceleration Algorithm

From Eq. (7) or (8), we can obtain the acceleration at a specific point in a trajectory by calculating the second derivative of position with respect to time,

[TABLE]

On the other hand, for each trajectory, we put the position and its first derivative, velocity, into the ephemeris framework to calculate the acceleration $\mathbf{a}_{\mathrm{eph}}$,

[TABLE]

where $\mathbf{a}_{\mathrm{Newton}}$ and $\mathbf{a}_{\mathrm{1PN}}$ are the Newtonian and first-order post-Newtonian acceleration from the major celestial bodies in the solar system considered as point mass, $\mathbf{a}_{\mathrm{fig}}$ is the acceleration due to the figure effects from the Sun, Earth and Moon, and $\mathbf{a}_{\mathrm{asteroid}}$ is the acceleration from the 340 asteroids’ Newtonian perturbation. The explicit interactions in our CGC ephemeris framework are fully described in reference \refcitewang&ni2012,wangni2015,wang2011.

Then the thruster acceleration to maintain the constant arm length trajectories is calculated by

[TABLE]

4 Thruster Requirement for Different Missions

4.1 B-DECIGO

The B-DECIGO proposed constant 100 km arm length mission orbit around the Earth[25, 26]. In our study about AMIGO mission which proposed the arm length with 10000 km arm length, we considered a multitude of possible orbital configurations[15]. In this work, the configurations AMIGO-E1, AMIGO-EML4 and AMIGO-S are selected to downscale by 100 times as the possible B-DECIGO orbital configurations. Therefore, we choose one geodesic orbit from the AMIGO configuration as fiducial, and find the other two S/Cs trajectories by scaling down the other two AMIGO arms in the constant-equal-arm formation. In this work, three B-DECIGO orbital configurations considered are as follow,

• •

B-DECIGO-AM-E1 around the Earth with semimajor axis of 100,000 km, this orbit configuration is downscaled from the AMIGO-E1 configuration[15].

• •

B-DECIGO-AM-EML4 near the Earth-Moon L4 point, this orbit configuration is shrunk from the AMIGO-EML4 configuration[15]

• •

B-DECIGO-AM-S which is LISA-like around the Sun with 10-degree trailing angle. In this case, the 100 km is downscaled from the AMIGO-S 10,000 km configuration[15].

The corresponding results are shown in Fig. 2. As we can see, the acceleration of the B-DECIGO-AM-E1 configuration is up to 25 $\mu$m/s2 in 180 days, the acceleration of the B-DECIGO-AM-EML1 configuration also could go up to 25 $\mu$m/s2 in 180 days, and the acceleration of the B-DECIGO-AM-S configuration is up to 0.15 nm/s2 in 600 days. The results clearly show that the heliocentric orbit B-DECIGO is much easier (4 orders of magnitude) to adjust as expected.

4.2 DECIGO

The DECIGO proposed constant 1000 km arm length in the heliocentric Earth trailing LISA-like orbit[26]. As the diagram shown in Fig. 3, six S/Cs form two concentric equilateral triangles at the same trailing angle position, and the orientation of two triangle have 180∘ difference in the formed plane. With respect to this position, there are another two triangle constellations with angle separation 120∘ and 240∘ in the ecliptic plane, respectively.

We assume the scientific observation of DECIGO starts from $t_{0}=$ JD2464694.0 (2036-Jan-1st 12:00:00) and calculation the mission orbits. In this work, we presume that two possible placements for the three constellations in the ecliptic plane which are

• •

three constellations at 30 deg (30-degree-a and 30-degree-b), 150 deg and 270 deg trailing angles with respect to position of the Earth, and the thruster acceleration compensations for each S/C at different constellations are shown in Fig. 4.

• •

three constellations at 60 deg (60-degree-a and 60-degree-b), 180 deg and 300 deg trailing angles with respect to position of the Earth, and the thruster acceleration compensations for each S/C at different constellations are shown in Fig. 5.

As we can read from Fig. 4 and 5, the thruster acceleration requirement is less than $\sim$200 pm/s2 in 1000 days. In the left panel of Fig. 4 and 5, the accelerations required for the five spacecraft change with the distance from the selected fiducial S/C1 in cluster a. An alternative way is to let one of the spacecraft of cluster b to follow a geodesic orbit and the other two spacecraft to follow this spacecraft; this may have the advantage that one spacecraft of cluster b will not be accelerated, and cluster a and cluster b will be more independent. Unlike other cases considered in this paper, we have not optimized the geodesic orbits given by the first selected orbit initial conditions; if optimization is implemented, the thruster requirement for DECIGO should be less than the result we estimated (although of the same order).

4.3 New LISA and TAIJI

The LISA proposed to use the drag-free technology to leave S/Cs orbiting with the gravitational field and the arm length change with time. The new LISA plans to use $2.5\times 10^{6}$ km nominal arm length[7]. The TDI is required to achieve the equivalent near-equal interference paths. TAIJI mission proposed a LISA-like orbit with nominal arm length $3\times 10^{6}$ km[8].

In our previous work[10], we worked out a set of the LISA-like mission orbits with the observation starting time at $t_{0}=$ JD2461853.0 (2028-Mar-22nd 12:00:00). Afterwards, we worked out another TAIJI mission orbit which is ahead of the Earth by 20∘ and form a $10^{8}$ km baseline with LISA. The large separation has high angular-resolution virtues for the short duration sources detection during joint LISA-TAIJI observation scenarios. By using the orbits we obtained, we explore the thruster requirement for the LISA and TAIJI missions to maintain a constant arm configuration. From the method we described in Section 3, the thruster acceleration compensations for the LISA and TAIJI mission are shown in Fig. 6. The acceleration requirement could be up to 2.5 $\mu$m/s2 for both LISA and TAIJI.

5 Fuel Requirement at Different Specific Impulses

In this section, we assume that the mass $M$ of a typical spacecraft weighs 1000 kg and derive propellant requirement at two different specific impulses 300 sec and 1000 sec. Typical liquid fuels have about 300 sec specific impulse. The colloid thruster fuel has the specific impulse in the range of 200-2000 sec. Let $I_{\mathrm{sp}}$ be the specific impulse in unit of sec. The thruster force $F_{\mathrm{thruster}}$ of a propeller is

[TABLE]

where $g_{0}$ (= 9.8 m/s2) is the gravity acceleration at surface of the Earth, and $\dot{m}$ (in unit of kg/s) is the rate of consumption of fuel. The consumed propellant $\Delta m$ is

[TABLE]

For small $\Delta m/M$ ratio, it can be approximated as

[TABLE]

where $T$ is the considered mission time. The yearly propellant in column 4 and 5 of Table 5 is averaged by using

[TABLE]

From this formula, we calculate the propellant required for 1000 kg spacecraft in Table 1 at specific impulses 300 sec and 1000 sec of thruster assumed. We compile our results on the required acceleration, thruster requirement, and fuel requirement at different specific impulse 300 sec and 1000 sec for AMIGO, AIGSO, B-DECIGO, DECIGO, LISA, and TAIJI in Table 5.

As we can read from Eq. (14), the propellant requirement is proportional to the $M$ of the S/C and inversely proportional to the specific impulse $I_{\mathrm{sp}}$ if $\Delta m/M$ is small. In these cases, the propellant requirements could be readily inferred from Table 5 for other values S/C mass and/or specific impulse. For the large $\Delta m/M$, the Eq. (13) should be applied.

6 Discussions and Conclusions

For the geocentric orbit options of B-DECIGO which we have explored, the propellant mass requirement is a concern. For the heliocentric orbit options of B-DECIGO and DECIGO, the propellant requirement could be easily satisfied.

From propellant requirement estimated in Table 5, this scheme is technologically feasible for both LISA and TAIJI as far as fuel is concerned. This means that, for LISA and TAIJI, both ordinary Michelson interferometry and first-generation Michelson interferometry would be good configurations for GW detection if other accompanying requirements can be met. The two accompanying requirements to fuel requirement of ordinary constant equal-arm Michelson interferometry implementation of LISA/TAIJI are (i) thruster noise requirement, (ii) and proof mass actuation requirement.

With large thruster power, the thruster noise increases. To minimize thruster noise, 2-stage thrusters/3-stage thrusters could be considered. With the pseudo-random code (PRC) ranging, sub-meter range accuracy can be achieved. A set of mN thrusters could provide the basic thruster force. The $\mu$N thruster could just provide the small residual precision adjustment to the acceleration needed to confirm to the PRC ranging measurement. The limiting noise would come from $\mu$N thruster noise and pseudo-random code range measurement noises.

The other requirement concerns proof mass actuation noise. In the application of $\mu$m/s2 acceleration to proof mass, the noise should not be greater than pm/s2 in some kind of average. A reference is needed for this measurement/monitoring of the actuation. Therefore an alternate proof mass is needed. Laser metrology has the required accuracy[47]. The two proof masses can alternate to become the reference masses. This way the required dynamical range can be achieved. However, the gap size limits the total range of one acceleration maneuver to about 2 mm. This limits the one acceleration maneuver time of AMIGO-S to about 500 sec. For AIGSO, B-DECIGO-AM-S, and DECIGO, the time limits of one acceleration maneuver are about 5000 sec, 1500 sec and 1500 sec, respectively. For LISA and TAIJI, this time limit is about 15 sec. The actuation induced Fourier spectral components needs to be subtracted. Whether the time limits for various mission implementation are enough need to be studied experimentally[48].

Constant equal-arm Michelson interferometry is preferred if these technical issue can be resolved, because it does not have the complication of the TDI. Moreover, constant equal-arm Michelson interferometry and the TDI could both be tested at the beginning of science mission and worked out in the same mission if pre-mission preparation was done.

Acknowledgments

We would like to thank Gerhard Heinzel (AEI), Hongbo Jin (NAOC), Jungang Lei (LSTT), Ziren Luo (IM, CAS), Liming Song(IHEP), Stefano Vitale (Trento U.), Bill Weber (Trento U.), Peng Xu (Lanzhou U.), and Zhen Yang (NSSC) for helpful discussions. This work was supported by National Key Research and Development Program of China under Grant Nos. 2016YFA0302002 and 2017YFC0601602, and Strategic Priority Research Program of the Chinese Academy of Sciences under grant No. XDB21010100.