Bound circumplanetary orbits under the influence of radiation pressure: Application to dust in directly imaged exoplanet systems

We examine the population of simply periodic orbits in the Hill problem with radiation pressure included, in order to understand the distribution of gravitationally bound dust in orbit around a planet. We study a wide range of radiation pressure strengths, which requires the inclusion of additional terms beyond those discussed in previous analyses of this problem. In particular, our solutions reveal two distinct populations of stable wide, retrograde, orbits, as opposed to the single family that exists in the purely gravitational problem. We use the result of these calculations to study the observational shape of dust populations bound to extrasolar planets, that might be observable in scattered or reradiated light. In particular, we find that such dusty clouds should be elongated along the star--planet axis, and that the elongation of the bound population increases with $\beta$, a measure of the strength of the radiation pressure. As an application of this model, we consider the properties of the Fomalhaut system. The unusual orbital properties of the object Fomalhaut~b can be explained if the observed light was scattered by dust that was released from an object in a quasi-satellite orbit about a planet located in, or near, the observed debris ring. Within the context of the model of Hayakawa \&Hansen (2023), we find that the dust cloud around such a planet is still approximately an order of magnitude fainter than the limits set by current JWST data.


INTRODUCTION
The late stages of planetary system assembly are expected to result in the production of copious amounts of dust, which can be observed due to its capacity for reprocessing the light from the central star.Imaging of this population of dust, either in scattered light or thermal emission, can provide information on the properties of the planetary system by virtue of the sensitivity of the dust to the gravitational influence of the planets in the system.
In Hayakawa & Hansen (2023), we present a model for the origin of thin dust rings, in which the dust is generated in irregular satellite systems.As an alternative to the common 'birth ring' model for dust populations -which posits an origin in a ring of colliding planetesimals -we invoke a population of irregular satellites that generates the dust in the local vicinity of the host planet.The effect of radiation pressure from the central star is to shift the pseudopotential so that the collisional cascade is halted when dust escapes via the L 2 point and goes into orbit exterior to the planet.We demonstrated that this dust, regulated by the contours of the pseudopotential, naturally yields a thin dust ring similar to those observed in some systems.
email hansen@astro.ucla.eduNot all of the dust immediately escapes to orbits exterior to the planet.The question we wish to address in this paper is whether other aspects of the model are amenable to observation.In particular, we wish to calculate the expected properties of dust trapped in stable orbits, and whether they can be observed.Perhaps the best known example of a thin ring system is that in orbit around the star Fomalhaut (Kalas et al. 2005).There have also been claims of a planetary object in this system (Kalas et al. 2008(Kalas et al. , 2013) ) although this has also been argued to be a dust cloud not associated with a planet (Currie et al. 2012;Galicher et al. 2013).The optical colors of this object suggest that the origin of the emission is scattered light from the primary, and Kennedy & Wyatt (2011) have suggested that the emission is the result of the collisional grinding down of a population of irregular satellites as posited above.
With the continued imaging of dust systems with HST, and the new capabilities of JWST becoming available, it is timely to reconsider the observability of dust bound, or in close orbit, around a planet, and to evaluate the expected signatures of this dust.Kennedy & Wyatt (2011) assume a distribution of dust that follows the orbits of the parent bodies about the planet, but particles in orbit around a planet will be subject to a variety of forces due to the radiation from the central star (Burns et al. 1979).
Thus, in § 2 we define the version of the Photogravitational Hill problem needed for our case.In § 3 we then generalise the families of periodic orbits discussed in the original Hill problem by Henon (1969Henon ( , 1970)).In § 4 we then discuss how these orbits will manifest themselves in the full restricted three-body co-ordinate system and how this might be observable.

