#### Abstract

In this paper, a modified targeting strategy is developed for missions on libration point orbits (LPOs) in the real Earth-Moon system. In order to simulate a station-keeping procedure in a dynamic model as realistic as possible, LPOs generated in the circular restricted three-body problem (CRTBP) are discretized and reconverged in a geocentric inertial system for later simulations. After that, based on the dynamic property of the state transition matrix, a modified strategy of selecting target points for station-keeping is presented to reduce maneuver costs. By considering both the solar gravity and radiation pressure in a nominal LPO design, station-keeping simulations about fuel consumption for real LPOs around both collinear and triangular libration points are performed in a high-fidelity ephemeris model. Results show the effectivity of the modified strategy with total maneuver costs reduced by greater than 10% for maintaining triangular LPOs.

#### 1. Introduction

Since the 1970s, libration points and bounded motions around them have been widely studied due to the unique dynamic properties near these equilibrium points in CRTBP. Until now, many missions have proven to be successful in placing a spacecraft in libration point orbits. On account of an unstable dynamic environment around libration points, effective strategies of station-keeping are also presented and examined to maintain a spacecraft within the reasonable vicinity of the desired trajectory or region. Moreover, the concept of using Earth-Moon libration point satellites for cislunar navigation has also been studied for decades.

Farquhar [1] presented a detailed analysis of the translation-control problem for satellites in the vicinity of collinear libration points and also proposed several possible missions around libration points. For example, the idea of using a lunar relay satellite operating around Earth-Moon to provide navigation capability for the far side of the Moon is described. Howell and Pernicka [2] and Pernicka [3] developed a patching technique incorporating with a differential correction scheme to construct Lissajous trajectories and “near-halo” orbits numerically of arbitrary size and length. A targeting method is presented by Howell and Pernicka [4] that uses maneuver executed at a discrete time interval for a spacecraft moving on halo or Lissajous orbits in the vicinity of of the Sun-Earth/Moon barycenter system. Gómez et al. [5] performed a target point approach and Floquet mode approach to estimate the station-keeping costs of a quasihalo orbit around in the Earth-Moon system. Folta and Vaughn [6] implemented numerical computations using high-fidelity models to simulate the various maneuver delta-V and temporal costs associated with orbits about the Earth-Moon libration points and also studied the transfer between libration orbit locations. In 2010, several station-keeping strategies are considered for application to ATREMIS by Folta et al. [7] and two approaches are examined to investigate the station-keeping problem in the Earth-Moon system in a full ephemeris model. Pavlak [8] introduced a numerical process of constructing nominal orbits around the Earth-Moon libration points in a high-fidelity model based on LPOs generated in the CRTBP model, and later Pavlak and Howell [9] proposed a near-optimal station-keeping method for long-term station-keeping in the Earth-Moon system. An adaptive robust controller for a spacecraft with unknown mass and boundary of perturbation forces is proposed by M. Xu and S. Xu [10] to deal with the station-keeping control for halo orbits around Sun-Earth . Folta et al. [11] summarized the results of ARTEMIS mission operations and used the optimal continuation strategy to obtain operational station-keeping results, which were also compared to orbit stability information. Qian et al. [12] investigated a numerical method that constructs perturbed halo/Lissajous orbits in a full solar system model, which is based on an improved multiple-shooting strategy that utilizes the initial conditions with consideration of the instantaneous state of libration points for each integral arc. Qian et al. [13] also proposed a new continuous low-thrust-based station-keeping strategy for quasiperiodic orbits around the Earth-Moon point in a high-fidelity model. The strategy was proven effective by simulation results with low-energy consumption. Zhang and Li [14] presented two complementary methods by exploiting the hyperbolic dynamics of the collinear libration points. Both methods need no information of the nominal orbits but utilize the manifolds to eliminate unstable components.

In Hill et al. [15], an autonomous orbit determination technique, known as “LiAISON Navigation,” was proposed for libration point satellites. Considering the strong asymmetry of the restricted three-body force field, the feasibility was demonstrated that the libration point navigation satellites can be determined autonomously in the case of only using scalar satellite-to-satellite tracking data. Based on the LiAISON technique, a study on the libration point satellite navigation system in the Earth-Moon system was presented by Zhang and Xu [16] to obtain candidate constellation architectures according to coverage ability and autonomous orbit determination performance. The current libration point satellite navigation system is mainly designed under the CRTBP, and further optimization of satellite constellations needs both studies of LPOs in a real Earth-Moon force model and station-keeping strategies for libration point satellites. Constructing the LPOs of variable amplitude and different orbital types in a high-fidelity ephemeris model is required for both navigation constellation design and station-keeping research.

