Abstract

The complexity of the GNSS and the several types of satellites in the MEO region turns the creation of a definitive strategy to dispose the satellites of this system into a hard task. Each constellation of the system adopts its own disposal strategy; for example, in the American GPS, the disposal strategy consists in changing the altitude of the nonoperational satellites to 500 km above or below their nominal orbits. In this work, we propose simple but efficient techniques to discard satellites of the GNSS by exploiting Hohmann type maneuvers combined with the use of the resonance to increase its orbital eccentricity, thus promoting atmospheric reentry. The results are shown in terms of the increment of velocity required to transfer the satellites to the new orbits. Some comparisons with direct disposal maneuvers (Hohmann type) are also presented. We use the exact equations of motion, considering the perturbations of the Sun, the Moon, and the solar radiation pressure. The geopotential model was considered up to order and degree eight. We showed the quantitative influence of the sun and the moon on the orbit of these satellites by using the method of the integral of the forces over the time.

1. Introduction

Global navigation satellite systems are a general denomination for constellations of navigation satellites, such as GPS (USA), GLONASS (Russia), Galileo (Europe), and Beidou (China), mainly placed in the medium earth orbit (MEO) region. Up to July 2013, the GPS had 31 active satellites in orbits with altitude near 20,000 km and 55 degrees of inclination. The Galileo system, at the same epoch, had four active satellites, with approximately 23,000 km of altitude and nominal inclination of 56 degrees. The Beidou system is composed by satellites in MEO, geosynchronous orbit (GEO), near circular, and inclined (55 degrees) geosynchronous orbit (IGSO). At the previously mentioned epoch this system had 14 satellites, four satellites being placed in MEO, five satellites in GEO, and five placed in IGSO [1]. The GLONASS system has 31 active satellites with inclinations ranging between 63° and 65° and 19,129 km of altitude. We will not consider the GLONASS system in this work, as well as the Beidou satellites with near circular orbits, because our technique is mainly devoted for satellites with inclinations around 56°, in which the 2 : 1 resonance between the argument of the perigee and the longitude of the ascending node is effective and causes a significant growth of eccentricity of the satellites.

In addition to the 2 : 1 perigee-ascending node resonance, the GNSS satellites are subject to resonances caused by the commensurability between the orbital period of these satellites and the rotation of the Earth (spin-orbit resonance). The GPS satellites are in a “deep” 2 : 1 resonance with the rotation of the Earth (the orbital period of these satellites is half of the rotation period of the Earth), and this resonance is dominated by the terms and of the geopotential [2]. The Galileo satellites are in 17 : 10 resonance, and the MEO satellites of the Beidou system are in 13 : 7 spin-orbit resonance [1]. The harmonic of the geopotential also plays an important role in the GNSS, because it causes short periodic variations in the semimajor axis [3, 4]. In general, the tesseral harmonics induce perturbations with short period and low amplitude in the orbital motion of the satellites. However, these terms may produce effects of high amplitude and long period [5, 6].

Due to the variety of orbits and resonances involved, the GNSS is a complex system, which turns the creation of a definitive strategy to dispose the satellites of this system into a hard, even impracticable, task. In this way, each constellation of the system adopts its own disposal strategy; for example, in the GPS system, the disposal strategy consists in changing the altitude of the nonoperational satellites to 500 km above or below their nominal orbits [7]. Due to the increase of debris around the Earth and the potential instability of orbits with inclination around 56 degrees [8, 9], several studies [1014] have been made to develop alternative strategies to dispose the satellites of the GNSS in a safer way. From the analysis of the initial conditions which lead to the previously mentioned instability caused by the 2 : 1 perigee-ascending node resonance, some works suggest moving the disposed satellites to regions such that the growth of the eccentricity does not allow the disposed satellites to invade the region of the operational satellites, at least for some acceptable time. This strategy may present some inconveniences, such as an accumulation of objects in the disposal region. On the other way, there are some works that suggest transferring the disposed satellites to appropriate orbits such that the new initial conditions lead to an increase of the eccentricity. The main goal is to lower continuously the periapsis radius of the disposed satellites, provoking their reentry in the atmosphere of the Earth. A positive characteristic of this approach is the decrease in the number of space debris, so reducing the collisional risk between inactive objects. Nowadays this number is not so high, but the continuous insertion of new satellites in this region certainly will change this scenario. Although the increase of the eccentricity may lead to a possible increase in the collisional risk between disposed and active objects, this risk can be minimized with a proper choice of initial conditions, which can accelerate the decay of these satellites. Therefore, the time that the satellite needs to reenter in the atmosphere is an important question related to the cost and success of this strategy. The effectiveness of disposal strategies for the GNSS satellites created so far can be found in [1].

