## 2 Mathematical Model  
The physical model adopted for the forces considered in the dynamics used here takes into account the Earth's flattening at the poles, the perturbations from the Sun and Moon, and the SRP. The equation of motion of the debris can be written as  
\[\ddot{\mathbf{r}}=\ddot{\mathbf{r}}_{G}+\ddot{\mathbf{r}}_{SRP}+\ddot{\mathbf{ r}}_{TBP_{Sun}}+\ddot{\mathbf{r}}_{TBP_{Mon}} \tag{1}\]  
where \(\ddot{\mathbf{r}}_{G}\) is the force due to the Earth's gravitational field, the terms \(\ddot{\mathbf{r}}_{TBP_{Sun}}\) and \(\ddot{\mathbf{r}}_{TBP_{Mon}}\) are the resulting forces due to the third-body gravitational attractions of the Sun and Moon, respectively, and the term \(\ddot{\mathbf{r}}_{SRP}\) represents the acceleration generated by SRP. The term \(\ddot{\mathbf{r}}_{G}\) is given by \(\ddot{\mathbf{r}}_{D}+\ddot{\mathbf{r}}_{I_{2}}\), where \(\ddot{\mathbf{r}}_{D}\) is the direct acceleration due to the gravitational central force of the Earth and \(\ddot{\mathbf{r}}_{I_{2}}\) is the disturbing acceleration given by the Earth's flattening at the poles. After several algebraic manipulations, Equation (1) can be defined as a function of the orbital elements of the space debris: \(a\)--semi-major axis, \(e\)--eccentricity, \(i\)--inclination, \(g\)--argument of perigee, \(h\)--longitude of the ascending node, and \(l\)--mean anomaly. In this work, we rewrite the equation of motion given by Equation (1) as a gradient of a potential and we present the disturbing potential as a function of the orbital elements. The algebraic development of the equations considered here is taken from references [27; 37]. After convenient algebraic manipulations, the equation to represent the oblateness of the Earth (\(\ddot{\mathbf{r}}_{I_{2}}\)), using the single-averaged method, can be written as [27; 37]  
\[<R_{J2}>=-\frac{1}{4}\frac{I_{2}R_{G}^{2}\,n^{2}(-2+3\sin^{2}(l))}{(1-e^{2})^{ 3/2}} \tag{2}\]where \(R_{e}\) is the equatorial radius of the Earth, \(n\) is the angular velocity of the satellite (debris) (mean motion, i.e., \(n=\sqrt{\mu/a^{3}}\)), and \(\mu\) is the Earth's gravitational parameter. The term due to the irregular shape of the Earth (\(J_{2}\)) considered in the present work represents the second term of the gravity field of the Earth, which represents the flattening at the poles. The disturbing potential due to SRP (\(\bar{\textbf{r}}_{SRP}\)) using a single-averaged model and assuming the Sun in elliptical and inclined orbits is written as [27]  
\[<R_{SRP}>=\frac{81}{64}\beta\frac{e{m_{\odot}}^{2}a_{\odot}}{(1-e _{\odot})^{2}}[-\frac{1}{27}(20\cos(i)-20)(\cos(i_{\odot})-1)({e_{\odot}}^{2}- 2/  
The last disturbing potential used in this present work is  
\[R=<R_{I_{2}}>+<R_{SRP}>+<R_{TBP_{\text{Sun}}}>+<R_{TBP_{\text{Mon}}}> \tag{5}\]  
Considering the perturbations previously commented, we replace Equation (5) into [PERSON]'s planetary equations and numerically integrate the set of nonlinear differential equations to obtain the variations of the orbital elements of the space debris. After the numerical integrations, we show the results in the next section. The idea is to move the satellite (debris) into an orbit with an altitude higher than the GEO region, where artificial satellites for communications, military use, meteorology, and other applications are located. In addition, this region presents deactivated satellites, parts of satellites after their fragmentation, etc. To avoid the risk of these non-operational objects colliding with operational satellites, it is necessary to clean the GEO region. Thus, we propose to elevate the orbits of the space debris to a region above this protected GEO \"ring\", which is commonly called a graveyard orbit, considering a solar sail such that the orbital perturbations included will help to amplify the growth of the eccentricity of the debris. As mentioned before and discussed in [27], the first-order SRP is the predominant perturbation due to the high \(A/m\) ratio of the satellite (debris) after it is coupled to a solar sail.