In this work, a discussion about nominal LPO design and a detailed process of constructing LPOs in an ephemeris model incorporating JPL DE405 ephemerides is proposed. In addition to the gravity of the Earth and the Moon, both solar gravity and radiation pressure are considered in the process of nominal LPO design in order to improve the accuracy of orbits and reduce control consumption in later simulations. Then, the targeting strategy is modified to fit nominal LPOs with divergence of orbital motions in the real Earth-Moon system. LPOs around both Earth-Moon collinear and triangular libration points are reconverged as nominal orbits for station-keeping simulations with a modified targeting strategy. The results of simulations show that the modified strategy performs well to select target points for real LPOs around triangular points.

The remainder of this paper is organized as follows. In Section 2, an introduction to the dynamic background of CRTBP and libration point orbits is provided, with a semianalytical computational method to construct the LPOs and a detailed CRTBP-to-ephemeris transition procedure for LPOs in the Earth-Moon system. After that, a typical description of the targeting strategy for station-keeping and a new modified method for optimizing the target point selection are given in Section 3, and Section 4 provides a discussion of nominal LPO design. Section 5 shows the simulation results of maneuver costs for real LPOs around both collinear and triangular points to examine the modified strategy. Some conclusions are drawn in Section 6.

#### 2. Dynamic Background

##### 2.1. CRTBP

The dynamic model used in this work is the circular restricted three-body problem (CRTBP). The equations of motions expressed in the barycentric synodic frame are written as [17] where with as the mass parameter of the system. For the Earth-Moon system, we have . and are the distances from the spacecraft to the Earth and the Moon, respectively. For convenience, the following system of units is introduced to normalize the variables used in practical computations: where is the mass of the Earth, is the mass of the Moon, and is the mean distance.

The equations of motion of CRTBP allow five well-known equilibrium points in the synodic coordinate, also known as the libration points located in the plane of motion of the primaries. Among these five points, three collinear points are located on the -axis and two triangular ones forming two equilateral triangles with the primaries (as shown in Figure 1).

##### 2.2. The LPOs

The abundant dynamical properties of five libration points offer different types of periodic or quasiperiodic orbits, so libration point orbits (LPOs) are represented to be used in many different missions. In this work, the LPOs generated from CRTBP in the Earth-Moon system are applied as initial orbits for later analysis in a more real dynamic model. Around the three collinear libration points , there are three types of LPOs: Lissajous orbits, Lyapunov orbits, and halo orbits. Specifically, Lyapunov orbits are specific forms of Lissajous orbits and include vertical Lyapunov orbits and planar Lyapunov orbits, which are caused by the in-plane or out-of-plane amplitude equaling to zero, respectively. For the two triangular libration points , there are three kinds of periodic orbits: the planar long and short periodic orbits and the vertical periodic orbits.

There are various methods in constructing LPOs. Both numerical differential correction and analytic construction are proven to be reliable. As a large number of LPOs of different type with different amplitudes need to be generated for station-keeping simulations, a computational method of high efficiency and precision is indispensable. Then, a semianalytical computational method is adopted to construct the LPOs. Using the Lindstedt-Poincare method, the analytical solutions of different LPOs are truncated at a relatively high order and can be obtained by series expansion. Through this semianalytical method, different types of LPOs with amplitudes of specific values could be generated easier than those numerical methods that need a series of differential corrections. Below are the resulting expressions without derivation. Construction details and definitions of variables can be found in [18–21].

For LPOs around collinear libration points [20],

For LPOs around triangular libration points [21],

##### 2.3. The Ephemeris Model and Nominal Orbits