In this work, we focus on the use of the strategy of eccentricity growth under the effect of the 2 : 1 perigee-ascending node resonance, which leads to the atmospheric reentry, and we compare this strategy with the direct discard using a Hohmann type transfer. This comparison is shown in terms of the velocity increment necessary to perform the atmospheric reentry. For the strategy that uses the resonance, grids of initial conditions were generated to evaluate the maximum eccentricity reached by the satellites after 250 years considering their nominal altitudes and also considering an apoapsis radius of 10,000 km above and below their nominal apoapsis radius. Testing different initial altitudes makes sense, since our purpose is to lower the periapsis of the orbit to drive the satellites to enter in the atmosphere within some prefixed time interval. We will see that these tests also show the “strength” of the resonance for different values of the semimajor axis of the orbit of the satellite, since, from the mathematical point of view, the 2 : 1 perigee-ascending node resonance always exists, no matter what is the altitude of the satellite.

As an additional study, we showed the quantitative influence of the sun and the moon on the orbit of these satellites by using the method of the integral of the forces over the time [15, 16]. The advantage of this method is that we can measure the total variation of the velocity caused by the moon and the sun without disregarding other perturbations; in this case, the perturbations of the geopotential and the radiation pressure of the sun which could interfere in the orbit of the satellite and, consequently, in the value of the perturbation of the sun and the moon. In the case of the sun, the total variation of the velocity is measured as a function of the perigee of the sun, and we show that this is an important element to take into account when planning the disposal of the GNSS satellites. For the moon, the most important element is its inclination.

2. The Resonance and Its Effects

In order to explain how the resonance 2 : 1 perigee-ascending node can affect the orbits of the GNSS satellites, initially we consider only the main disturbers, namely, the oblateness of the earth and the gravity of the sun. As usually, the oblateness is the dominant part of the disturbing function of a satellite in the MEO, and then we can use the single averaged form of the oblateness [13] to define the resonances which involve the perigee and the ascending node of the satellites in MEO. The expression for the single averaged oblateness is given by where is the mean motion of the satellite and is the mean equatorial radius of the Earth.

In this case, the main frequencies of the system are given by The ratio between and is For integer we have the special resonances which do not depend on the semi-major axis. These resonances usually affect the eccentricity [17]. For we have for or . Another classical resonance occurs when , so that , and this resonance affects the GLONASS satellites, but those are not considered in the present work.

The closer the inclination of the GNSS satellite is to 56.06°, the stronger will be the effect of the resonance on the eccentricity [13]. In order to see the effects of this resonance, the osculating equations of motion of a satellite will be integrated. As we mentioned before, only for this part of this work, as perturbations, we consider just the sun and the oblateness of the earth (the complete Cartesian equations involving the remaining perturbations will be given in the next section). Figure 1 shows the effects of the resonance on the eccentricity and on the resonant angle. Note that an initial small eccentricity reaches a significant increase.

For the inclination °, the dominant resonant term in the double averaged solar disturbing function is . Neglecting the remaining terms of (a full development of the single and double averaged disturbing function of the Sun is found in [13]), the Hamiltonian of the problem is where is the gravitational constant, is the mass of the sun, and , , , are, respectively, the magnitude of the position vector and the inclination and the ascending node of the sun: and .

Taking , , , mean anomaly, , , the set of the Delaunay variables and performing a Mathieu canonical transformation [18], we have where , , , and .

The orbit of the Sun is assumed to be Keplerian and we also considered , . In variables, our Hamiltonian is a one-degree of freedom problem, whose dynamics is very similar to the Lidov-Kozai resonance. In Figure 2 we consider an initial eccentricity and semimajor axis  km. This figure is very instructive: note that in the bottom part there is a large region where the satellite remains some finite time with very small eccentricity. It corresponds to the region where , which is the region of interest for the minimum eccentricity strategy study shown in [13]. On the other hand, we have the counterpart of this situation at the top of the figure: very high eccentricity, which occurs again for . This is our region of interest, because we intent to use this increase of the eccentricity to discard the satellites. It is worth noting that if the semimajor axis is high, the effect of the moon cannot be neglected, so that the problem is no more a one degree of freedom problem. In this case the search of the () pair such that the eccentricity increases must be done by integrating the complete equations of motion. A full characterization of this resonance is found in [13].

