## 1 Introduction  
The analysis of orbit uncertainty propagation is a key step in space mission design [1, 2]. Because of the existence of uncertainty, e.g., injection errors and measurement errors, the spacecraft will deviate from their nominal orbits [3, 4], which may change the orbital stability and impact the performance of exploration mission [5, 6]. Therefore, it is necessary to investigate the evolution of orbital uncertainty under multi-source errors [7], which provides feedback coupling inputs for control and navigation subsystem design [8, 9]. Generally, uncertainty propagation is used to characterize the spacecraft's state moment (usually represented by mean and covariance matrix) or its probability density function (PDF) [10, 11]. Theoretically, the uncertainty description of the orbital state can only be accurately obtained by solving Fokker-Plank Equation (FPE) [12, 13] or carrying out Monte Carlo (MC) simulations [9]. Though the FPE provides accurate solutions for the propagation of the PDF, it is usually difficult to solve the equation as it is a high-dimensionally defined and nonlinear perturbed Keplerian-driven partial differential equation. The MC simulation is always considered a reliable method whose result approaches the true probability distribution when the number of MC samples tends to be infinity [14]. As a result, either FPE or MC simulations suffer from drawbacks related to a heavy computation burden [9]. To improve efficiency, some approximations are made, which can be categorized into gradient-enhanced methods and gradient-free methods [15, 16]. Gradient-enhanced methods, such as state transfer matrix (STM) [17; 18] and state transition tensors (STT) [19; 20; 21], analyze the orbital uncertainty by simplifying the relative nonlinear motion with the Taylor series terms [22; 23]. For gradient-enhanced methods, the gradients or higher-order derivatives are incorporated for approximation [20], and analytical or semi-analytical solutions are obtained [24; 25]. However, gradient-enhanced methods are only suitable for cases whose gradient information is easy to obtain. For some of the practical problems, the dynamic motion is complex, discontinuous, and nondifferentiable, and the derivatives are unavailable [16]. The gradient-free methods (also called non-intrusive methods [26]) allow for the usage of already existing orbit propagation tools as the black-box and avoid any simplifications on the fully perturbed orbital dynamics [27; 28], making them conveniently extend into orbit uncertainty propagation problems [29; 30]. Recently, some sample-based methods, for example, unscented transformation (UT) [31; 32], Conjugate Unscented Transform (CUT) [33], and Cubature rules (CRs) [34; 35], have been well investigated. These methods approximate the probability distribution at a future time by nonlinearly integrating a few sigma samples, which are deterministically chosen according to the given initial distribution [9]. The drawbacks of the above sample-based methods are that they can only deliver a second-order approximation of the distribution, which is inadequate for long-term propagation whose ultimate uncertainty distribution inevitably turns out to be highly non-Gaussian [36; 37]. In order to fully characterize the non-Gaussian orbital uncertainty, the Gaussian mixture model (GMM) is commonly employed [38; 39]. The basic idea of GMM is that an arbitrary PDF can be approximated, as closely as desired, by a finite sum of weighted Gaussian components (also called Gaussian PDFs). For each Gaussian component of GMM, only the mean and covariance need to be propagated because a Gaussian PDF can be completely determined with its first two moments [40; 41]. Therefore, the above methods, such as QRs, UT, CUT, and CRs, can be used to propagate each Gaussian component. The GMM can better represent non-Gaussian PDFs and results in a more accurate prediction [42]. However, the GMM suffers from a heavy computation burden.