The LPOs generated under the CRTBP have proven to be useful and effective as reference orbits for deep space mission design, satellite constellation analysis, and station-keeping research since LPOs are generally periodic and available with specific amplitude and different orbital types. In the real Earth-Moon system, lunar eccentricity and perturbing effects due to additional gravitational bodies or solar radiation pressure should be taken into account for accurate mission design. When the dynamic model is transformed to a more realistic one than CRTBP, the LPOs no longer keep complete periodic motion due to the influence of lunar eccentricity and other perturbations. For collinear libration points, perturbations can cause the trajectories to diverge significantly due to the nature of the saddle-type libration point. For the triangular case, the special dynamic environment of restricted three-body problem still maintains the orbits operating within a close torus around the original orbits over a certain period of time. The libration point orbits constructed under the ephemeris model after a shooting method could be continuous in limited areas for years, so it is feasible to use such orbits as reference orbits for short-term missions in the real Earth-Moon system. These orbits are more realistic for precise satellite constellation design in the libration point navigation system and could also serve as nominal orbits in station-keeping simulations to effectively reduce the total maneuver costs.

Before referring to the targeting strategy, the definition and construction process of nominal orbits in this work ought to be elaborated in detail. During a station-keeping simulation of LPOs in the targeting strategy, it is desired to maintain the actual spacecraft trajectory within some limited regions around the reference orbit, which is predesigned and called “nominal orbits.” The aim of performing maneuvers at specified times is to drive the spacecraft trajectory back to the acceptable region around the nominal orbits. Maneuver costs computed in the targeting strategy may be significantly influenced by the accuracy of the nominal orbits. Therefore, selection of the nominal orbits can be critically important to avoid unnecessary consumption. However, once the nominal orbits are constructed, it is required to keep the balance between the fidelity of the dynamic model and the computational cost. Considering much more perturbations, especially the dissipative ones, will result in an overlong orbit designing process and even failures in the continuity of trajectory.

To construct nominal orbits in a higher-fidelity model, a general method for transitioning LPOs from the CRTBP to an Earth-Moon-Sun ephemeris model is presented by Pavlak [8]. Before introducing the transition procedure, the equation of motion for the ephemeris model, also known as -body relative equations of motion, is given as where subscripts , , and are for the body of interest, central body, and perturbing bodies, respectively. The vector represents the position of each perturbing body with respect to the central body, and represents the gravitational constant in dimensional units. The vector is obtained from vector subtraction, i.e., , and represents the position of each perturbing body relative to the body of interest, .

By incorporating JPL DE405 ephemerides and multiple-shooting algorithms, a converged trajectory in the CRTBP could be numerically transited to an inertial ephemeris model, considering the gravity of the Earth, the Moon, and the Sun. A compendious CRTBP-to-ephemeris transition procedure is summarized as follows: (1)Compute the desired orbit in the CRTBP(2)Discretize the orbit into patch points(3)Shift the CRTBP rotating states to primary-centered rotating states(4)Dimensionalize the states using instantaneously defined characteristic quantities(5)Apply the transformation to the ephemeris model(6)Nondimensionalize the inertial J2000 states using the standard characteristic quantities(7)Reconverge the solution in the ephemeris model

All the related frame transforming algorithm and multiple-shooting techniques can be found in the reference, and some modification was added corresponding to the subsequent station-keeping research in this work. Since the transition is only applied to LPOs around collinear points in the Earth-Moon system in the previous study, the feasibility of such transition procedure still remains to be assessed for those orbits around triangular points.

Based on this CRTBP-to-ephemeris transition for LPOs, some typical orbits around both collinear and triangular points are reconverged and constructed in a high-fidelity dynamic model that considers the gravity of the Earth, the Moon, and the Sun. As shown in Figure 2, the coordinate system of reconverged orbits is transformed back to the synodic one in order to compare these newly constructed orbits with original patch points generated in the CRTBP. To convert LPO trajectories in the ephemeris model back into the synodic one is almost the reverse process of the CRTBP-to-ephemeris transition and would not be explained in detail.

**(a)**

**(b)**

**(c)**

**(d)**

For LPOs around collinear points, the reconverged orbits are very close to the patch points. For LPOs around triangular points, the orbits could also be converged in the ephemeris model but the trajectories are no longer within limited regions around the original patch points. Both planar long periodic (LP) orbits and planar short periodic (SP) orbits enlarge their amplitudes to about over 20,000 km, and some special orbital structures like torus also appear. This difference could be explained by the complicated dynamic environment of triangular points, such as the great influence of solar gravity, but the main focus of the later work is to modify the traditional targeting strategy so that it could fit the real LPOs around triangular points with specific orbital structures.

