Initial Semimajor Axis Research Articles

Influence of planets on debris discs in star clusters – II. The impact of stellar density

ABSTRACT We present numerical simulations of planetary systems in star clusters with different initial stellar densities, to investigate the impact of the density on debris disc dynamics. We use lps+ to combine N-body codes nbody6++gpu and rebound for simulations. We simulate debris discs with and without a Jupiter-mass planet at 50 au, in star clusters with $N=$ 1k–64k stars. The spatial range of the remaining planetary systems decreases with increasing N. As cluster density increases, the planet’s influence range first increases and then decreases. For debris particles escaping from planetary systems, the probability of their direct ejection from the star cluster decreases as their initial semimajor axis ($a_0$) or the cluster density increases. The eccentricity and inclination of surviving particles increase as cluster density increases. The presence of a planet leads to lower eccentricities and inclinations of surviving particles. The radial density distribution of the remaining discs decays exponentially in sparse clusters. We derive a general expression of the gravitational encounter rate. Our results are unable to directly explain the scarcity of debris discs in star clusters. Nevertheless, given that many planetary systems have multiple planets, the mechanism of the planet-cluster combined gravitational influence on the disc remains appealing as a potential explanation.

Triton’s Captured Youth: Tidal Heating Kept Triton Warm and Active for Billions of Years

Neptune’s moon Triton has two remarkable attributes: its retrograde orbit suggests that it was captured from the Kuiper Belt, and Triton has one of the youngest surfaces of all the icy satellites in the solar system. Soon after capture, Triton experienced strong diurnal tides raised by Neptune, which caused intense deformation, heating, and melting of its ice shell as its highly eccentric initial orbit was circularized. While previous studies have suggested that Triton’s orbit would have circularized early in solar system history, we show that internal feedbacks between tidal heating and ice shell melting significantly reduce the orbital evolution rate, causing strong tidal heating to persist for billions of years. We simulate Triton’s post-capture evolution over a range of initial semimajor axes and ice shell properties. We find that Triton’s ice shell would have been extremely thin (1–10 km) for a period of 1–4 billion years, with tidal stresses strong enough to fracture the entire ice shell down to the subsurface ocean. A final phase of intense geologic activity may have occurred after tidal dissipation waned, in which late-stage ice shell thickening caused ocean pressurization potentially sufficient to refracture the ice shell and push water to the surface. Such overpressurization could have caused recent massive cryovolcanic resurfacing, perhaps explaining Triton’s geologically young surface. It is therefore possible that Triton’s youthful surface and its origin as a captured satellite may in fact be related. A long-lived subsurface ocean and extended thin ice shell period also greatly increase Triton’s astrobiological potential.

Free-floating Planets, Survivor Planets, Captured Planets, and Binary Planets from Stellar Flybys

In star clusters, close stellar encounters can strongly impact the architecture of a planetary system or even destroy it. We present a systematic study of the effects of stellar flybys on two-planet systems. When such a system experiences flybys, one or both planets can be ejected, forming free-floating planets (FFPs), captured planets (CPs) around the flyby star, and free-floating binary planets (BPs); the remaining single-surviving planets (SSPs) can have their orbital radii and eccentricities greatly changed. Through numerical experiments, we calculate the formation fractions (or branching ratios) of FFPs, SSPs, CPs, and BPs as a function of the pericenter distance of the flyby, and use them to derive analytical expressions for the formation rates of FFPs, SSPs, CPs and BPs in general cluster environments. We find that the production rates of FFPs and SSPs are similar (for the initial planet semimajor axis ratio a 1/a 2 = 0.6–0.8), while the rate for CPs is a few times smaller. The formation fraction of BPs depends strongly on a 1/a 2 and on the planet masses. For Jupiter-mass planets, the formation fraction of BPs is always less than 1% (for a 1/a 2 = 0.8) and typically much smaller (≲0.2% for a 1/a 2 ≲ 0.7). The fraction remains less than 1% when considering 4M J planets. Overall, when averaging over all flybys, the production rate of BPs is less than 0.1% of that of FFPs. We also derive the velocity distribution of FFPs produced by stellar flybys, and the orbital parameter distributions of SSPs, CPs, and BPs. These results can be used in future studies of exotic planets (including FFPs) and planetary systems.

Scattering of Giant Planets and Implications for the Origin of the Hierarchical and Eccentric Two-planet System GJ 1148