For a feasible disposal strategy, the time that the eccentricity demands to reach high values is a factor to take into account. In the next section we will consider this criterion in the search of initial conditions that causes an increase in the eccentricity.

3. Initial Conditions for the Strategy of Increasing the Eccentricity

In the previous section, we considered only the effects of the sun and the oblateness of the earth. Moreover, in the presence of the resonance, the main effects are now governed by the long term variations, since we considered the averaged problem. From a theoretical point of view, this averaged system is quite efficient to highlight the basic dynamics that affects the eccentricity of the GNSS satellites. However, for a more complete and realistic study, from now on, we need to include more perturbations.

In this section we want to find some special initial conditions such that the eccentricity of the satellites reaches high values within an interval of 250 years. The strategy to search for these particular initial conditions is guided from the theoretical approach described in the previous section.

We integrate the exact equations of motion of the GNSS satellites to be disposed (with inclination around 56°), under the effects of the sun, the moon, the geopotential with degree and order eight, and the radiation pressure coming from the Sun with eclipses. We use the RADAU integrator, a fast and precise numerical code [19], written in FORTRAN language and compiled using a Linux system.

The acceleration of the satellite using the exact system is where is the gravitational constant, , , are the masses of the earth, the sun, and the moon, respectively. , , are, in order, the position vector of the satellite, the sun, and the Moon. is the acceleration due to the disturbing geopotential up to order and degree eight, which is enough to ensure sufficient precision for this type of study [20]. We use a recursive geopotential model [2023], with the coefficients of the EGM2008 (Earth Gravity Model) [24]. is the acceleration due to the radiation pressure of the sun with eclipses [2527].

Figures 3 and 4 show the initial pair () and the maximum eccentricity reached in 250 years. The other initial elements are  km, , and °. We consider two values for the inclination of the moon, ° and 28.58°, with respect to the mean equator of the earth, and we can see that the initial value of the inclination of the moon is important. The two straight lines represent the exact condition . Note that now we have added the moon to our system, so that its effects can enhance the time variation of the eccentricity. The very smooth behavior (eccentricity-resonant angle) seen in Figure 2 certainly will be lost as well as the position and the magnitude of the maximum and minimum values of the eccentricity. This is expected since now the satellite will be under the action of the disturbance of the moon whose effect is about two times larger than the sun, as we will show in Section 5. We also should point out that, in most of the cases shown in Figures 3 and 4, the eccentricity does not reach high values in less than 200 years. Although this time it is below 250 years, this is not acceptable in practical terms, because it can increase the risk of collision between active and disposed satellites.

Figure 5 shows the results for a complete set of initial conditions of a GPS satellite (US catalog 22108) for the epoch April 4, 2012. The initial conditions of the moon and the sun are adjusted for the same epoch, because the initial conditions of the perigee of the sun are significant (the quantitative impact of the inclination of the moon and the perigee of the sun will be shown in Section 5). We note the growth of the regions with high eccentricities. However, again, in all cases the time for the eccentricity growth is very large (200 years). A more realistic condition requires not only the increase of the eccentricity, but the decay of the periapsis radius to values below than 6,578 km (periapsis altitude of 200 km). This value depends on the area to mass ratio of the satellite [28], and we consider that after the periapsis radius reach this value, the residual atmospheric drag will make the satellite reenter in the atmosphere. This value of periapsis is successfully used in studies of deorbiting of the GNSS satellites using low-thrust propulsion and solar radiation pressure [1]. Therefore, in our simulations we are interested in obtaining sets of initial conditions that satisfy this precise concept, so that after some prefixed time, we can ensure that the satellites really end up reentering into the atmosphere of the earth.