#### 3. Targeting Strategy and Modification

##### 3.1. Targeting Strategy

A targeting strategy, also known as the target point approach, is a common method for station-keeping missions. The approach was first presented to be used for LPOs in the Earth-Moon system in 1993 by Howell and Pernicka [4].

Unmodeled perturbations and orbit injection error will result in a drift of the spacecraft from the nominal path, and the unstable nature of libration point orbits will soon amplify the deviation. To maintain the actual spacecraft trajectory within some torus centered about the predesigned nominal orbits, timely adjustment such as impulsive maneuvers computed by the targeting strategy is required. The targeting strategy is a control method based on linear feedback and could efficiently compute the maneuvers to be performed at specified times to drive the spacecraft trajectory back to the acceptable region near the predesigned nominal orbits.

The targeting strategy computes maneuvers by minimizing a weighted cost function that is defined in terms of delta-v with the position and velocity deviations from a nominal orbit at specified times. Therefore, the state transition matrix of the nominal orbits must be available at any time along the orbit. Inaccurate representations may contribute to higher maneuver costs. For convenience, the state transition matrix of the nominal orbit is partitioned as where , , , and are submatrices. Let be the vector representing the deviation of the actual trajectory from the nominal orbit at time due to a velocity perturbation at time , a position perturbation at time , and a corrective maneuver performed at time . Then, is the future time at which the predicted actual path will be compared to the nominal one. The linear approximate for the position deviation at target may be expressed as

Then, the cost function is then defined by Howell and Pernicka [4]: where superscript denotes the transpose. The weight matrix is the symmetric positive definite, while and are the positive semidefinite. To minimize the cost function shown in equation (9), the result produces the optimal maneuver as [4]

For equation (10), maneuvers were assumed to occur at and based on the use of two target times which are chosen arbitrarily. The weighting matrices , , and are assumed as constant during the later station-keeping simulation.

##### 3.2. Modified Method

The traditional targeting strategy for LPOs is elaborated in the above section. Based on linear feedback, maneuvers at specific times are computed efficiently under the target strategy, so the spacecraft trajectory is maintained within the limited region around the nominal orbit. For LPOs generated under the CRTBP model, the periodicity and orbital shapes make it simple to select target points for the targeting strategy, such as the state points after the half period or one period in the nominal orbit or the state points with specific position features. These methods of determining target points perform well for periodic orbits but still could be improved for those reconverged orbits constructed in the ephemeris model. With a view to the dynamic properties of the state transition matrix, a modified strategy is presented to choose the proper target points in a relatively reasonable and concise way for LPOs constructed in the ephemeris model, which even own chaotic orbital structures.

Notice the upper-right section of a state transition matrix as follows:

This submatrix contains the linear effect of little variation of the initial velocity to the final position in a specific trajectory. Once the actual position of the spacecraft exceeds the preset positional limit, the position deviation is assumed to be known at that point and the preconstructed nominal orbits contain dynamic information in the following time. By solving the linear equations as follows,

We can obtain a pseudo maneuver vector corresponding to initial time and end time ; then, we define a ratio between the modulus of this pseudo maneuver vector and the time interval as

The ratio partly reveals the size of maneuver cost of each state point as a target point in the following time. To select a relatively capable state point for the computation of the targeting strategy may result in lower maneuver consumption.

In later station-keeping simulations for LPOs, both the traditional targeting strategy and modified method are performed to get average maneuver costs as a comparison. Whenever the position deviation between the spacecraft and reference orbit is out of limit, a maneuver needs to be carried out to maintain the spacecraft operating in a reasonable region. According to the modified method, we may choose the state point with minimum ratio during the first following period as the first targeting point, while the period here refers to the original period of the corresponding orbit in CRTBP. The second targeting point is also determined this way in the second period, and a restriction is also added that the time interval between the two points must be longer than a quarter of the corresponding period. This modified strategy tries to pick out relatively reliable state points as target points in calculating maneuvers so that the total control consumption for a certain period of time may be reduced. Detailed simulation results for LPOs are given in the following sections.

#### 4. Nominal LPO Design

##### 4.1. Dynamic Models

The construction of nominal orbits in the ephemeris model has been thoroughly introduced in Section 2. By incorporating JPL DE405 ephemerides, the nominal orbits generated in this method will be finally converted into an ephemeris model, the geocentric celestial reference frame (GCRF). The GCRF could be regarded as a local inertial dynamic frame for the libration point orbits in the Earth-Moon system.