The GJ 1148 system has two Saturn-mass planets orbiting around an M dwarf star on hierarchical and eccentric orbits, with orbital period ratio of 13 and eccentricities of both planets of 0.375. The inner planet is in the regime of eccentric warm Jupiters. We perform numerical experiments to study the planet–planet scattering scenario for the origin of this orbital architecture. We consider a third planet of 0.1M J (Jupiter's mass) in the initial GJ 1148 system with initial orbital separations of 3.5, 4, and 4.5 mutual Hill radii and initial semimajor axis of the innermost planet in the range of 0.10–0.50 au. The majority of scattering results in planet–planet collisions, followed by planet ejections, and planet–star close approaches. Among them, only planet ejections produce eccentric and widely separated two-planet systems, with some having similar orbital properties to the GJ 1148 system. We also examine the effects of general relativistic apsidal precession and a higher mass of 0.227M J for the third planet. The simulation results suggest that the GJ 1148 system may have lost a giant planet. We also perform simulations of the general problem of the origin of warm Jupiters by planet–planet scattering. As in the GJ 1148 simulations, a nontrivial number of stable two-planet systems are produced by ejection, which disagrees with the result from a previous study showing that two-planet systems arise exclusively through planet–planet collisions.

The chaotic four-body problem in Newtonian gravity – II. An ansatz-based approach to analytical solutions

ABSTRACT In this paper, we continue our analysis of the chaotic four-body problem and our study of binary–binary interactions in star clusters. We present a general ansatz-based analytical treatment using statistical mechanics, where each outcome of the four-body problem is regarded as some variation of the three-body problem. For example, when two single stars are produced (the 2 + 1 + 1 outcome), each ejection event is modelled as its own three-body interaction by assuming that the ejections are well separated in time. This is a generalization of the approach adopted in Paper I, based on the density-of-states formalism. There are three possible outcomes for the four-body problem with negative total energies: 2 + 2, 2 + 1 + 1, and 3 + 1. For each outcome, we apply an ansatz-based approach to deriving analytical distribution functions that describe the properties of the products of chaotic four-body interactions involving point particles. To test our theoretical distributions, we perform a set of scattering simulations in the equal-mass point-particle limit using FEWBODY, where we vary the initial ratio of binary semimajor axes. We compare our final theoretical distributions to the simulations for each particular scenario, finding consistently good agreement between the two. The highlights of our results include that binary–binary scatterings act to systematically destroy binaries producing instead a binary and two ejected stars (when the initial binary semimajor axes are similar) or a stable triple (when the initial semimajor axes are very different). The 2 + 2 outcome produces the widest binaries, and the 2 + 1 + 1 outcome produces the most compact binaries.

Stability of coorbital planets around binaries

In previous hydrodynamical simulations, we found a mechanism for nearly circular binary stars, such as Kepler-413, to trap two planets in a stable 1:1 resonance. Therefore, the stability of coorbital configurations becomes a relevant question for planet formation around binary stars. For this work, we investigated the coorbital planet stability using a Kepler-413 analogue as an example and then expanded the parameters to study a general n-body stability of planet pairs in eccentric horseshoe orbits around binaries. The stability was tested by evolving the planet orbits for 105 binary periods with varying initial semi-major axes and planet eccentricities. The unstable region of a single circumbinary planet is used as a comparison to the investigated coorbital configurations in this work. We confirm previous findings on the stability of single planets and find a first order linear relation between the orbit eccentricity ep and pericentre to identify stable orbits for various binary configurations. Such a linear relation is also found for the stability of 1:1 resonant planets around binaries. Stable orbits for eccentric horseshoe configurations exist with a pericentre closer than seven binary separations and, in the case of Kepler-413, the pericentre of the first stable orbit can be approximated by rc,peri = (2.90 ep + 2.46) abin.

Orbit injection of planet-crossing asteroids

ABSTRACT Solar system Centaurs originate in trans-Neptunian space from where planet orbit crossing events inject their orbits inside the giant planets’ domain. Here, we examine this injection process in the three-body problem by studying the orbital evolution of trans-Neptunian asteroids located at Neptune’s collision singularity as a function of the Tisserand invariant, T. Two injection modes are found, one for T > 0.1, or equivalently prograde inclinations far from the planet, where unstable motion dominates injection, and another for T ≤ 0.1, or equivalently polar and retrograde inclinations far from the planet, where stable motion dominates injection. The injection modes are independent of the initial semimajor axis and the dynamical time at the collision singularity. The simulations uncovered a region in the polar corridor where the dynamical time exceeds the Solar system’s age suggesting the possibility of long-lived primordial polar trans-Neptunian reservoirs that supply Centaurs to the giant planets’ domain.

A dependence of binary and planetary system destruction on subtle variations in the substructure in young star-forming regions

ABSTRACT Simulations of the effects of stellar fly-bys on planetary systems in star-forming regions show a strong dependence on subtle variations in the initial spatial and kinematic substructure of the regions. For similar stellar densities, the more substructured star-forming regions disrupt up to a factor of 2 more planetary systems. We extend this work to look at the effects of substructure on stellar binary populations. We present N-body simulations of substructured, and non-substructured (smooth) star-forming regions in which we place different populations of stellar binaries. We find that for binary populations that are dominated by close (<100 au) systems, a higher proportion are destroyed in substructured regions. However, for wider systems (>100 au), a higher proportion are destroyed in smooth regions. The difference is likely due to the hard–soft or fast–slow boundary for binary destruction. Hard (fast/close) binaries are more likely to be destroyed in environments with a small velocity dispersion (kinematically substructured regions), whereas soft (slow/wide) binaries are more likely to be destroyed in environments with higher velocity dispersions (non-kinematically substructured regions). Due to the vast range of stellar binary semimajor axes in star-forming regions (10−2 to 104 au), these differences are small and hence unlikely to be observable. However, planetary systems have a much smaller initial semimajor axis range (likely ∼1–100 au for gas giants) and here the difference in the fraction of companions due to substructure could be observed if the star-forming regions that disrupt planetary systems formed with similar stellar densities.

The Impact of Tidal Migration of Hot Jupiters on the Rotation of Sun-like Main-sequence Stars

The tidal interactions of planets affect the stellar evolutionary status and the constraint of their physical parameters by gyrochronology. In this work, we incorporate the tidal interaction and magnetic braking of the stellar wind into MESA and calculate a large grid of 25,000 models, covering planets with masses of 0.1–13.0 M J with different orbital distances that orbit late-type stars of different metallicities. We also explore the effect of different stellar initial rotations on the tidal interactions. Our results show that in the case of tidal inward migration, the stellar rotation periods are always lower than that of the star without planet before the planet is engulfed and the difference in the rotation period of its host star always increases with time. After the planet is engulfed, the stellar rotation periods are still lower than that of star without planet, but the difference of periods can be quickly eliminated if the star has a thick convective envelope (smaller mass and larger metallicity), regardless of the mass of the planet and the initial rotation period of the star. In the case of stars with thinner convective envelopes (larger mass and smaller metallicity), the stars will be spun up and remain the faster rotation in a long time. Meanwhile, the planet is easily swallowed and the period differences are large if the initial rotation period of its host star is higher. Finally, we also study the evolution of WASP-19 and estimate the range of tidal quality parameter and the initial semimajor axis as (0.035 ± 0.004) au.

Modelling uncertainties in wide binary constraints on primordial black holes

ABSTRACT Dark matter in the form of compact objects with mass Mco ≳ 10 M⊙ can be constrained by its dynamical effects on wide binary stars. Motivated by the recent interest in primordial black hole dark matter, we revisit the theoretical modelling involved in these constraints. We improve on previous studies in several ways. Specifically, we (i) implement a physically motivated model for the initial wide-binary semimajor axis distribution, (ii) include unbound binaries, and (iii) take into account the uncertainty in the relationship between semimajor axis and observed angular separation. These effects all tend to increase the predicted number of wide binaries (for a given compact object population). Therefore, the constraints on the halo fraction in compact objects, fco, are significantly weakened. For the wide binary sample used in the most recent calculation of the constraints, we find the fraction of halo dark matter in compact objects is fco < 1 for $M_{\rm co} \approx 300 \, \mathrm{ M}_{\odot }$, tightening with increasing Mco to fco < 0.26 for $M_{\rm co} \gtrsim 1000 \, \mathrm{ M}_{\odot }$.

Residual eccentricity of an Earth-like planet orbiting a red giant Sun

Context. The late phases of the orbital evolution of an Earth-like planet around a Sun-like star are revisited in order to consider the effect of density fluctuations associated with convective motions inside the star. Aims. Such fluctuations produce a random perturbation of the stellar outer gravitational field that excites a small residual eccentricity in the orbit of the planet. This counteracts the effects of tides, which tend to circularize the orbit. Methods. We computed the power spectrum of the outer gravitational field fluctuations of the star in the quadrupole approximation and studied their effects on the orbit of the planet using a perturbative approach. The residual eccentricity is found to be a stochastic variable showing a Gaussian distribution. Results. Adopting a model of the stellar evolution of our Sun computed with Modules for Experiments in Stellar Astrophysics (MESA), we find that the Earth will be engulfed by the Sun when it is close to the tip of the red giant branch phase of evolution. We find a maximum mean value of the residual eccentricity of ~0.026 immediately before engulfment. Considering an Earth-mass planet with an initial orbital semimajor axis sufficiently large to escape engulfment, we find that the mean value of the residual eccentricity is greater than 0.01 for an initial separation of up to ~l.4 au. Conclusions. The engulfment of the Earth by the red giant Sun is found to be a stochastic process instead of being deterministic as assumed in previous studies. If an Earth-like planet escapes engulfment, its orbit around its remnant white dwarf (WD) star will be moderately eccentric. Such a residual eccentricity of on the order of a few hundredths can play a relevant role in sustaining the pollution of the WD atmosphere by asteroids and comets, as observed in several objects.

Spatial distribution of isotopes and compositional mixing in the inner protoplanetary disk

The mass-independent isotopic signatures of planetary bodies have been widely used to trace the mixing process that occurred during planet formation. The observed isotopic variations among meteorite parent bodies have been further linked to the modeled mass-weighted mean initial semimajor axes in N-body simulations, assuming a spatial isotopic gradient in the inner protoplanetary disk. However, nucleosynthetic isotopic anomalies of nonvolatile elements and mass-independent oxygen isotopic variation (∆17O) show different relationships with distance from the Sun. Therefore, it is crucial to know whether isotopes were distributed systematically with heliocentric distance in the inner protoplanetary disk. In this study, we performed N-body simulations on compositional mixing during the collisional accretion and migration of planetary bodies to investigate the spatial distributions of Cr and O isotopes in the inner protoplanetary disk. The modeled mass-weighted mean initial semimajor axes of the parent bodies of noncarbonaceous (NC) meteorites and terrestrial planets were used to calculate the isotopic compositions of these bodies. Our simulations successfully reproduced the observed nucleosynthetic Cr isotopic anomaly among Earth, Mars, and the NC meteorite parent bodies, consistent with a spatial gradient of isotopic anomalies in the inner disk. Asteroids originating from different regions in the inner disk were transported to the main belt in our simulations, resulting in the Cr isotopic anomaly variation of the NC meteorite parent bodies. However, the ∆17O distribution among the terrestrial planets and the NC meteorite parent bodies could not be reproduced assuming a ∆17O gradient in the inner protoplanetary disk. The spatial gradient of the nucleosynthetic isotopic anomaly may be a result of changing isotopic compositions in the infalling materials, or reflect the progressive thermal processing of presolar materials. In contrast, the absence of a ∆17O gradient reflects that the oxygen isotopic mass-independent fractionation might have altered the spatial distribution of the nucleosynthetic ∆17O variation in the inner protoplanetary disk before protoplanets formed.

Production of hot Jupiter candidates from high-eccentricity mechanisms for different initial planetary mass configurations

ABSTRACT Hot Jupiters (HJs) are giant planets with orbital periods of the order of a few days with semimajor axis within ∼0.1 au. Several theories have been invoked in order to explain the origin of this type of planets, one of them being the high-eccentricity migration. This migration can occur through different high-eccentricity mechanisms. Our investigation focused on six different kinds of high-eccentricity mechanisms, namely, direct dispersion, coplanar, Kozai–Lidov, secular chaos, E1 and E2 mechanisms. We investigated the efficiency of these mechanisms for the production of HJ candidates in multiplanet systems initially tightly-packed in the semimajor axis, considering a large set of numerical simulations of the exact equations of motion in the context of the N-body problem. In particular, we analyzed the sensitivity of our results to the initial number of planets, the initial semimajor axis of the innermost planetary orbit, the initial configuration of planetary masses, and to the inclusion of general relativity (GR) effects. We found that the E1 mechanism is the most efficient in producing HJ candidates both in simulations with and without the contribution of GR, followed by the Kozai–Lidov and E2 mechanisms. Our results also revealed that, except for the initial equal planetary mass configuration, the E1 mechanism was notably efficient in the other initial planetary mass configurations considered in this work. Finally, we investigated the production of HJ candidates with prograde, retrograde, and alternating orbits. According to our statistical analysis, the Kozai–Lidov mechanism has the highest probability of significantly exciting the orbital inclinations of the HJ candidates.

Analysis of ODMSP-compliant near-circular GPS disposal orbits and resulting long-term collision risk

The U.S. Government Orbital Debris Mitigation Standard Practices (ODMSP) released in December 2019 include several disposal options that are applicable to GPS satellites. One option is to leave a satellite on a near-circular disposal orbit with limited long-term eccentricity growth. Many currently operational GPS satellites are planned to move to near-circular disposal orbits above the GPS operational constellation and below the BeiDou Satellite System (BDS). The conditions in the ODMSP specify that GPS disposal orbits avoid crossing the semi-synchronous zone (altitude range 20,182 +/- 300 km, occupied by GPS) or the BDS operational altitude for 100 years. An analysis was performed to determine ways to achieve ODMSP-compliance of near-circular GPS disposal orbits and limit associated collision risk. Disposal orbit propagations over 100 years were used to determine the available initial semi-major axis (SMA) range for compliant disposal orbits. Next, the GPS constellation and satellite disposals were simulated over a 200-year time frame. The disposal orbits were propagated and the collision probability between disposed satellites was then evaluated. Results show an available initial SMA range of 363 km for an initial eccentricity of 0.0003 and a range of only 3 km for an initial eccentricity of 0.001, which cover worst-case eccentricity growth scenarios. Results also show that a disposal strategy with initial eccentricity of 0.0003 and initial SMA uniformly spread over 363 km has the lowest collision probability for future disposed GPS satellites.

Forming short period sub-stellar companions in 47 Tucanae – II. Analytical expressions for the orbital evolution of planets in dense environments

ABSTRACT Short period, massive planets, known as hot Jupiters (HJs), have been discovered around ∼1 per cent of local field stars. The inward migration necessary to produce HJs may be ‘low eccentricity’, due to torques in the primordial disc, or ‘high eccentricity’ (HEM). The latter involves exciting high orbital eccentricity, allowing sufficiently close passages with the host star to raise circularizing tides in the planet. We present an analytical framework for quantifying the role of dynamical encounters in high density environments during HEM. We show that encounters can enhance or suppress HEM, depending on the local stellar density and the initial semimajor axis a0. For moderate densities, external perturbations can excite large eccentricities that allow a planet to circularize over the stellar lifetime. At extremely high densities, these perturbations can instead result in tidal disruption of the planet, thus yielding no HJ. This may explain the apparent excess of HJs in M67 compared with their local field star abundance versus their apparent deficit in 47 Tuc. Applying our analytical framework, we demonstrate that for an initial massive planet population similar to the field, the expected HJ occurrence rate in 47 Tuc is fHJ = 2.2 × 10−3, which remains consistent with present constraints. Future large (sample sizes ≳105) or sensitive transit surveys of stars in globular clusters are required to refute the hypothesis that the initial planet population is similar to the solar neighbourhood average. Non-detection in such surveys would have broad consequences for planet formation theory, implying planet formation rates in globular clusters must be suppressed across a wide range of a0.

The growth and migration of massive planets under the influence of external photoevaporation

ABSTRACT The formation of gas giant planets must occur during the first few Myr of a star’s lifetime, when the protoplanetary disc still contains sufficient gas to be accreted on to the planetary core. The majority of protoplanetary discs are exposed to strong ultraviolet irradiation from nearby massive stars, which drives winds and depletes the mass budget for planet formation. It remains unclear to what degree external photoevaporation affects the formation of massive planets. In this work, we present a simple one dimensional model for the growth and migration of a massive planet under the influence of external FUV fields. We find that even moderate FUV fluxes $F_\mathrm{FUV}\gtrsim 100 \, G_0$ have a strong influence on planet mass and migration. By decreasing the local surface density and shutting off accretion on to the planet, external irradiation suppresses planet masses and halts migration early. The distribution of typical stellar birth environments can therefore produce an anticorrelation between semi-major axis and planet mass, which may explain the apparent decrease in planet occurrence rates at orbital periods Porb ≳ 103 d. Even moderate fluxes FFUV strongly suppress giant planet formation and inward migration for any initial semi-major axis if the stellar host mass $M_*\lesssim 0.5\, {\rm M}_\odot$, consistent with findings that massive planet occurrence is much lower around such stars. The outcomes of our prescription for external disc depletion show significant differences to the current approximation adopted in state-of-the-art population synthesis models, motivating future careful treatment of this important process.

Probing Planets with Exomoons: The Cases of Kepler-1708 b and Kepler-1625 b

Abstract The tidal interactions between a planet and moon can provide insight into the properties of the host planet. The recent exomoon candidates Kepler-1708 b-i and Kepler-1625 b-i are Neptune-sized satellites orbiting Jupiter-like planets and provide an opportunity to apply such methods. We show that if the tidal migration time is roughly equal to the age of these systems, then the tidal dissipation factor Q for the planets Kepler-1708 b and Kepler-1625 b have values of ∼3 × 105–3 × 106 and ∼1.5 × 105–4 × 105, respectively. In each case, these are consistent with estimates for gas-giant planets. Even though some work suggests an especially large semimajor axis for Kepler-1625 b-i, we find that this would imply a surprisingly low Q ∼ 2000 for a gas giant unless the moon formed at essentially its current position. More detailed predictions for the moons’ initial semimajor axis could provide even better constraints on Q, and we discuss the formation scenarios for a moon in this context. Similar arguments can be used as more exomoons are discovered in the future to constrain exoplanet interior properties. This could be especially useful for exoplanets near the sub-Neptune/super-Earth radius gap where the planet structure is uncertain.

Planet–planet scattering in presence of a companion star

ABSTRACT Planet–planet (P–P) scattering is a leading dynamical mechanism invoked to explain the present orbital distribution of exoplanets. Many stars belong to binary systems; therefore, it is important to understand how this mechanism works in the presence of a companion star. We focus on systems of three planets orbiting the primary star and estimate the time-scale for instability, finding that it scales with the Keplerian period for systems that have the same ratio between inner planet and binary semimajor axes. An empirical formula is also derived from simulations to estimate how the binary eccentricity affects the extent of the stability region. The presence of the secondary star affects the P–P scattering outcomes, causing a broadening of the final distribution in semimajor axis of the inner planet as some of the orbital energy of the planets is absorbed by the companion star. Repeated approaches to the secondary star also cause a significant reduction in the frequency of surviving two-planet systems in particular for larger values of the inner planet semimajor axis. The formation of Kozai states with the companion star increases the number of planets that may be tidally circularized. To predict the possible final distribution of planets in binaries, we have performed a large number of simulations where the initial semimajor axis of the inner planets is chosen randomly. For small values of the binary semimajor axis, the higher frequency of collision alters the final planet orbital distributions that, however, beyond 50 au appear to be scalable to wider binary separations.

Symmetry Breaking in Dynamical Encounters in the Disks of Active Galactic Nuclei

Abstract The disks of active galactic nuclei (AGNs) may be important sites of binary black hole (BBH) mergers. Here we show via numerical experiments with the high-accuracy, high-precision code SpaceHub that broken symmetry in dynamical encounters in AGN disks can lead to asymmetry between prograde and retrograde BBH mergers. The direction of the hardening asymmetry depends on the initial binary semimajor axis. Under the assumption that the spin of the BHs becomes aligned with the angular momentum of the disk on a short timescale compared with the encounter timescale, an asymmetric distribution of mass-weighted projected spin χ eff is predicted in LIGO–Virgo detections of BBH mergers from AGN disks. In particular, this model predicts that positive χ eff BBH mergers are most likely for encounters with massive tertiaries in migration traps at radial distances ≳500–600 gravitational radii.

Air-breathing electric propulsion: Flight envelope identification and development of control for long-term orbital stability

Air-breathing electric propulsion (ABEP) enables long duration missions at very low orbital altitudes through the use of drag compensation. A system-level spacecraft model is developed, using the interaction between thruster, intake and solar arrays, and coupled to a calculation of the drag. A quadratic solution is found for specific impulse and evaluated to identify the thruster performance required for drag-compensation at varying altitudes. An upper altitude limit around 190 km is based on a minimum thruster propellant density, resulting in required thruster performance values of Isp>3000s and T/P>8mN/kW for a realistic ABEP spacecraft. The orbit of an air-breathing spacecraft is propagated with time, which highlights the prescribed orbit eccentricity due to non-spherical gravity and therefore an increased variability in the atmospheric conditions. A thruster control law is introduced which avoids a divergent altitude behaviour by preventing thruster firings around the orbit periapsis, as well as adding robustness against atmospheric changes due to season and solar activity. Through the use of an initial frozen orbit, thruster control and an augmented T/P, a stable long-term profile is demonstrated based on the performance data of a gridded-ion thruster tested with atmospheric propellants. An initial mean semi-major axis altitude of 200 km relative to the equatorial Earth radius, a spacecraft mass of 200 kg, Isp=5455s and T/P=23mN/kW, results in an altitude range of around 10 km at altitudes of 160–183km during a period of medium to high solar activity.