THE PHOTOGRAVITATIONAL HILL PROBLEM AND HENON ORBIT FAMILIES
The motion of dust in stellar or planetary systems is well suited to description in the limit of the restricted three-body problem, as the mass of dust particles is so small as to be accurately described as a test particle.In the case where the two massive bodies have a circular orbit, the dynamics is well-known and regulated by a conserved integral, the Jacobi integral (e.g.Murray & Dermott 2000).
In the limit where the mass ratio µ between the two massive bodies is small, the dynamics in and around the sphere of influence of the less massive body can be fruitfully described in the context of the Hill problem (Hill 1878a,b,c), wherein the dynamics is described in a coordinate frame centered on the smaller body (the planet).Henon (1969) presented an elegant analysis of the simply periodic orbits in a rescaled version of Hill's problem in the limit where µ → 0. We wish here to revisit this question when the test particles are subject to a radiation pressure force from the central star, in order to describe the kinds of structures we might observe in a population of dust particles generated in close proximity to a planet.
In the original restricted three body problem, the massive bodies are separated by a semi-major axis of unity, and the center of mass in the co-rotating frame lies at (0, 0, 0).If we wish to describe the motion about the planet, it is helpful to move the origin of the co-ordinate system to the position of the planet (at x = 1 − µ and y = 0).Henon (1969Henon ( , 1970) ) noticed that one could describe the dynamics quite generally by scaling out the mass ratio in the definition of the circumplanetary coordinates ξ, η and ζ.We follow this approach, resulting in the circumstellar-to-circumplanetary coordinate transformation x = 1 − µ + µ 1/3 ξ, y = µ 1/3 η and z = µ 1/3 ζ.The effect of the radiation pressure is to reduce the effective gravity of the central, luminous, object by a factor (1 − β), where β parameterises the strength of the radiation pressure.The resulting dynamics of the test particle is given by the equations where ∆ 2 = ξ 2 + η 2 + ζ 2 and we have dropped terms of higher order than µ 1/3 .We have also restricted ourselves to the radial component of the radiation pressure force.The scattering of light by an object in motion around the source also results in an azimuthal component, which gives rise to the Poynting-Robertson drag.However, this force is smaller than the radial component (Burns et al. 1979) by a factor ∼ v/c (where v is the particle's orbital motion).For objects on scales ∼ 10-100 AU, this is < 3 × 10 −5 .This is much smaller than terms of order µ 1/3 and thus is neglected here.The dynamics of small particles in the restricted three body problem with radiation pressure has been described in many prior studies (Radzievskii 1950(Radzievskii , 1953;;Schuerman 1980;Simmons et al. 1985;Kushvah 2008;Zotos 2015, e.g).These analyses describe how the nature and stability of the Lagrangian equilibria evolve as the balance between gravity and radiation changes.Of particular interest to us are studies such as those by Markellos et al. (2000); Kanavos et al. (2002); Perdiou et al. (2012), which seek to describe the dynamics of particles near a planet while subject to radiation pressure, and which catalog the kinds of simply periodic orbits that arise.
However, these prior analyses of the 'photogravitational problem' make an approximation that is not appropriate for the application we seek.After the derivation of the above equations, the authors next step is to set β = µ 1/3 Q 1 and then take the limit µ → 0. This achieves an elegant rescaling that removes the mass dependance, in the same manner as the original classic work of Henon.However, it comes at a cost -in order for Q 1 to be a constant, β → 0 in lockstep with µ, so this only applies in the limit of small radiation pressure.This may be appropriate to applications such as describing the slow drift of spacecraft under the influence of stellar radiation (e.g.Giancotti et al. 2014;García Yárnoz et al. 2015), but it is not appropriate for describing the dynamics of dust, where the value of β need not be infinitesimal.
Thus, we opt to keep β finite, so that we may retain our ability to describe the dynamics of dust experiencing significant perturbations due to radiation pressure.This comes at the price that our description is no longer scale free -for a given β we must also specify the value of µ.We will, however, denote Q ′ = β/µ 1/3 , as prior authors have done, to enable direct comparison between their results and ours.
The equations ( 1) and (2) still admit an equivalent of the Jacobi integral, now given as (4) The difference between these equations and those of the traditional 'photogravitational Hill problem' (Markellos et al. 2000;Kanavos et al. 2002) leads to some qualitative differences, as we will see below (although they represent a subset of the equilibria discussed in the more general restricted three body problem - Simmons et al. (1985)).

The Planar Equilibrium Points
The equilibrium points of the equations for the restricted three body problem define the transitions between different classes of dynamics.In the classic reduced three body problem, there are five equilibrium pointsthe Lagrange equilibria.Three of them -the colinear points L 1 , L 2 and L 3 -lie on the line between the two gravitating masses in the co-rotating frame.The other two -the triangular points L 4 and L 5 -lie on either side of the secondary, at angles ±60 • .
The introduction of the radiation force shifts the positions of the equilibria and even introduces a new set of out-of-plane equilibria in the case of very strong radiation pressure (Radzievskii 1950(Radzievskii , 1953;;Schuerman 1980;Simmons et al. 1985).In order to understand the dynamics of dust particles, we wish to understand the implications of these changes in the Hill limit (i.e. when µ ≪ 1).

The Colinear Points
In the purely gravitational case (β = 0), the equations in the Hill limit require η = 0 for equilibria to exist, which implies ∆ = |ξ|.When we feed this condition into equation (1), we find that the colinear equilibria must satisfy which can be written as two different cubic equations for positive and negative ξ, with different signs for the second term.
In the case of β = Q ′ = 0, this yields the usual criterion ξ = ±3 −1/3 .Thus, the L 1 (negative ξ) and L 2 (positive ξ) equilibria lie at the same distance from the planet, on either side.The distant L 3 equilibrium vanishes in the limit.For finite β, there are two different cubic equations, depending on the sign of the Q ′ term in equation ( 5), which will yield different values of |ξ| on either side of the planet.Although cubic equations can yield multiple solutions, we note that not all of the solutions to equation ( 5) are valid, because valid solutions must have the correct relationship between ξ and |ξ| as determined by the sign of the Q ′ term.The result of this is that the L 1 and L 2 points now lie at different distances from the planet.As we increase β, the L 1 point moves closer to the star and the L 2 point moves closer to the planet.The L 3 point also remains absent in this case.The character of the equipotentials also changes with finite β.Interior to L 2 , the equipotentials remain centered on the planet, and so the zero-velocity contours restrict particles to circumplanetary motion.However, between the contour passing through L 2 and that passing through L 1 , the zero-velocity curves open up and allow particles to escape to circumstellar orbits (this is the essential element that allows narrow dust rings to form in the model of Hayakawa & Hansen (2023)).As a result, the range of circumplanetary orbits is not bounded by L 1 and L 2 but by L 2 and a location L ′ 2 , representing the innermost edge of the pseudopotential contour that passes through L 2 .
The nature of the equilibria will depend on the second derivatives of the pseudopotential.At low values of β, both the L 1 and L 2 points remain saddle points, as in the purely gravitational case.However, as β increases, the L 1 point becomes a potential minimum when ξ 1 = −β −1/3 (this is when ∂ 2 C H /∂η 2 changes sign.) Figure 1 shows how the locations of L 1 , L 2 and L ′ 2 vary with the value of β (for the case µ = 0.001).

Triangular Points in the Hill Problem
The presence of the β dependant term in equation ( 2) enables the existence of additional equilibria beyond the colinear points.The new solution condition is ∆ = β −1/3 , which amounts to a condition on a specified distance from the planet.If we put this into equation ( 1) we must satisfy This yields something very like the L 4 and L 5 points, which don't appear in the usual Hill or photogravitational Hill problems.This is a consequence of the β dependant terms that remain even after the Taylor expansion of the potential about the location of the secondary.These terms break the usual symmetry that exists in the other versions of the Hill problem.As such, these points are not really a qualitatively new equilibrium, but they will have an influence on the kinds of periodic orbits we seek to describe below.Figure 2 shows an example for β = 0.1 and Q = 1 (so µ = 10 −3 ) compared to the β = 0 case for the same mass (in red).The red contour shows the seperatrix for the β = 0 case, delineating the boundary between circumplanetary orbits, interior and exterior heliocentric orbits, and the forbidden region.The inclusion of finite radiation pressure introduces a new set of solutions, defined by the condition ∆ = β −1/3 and shown by the dotted circle.This leads to the solution defined by equation ( 6), resulting in the new L 4 and L 5 equilibria as shown by the new extrema in the solid black contours.We see that the circumplanetary region is now compressed and there is a region of mostly-exterior heliocentric orbits that now pass interior to the planet briefly, between the L 1 point and the circumplanetary region.
We dub these points the Triangular points because of their obvious connection to the L 4 and L 5 equilibria in the traditional restricted three body problem, but we must also note the features that do not exhibit an exact correspondence.In particular the angle these equilibria make with the primary-secondary axis is not 60 • , and depends on β.As β increases, these points move closer to the axis, and will converge when Since Q ′ = β/µ 1/3 , we can expand in the small parameter µ to derive an approximate criterion For µ = 10 −3 , this criterion yields β con = 0.307 (a direct numerical solution yields the root β con = 0.306).We also note that the criterion for L 4 and L 5 to merge on the axis also corresponds to the condition for the change of L 1 from a saddle point to a potential minimum.Note that, with our sign convention, this does not correspond to a fixed point for stable orbits.As we shall see, L 1 remains a locus for unstable equilibria.

Equilibria out of the plane
In the case of very strong radiation pressure, a new set of equilibria is possible that simply do not appear in the purely gravitational problem (Schuerman 1980;Simmons et al. 1985).This occurs if the radiation pressure from one of the objects is sufficient to completely overwhelm the gravitational influence of that body.In this case it is possible to achieve equilibria that lie out of the orbital plane.
In the context of our problem, this will only occur if β > 1, in which case the right hand side of equation 3 can be set to zero if ∆ = (β − 1) −1/3 , even if ζ ̸ = 0.In the usual context of dust in orbit around a star, β > 0.5 leads to unbound trajectories, but the presence of the planetary gravity may open up the possibility of equilibria featuring very small dust that feels strong radiation pressure.
Fig. 3.-The contours show equipotentials in the η = 0 plane for the case µ = 0.001 and β = 1.0008.We see that there is a narrow libration region about this equilibrium point, potentially available to very small dust whose radiation pressure can outweight the gravity of the host star.
Further conditions on these equilibria are that η = 0 (from equation 2) -so that these equilibria lie in the plane passing through the planet and the star -and that We also have that ∆ 2 = ξ 2 + ζ 2 in this case, so that the requirement that ζ 2 > 0 places some restrictions on the solution, namely that Since Q ′ = β/µ 1/3 , this places a restriction on the mass, namely Thus, for planetary mass ratios, such equilibria are possible, but only for particles whose β is infinitesimaly above unity.As an example, Figure 3 shows the contours of the potential about the equilibrium point at ξ = −10.032and ζ = 3.924 for the case of µ = 0.001 and β = 1.0008.Note that the range on the horizontal axis is quite small -the libration region is very narrow in ξ, although it can extend over several scale heights in ζ, as demonstrated by the large range on the vertical axis.Thus, these equilibria are narrowly confined to a essentially a line-segment above the plane, on the line between the planet and the star Although these equilibria exist in the Hill expansion version, they occur only close to the primary (as the radiating body) and, as such, are not relevant to the question of dust structures in the vicinity of the secondary.
Our principal goal is to understand the nature of longlived dust particles orbiting a planet, such that they might be imaged in scattered light.In order to better understand this, we wish to understand the orbits that are stable around the planet, subject to the combined effects of gravity and radiation pressure.Henon (1969Henon ( , 1970) ) presented a classic analysis of the different families of periodic orbits in the Hill problem.We wish to understand how the effects of radiation pressure change this.Before we do that, it is worth briefly reviewing the simply periodic orbit families as classified by Henon.
Families a and c represent the (unstable) extensions of the libration orbits about the Lagrange equilibria at L 1 and L 2 .The families g and f represent the prograde and retrograde satellite orbits of the planetary body.The prograde orbit family splits, at critical point g1, with the introduction of a second orbital family g ′ (which contains both a stable and unstable equilibrium for some ranges of C H ). Family g also becomes doubly periodic at a critical point g2.

Prograde Orbits
A detailed survey of the simple periodic orbits requires direct numerical integration of the equations of motion ( 1),( 2) and ( 3).Here we will restrict ourselves to the orbital plane, as the out-of-plane equilibria are too far from the planet to be realistically included in the Hill limit.We adopt a similar orbital classification procedure as discussed in Henon (1969).We begin integrations with η = 0 and ξ = 0 and choose an initial value for ξ, and a value for η based on an assumed value of C H .The positive sign of η is chosen and we record the ξ and ξ every time the orbit crosses the η = 0 plane moving in the positive direction.This is used to construct a Poincare plot.Equilibria are defined as those orbits for which ξ does not change from the initial value on the first complete orbit (i.e.simple periodic orbits in the definition of Henon).The stability of each equilibrium is evaluated by constructing the linear mapping of orbits surrounding the equilibrium and evaluating the eigenvalues of the resulting matrix (Tremaine 2023).We map out the families of simply periodic orbits and will also describe a subset of the doubly periodic orbits which can impact the stability of some of the simply periodic families.
Figure 4 shows the resulting prograde equilibria, for the mass ratio µ = 10 −3 .In the upper panel, the red curves represent the original orbital families for β = 0 from Henon (1969).The black curves show the case for weak radiation pressure (β = 0.001, so Q = 0.02).The red curves show the original Henon orbit families of g and g ′ (and the unstable family a).In the case of β = 0, the appearance of family g ′ occurs at C H = 4.4999 and ξ 0 = 0.2835 and all three equilibria move away from this point as C H decreases.In the case of non-zero β we find that, instead of a point of intersection, the new family G ′ emerges at lower values of ξ than the equilibrium on the G family at the corresponding C H .The evolution, in this case, is more akin to an avoided crossing of two families, G and G ′ , each of which contains a portion of the original g and g ′ .A similar behaviour occurs in the prior definition of the photogravitational Hill problem (Markellos et al. 2000;Kanavos et al. 2002), albeit with different G and G ′ (corresponding to linking different combinations of the sections of g and g ′ ).In our case, the family G extends from high C H all the way to the Lagrange point, with the equilibria distorting from quasi-circular shapes (deep in the potential well) to more elongated shapes as they approach the Lagrange point.The family G ′ , on the other hand, only exists at moderate C H , linking a stable branch of compact eccentric orbits to the unstable branch that was originally the extension of the g family.
The blue curves in the upper panel of Figure 4 show the equilibria for the case of β = 0.01.The G family shifts to larger C H at fixed ξ, while the G ′ family shifts to the left (lower C H at fixed ξ).Note also that the A family (the analog of the a family) shifts down (to lower ξ at fixed C H ), so that the A and G ′ families start to approach each other as β increases.The lower panel of Figure 4 shows the evolution of the orbital families as β increases.The magenta curves show the case of β = 0.05.As the radiation pressure grows, the forbidden region moves to large C H and G family shifts in the same direction.Conversely, the G ′ family shifts down to lower C H .The black curves show the case for β = 0.1.In this case, the family G ′ no longer exists.
In the case of β = 0, the association of the two branches of the g ′ family are natural, as their orbit shapes are the same except for a simple rotation of 180 • about the η axis.Once β > 0, the symmetry of the potential is broken and the equivalent solutions are no longer as similar in shape.Indeed, it is now the two solutions of the G ′ branch that retain similar shape (initially, since one is unstable).Nevertheless, the existence of two stable equi- In the upper panels we show the Poincare plots.Each shows two clear stable equilibria and one saddle point between them.However, the β = 0 case is much more symmetric.The differences can also be seen in the lower panels, which show the orbits corresponding to the equilibria (red is unstable in the long term).On the left, the two stable equilibria are close to reflectionsymmetric in ξ, but the plot on the right shows that this symmetry has been broken by the presence of radiation pressure.It is now the two G ′ solutions that share a similar shape -because they emerge together as C H drops below the critical value.
libria and one unstable equilibrium remains even with finite β, as demonstrated by Hamilton & Krivov (1996).
Figure 5 demonstrates how the orbital shapes are affected by non-zero β.The two upper panels show Poincare plots of orbits crossing the η = 0 plane with η > 0, for C H = 4.4 in both the β = 0 and β = 10 −3 cases.We see that the stable orbits are centered around the two equilibria, tightly confined within the separatrix that passes through the unstable saddle point.The two lower panels show the shapes of thie corresponding three equilibria in the ξ-η plane, with the unstable orbit shown in red.The similarity of the two G ′ orbits to each other, for finite β, is clear.
The G family solution continues to evolve to ever more eccentric shapes as C H decreases.As C H decreases, the shape of the inner stable solution on the G ′ branch also evolves to a more eccentric shape and approaches the same dichotomy with the G solution as in the β = 0 case.Figure 6 shows the case for C H = 4.3.We see that the two equilibria are now islands of stability with an intervening unstable region, and that the remaining stable equilibria are now more similar (and more eccentric).
As β increases, the gap in C H between the G and G ′ families increases, and the G ′ family eventually becomes completely unstable for β > 0.01.For β > 0.095, the G ′ family of equilibria disappears.

Retrograde Orbits
A striking feature of Henon's original analysis was that the family f of retrograde simply periodic orbits remained stable to arbitrarily large distances from the planet, merging into the class of orbits referred to as quasi-satellites.Figure 7 shows how this family adjusts to the strength of the radiation pressure.
For small values of β the nature of the retrograde family (called F here) remains similar to the β = 0 case.However, for β = 0.1, the curve shifts to larger separations more quickly, and begins to exhibit qualitatively new features for β = 0.2.In the β = 0.2 case, the F family does not extend all the way to arbitrarily small C H but loops back to larger C H at larger |ξ|, forming an unstable branch.Figure 17 in Appendix A demonstrates this evolution.
The middle panel of Figure 7 also shows a second solution branch, called F ′ , which does not appear in the β=0 case.Appendix B demonstrates that this can be understood to be the result of the fact that the asymptotic solution at large distances from the planet now splits into two branches for β > 0. The solution branch F ′ corresponds to the second branch of this asymptotic solution.At large distances, this forms a second stable branch, shown as the blue curves in Figure 17 in Appendix A. The co-existence of these two stable branches, F and F ′ is shown in Figure 8 for the value C H = 4 and β = 0.1.
Another striking feature of these equilibria is that there appears a 'vortex-like' feature around C H ∼ 1 and ξ ∼ −1.9, where several different equilibrium curves converge.As discussed in Appendix A, and Figure 20, this appears to be the result of the interaction of the orbits with the triangular equilibrium points.
We also show the evolution of the G orbits (equivalent to the g ′ from Henon (1969)).These result from the part of the G family that becomes doubly periodic, where the condition that the path cross the η = 0 axis in the positive direction is satisfied for both positive and The middle panel shows the orbital families for the case β = 0.1.The family F is the equivalent of f and is also stable except for a brief crossing with the unstable G ′′ family.We also find the family F ′ , associated with the second asymptotic solution, and which is also stable at large enough distances.The shapes of the unstable families C,G and G ′′ are affected by the presence of the triangular equilibria, as discussed in the appendix.The bottom panel shows the case for β = 0.2.We see here that the family F no longer continues to arbitary distances.In all panels, the shaded region represents the forbidden region for that particular β. negative ξ.Thus, these curves form a pair with curves in the prograde case (Figure 4).In addition, we find another family of orbits, dubbed G ′′ , which is also doubly periodic, but both crossings of the η = 0 axis occur at ξ < 0. These orbital families are outlined in Figure 14 and Figure 18 respectively.Finally, we note that the F and F ′ families are the only ones that extend to large distances for β = 0.1, namely that C and the G, G ′′ families have only a finite extent.This is consistent with the absence of any corresponding asymptotic solutions ( § B.2).

Available Stable orbits
Our goal here is to determine the parameter space available for stable orbits under the influence of both planetary gravity and stellar radiation pressure.The interest in the simply periodic orbits is that stable regions are surround these orbits, so that they form the 'scaffolding' around which the structure of long-term stable orbits is built.Figures 5, 6 and 8 all demonstrate that the long-term stable orbits represent families that surround the different equilibria.
To properly explore this parameter space, we integrate the equations of motion for a range of initial conditions specified by a choice of ξ 0 and C H , with η = 0 and ξ = 0, for a time=200 (where the period of the planetary orbit is 2π). Figure 9 shows the outcomes of this parameter scan for µ = 0.001 and three choices of β = 0, β = 0.1 and β = 0.2.In this figure, a point is plotted at the corresponding value of ξ 0 and C H if the trajectory  remains within ∆ = 10 of the planet at the end of the simulation.
The β = 0 case closely resembles Figure 12 of Henon (1970), demonstrating the existence of stable orbits, both prograde and retrograde, for C H > 4.27, and an extension of the stable retrograde family to arbitrarily large distances (this family is characterised as the quasisatellites).The middle panel of Figure 9 shows the case of moderate radiation pressure (β = 0.1).We see that the prograde and retrograde stable regions at larger C H are now separated by a region where the orbits are unstable.This is mostly the consequence of orbits getting excited to essentially radial orbits and hitting the planet (Burns et al. 1979;Hamilton & Krivov 1996;Zotos 2015).The retrograde family still extends to larger radii and actually broadens, along with the presence of a second family of stable orbits at larger radii and positive C H .This represents the second family of asymptotic solutions, F ′ , discussed in Appendix B. Much of the other structure observed in Figure 7 does not appear in these figures because those additional equilibria are unstable.However, the F ′ family is robust and also appears in the full restricted three-body problem with radiation (as can be seen from Figure 10 of Zotos (2015), for example -albeit for the µ = 0.5 case).
Finally, the lower panel shows the case of strong radiation pressure β = 0.2.Inside the Hill sphere the prograde and retrograde stable regions are again separated by a large region of mixed stable and unstable orbits.The retrograde family also does not now extend to arbitrarily large distances, as expected based on our asymptotic solutions.

Mass Dependance
Our discussion thus far has focussed on the case of a mass ratio µ = 10 −3 , appropriate for a Jupiter-mass planet orbiting a Solar mass star.However, as we noted in our original derivation, our desire to retain the effects of non-infinitesimal β means that our equations are not scale free as in the case of the original Hill problem.Thus, we must investigate the effect of mass.
Fortunately, most of the structure in the solutions is determined by the β parameter itself, and the role of µ is primarily to determine a shift along the ξ axis, through the effect of the Q ′ parameter.This is demonstrated in Figure 10, which shows the orbit families for the case of β = 0.1 and three different mass ratios (µ = 10 −4 , µ = 10 −3 and µ = 10 −2 ).We see that most of the structure is maintained, but that the lower mass families are shifted to lower C H .If we examine the asymptotic solutions, we note that the frequencies ω depend only on β, not Q ′ , and so the existence of the asymptotes should be independent of µ.
At the higher mass ratio, there are some differences in the details, as the G ′ prograde family returns, even for β = 0.1.

THE FULL RESTRICTED THREE-BODY PROBLEM AND OBSERVATIONAL SIGNATURES
The direct observation of a self-luminous exoplanet requires the detection of an unresolved point source, usually at the infrared wavelengths, where self-powered emission is expected to peak.However, if there is a substantial population of circumplanetary dust, the scattering of stellar light by the dust may provide an observable signature at shorter wavelengths more characteristic of Fig. 10.-All three panels show the orbital families for the case β = 0.1.In the upper panel, the secondary is only µ = 10 −4 of the total system mass.In the middle panel, µ = 10 −3 (so this is equivalent to the middle panel of Figure 7).In the lower panel the mass ratio is µ = 0.01.Once again, stable equilibria are maked by solid lines and unstable equilibria are marked by dotted lines.
the stellar emission.Similarly, with a large enough surface area, the thermal emission from the dust may overwhelm that of the planet, although this is a function of wavelength (Kennedy & Wyatt 2011).Furthermore, such populations are more extended and thus potentially resolvable.Therefore, we wish to now place our results of the previous section within the context of the full restricted three body problem, and to identify potentially observable signatures.We will focus primarily here on the stable orbits, which could be populated by long-lived dust populations, and on those families within or near the planetary Hill sphere.There will also be long-lived families associated with the L 4 and L 5 points in the nonzero β case, but we will not treat them here.

The shape of a trapped dust population
To model the potential signature of a circumplanetary dust population, we repeat the orbital integrations of the prior sections, but now within the context of the full restricted three-body problem, so that X = 1 − µ + µ 1/3 ξ and Y = µ 1/3 η. Figure 11 shows three examples of the resulting orbits, chosen to represent the different classes of bound orbit, in the case µ = 10 −3 and β = 0.1.In each case, the choice of initial parameters was chosen from the middle of the corresponding stability region in Figure 9, albeit integrated now in the full restricted three body potential.So, these do not correspond exactly to the relevant simply periodic orbits, but are representative of the broader family of stable orbits that precess about the periodic orbit.The three examples shown are for one of the widest stable prograde orbits (black curve), a retrograde orbit of similar width (red curve) and a quasisatellite retrograde orbit (green curve).
To properly represent the observations, we must go beyond the two dimensional orbits.Thus we now evolve the full three dimensional set of equations, and choose our initial particle trajectories with positions initially drawn from a singular isothermal sphere radius distribution, spherically distributed except that we exclude starting positions with inclinations between 60 • -130 • of the orbital plane.This is based on the assumption that the dust initial distribution should follow that of the parent bodies.The observed irregular satellites in the Solar system show this deficit (Jewitt & Haghighipour 2007), which is believed to be a consequence of the dynamical instability of such orbits (Hamilton & Burns 1991;Nesvorný et al. 2003).We choose initial velocities based on the circular velocity about the planet and random orientations.We integrate these forward for t = 200 in our time units (where t = 2π is one orbital period for the secondary about the primary).
Figure 12 shows contour plots of the resulting density of particles for the cases β = 0.1 (left panels) and β = 0.2 (right panels), in both the X-Y plane (at Z=0) and in the X-Z plane (at Y=0).The distribution of particles is primarily confined within the seperatrices that pass through L 2 and it is this surface which is responsible for the elongation of the observed distribution.The half-max and quarter-max contours are marked in red, demonstrating this elongation, in both the X-Y and X-Z planes.The principal feature is the elongation of the observable population along the axis between the secondary and the primary.This is to be expected given the ellipticity of the orbits shown in Figure 11, and also discussed in Hamilton & Krivov (1996).In particular, we see that the combination of prograde and retrograde orbits are expected to give a tapered appearance to the contours, with the flattening increasing away from the direction of the primary.
If we characterise the shape of the original parent population, by following the above model with β = 0, we find that the Full-Width Half Max (FWHM) of the surviving particles is 0.28 of the Hill diameter in the X direction, 0.27 in the Y direction, and 0.20 in the Z direction.These proportions increase to 0.39, 0.37 and 0.30 respectively, if we take the widths at quarter of the maximum (FWQM).Thus, the underlying population is spherically symmetric in projection, despite the Kozai-induced holes at the pole.
The results of the β = 0.1 integration show FWHM of 0.30 in the X direction, 0.24 in the Y direction and 0.20 in the Z direction (with 0.40, 0.31 and 0.29 at FWQM).Thus, the distribution of surviving dust is elongated with an aspect ratio ∼ 1.3 : 1.For β = 0.2 the elongation becomes even more extreme, with FWHM of 0.31 in X, 0.21 in Y and 0.17 in Z (with 0.43, 0.29 and 0.33 respectively at FWQM).Thus, the elongation of the dust population increases with increasing β.
This suggests a potential signature of the influence of radiation pressure -the observation of an extended distribution of dust yields a potentially resolvable target in thermal or scattered emission, and one that should become more elongated as one approaches shorter wavelengths that probe smaller particles and larger β.

Fomalhaut: A case study
Although the above model can be applied to any star/planet system imaged in scattered or thermal light, the star Fomalhaut and its attendant dust structures offers the most complete application to date.The narrow dust ring imaged by Kalas et al. (2005) has long stimulated discussion about the origin of this dust and how planets might sculpt it (Quillen 2006;Chiang et al. 2009), and the detection of a candidate planet Fomalhaut b (Kalas et al. 2008(Kalas et al. , 2013) ) has further amplified this discussion, including alternative interpretations (Currie et al. 2012).
The narrowness of the dust ring is another feature that needs explanation.Hayakawa & Hansen (2023) offer an alternative to the standard birth ring of planetesimals -possibly sculpted by additional planets (Boley et al. 2012), in which the dust is generated in a circumplanetary cloud of irregular satellites and then spirals out through the L 2 point, sherpherded into a thin ring by the restrictions imposed by the relevant Jacobi constant.The discussion of Hayakawa & Hansen (2023) focusses on the ring morphology, but the potential observability of the dust bound to the planet is a further interesting consideration.Kennedy & Wyatt (2011) discuss the application of this model to the observations of Fomalhaut b, but can be considered relevant to any putative irregular satellite cloud, either at the location of Fomalhaut b or elsewhere in the system.Kennedy & Wyatt (2011) discuss a dust cloud whose spatial distribution follows that of the parent bodies, but our calculations above show that the dust dynamics under the influence of radiation pressure result in an elongated structure, with the long axis pointing along the star-planet line.Interestingly, such a signature was discussed by Kalas et al. (2013) for Fomalhaut b, although it could not be discounted that this was a residual signature of speckle noise.
An additional signature of the Hayakawa & Hansen (2023) model is that the planet should lie slightly interior to the inner edge of the dust ring, by virtue of the geometry of the equipotential that passes through the L 2 point (through which the unbound dust escapes).This implies that the planet should lie at ∼ 0.92 the semimajor axis as measured by the shape of the ring.Once again, Fomalhaut b lies at about the correct distance from the ring, although the proper motion of the planet does not conform to this model any more than it does for the gravitational sculpting model of Quillen (2006); Chiang et al. (2009). .The fading and eventual disappearance of Fomalhaut b (Kalas et al. 2013;Gaspar & Rieke 2020) in the optical suggests that this object is not a planet but a dust cloud undergoing dispersal over a timescale of decades.The origin of the dust cloud has been postulated to be the result of a collision between planetesimals (Currie et al. 2012;Kalas et al. 2013;Gaspar & Rieke 2020) although the origin of the parent population is still somewhat in question (since it is not located in the nominal birth ring associated with the dust disk).The discovery of interior belts (Gáspár et al. 2023) is a possible source, although the path from the quasi-circular interior belts to a coherent, highly elliptical, dispersing orbit has yet to be described in detail.
The population of quasi-satellites associated with Henon's f family offers a potential parent population.If the planet associated with the dust disk is accompanied by a population of planetesimals in quasi-satellite orbits, then these objects would have heliocentric orbits with semi-major axes within ∼ 20% of the planetary orbit, and with eccentricities ∼ 0.1-0.3.However, if such an object were to release dust that is subject to radiation  2013).The red curve shows the evolution of the heliocentric parameters of quasisatellite of a Jupiter-mass planet orbiting at the distance of the primary Fomalhaut dust belt.The blue curve shows the orbital evolution if we start with the same initial conditions as the red curve, but assume radiation pressure with β = 0.1.The green curve corresponds to β = 0.2.Semi-major axes are normalised by the observed semi-major axis of the debris belt, with the planet offset relative to that according to the model of Hayakawa & Hansen (2023).
pressure, the orbit of the dust would deviate from that of the parent body, because the F and F ′ orbits in Figures 7  and 9 are not the same as for the β = 0 case.Indeed, the inferred heliocentric parameters for particles on these orbits can deviate substantially from those characteristic of the planet itself, even reaching apparently unbound values for part of the orbit (see appendix C for an example).
Figure 13 demonstrates how radiation pressure would cause the dust to deviate.Three orbits are shown.The red trajectory shows the variation in semi-major axis (scaled relative to the parent planet) and eccentricity of a single quasi-satellite trajectory in the case β = 0.If one takes the same initial position and velocity, but calculates the trajectory with β = 0.1 (blue curve) and β = 0.2 (green curve), one finds that the dust orbit makes several loops in semi-major axis and eccentricity, which pass through the region occupied by the measured orbital parameters of Fomalhaut b.For β = 0.2, the dust does pass through the observed region of parameter space, but is ejected pretty quickly.However, for β = 0.1, the orbit lies within the stable part associated with the F ′ family, and makes repeated passages through the observed region.These results suggest that the dust involved in the cloud should have β ∼ 0.1 to increase the chances of observation.For values close to β ∼ 0, the orbital parameters do not approach those observed, and for β ∼ 0.2, the dust leaves the system rapidly.For β ∼ 0.1, orbits of this type complete many orbits around the primary, spending ∼ 4% of their time within the inferred eccentricity range.This is consistent with the fact that the observa-tions from 2008-2020 span < 1% of the estimated orbital period at the inferred separation of Fom b, and can be described by a single set of osculating orbital parameters.
If Fom b is a dissipating dust cloud, the rate at which such objects are generated is very uncertain.It takes ∼ 700 years for the eccentricity to grow to the observed value from the original β = 0 orbit, so that any individual event has a low probability of detection within a 10 year window.The fact that one has been observed over the course of our brief observations suggest that the rate of dust cloud generation is non-negligible (or we have been very lucky).Although it may be challenging to justify a high rate of collision in the system, we should note that the Solar system contains populations of 'active asteroids' that need not experience collisions to eject dust (e.g.Jewitt 2012).If more such dust clouds are discovered in the future, their kinematics should be similar to that described in Figure 13 in this model, whereas they should exhibit a much wider range of parameters if they are being seeded by the interior belts observed in the system.

Fom b as an L6 point?
Much has been made of the curious orbit that Fom b appears to trace.Most of the attempted physical explanations assume an orbit close to coplanar with the disk.If we were to hypothesize that this feature was, instead, associated with the out-of-plane equilibria due to radiation pressure, could this very different geometry provide a better explanation for Fom b?The short answer is no.One can plot the locus of possible positions of the L 6 and L 7 points, assuming that the Fomalhaut ring indicates the orbital plane for any secondary.For different values of β, these points lie at different heights above the plane which, in turn, project to different locations on the sky plane.However, none of them come close to Fomalhaut b.
Somewhat more strikingly, the resulting locus does pass close to the location of the feature dubbed the 'Great Dust Cloud' in the JWST images of this system (Gáspár et al. 2023).However, it requires very small particles to achieve β > 1, and so the 23.5 µm bandpass seems like an unlikely wavelength to identify such features.Furthermore, ground-based imaging suggests this feature may be a background feature (Kennedy et al. 2023) and so this agreement appears to be a coincidence.

Searching for the True Planet
We have thus far focussed on the potential shape of a population of trapped dust.To estimate the brightness, we must account for the finite lifetime of dust trapped within the Hill sphere, subject to both erosion due to collisional evolution, and to the loss of material through the L 2 point to the thin disk.In Hayakawa & Hansen (2023) we estimated the average time for escape to be ∼ 10 7 years for grains with β ∼ 0.1-0.2 in orbit around a Jupiter mass planet at 140 AU.In similar fashion as Kennedy & Wyatt (2011), we can estimate the timescale in a population of irregular satellites to generate dust by collisions.In a collisional cascade, the mass and mass loss rate is dominated by the most massive bodies, denoted by radius s max , but the surface area for emission and scattering is dominated by the smaller bodies, of radius s.The rate of collisions amongst the larger bodies is assuming a total mass M irr , scaled here relative to a Lunar mass, confined to a volume ∼ 0.3 Hill radii.We truncate the collisional cascade at dust radii s ∼ 10µm, because β ∼ 0.18(10µm/s) for the Fomalhaut system and this is approximately the β value where dust escapes before being ground down further (Hayakawa & Hansen 2023).Thus, an irregular satellite cloud of this size would be in approximate steady state.We note that this age is younger than the overall age of the Fomalhaut system, but we expect irregular satellites to be captured during three-body interactions during epochs of dynamical instability in a planetary system, which can occur at ages of hundreds of millions of years (as has been hypothesized in our own Solar system).Assuming a Dohnanyi (1969) power law, extending down from 100 km to 10µm, we estimate a total cross-sectional area in dust ∼ 1.5×10 23 cm 2 ( similar to that estimated by Kennedy & Wyatt (2011) for this system with slightly different assumptions).
To assess whether this is observable in the recent JWST images (Gáspár et al. 2023), we estimate that the dust should be at a temperature ∼ 47K (a/140AU ) −1/2 .At 23µm, and at a distance of 7.7 pc, this implies an object with total observable flux ∼ 0.4µJy, or ∼ 1µJy per square arcsecond, if spread over an angular width of ∼ 0.4 ′′ (corresponding to 0.3 Hill radii).This is ∼ 150 times fainter than the 'Great Dust Cloud' and about ten times fainter than the 1σ noise level estimated by Gáspár et al. (2023).
The above estimate assumes a blind search, but one might also help to narrow the field if one associates the Fom b orbit with an associated quasi-satellite.However, quasi-satellites can exhibit large excursions relative to their guiding center planets, so the identification of Fom b with such a population does not narrow the location of any particular planet (other than locating it in the western half of the orbit).
We conclude that the nominal model is not yet easily detectable with the extant JWST data, but does lie within an order of magnitude of the current noise levels.Thus, either deeper observations or a more optimistic estimate for the size of the dust population could bring observations and theory into comparable brightness levels.

CONCLUSIONS
The calculations in this paper represent a follow-on from the calculations of Hayakawa & Hansen (2023).In particular, if we posit that geometrically thin dust rings are the product of the grinding down of an irregular satellite population trapped within the Hill sphere of an extrasolar planet, then we seek to identify observational signatures of a dust population that remains bound to the planet.
To that end, we have studied the orbital properties of dust affected by radiation pressure and the gravity of the secondary body, and classified the kinds of periodic solutions that may occur.We have extended our study to values of β that exceed the infinitesimal limit treated in many previous studies.We find that the orbital families of the β = 0 case are reproduced for small β but some of them disappear or are substantially modified at larger values of β.
The modified orbital properties at finite β suggest that a bound dust cloud should be elongated along the primary-secondary axis, with the level of elongation increasing with β.This suggests that resolved sources whose shapes change with the observed wavelength may be the signature of radiation pressure at work.
As an application of our model we consider the Fomalhaut system and its enigmatic planet/dust cloud Foma-lhaut b.We find that the curious orbital properties of Fom b could be explained if the parent body of the dust cloud was following a quasi-satellite trajectory relative to a yet-undiscovered planet orbiting in the narrow dust disk.The effect of radiation pressure on dust released from such bodies can produce orbital elements consistent with those observed.
Data availability: The data underlying this article will be shared on reasonable request to the corresponding author.
The author thanks the referees for constructive comments.This research was supported by NASA grant 80NSSC20K0266.This research has made use of NASA's Astrophysics Data System Bibliographic Services.
The primary family of stable prograde orbits in Henon's analysis is g, intersecting with a secondary family g ′ that emerges below a critical value of C H .The equivalent families G and G ′ for non-zero β do not map directly onto g and g ′ because they do not intersect.
Figure 14 shows the generalisation G of this family (g) for µ = 0.001 and β = 0.001 (black), β = 0.05 (red) and β = 0.1 (blue).Each panel is chosen to have the same ξ(0), which means that the C H will be different for different β.The progression from the upper left to bottom right follows the decrease of C H (for fixed β).The principal qualitative effect of increasing β is compression of the orbit shape in the η direction, i.e. the effect of radiation pressure is to elongate the orbit along the primary-secondary axis.
We have designated this family as G because we wish to use this for the main family of stable prograde orbits.However, the comparison with Figure 14 and the corresponding Figures 4, 5 and 6 of Henon (1969) shows that the orbits for lower C H represent the generalisation of g ′ rather than g.Note also that, in the lower two panels, the solutions are double periodic and include an intersection moving in the positive direction at ξ < 0. These are the origin of the G curves in Figure 7.
Figure 15 shows the orbits of the G ′ family, as a function of C H .This is composed of a combination of g and g ′ from the β = 0 case, although most of them correspond to the g family.The black curves show β = 0.001 and the red curves show the case for β = 0.05.This family shrinks in extent as β increases, and so there are no red curves in the upper right and lower left panels, because there are no solutions for these values of ξ(0).Furthermore, there are no blue curves because this family is absent for β = 0.1.We show green curves in those cases where there is no solution at β = 0.05.In these cases, C H is chosen to have the same value for β = 0.001 and β = 0.05.To complete the discussion of prograde orbital families, Figure 16 shows the evolution of family A (the orbits of libration about the L 2 point), for different C H and β = 0.05 and β = 0.1.Once again, the comparison in each panel is made with fixed ξ(0).We see that increasing β makes the orbit more compact.The general shape trend is similar to that in the β = 0 case, with the orbits tending towards the consecutive collision shape as C H becomes more negative.
The most important retrograde orbital family is F , the generalisation of Henon's f family.Figure 17 shows the evolution of this family as C H decreases, for β = 0.1 and β = 0.2.As seen previously in Figure 7, the extension of F to arbitrarily large distances remains for β < 0.11, but fails for larger values.This is why there is no red orbit in the bottom right panel of Figure 17.
Figure 17 also demonstrates the appearance of the second branch of the F family, which we have called F ′ .This is represented by the second asymptotic solution discussed in § B, and is shown by the blue curves.The orbital shape is similar but the extension is more pronounced in the η direction.
The G curves in Figure 7 are simply the second axis crossing of those shown in Figure 15, but the β > 0 cases also show a doubly periodic family G ′′ with both crossings at ξ < 0. The shape of this family is shown in Figure 18.This family is responsible for a lot of the structure in the middle and lower panels of Figure 7.The β = 0.1 case shows how the two loops converge as the solution approaches the turning point in C H , and also how a second branch emerges to track the F ′ family as the doubly periodic version.The β = 0.2 case does not extend over as large a range.Once again, we compare orbits with the same initial ξ(0) and choose C H accordingly.
Figure 19 rounds out the census of orbit families, showing the Family C of libration orbits about the L 1 point.As can be seen in Figure 7, for β = 0.1 the family does not extend to arbitrarily large distances, and so some of the panels show only the β = 0.05 case.
One of the more striking features of the β = 0.1 and β = 0.2 retrograde orbits is the 'vortex-like' structure apparent in the middle and lower panels of Figure 7.To better understand the nature of this feature, Figure 20 shows the full set of retrograde equilibria that pass through ξ(0) = −1.9, for the case of β = 0.1.There are a total of eight different curves shown between the two panels.Some of the orbits correspond to the standard incarnations to be expected (the F orbit and the leftmost G ′′ orbit) but several of the orbital families exhibit a doubling of the number of potential equilibria (this is responsible for the vortex-like structure in Figure 7).
We can trace the extra orbits to the appearance of the L 4 /L 5 analogues (shown as crosses in Figure 20) because we note that several of the new orbits retain the overall shape of the parent family, but feature an extra loop around the equilibrium point (compare the blue and cyan orbits in the right hand panel, or the red and magenta orbits).So, we can trace the extra features in the orbital distribution to the presence of these potential minima, which affect the orbital shape if the orbit passes too close to them.It is important to note that these particular solutions within the context of the Hill problem will likely experience significant modification when translated to the full restricted three-body problem, as the positions of the L 4 and L 5 equilibria will be affected by the inclusion of higher order terms truncated in the expansions used to define the Hill problem.We have included them here for completeness, and note that they have little role in determining the potentially observable features, as these orbits are unstable.

B: ASYMPTOTIC SOLUTIONS
We can gain some preliminary insight into the influence of radiation pressure by considering the asymptotic limits of the periodic orbits in the limit ∆ → ∞, so we are looking for simply periodic solutions to the equations We will adopt a trial solution of ξ = ξ 0 + K 1 cos ωt and η = K 2 sin ωt (assuming a judiciously chosen origin for time).Note that there is no term linear in t as in Henon's version, because of the βη term on the right hand side of equation (B2) and because the presence of the Q ′ term is sufficient to determine ξ 0 (which would require a contribution from the η term otherwise).
Inserting this trial solution into equations (B1) and (B2) implies a requirement on the frequency, namely which yields the solution For β = 0, the solutions are either ω 2 = 0 or ω 2 = 1 (Henon's solution).For finite β, we have real solutions for β < β crit = 0.11095.Larger values yield complex frequencies, suggesting that there are no stable asymptotic solutions for larger β. Figure 21 demonstrates the form of these solutions.
The quasi-satellites To complete the solution, ξ 0 = −Q ′ /(3 − 2β).The expression for the constant, in this limit, is This recovers C H = −K 2 1 in the β → 0 limit, demonstrating that this is the asymptotic form of Henon's family f .We note that the sign of the coefficient to K 1 in equation (B5) changes when ω 2 = 2β(3 − 2β)/(1 + β), which also happens to mark the maximum β at which ω is real.
In the purely gravitational case, Henon showed that the retrograde family f show stable equilibria in the asymptotic limit of large distances, and we now understand that these are associated with the class of orbits termed 'quasisatellites' (Mikkola & Innanen 1997;Wiegert et al. 2000;Giuppone et al. 2010).These are orbits whose heliocentric Fig. 17.-These panels correspond to the famiily F of retrograde orbits.In each panel, the black orbits refer to the case of β = 0.1 and the red orbits for β = 0.2.The crosses show the locations of the respective fixed points.The blue curves represent the F ′ family of the β = 0.1 case.The magenta curve represents the F ′ solution in the case β = 0.2.Once again, the overall trend with increasing β is to make the orbits more compact.
parameters share the same semi-major axis as the planet, but which exhibit finite eccentricities, and which circulate about the planet, on scales larger than the Hill sphere, when viewed from the planetary rest frame -as illustrated in Figure 11.Our results demonstrate that the existence of this class of orbits is not universal under the influence of radiation pressure.For strong enough β there is no longer such a stable asymptotic solution.
Furthermore, even for β < β crit , we find two solutions for ω 2 .One applies in the limit C H → −∞ (as occurs in the purely gravitational case), and the other in the limit C H → +∞.This suggests that finite β below the critical value should actually yield two equilibrium families at fixed K 1 , rather than the single one in the β = 0 limit.The solutions are given by Thus, both solutions correspond to ellipses, although the amplitude in the η direction will be different at the same K 1 (also corresponding to a different C H ).This second asymptotic solution does not occur in the Photogravitational Hill problem (Markellos et al. 2000;Kanavos et al. 2002;Perdiou et al. 2012), where only the shift due to the Q ′ is introduced.However, stable orbits corresponding to this equilibrium can be observed in numerical orbital integrations of the full restricted three body problem with radiation (e.g.Zotos 2015).As we shall see below, this solution may also have observational implications for the properties of the enigmatic object Fomalhaut b.

Orbits with Consecutive Collisions
Henon also discussed a class of solutions wherein the orbits exhibited excursions to arbitrarily large distances, but also returned periodically to ξ = 0, η = 0 to scatter off the secondary, terming these the "orbits with consecutive collisions".These proved to be the asymptotic forms of the families a, c, g and g ′ .However, an integral part of these asymptotic solutions were terms that provided the η component of the solution with either an additive constant or a term linear in time.As noted above, equation (B2) no longer admits such solutions, and so we find no equivalents to these asymptotic solutions in this case.

C: HELIOCENTRIC PARAMETERS
In order to determine whether a directly imaged feature is bound to the host star, one measures the proper motion (and radial velocity, if one is lucky) relative to the star and converts the results into a set of Keplerian orbital parameters.If the object is experiencing additional accelerations due to an unseen companion, the resulting parameters may take values different from those expected given the position relative to the star.
Figure 22 shows the calculated heliocentric Keplerian semi-major axis and eccentricity for two example quasi-satellite orbits -one from the F family (shown in red) and one from the F ′ (shown in black).These are both examples of stable quasi-satellite orbits for the case of µ = 0.001 and β = 0.1.We see that, during the wide excursions relative to the planet, the inferred heliocentric parameters can vary wildly, including reaching apparently unbound values in the case of the F ′ orbit.We also compare this to the values inferred for the dusty object Fomalhaut b Kalas et al. (2013); Gaspar & Rieke (2020).The inferred orbital parameters of this object were the result of much discussion in the past, but are quite consistent with dust in a quasi-satellite orbit.We discuss this more directly in § 4.2.1.This paper was built using the Open Journal of Astrophysics L A T E X template.The OJA is a journal which provides fast and easy peer review for new papers in the astro-ph section of the arXiv, making the reviewing process simpler for authors and referees alike.Learn more at http://astro.theoj.org.Gaspar & Riecke (2020) and the open point those from Kalas et al. (2013).Semi-major axis is measured relative to the planet (which is taken to be the dust ring median for the observations).We note again the black curve is a bound equilibrium simply periodic orbit -the unbound heliocentric parameters are the result of accelerations due to the planet and not modelled in the Keplerian fit.

Fig. 1 .
Fig. 1.-The dotted lines indicate the position of the saddle points L 1 and L 2 , for the stated value of β.The L 1 curve turns solid at the point that it transitions from a saddle point to a potential minimum.The dashed contour indicates the location L ′ 2 that marks the inner edge of the pseudopotential that passes through the L 2 point.The blue contours are for a mass ratio of µ = 10 −4 , the black contours are for µ = 10 −3 and the red contours are for µ = 10 −2 .

Fig. 2 .
Fig. 2.-The black contours show equipotentials for the case of µ = 10 −3 and β = 0.1 (so Q ′ = 1).The red contour shows the separatrix for the case of β = 0 (and the same mass).The dotted line represents the circle upon which the Triangular equilibria are expected, and we can see the appearance of these points in the contours.

Fig. 4 .
Fig. 4.-The dark shaded region indicates the region in which orbits are disallowed.In the upper panel this is shown for β = 0.The red curves show the prograde equilibrium orbit families for the case β = 0.The black curves show how these shift once a small radiation pressure (β = 0.001, Q ′ = 0.02) is included.The blue curve shows the case for β = 0.01 (Q ′ = 0.2).In the lower panel, the magenta curves show the case for β = 0.05 and the black curves show the case for β = 0.1.The shaded forbidden region shown here is for β = 0.1 and the magenta long dashed line shows the boundary of the forbidden region in the case β = 0.05.In both panels, the solid lines indicate those equilibria that are stable, and the dotted lines indicate the unstable portions of the G/G ′ families.The dashed lines indicate the unstable A family.

Fig
Fig.5.-The two panels on the left are for the case of β = 0 and the two panels on the right are for β = 10 −3 .All integrations are for the case C H = 4.4.In the upper panels we show the Poincare plots.Each shows two clear stable equilibria and one saddle point between them.However, the β = 0 case is much more symmetric.The differences can also be seen in the lower panels, which show the orbits corresponding to the equilibria (red is unstable in the long term).On the left, the two stable equilibria are close to reflectionsymmetric in ξ, but the plot on the right shows that this symmetry has been broken by the presence of radiation pressure.It is now the two G ′ solutions that share a similar shape -because they emerge together as C H drops below the critical value.

Fig
Fig.6.-Thetwo panels on the left are for the case of β = 0 and the two panels on the right are for β = 10 −3 .All integrations are for the case C H = 4.3.In the upper panels we show the Poincare plots.Each again shows two clear stable equilibria and one saddle point between them.The lower panels show the orbits corresponding to the equilibria in each case.The comparison between the two lower panels shows that the shapes of the G and G ′ orbits evolve towards similar shapes as C H drops, approaching a similar structure as the original β = 0 case.

Fig
Fig. 7.-The upper panel shows the retrograde orbital families in the original β = 0 case.The black solid curve shows the original stable family and the black dotted curves shows the unstable equilibrium family c and the negative intersection of the g ′ family.The middle panel shows the orbital families for the case β = 0.1.The family F is the equivalent of f and is also stable except for a brief crossing with the unstable G ′′ family.We also find the family F ′ , associated with the second asymptotic solution, and which is also stable at large enough distances.The shapes of the unstable families C,G and G ′′ are affected by the presence of the triangular equilibria, as discussed in the appendix.The bottom panel shows the case for β = 0.2.We see here that the family F no longer continues to arbitary distances.In all panels, the shaded region represents the forbidden region for that particular β.

Fig. 8 .
Fig.8.-The two panels show the Poincare plot for the retrograde orbits in the case µ = 10 −3 , β = 0.1 and C H = 4.The red curves delineate the edges of the forbidden region, which separates the F region of stability (on the right) from the F ′ region (on the left).In the left panel, the diffuse band of points near the red contour correspond to the quasi-circular orbits about the primary in the full restricted three-body problem.This is a diffuse band because these do not correspond to a fixed point family in the Hill problem.

Fig. 9 .
Fig.9.-In each panel the red region is the range of forbidden orbits.Black points indicate initial conditions that produce orbits which remain within ∆ = 10 after t = 200 and so are considered stable.The strength of the radiation pressure increases from the top (β = 0) to the botton (β = 0.2).We see that the radiation pressure initially increases the range of stable distance orbits (see middle panel) but eventually truncates the stable branch of retrograde orbits (lower panel).

Fig. 11
Fig. 11.-The black orbit represents one of the widest prograde orbits that is stable, in the case of µ = 10 −3 and β = 0.1.The red curve represents a stable retrograde orbit and the green curve represents a retrograde quasi-satellite orbit.The dotted curve represents the Hill sphere radius, centered on the secondary at (0.999,0).All orbits are shown in the reference frame co-rotating with the secondary.

Fig. 12
Fig. 12.-The contours represent surface density of points for particles orbiting about a secondary with mass ratio µ = 10 −3 and radiation pressure β = 0.1 (left panels) and β = 0.2 (right panels).The blue solid point represents the position of the secondary.The red solid contours represent the value corresponding to half the maximum and a value a quarter of the maximum.The dotted contours represents the equipotential that pass through the L 1 and L 2 points.It is the confinement within the L 2 contour that is responsible for the elongated shape of the distribution.The upper panel represents the X-Y plane, and the lower panel the X-Z plane.

Fig. 13
Fig. 13.-The solid point shows the estimated orbital parameters of Fom b fromGaspar & Rieke (2020), while the open point represents the earlier measurement ofKalas et al. (2013).The red curve shows the evolution of the heliocentric parameters of quasisatellite of a Jupiter-mass planet orbiting at the distance of the primary Fomalhaut dust belt.The blue curve shows the orbital evolution if we start with the same initial conditions as the red curve, but assume radiation pressure with β = 0.1.The green curve corresponds to β = 0.2.Semi-major axes are normalised by the observed semi-major axis of the debris belt, with the planet offset relative to that according to the model ofHayakawa & Hansen (2023).

Fig. 14 .
Fig. 14.-These panels show how the shape of the G family of orbits changes with C H .In the upper left panel, we show the equilibria with ξ(0) = 0.24 for β = 0.001(black), β = 0.05 (red) and β = 0.1 (blue).Also shown, as crosses, are the corresponding L 1 ,L 2 and L ′ 2 points for each β.The upper right panel shows the equivalent for orbits close to the widest point of the stable branch in Figure 4.The lower left panel shows the shapes of the orbits after the family has become doubly periodic, and the lower right panel shows the continued evolution towards lower C H .

Fig. 15 .
Fig.15.-These panels show the evolution of the G ′ family of orbits.In the upper left panel, we show orbits for β = 0.001 (black) and β = 0.05 (red).The linear equilibrium points are shown as crosses.The two solutions have the same starting value ξ(0).In the upper right and lower left, there are no solutions for β = 0.05 with equivalent ξ(0).In these cases, orbits with the same C H are shown as green curves.Finally, in the lower right, we show solutions for both β and the same ξ(0).

Fig. 16
Fig. 16.-In each panel, the black orbits refer to the case of β = 0.05 and the red orbits for β = 0.1.The orbits correspond to the family A of equilibria.The crosses show the locations of the respective fixed points.In particular, this family corresponds to the librations about the L 2 point.As C H decreases (from upper left towards lower right) the orbits acquire larger amplitudes and become more distorted.

Fig. 18
Fig.18.-In each panel, the black orbits refer to the case of β = 0.1 and the red orbits for β = 0.2.The crosses show the locations of the respective fixed points.The blue curves represent the F ′ family of the β = 0.1 case.Once again, the overall trend with increasing β is to make the orbits more compact.

Fig. 19
Fig. 19.-The equilibria shown here represent the family C of librations about the L 1 point.In each panel, the black orbits refer to the case of β = 0.05 and the red orbits for β = 0.1.The crosses show the locations of the respective fixed points.

Fig. 20
Fig. 20.-All orbits in this figure are equilibria for the case µ = 0.001, β = 0.1 and with starting condition ξ(0) = −1.9 and ξ = 0.The value of C H is different, however.In the left hand panel, the black curve represents the F family, and the red curve indicates the leftmost intersection of the G ′′ family.Finally, the blue and cyan curves represent two incarnations of the C family.The right-hand panel shows two incarnations each of the G ′′ family (blue and cyan) and the F ′ family (red and magenta).Also show, as crosses, are the stationary points of the potential.The equivalents of the L 4 and L 5 points are responsible for the doubling of the orbital families, as can be seen by the loops induced in the red and blue curves.

Fig. 22
Fig. 22.-The red curve shows the inferred heliocentric parameters for a quasi-satellite orbit drawn from the F family, in the case µ = 0.001 and β = 0.1.The initial conditions were C H = −5 and ξ(0) = −4.1455.The black curve shows the equivalent for a member of the F ′ orbital family, with C H = 5.85 and ξ(0) = −5.The solid point represents the Fom b parameters estimated by Gaspar & Riecke (2020) and the open point those fromKalas et al. (2013).Semi-major axis is measured relative to the planet (which is taken to be the dust ring median for the observations).We note again the black curve is a bound equilibrium simply periodic orbit -the unbound heliocentric parameters are the result of accelerations due to the planet and not modelled in the Keplerian fit.