In this work, the nominal LPOs are constructed in a dynamic model not only including the gravity of the Earth, the Moon, and the Sun but also containing the solar radiation pressure (SRP) perturbation. Because the force direction of SRP is almost the opposite of solar gravity under the general simplified model, the construction of nominal orbits will be hardly influenced by taking the SRP perturbation into account and the maneuver costs can be reduced due to using nominal orbits generated in a higher-fidelity model. Once the nominal orbits are determined, the other perturbations can also be enumerated for station-keeping demand. In Table 1, we can contrast the dynamic models between the nominal orbit construction and the actual station-keeping simulation environment with the perturbations included. It should also be noted that the area-to-mass ratio of SRP perturbation is set to 10^{9}, as a typical value of a spacecraft.

##### 4.2. Discussion of Nominal LPOs

As mentioned in the dynamic background, there are commonly three kinds of periodic orbit around the collinear libration points: vertical Lyapunov (VT) orbits, planar Lyapunov (PL) orbits, and halo orbit. For triangular libration points, planar long periodic (LP) orbits, planar short periodic (SP) orbits, and vertical periodic (VP) orbits are three kinds of general periodic orbits. Although the nominal LPOs may no longer keep complete periodicity after the CRTBP-to-ephemeris transition, we still label each nominal LPO based on the type and amplitude parameter of its original periodic orbit in CRTBP.

Generally, the nominal LPOs reconverged in the ephemeris model are similar to their original orbits in CRTBP for all types of LPOs, but we find that the amplitude enlargement exists in triangular LPOs with relatively small amplitudes after the CRTBP-to-ephemeris transition. A planar short periodic orbit of original amplitude of 4,000 km would be enlarged into an orbit with an actual trajectory range of over 20,000 km around the point. As shown in Figure 3, the reconverged nominal LPOs are operating far beyond the range of their corresponding patch points discretized from initial LPOs in the CRTBP model. The main cause of such enlargement may be the considerable solar gravity, which is not included in CRTBP and leads to divergence of orbital motions around the triangular libration point.

**(a)**

**(b)**

After constructing nominal triangular LPOs of various amplitudes, the effects of amplitude enlargement are found to be weakened with the increasing amplitudes of LP and SP orbits and almost disappear when the amplitudes are increased to a certain size. For both the SP and LP orbits around the point, the enlargement almost no longer exists when the amplitude parameter reaches about 20,000 km. The relation between the final amplitude and the original amplitude is summarized in Table 2. On this condition of the nominal orbits around the triangular libration point , we may begin the station-keeping simulations for those orbits whose amplitudes cover a wide range from 4,000 km to more than 20,000 km.

Also mentioned are the vertical periodic (VP) orbits, which are also influenced by solar gravity and solar radiation pressure. The construction of small-amplitude VP orbits shows that VP orbits around the triangular point generate planar motion components in the - plane after the reconverging process. The trajectories of real VP orbits are just similar to those of real SP orbits around , so VP orbits would not be studied in later simulations.

#### 5. Station-Keeping Simulations

Before we start the station-keeping simulation, a parameter setting must be provided in advance to ensure that the simulation process is entirely clear and specific. Based on the LiAISON Navigation technique, a systematic study on the libration point satellite navigation system was presented by Zhang and Xu [16], where the Earth-Moon four-satellite constellation architectures were obtained from a joint consideration of autonomous orbit determination performance. In a 180-day period, the maximum position determination error can be within a few meters, which is considered a conservative reference of the tracking error in this work.

##### 5.1. Parameter Setting

For both LPOs around the collinear libration points and the triangular libration points, it is common to set some error parameters before the station-keeping simulation, including orbit injection error, tracking error, and maneuver error. By incorporating with other necessary baseline values, the parameter setting of the targeting station-keeping strategy in this work is shown in Table 3.

##### 5.2. Results of Collinear LPOs

Using the targeting strategy with parameters set before, Monte Carlo simulations based on 100 trials are undertaken for nominal LPOs of different types and amplitudes around collinear libration points in the real Earth-Moon system. The simulation results, average maneuver costs for 180 days, are shown in Tables 4 and 5. As a comparison, the mean costs computed by the modified strategy and the mean costs computed by the general targeting strategy that selects every two target points in fixed time intervals are both given to examine the efficiency of the modified method.