In an attempt to reduce the time that the eccentricity demands to reach high values, we simulate orbits with a pair of initial semimajor axis and eccentricity so that the initial value of the apoapsis radius of the satellites is 10,000 km below the value of the apoapsis radius of the orbits presented in Figure 5. The result is shown in Figure 6. Contrary to what is desired, note that the region with high values of eccentricity almost vanished. This can be explained by our averaged model (4): the smaller is the semimajor axis (altitude), the weaker is the effect of the resonance; that is, for satellites in lower orbits (small semi-major axis) the net effect of this resonance is negligible, in spite of its existence. Figure 6 clearly shows this fact. In this sense, we take the opposite scenario of Figure 6, an initial pair of semimajor axis and eccentricity so that the apoapsis radius is 10,000 km above the value used in Figure 5. Figure 7 shows these simulations and we can see that the regions with high eccentricity are larger than in the previous figures, as expected.

As explained before, an analysis of which initial conditions in Figure 7 drive the periapsis altitude to less than 200 km is necessary. Figures 8 and 9 show the complementary information for that. Figure 8 shows the regions where the periapsis altitude decreases to values less than or equal to 200 km within 250 years. Comparing Figure 8 with Figure 7, we see that although the eccentricity can reach high values (0.8 for example), this is not enough to lower the periapsis altitude to penetrate into the atmosphere (in our case 200 km). Moreover, to decide the best initial pair (perigee and longitude of ascending node) to discard a satellite, we also need to consider the amount of time which will be spent by the satellite to reach 200 km of periapsis altitude. The search of these conditions is easily done directly from the simulations. In fact, fortunately, in the plane, the region covered by the conditions we are looking for is significant. In Figure 9, in blue color, we have an interesting view of the main initial conditions. In Figure 10, the time evolution of one specific orbit is given. Thus, increasing the semimajor axis of the satellite (via orbital transfer) so that the apoapsis radius raises, for example, 10,000 km above its nominal value, is a viable way to exploit the resonance to disposal the GNSS satellites. In the next section we make a comparison of the cost (in terms of the increment of velocity) between the disposal via resonance and the direct deorbiting.

In order to ensure that the method of increasing the apoapsis radius to enhance the effect of the 2 : 1 perigee-ascending node resonance works for all the satellites of the GNSS with inclination around 56°, we present Figure 11, which show a satellite of the Galileo system (US catalog 28922) with a semimajor axis so that its apoapsis radius has 10,000 km above the nominal value. Analogously, Figure 12 shows a case for a satellite of the Beidou system (US catalog 31115) also with apoapsis radius 10,000 km above the nominal one.

Figures 13 and 14 show, respectively, the orbits of Figure 11 which achieve 200 km of periapsis altitude and the time spent for these satellites to achieve 200 km of periapsis altitude. Figures 15 and 16 present the same information for Figure 12.

4. The Direct Deorbiting Compared with the Disposal via Resonance

The transfers chosen for the direct deorbiting and for the increase of the semimajor axis in the technique of disposal via resonance are both of the Hohmann type [29], because, in both cases, they are the simplest and fastest solution.

Classically, the Hohmann transfer is coplanar and consists in the application of two increments of velocity : the first one () at the periapsis of the initial orbit (point A, in Figure 17), that puts the satellite in the transfer orbit (dashed), and the second one () applied at the apoapsis of the transfer orbit (point B) that puts the satellite in the final orbit, which has the apoapsis coincident with that of the transfer orbit.

Some differences appear between the original Hohmann transfer and the transfer that we use. Both orbits (initial and final) are elliptical in our transfer, including the one for the direct deorbiting. The transfer occurs from point B to point A, and the second does not exist, because we consider that the satellite has already entered in the atmosphere. is in the direction opposite to the motion of the satellite. For the disposal via resonance, the transfer occurs from the point A to point B, but is not required, because this transfer is performed only to increase the apoapsis of the initial orbit.

For the direct deorbiting, the apoapsis radius and the velocity at the point B of the orbit 3 are given by where and are the semimajor axis and the eccentricity of the orbit 3 and is the gravitational parameter of the Earth.

The periapsis radius is the equatorial radius of the earth and the apoapsis radius coincides with the apoapsis of the orbit 3. Thus, we can calculate the eccentricity , the angular momentum , and the velocity in the apoapsis of the orbit 2 [30]: The increment of velocity in B (and consequently the increment to perform the direct deorbiting), is