In general, the results show that the modified strategy slightly reduces the maneuver costs for collinear LPOs, even though the reductions in consumption are generally less than 2%. For further discussion, some detailed maneuvering results of halo orbits around are given in Table 6, including the average time interval between each maneuver and average single consumption for both original and modified strategies. The detailed results show that the average consumption of both methods for each specific orbit is almost the same, but the modified strategy slightly lengthens the time interval between each maneuver. Thus, we find that the modified strategy reduces the total consumption by making each maneuver more effective though the effect for LPOs around collinear points is much limited.

The results are not satisfactory, but we still keep the data for further discussion and later comparison with the simulation results of triangular LPOs. The collinear LPOs keep very similar orbital shapes to original orbits after the CRTBP-to-ephemeris transition, so the modified method shows little advantage over the traditional strategy. Besides, the costs per 180 days are becoming slightly larger with the increasing amplitudes for each type of the nominal orbits and the costs of real LPOs around are slightly less than those around .

##### 5.3. Results of Triangular LPOs

It is general to repeat the same station-keeping simulation procedure as before for triangular LPOs. It should be noted first that the nominal orbits constructed in the same CRTBP-to-ephemeris transition from the typical LP orbits and SP orbits show some different orbital characteristics from those around collinear points. The modified method is presented to select target points for these nominal orbits, and its performance will be examined in the following station-keeping simulations for triangular LPOs in the real Earth-Moon system.

The results of station-keeping simulations by both original and modified strategies are shown in Table 7. The percentage reduction of consumption costs between two strategies is also provided. In addition, some detailed maneuvering results are also given in Table 8 and only several specific orbits are taken as examples. The average single consumption is still almost the same for two methods, and the average time interval between each maneuver is lengthened by the modified strategy, which causes reductions in total consumption.

From the above results and discussions, it can be found that the performance of the modified strategy is relatively better and the percentage reductions of total maneuver costs by using the modified strategy are all greater than ten percent; thus, the modified method of selecting target points proves to be effective for these reconverged orbits around Earth-Moon triangular points. Besides, the average maneuver costs of real triangular LPOs are relatively larger than the costs of collinear LPOs for the same trajectory duration of 180 days.

#### 6. Conclusions

This paper presents a detailed study on the station-keeping simulations of libration point orbits in the real Earth-Moon system by performing a targeting strategy. In order to construct nominal orbits in a high-fidelity dynamic model, a CRTBP-to-ephemeris transition procedure is provided for different types of LPOs generated in CRTBP. After that, the targeting strategy is modified for LPOs in the real Earth-Moon system because of the diverging orbital structure of various LPOs after the CRTBP-to-ephemeris transition. The modified method provides a way to select proper target points on the basis of the dynamic properties of the state transition matrix. Based on high-precision nominal orbits and the modified targeting strategy, station-keeping simulations are carried out for real LPOs of three different types around and and two types around . The results of collinear LPOs show that the modified strategy is not so satisfactory for LPOs around and and the maneuver costs are all less than 15 m/s in the trajectory duration of 180 days. In station-keeping simulations of triangular LPOs, the modified method for optimized selection of target points proves to be effective, with average maneuver costs reduced by more than 10%. The simulation results show that the maneuver costs are not greater than 25 m/s in the trajectory duration of 180 days for the orbits around . In general, all the results show that the maneuver costs are acceptable for a deep space mission and libration point navigation satellite system.

This paper applies a targeting strategy with a modified method of selecting target points to both collinear and triangular LPOs. The construction process of nominal LPOs and the results of maneuver costs could serve as a reference for future study of libration point navigation system design and targeting strategy application.

#### Data Availability

No data were used to support this study.

#### Conflicts of Interest

All the authors declare that there is no conflict of interest regarding the publication of this paper.

#### Acknowledgments

This work was carried out with financial support from the National Basic Research Program 973 of China (No. 2015CB857100), the Base of National Defence Scientific Research Fund (No. 2016110C019), the National Natural Science Foundation of China (Nos. 11603011 and 61803027), the Shanghai Space Science and Technology Innovation Fund (No. SAST2017088), and the National Natural Science Foundation (No. 11705088).