Table 1 shows the results of the direct deorbiting transfers applied on some elements of the GNSS that have inclinations around 56° and that resides in MEO. The initial conditions were collected in the Celestrak Catalog [31] and correspond to the epoch April 18, 2012 and are shown in Table A1 in Supplementary Material available online at http://dx.doi.org/10.1155/2014/382340. The initial conditions of the sun and the moon also correspond to the same epoch. The US catalog number is the satellite identification, , , are the semimajor axis, inclination, and eccentricity of the initial orbit, is the increment of velocity to execute the transfer, is the period of transfer, and the eccentricity of the final orbit. Figure 18 presents the tridimensional representation of the transfers for the satellites 20959 and 22014, where , , and are the components of the radius vector of these satellites.

The procedure to calculate the to increase the apoapsis in the disposal via resonance is quite similar to the previous calculation and is omitted here. The final orbit in this case is more eccentric and this fact increases the solar perturbation. We impose a limit to the augment of the apoapsis radius, which is equal to the one used to deorbit a satellite at the same altitude. We also limit the apoapsis height in order to avoid the crossing orbit with the GEO satellites. Table 2 shows some results of this transfer for the GNSS satellites, where is the increment in the apoapsis radius, and are the final semimajor axis and eccentricity, is the transfer time, and is the velocity increment that has to be applied to the satellite in order to increase the apoapsis radius in . The time of transfer does not include the time required to adjust the values of the pair , because this adjust occurs without any cost, by using the natural variation of these elements of the satellite, so, the time spent in this adjustment depends on how much the pair chosen is near the real values of these elements at the moment of the desired disposal.

Comparing the required by the direct deorbiting with the required by the disposal via resonance, we can note that the most economical strategy is the disposal via resonance, even for high values of .

In order to demonstrate the entire process of the disposal via resonance, we take, as an example, the case of the satellite 22108 using 10,000 km of to make the disposal via resonance. From Figures 8 and 9 we can choose the initial pair that makes the satellite reach 200 km of periapsis altitude in the shortest time. This pair must be chosen as close as possible to the current value of the satellite, so, we save time waiting for the satellite to reach the desirable value of , due to the natural precession of and . As soon as these values are reached, the maneuver to increase the apoapsis radius is performed, raising the apoapsis radius by 10,000 km. The consequence of this process is to place a satellite in a more eccentric orbit (natural consequence of the maneuver), which reaches 200 km of periapsis altitude in 50 years (in this case). This time is comparable with the time that the upper stages of the rockets that send the satellites Beidou into orbit takes to reenter in the atmosphere of the earth [32].

We can decrease the time of discard by taking other values of and keeping the same values of and from Figure 9. Figure 19 shows the time in which the periapsis radius reaches the value of 200 km of altitude during the numerical integration for several values of . We can verify that if we take a pair in the unstable region of Figure 8 (and consequently Figure 9) and increase the apoapsis radius, the time of disposal is reduced. However, from equal to 18,000, the time has an erratic distribution, which leads to an optimum value of around 10,000 km.

5. Quantitative Analysis of the Perturbation of the Moon and the Sun

The integration of the equations of motion (6) can give us interesting information about the influence of the initial conditions of the moon and the sun on the evolution of the orbit of the satellite. We have seen that the inclination of the moon changes the dynamics of the resonance studied, and the value of the initial argument of the perigee of the sun also can affect the orbits of these satellites. By using the method of the integral of the forces over the time that was recently proposed by Prado [15] that has many applications [16, 20], we can measure the magnitude of the acceleration due to the moon (as a function of its initial inclination) and to the sun (as a function of its initial argument of perigee) for several periods of time.

After the integration of the system, for an arbitrary period of time , by using the RADAU integrator, the magnitude of the accelerations due to the perturbation of the moon and to the sun is stored. We integrate those magnitude of the accelerations with the method of Simpson 1/3 [32] with respect to the period of time , obtaining the total velocity contribution of the disturbing body (Moon and Sun) during the time interval .

Figure 20 shows the total velocity applied to a satellite of the GNSS constellation due to the perturbation of the sun as a function of the initial perigee of the sun for 1 day (a), 100 days (b), 1 year (c), and 100 years (d). The predictability of the effects ends in a period of around 70 years. Figure 21 shows the specific total velocity due to the perturbation of the moon as a function of the initial inclination of the moon over 1 day (a), 100 days (b), 1 year (c), and 100 years (d). The predictability of the effects ends in a period of about 50 years. We can verify, comparing Figures 20 and 21, that the effect of the Moon is twice the effects of the Sun for the GNSS, as expected. Trough Figure 21 we also verify that the higher is the inclination of the Moon, the greater is the perturbation on the GNSS satellite, which confirms what we mentioned in the previous sections. Also, the values of the total velocity due to the sun and the moon, can be useful in the creation of laws of control for the GNSS satellites, in order to attenuate the perturbations of the Sun and the Moon.

6. The Effect of the Resonance over Beidou IGSO Satellites

The Beidou satellites are composed, as we mentioned before, by satellites in MEO, with inclination near 56°, GEO satellites, and satellites in inclined geosynchronous orbits (IGSO), with inclinations near 56°. As we saw in the previous sections, in this inclination, the resonance 2 : 1 perigee-ascending node can cause great increase in the eccentricity of the orbits of the satellites, even though these satellites are in near circular orbits. In the IGSO Beidou satellites, this increase in the eccentricity can amplify the risk of collision among functional satellites. Figures 22 and 23 show the maximum eccentricity achieved in 250 years (by integrating (6) with the RADAU integrator) for the Beidou satellite 37210 (0.7716° of inclination) and 37384 (55.6720° of inclination), respectively, both with initial conditions corresponding to the epoch April 18, 2012. We can see that, for the satellite 37210, there is no significant increase in the inclination over the 250 years of integration, whereas the eccentricity of the satellite 37384 reached high values. The comparison between Figures 22 and 23 shows us that the initial position of the pair may be taken into account when the IGSO satellites are positioned in their nominal orbits in order to avoid the effects of the 2 : 1 perigee-ascending node resonance.

7. Conclusions

We found an alternative method to dispose the satellites of the GNSS which are in MEO and have inclinations around 56°, by using the resonance associated to an orbital transfer to increase the apoapsis radius of the satellite to be disposable. We found sets of initial conditions such that the eccentricity of the satellite reaches high values, usually in less than 250 years. In comparison with the direct deorbiting, the alternative method has advantage in terms of fuel consumption, which is shown in terms of the increment of velocity that has to be applied to the satellite. The time required for the disposal is clearly shorter when compared to the time required to dispose the upper stages of the rockets used to place the satellites of the GNSS in orbit. Our technique provides the initial disposal conditions which lead to high values of eccentricity much faster than those presented in previous work. The additional information of Figures 8 and 9 guarantee the reentry in the atmosphere in the shortest possible time.

Although the focus of our work is the MEO satellites of the GNSS, and since the eccentricity grows faster when the satellite is closer to the sun, we wonder on the long term stability of the IGSO members of the Beidou constellation. Their inclinations are close to 56° and we show that the effects of the resonance are more intense for these satellites. Based on the results, we suggest that the effects of the resonance may be taken into account when the IGSO satellites are positioned in their nominal orbits in order to avoid the effects of the 2 : 1 perigee-ascending node resonance.

In this work we also used the method of the integration of the forces to measure the total velocity applied to a satellites of the GNSS due to the sun and the moon, taking into account different values for the perigee of the sun and for the inclination of the moon. This information may be useful, in the future, for creating laws of control for the MEO satellites, in order to avoid seasonality in the perturbation of the sun and the moon.

Conflict of Interests

The authors declare that there is no conflict of interests regarding to the publication of this paper.

Acknowledgments

The authors wish to express their appreciation for the support provided by Grants nos. 473387/2012-3, 304700/2009-6, and 305834/2013-4, from the National Council for Scientific and Technological Development (CNPq); Grants nos. 2012/21023-6 and 2014/06688-7, from São Paulo Research Foundation (FAPESP), and the financial support from the National Council for the Improvement of Higher Education (CAPES). This research was supported by resources supplied by the Center for Scientific Computing (NCC/GridUNESP) of the São Paulo State University (UNESP).

Supplementary Materials

The initial conditions of the active GNSS satellites for the epoch April 18, 2012, presenting the identification of the satellite (US catalog number and name), semimajor axis (𝑎), eccentricity (𝑒), inclination (𝐼), longitude of the ascending node (Ω), argument of perigee (𝜔), and mean anomaly (𝑙).The US catalog number is the same number of the NORAD catalog.

  1. Supplementary Material