Different perspectives on the b-plane: perturbation effects and use for resonant flyby design

Orbital resonances can be leveraged in the mission design phase to target planets at different energy levels. On the other side, precise models are needed to predict possible threatening returns of natural and artificial objects closely approaching a target planet. To this aim, we propose a semi-analytic extension of the b-plane resonance model to account for perturbing effects inside the planet’s sphere of influence. We compute the actual values of the perturbing coefficients by means of precise numerical simulations, whereas their expression stems from the properties of hyperbolic trajectories and asymptotic planetocentric velocity vectors. We apply the proposed b-plane model to design ballistic resonant flybys by solving a multilevel mixed-integer nonlinear optimization problem.


Introduction
The condition for two orbiting bodies to be in orbital resonance with each other is conceptually simple and is met when the two orbital periods are in a ratio that equals that of two integers. Such a physical phenomenon has been widely explored and exploited in the past two decades.
Orbital resonances can be used to further characterize the planetary protection analysis, which in general poses a constraint on the design of end of life and disposal maneuvers. The impact probability of human-crafted objects with specified planets or celestial bodies must remain below a certain threshold, under a certain confidence level up to some epoch forward in time. For instance, European Space Agency (ESA)'s planetary protection requirements (Kminek 2012) impose a planet-dependent maximum impact probability for the 100 upcoming years, accurate at 95% confidence level. Any of the mission disposal objects must be considered, as well as any perturbing effects on the two-body dynamics, uncertainties on maneuvers or initial conditions and possible failures of the propulsion systems must also be accounted for. An uncontrolled close approach of a disposal object with a given body is therefore a potential impact threat; studying orbital resonances allows the extension of the planetary protection analysis considering also future returns at the moment of first encounter. Colombo et al. (2016) made some work in this direction developing the SNAPPshot suite, a tool that performs the planetary protection compliance analysis according to the just mentioned requirements (Kminek 2012). SNAPPshot's original version estimates the impact probabilities performing a full Monte Carlo simulation, by sampling the state probability distribution that stems from either the uncertainty on the initial condition, the interplanetary injection maneuver or the failure of the propulsion system. SNAPPshot also represents possible resonant returns that may arise from the uncertainty sampling, using the b-plane formalism (Valsecchi et al. 2003). Romano (2020) improved the SNAPPshot suite by extending the pool of available numerical integration schemes and implementing a line sampling procedure, a Monte Carlo-based technique that improves the accuracy of the impact probability estimation and lowers the overall computational burden of the analysis.
Planetary defense is another direct application of orbital resonance theory, aiming to understand whether a near-Earth asteroid returns and threatens the planet again. Valsecchi et al. (2003) showed the b-plane potential for the resonance analysis, using the asteroid orbit characterization as principal application in their work. They also modeled the so-called keyholes (Chodas 1999), which are the counter-image of the Earth at the consequent resonant return mapped onto the b-plane of the current close encounter (Valsecchi et al. 2003). Keyholes identify the regions of the current b-plane that lead to future impacts, not simply to new crossings of the sphere of influence. Petit (2018) exploited the keyhole definition to design asteroid deflection maneuvers so that a consequent impact and resonance is avoided. The variety of works on the asteroid identification and characterization is extremely wide on its own, and it is beyond the scopes of this work to deepen the topic. Here, we focus on the b-plane model and its planetary protection and flyby design applications.
Orbital resonances are exploited in mission design applications, and repeated flybys over the same body can allow complex and composed trajectory deflections with low fuel consumption that could not be achieved with a single gravity assist maneuver. Most recent works explore the concept in the three-body problem and different moon systems.
Within the Earth-Moon system, Topputo et al. (2005) showed that a low-energy transfer to the L1 libration point is similar to a 5/2 resonance with the Moon. Topputo et al. (2008) studied a peculiar Earth-Moon resonance with a Poincare map approach to show its close link with weak capture into Moon orbits. Then, they also showed that such resonant orbits may lead to escape toward heliocentric trajectories, accounting for solar gravity within the four-body environment. Ceriotti and McInnes (2014) investigated the resonance exploitation for continuous polar observation missions, using the resonance with the Moon particularly for a ballistic flip of line of apsides. Vaquero and Howell (2014b) exploited the design of resonant arcs to target the libration points in the Earth-Moon system in the circular restricted three-body problem. Vaquero and Howell (2014a) also proposed a solution for targeting the libration point L5 exploiting dynamical systems theory and unstable invariant manifolds together with Moon-resonant orbits. Yárnoz et al. (2016) showed that the Sun gravitational perturbation can be used as free acceleration within multiple gravity assist trajectories in the Earth-Moon system, and applied the concept to design transfers to the Earth-Moon L2 point and toward the asteroid Phaethon. Oshima et al. (2017) searched for low-thrust transfers to the Moon in the circular restricted three-body problem, showing that repeated high-altitude lunar flybys reduce the total fuel consumption.
The other major set of works on orbital resonance exploitation regards the design of missions within Jupiter's and Saturn's moon systems. A proper trajectory choice is crucial, as major constraints on payload and fuel consumption are present. Johannesen and D'Amario (2000) designed the endgame of the Europa Orbiter mission within the Jovian system as a series of nearly resonant orbits with the Jovian moon Europa itself. Within the framework of the circular restricted three-body problem, Ross and Scheeres (2007) studied possible multiple gravity assist trajectories for quasi-ballistic captures and escapes in low-energy orbits. Campagnola and Russell (2010) introduced the Tisserand-Poincare graph for the study of ballistic endgames, showing the important role of resonant orbits patched with high-altitude flybys in the low fuel consumption and applying it to transfers between the Jupiter's moons Europa and Ganymede and between halo orbits from the Jupiter-Ganymede to the Jupiter-Europa system. Woolley and Scheeres (2011) showed that, albeit the most efficient strategy for a cheap capture would require infinite flybys, fuel savings up to 50% can be achieved performing repeated gravity assists for v-infinity leveraging, without requiring excessive times of flight. Lantoine et al. (2011) used the Tisserand-Poincare graph to target trajectories in the patched three-body problem to optimize the fuel consumption in a transfer between a close resonant orbit of Ganymede to a close resonant orbit of Europa. Campagnola and Kawakatsu (2012) studied 3-D resonant hopping strategies to connect resonant orbits over Jupiter's moons, for JAXA's Jupiter Magnetospheric Orbiter mission. Vaquero and Howell (2013) exploited dynamical systems theory, invariant manifolds and Poincare maps to design ballistic transfers between libration point and resonant trajectories in the Saturn-Titan system to target the moon Hyperion. Campagnola et al. (2014) introduced the Tisserand-leveraging transfers to design low energy and applied them to ESA's JUICE mission to Ganymede. They also used low-thrust models within Tisserand-leveraging transfers to explain the lunar resonances of the Moon orbiter SMART-1.
In this work, we use orbital resonances for trajectory design purposes in a different framework compared to what mentioned above, sticking to interplanetary resonant transfers and adopting the b-plane formalism, that allows the study of orbital resonances at the moment of first close encounter (Valsecchi et al. 2003). The theory is based on Opik's variables (Opik 1976) and continues the work on the geometry and the modeling of close encounters initiated by Carusi et al. (1990) and . The formalism has been further developed in the past years by Valsecchi (2006) and Valsecchi et al. (2015) providing better insight on the geometry of close encounters and quasi-collision conditions (Valsecchi 2006). Valsecchi et al. (2015) obtained a fully analytic solution for the post-encounter orbital parameters, for a given entry condition to the sphere of influence (Valsecchi et al. 2015), under the assumptions of two-body zero-radius sphere of influence theory and conservation of the Tisserand parameter T . A specified post-encounter semi-major axis is described by a circle in the b-plane (Valsecchi et al. 2003), tightly connected to the resonance definition because of its relation with the orbital period.
This work refines the b-plane model for a given post-encounter semi-major axis accounting for general perturbing effects. This is a direct consequence of better modeling the rotation of the planetocentric velocity vector experienced during gravity assists. The standard b-plane theory is summarized in Sect. 2. Section 3 presents a perturbed model for the b-plane circles, the loci of points of common post-encounter semi-major axis. We obtain the key physical quantities to extend the b-plane circle model to the general perturbed case with numerical simulations. The circular locus of point is extended to a belt-shaped one, also considering almost perfectly phased returns among the possibly threatening ones for planetary protection applications. We validate the model against a Monte Carlo simulation of Solar Orbiter's upper stage of launcher, visualizing which subset of the simulated cloud on the b-plane is in actual, i.e., simulated, resonance with Venus. We finally apply the proposed b-plane model to the design of interplanetary resonant orbits in Sect. 4. We study the case of Solar Orbiter (SolO) (EADS-Astrium 2011; European Space Agency (ESA) 2011) with its repeated flybys of Venus aimed at rising the orbit inclination. An optimization problem is formulated and solved, designing resonant interplanetary trajectories with ballistic gravity assist maneuvers. We obtain a solution already close to ESA's optimized resonant trajectory in a few seconds only.

Deflection and b-plane reference frame
By assuming that the planet follows a circular orbit around the Sun, the reference frame of analysis ( Fig. 1) was first introduced by Carusi et al. (1990). The system is centered on the planet's center of mass, and the (x, y, z) axes are directed as the heliocentric position, velocity v p and orbital angular momentum of the planet, respectively.
The remaining quantities stand for the planetocentric velocities U and U before and after the close approach, the spherical angles θ and θ with respect to the planet's velocity. The angles φ and φ locate the maximum circles identified by the y-axis passing through the U and U directions, respectively, measured from the maximum circle of the plane (y, z). γ identifies the deflection magnitude (i.e., the turn angle of the flyby) and ψ its direction, given as the internal angle between the sides γ and θ of the spherical triangle θγ θ , formed by the vectors U, U and v p (Carusi et al. 1990;Valsecchi et al. 2003). The flyby effect is modeled as an instantaneous rotation of the vector U without magnitude change.
First defined in Appendix of Greenberg et al. (1988) as (x , y , z ) in the framework of Öpik's theory (Opik 1976) and corresponding to the (ξ, η, ζ ) frame of Valsecchi and Manara (1997) and Valsecchi et al. (2003), the axes (ξ,η,ζ ) of the b-plane reference frame are defined as: The impact parameter b is defined in the b-plane as (Milani et al. 2002) (2) Fig. 1 Graphical representation of the reference frame of analysis. Picture from Carusi et al. (1990) Recalling Fig. 1, from an interplanetary point of view the flyby can be modeled as an instantaneous rotation of U into U . The resulting velocity deflection can be determined by solving the spherical triangle (θ, γ, θ ), with ψ acting as internal angle, obtaining the new orientation θ . The cosine law for spherical geometry gives (Valsecchi et al. 2003): cos θ = cos θ cos γ + sin θ sin γ cos ψ. (3)

Resonant circles
The b-plane coordinates and the impact parameter provide a straightforward definition of the deflection direction angle ψ (Carusi et al. 1990;Valsecchi et al. 2003): and introducing the quantity c = m/U 2 , where m is the planet's mass expressed in solar masses, γ can be identified as (Valsecchi et al. 2003): We introduce the resonance definition, i.e., that a flyby leads to a resonant return when the new small object's interplanetary orbit (with superscript ) satisfies the condition (Valsecchi et al. 2003) hT obj = kT pl with T to identify the orbital periods for the planet ( pl) and of the small object (obj) and (k, h) positive integers, respectively, the number of orbits of planet and small object until the next close approach. In the adopted unit system (Valsecchi et al. 2003), the planet's orbital period is T pl = 2π and the small object's one is linked to the new interplanetary orbit semimajor axis by T obj = 2πa 3/2 . Therefore, a resonant post-encounter semi-major axis can be determined just by k and h (Valsecchi et al. 2003). The subscript R is added to remark the resonance feature of this new orbit: The resonant condition is fully determined by the angle θ R (Valsecchi et al. 2003): This post-encounter angle must satisfy the deflection equation (Equation (3)). The resonant circle equation proposed in Valsecchi et al. (2003) is then: which is equivalent to

Perturbed semi-analytic extension
We semi-analytically extend the b-plane resonance model, considering also non-perfectly phased returns as resonant ones (Sect. 3.1) and accounting for perturbation effects inside the sphere of influence (Sect. 3.2).

Resonant belts
Recalling Equations (7) and (8), as in Colombo et al. (2016) we introduce a quasi-resonance definition, i.e., an object is considered in resonance with a given planet whether the condition for a generic resonance k/h is satisfied, with the quantity * to be an arbitrary value identifying the quasi-resonance threshold. We can be conceptually revert the definition, in order to obtain the values of k/h that correspond to the quasi-resonance boundaries: Re-applying Equation (8), we obtain two new quasi-resonant post-encounter semi-major axes: which in turn result in two new values for cos θ R (as in Eq. (9)) and the parameters of two new resonant circles (as in Eq. (12)).

Perturbed deflection
Some work to model perturbing effects in hyperbolic trajectories has already been made. For example, Anderson and Giampieri (1999) present an analytic solution for perturbing angles in a formalism close to Öpik's variables (Opik 1976), by means of the Born approximation (Fowler et al. 1927). A few of the geometrical considerations made in Anderson and Giampieri (1999) remain valid in this work, even though we assume no specific expression of the perturbing force and do not apply superposition of effects. Starting with the perturbing effects, inside the sphere of influence any perturbation acts by modifying the classical spherical triangle introduced by Carusi et al. (1990) and Valsecchi et al. (2003) for the angles: with the superscript * to identify the real quantities after the flyby and the clean symbols to denote the analytic ones, from the usual resonance theory (Valsecchi et al. 2003). Such angles are relative to the single planetocentric trajectory, i.e., in general different variations appear for different trajectories undergoing the same perturbations. Then, we model perturbing effects as the angular displacements given by ( γ , ψ, θ ), which modify the two-body asymptotic velocity deflection given by (γ, ψ, θ ).
Modeling the variation in the deflection γ * and θ * is straightforward: Re-arranging Equation (16) and with the definition of γ from Equation (6), γ becomes and θ , with the expression for θ as in Equation (9 with the symbol a * to denote the actual post-encounter semi-major axis, without the resonance subscript to preserve for now the generality of the definition. Following the definition in Equation (19), possible variations of the asymptotic velocity magnitude U can be included in θ . The computation of the perturbation ψ requires some observations on the encounter geometry. By definition, the b-plane is perpendicular to the vector U, which identifies theη direction (Equation (1)) and is common to both the simulated and the theoretical deflections. Particularly, we observe that: • The angle ψ is measured clockwise on the b-plane and counterclockwise on the threedimensional reference frame of Fig. 1 (Carusi et al. 1990;Valsecchi et al. 2003). • The vector U × U * identifies the principal rotation direction of U into U * and must lie on the b-plane, because of the cross-product properties.
We denote with ψ and ψ * the deflection direction angles, theoretical and actual, respectively, according to the definition by Carusi et al. (1990) and Valsecchi et al. (2003). ψ * is related but not equal to the angle identified by the vector U × U * and the +ζ direction. If we definẽ ψ * as the counterclockwise angle measured from the −ζ direction to the vector U × U * , we have that |ψ * | = |ψ * | + π/2.
Similarly, switching to a counterclockwise measure from the −ζ direction also for the angle ψ allows us to introduce the related angleψ. Its definition is similar to Equation (4) and also stems from the b-plane coordinates ξ and ζ : (21) Figure 2 shows a graphical representation of these handle anglesψ andψ * . These handle angles allow us to retrieve the value of ψ because they still refer to the rotation direction vector and the chosen b-plane point (ξ, ζ ). To obtainψ * , we first express the vector U × U * in the b-plane reference frame with 1 and then, observing the geometry of Fig. 2, we retrieveψ * using its cosine: Since both the anglesψ andψ * are measured counterclockwise on the b-plane, we add a change of sign to finally comply with the original b-plane description (Carusi et al. 1990;Valsecchi et al. 2003), obtaining:

Perturbed circle parameters
The perturbing angles ( γ , ψ, θ ) modify the circle parameters presented in Equation (12), since their definition has a direct impact on the variables that mathematically define the circles. Consequently, the resonant circle equations need to be re-derived. The spherical triangle equation highlighting the perturbations is then Since the perturbing effects must be small compared to the main two-body effect, the overall shape should be still identified by circles centered on theζ axis. Then, we find the two intersections of such circles with the b-planeζ axis, linked to the new circle parameters as: which correspond to the solution of the quadratic equation that is obtained by setting ξ = 0 and exploiting b 2 = ζ 2 + ξ 2 in Equation (25), after introducing the b-plane definition of the trigonometry functions of γ and ψ. The full derivation is reported in Appendix A.
The trigonometry function properties lead to the definition of the following perturbed circle parameters: If the angles γ , ψ and θ were all zero, the original solution of Equation (12) from Valsecchi et al. (2003) would be obtained.

Validation: SolO's upper stage of launcher planetary protection
A test case taken from a study conducted in Colombo et al. (2016) over Solar Orbiter's (SolO) upper stage of launcher is shown in Fig. 3, the same presented for the planetary protection analysis in Fig. 5. The original analysis and data reported in Colombo et al. (2016) were referred to a mission profile with launch on October 2018 (European Space Agency (ESA) 2011), while the actual launch eventually took place in February 2020. The nominal initial condition is (J2000 reference frame, centered on the Sun) 2 :  Fig. 3, where we held the entry velocity constant and equal to: We multiplied the same initial velocity by a factor 0.1 to generate the low-velocity case presented in Fig. 4.

B-plane deflection model
We have obtained the b-plane region presented in Figs. 3 and 4, simulating with Matlab ® a fine grid of points, all with the planetocentric velocity of the respective nominal conditions but with different initial positions. We generated the samples in b-plane polar coordinates and divided the whole domain in several regions; for the sake of conciseness, we present only the low-and high-altitude results, for a limited b-plane region. We compute the deflection as the rotation of the planetocentric velocity due to the angles γ and ψ, and we then compare it to the results from numerical simulations in the relativistic N-body environment. In the corrected b-plane model, the angles γ and ψ introduce the perturbing effects, computed from the simulation of the nominal sample. We then apply the same numerical values of γ and ψ to the whole cloud of simulated b-plane points. The perturbation on θ affects the circles for a specified post-encounter semi-major axis only; thus, we show it directly on the planetary protection analysis.
In the SolO case, the implemented correction ( Fig. 3b and d) improves the deflection model, obtaining exact deflections at the point where the perturbations are computed and highly accurate results nearby, although the error obtained with the standard deflection is already small (Fig. 3a and c). High-altitude flybys are less affected by errors, because of the lower time spent within the planet's sphere of influence and the lower nominal deflection, which introduces some error amplification itself.
The low-speed case shown in Figure 4 partly exposes the singularity of the b-plane formalism visible in the definition of cos θ (Equation (9)). The error becomes large at both low and high altitudes (Fig. 4a and c), and despite enormously improved by the perturbed model ( Fig. 4b and d), the reliability of the proposed approach becomes extremely localized. This is Fig. 3 Relative errors against full numerical simulation, SolO-like close approach due to the b-plane theory itself, not suited to handle low-energy cases where the three-body problem would better describe the dynamical regime. Still, the simulation-based computation of the perturbing angles keeps the deflection error small nearby the reference where such angles are computed.

Perturbed resonant belts
We obtained the results presented in this subsection using SNAPPshot to perform a planetary protection analysis for SolO's upper stage of launcher. We generated the cloud of points shown in Fig. 5a and b using the Monte Carlo analysis and consequent b-plane representation of the resonances already available in SNAPPshot. We computed the perturbing coefficients ( γ , ψ, θ ) only at the barycenter of the simulated cloud, expecting accurate results near it. This is mainly due to SNAPPshot's current features that allow the user to retrieve the full propagation data only of a single specified sample of the cloud (Colombo et al. 2016). We present a visual check in Fig. 5a and b, whereas the numerical values for the relative error are available in Table 2. We show directly the resonant belts, both for a better visualization A. Masat et al.  Table 2 show. Figure 5a and b also offers a visual feeling of the improvement reached: The yellow dotted loci of points, being the b-plane representation of numerically simulated resonances, are now almost perfectly predicted by the analytical belts, drawn in black. Furthermore, the neat separation from the non-resonant close approaches (dark gray dots) is also precisely predicted.

Ballistic resonant flyby design
In a planetary protection application, one wishes to identify orbital resonances in order to minimize the probability of a resonant return. Here, instead, we exploit the resonant circle concept in the opposite way, using the degree of freedom left in the patched conics  approximation 3 to force resonant flybys, trying to achieve a total ballistic deflection unfeasible for a single flyby. About the perturbed b-plane model proposed in this work, we are not aware of any application of this kind in the literature, particularly dealing with open interplanetary resonances. The model is yet to be completed to the most general extent, i.e., it still relies on the patched conics approximation for the definition of resonance. In other words, even though there we can already semi-analytically refine the deflection of the planetocentric velocity vector, constraining the corresponding modified resonant circles does not provide the same v of their two-body parents and therefore does not correspond to the v required for that specified resonance to happen. In any case, the innermost structure of the algorithm we present remains unchanged, since only the interfaces with the interplanetary environment will be adapted in future works.
The proposed algorithmic method features a modular structure, where at the highest level the resonant flybys are designed by solving a global optimization problem, under a fully two-body patched conics approximation. At a lower level, the perturbed resonant deflections are designed to minimize their difference against each correspondent optimal two-body v. Some work in a similar direction, even though not specifically meant for the optimal trajectory design purposes of this work, was already done in Valsecchi (2006) and Valsecchi et al. (2015), where the authors developed an analytic solution for the post-encounter Keplerian parameters implementing the b-plane deflection formalism. In Valsecchi (2006) and Valsecchi et al. (2015), a fully b-plane and Keplerian description was enough to achieve a fully analytic solution. In this work, we propose instead a hybrid Cartesian, b-plane and Keplerian approach,

Assumptions
First of all, a major limitation appears when using the b-plane formalism for design purposes: The deflections can be described in this framework only if fully natural. An artificial maneuver modifies the outgoing asymptote, the hyperbolic excess velocity magnitude and consequently the geometry described in Carusi et al. (1990). The design is limited to initial and final orbits featuring the same Tisserand parameter T . Recalling the relation in Carusi et al. (1990) and Valsecchi et al. (2003): The two-body patched conics definition of orbital resonance naturally embeds another assumption: Under a Cartesian description of the two orbits, the position remains fixed for all the resonant flybys designed, with only the heliocentric velocity v to change according to some v to have consequent resonant returns. This new velocity v + v becomes the incoming condition of the consequent flyby. In this way, at least part of the combinatorial nature of the multi-flyby design problem can be neglected; having the encounter position always fixed in time and space reduces the problem to identifying the velocities only.
Even when dealing with resonances some combinatorial issues exist, i.e., what actual resonance to choose at each flyby. Once a specified resonance is selected, there are again infinite possible solutions that satisfy such a resonance condition, which represent the degree of freedom exploited in the presented algorithm.
Apart from the gravitational constants, the reference frame and the ephemerides data, the algorithm requires: • A set of admissible resonances S = k/h| i . • An initial orbit (r 0 , v 0 ). • An ultimate target orbit (r f , v f ).
• The number of close approaches N to be designed to reach such a final trajectory. Figure 6 shows a block scheme diagram of the overall algorithm, described in detail in the following sections.

Highest level
A sub-v target is assigned to each flyby, and the planetocentric trajectory is selected in order to minimize the difference between its own deflection and the specified sub-target. The highest optimization level determines the best set of sub v targets that eventually minimizes the artificial contribution needed to match the final orbit. We can formulate the optimization problem as: Three remarks arise from the possible ways of computing | v residual |: • The summation of the sub-targets v (i) target must lead to the global v target , but the optimal solution could still feature some residual at each flyby. The resonance condition is a constraint on the v to ensure consequent return; matching the global deflection does not require anything specific on the single flybys.
• We should expect some residuals arising because of the resonance constraint; in spite of this, the summation of all the actual v-s given must still lead to the final target. • The last flyby may not be resonant and should match the final orbit, allowing a much more flexible design of the last deflection.
A possible design choice is then splitting evenly the residual v of the current flyby over all the next sub-targets. This allows the design of resonant flybys up to and including the penultimate one, without any artificial v, optimizing the last and more flexible deflection to reach the final target.
For N flybys, we can write the sub-target update as: The sub-targets are a sort of guidance for the design of each flyby. Here lies the drastic reduction in the combinatorial nature of the problem, replaced by the solution of a nonsmooth optimization leading to the optimal set of brute-force designed resonant deflections. Figure 7 shows a summary of the algorithm flow for the presented optimization level.

Resonant flyby optimization
Given a set S of resonances, for a specified sub v the algorithm finds the best resonant trajectory whose deflection gets the closest to the sub-target. We perform this with another two-level local optimization algorithm: 1. Select one of the resonances in S.  We still need to introduce some constraints to bound the trajectories in a feasible region (i.e., the impact parameter must be at the same time not too low and not too high, b min ≤ b ≤ b max ). Figure 8 summarizes the algorithm flow.

Unperturbed case
We model each deflection using the b-plane formalism. A distinction between resonant and free flybys is required for the first step: A specific resonance acts similarly to imposing a constraint on the position on the b-plane; the point must belong to a circle of center'sζ coordinate D and radius R. Identifying with α a counterclockwise angle measured from thê ξ direction (exactly as the angle in simple polar coordinates), given D and R a generic point is identified by: and the impact parameter is defined as whereas for a non-resonant case we obtain ξ and ζ by where b acts directly as the radius variable in b-plane polar coordinates. The non-resonant case removes the resonant circle constraint, introducing one more degree of freedom on b.
We compute the deflection angle γ as in Valsecchi et al. (2003): The rotation direction lies on the b-plane, and its orientation is strictly linked to the angle ψ, defined as: Particularly, the planetocentric velocity points toward the b-plane and the angle ψ measures a clockwise rotation from theζ direction to the b-plane point itself. In addition, we observe that the velocity vector is always rotated toward the center of the plane. Consequently, a possible deflection direction computation strategy rotates the direction −ζ counterclockwise of an angle ψ + π/2. Then, the rotation of the Cartesian incoming planetocentric velocity U i into the outgoing one U i at the ith flyby follows these steps: 1. Rotate −ζ counterclockwise of ψ + π/2. 2. Represent this new vector in the Cartesian reference frame, applying a vector rotation whose matrix arises from the b-plane axes definition (Equation (1)). 3. Include γ to build the principal rotation vector in the Cartesian reference frame. 4. Apply the rotation to U i .

Defining the residual of flyby i as
at the lowest level, the optimization process finds a solution where the quantity | v (i) residual | 2 is minimized. For the last free flyby, the algorithm selects the best coordinates (b, α), whereas for a resonant flyby the choice is limited to the coordinate α and the correspondent circle parameters.

Perturbed case
We rely on a fine b-plane mapping of γ , ψ and θ to account for perturbing effects. We report the interpolation and mapping strategies we adopted in this concept validation in Appendix B, as well as the way to compute the perturbed resonant circles from a map of perturbing angles.
Since the resonance definition is still strictly two-body patched conics, the current implementation of the algorithm targets these resonant orbits. We expect some residuals given the circle modification, i.e., because we compute the perturbed vs so that their difference with respect to the already optimized two-body deflections is minimized. We do not mean these vs to be maneuvers to be implemented, and our current purpose is only to show the difference between the perturbed and the unperturbed cases.
For the resonant case, we apply the same optimization process, constraints and deflection presented in the unperturbed application. The only difference lies on including the perturbing angles γ and ψ in the deflection computation. The perturbed last free flyby is again free of any constraint, and solution is free to comply with the last deflection required. Yet again the deflection model remains the same, but having both b and α as optimization variables and including the perturbing effects.

Validation:SolO-like resonant orbit design
We performed the optimization with Matlab ® and JPL ephemerides through the SPICE toolkit (Acton 1996), using the function patternsearch.m at the highest level (optimal v target search) and a multi-start fmincon.m for the unperturbed b-plane deflections, since the formulation is non-convex albeit smooth enough to allow the convergence of descent methods if starting from a good enough initial guess. In general, for the perturbed model the interpolation of the perturbing parameters makes the problem non-smooth even at the planetocentric design; thus, we implemented a patternsearch.m optimization at this level as well, with an fmincon.m refinement performed on the direct search optimal solution.
We ask the resonant trajectory design algorithm to search for an optimal solution similar to SolO's flybys of Venus (EADS-Astrium 2011), where resonant gravity assists raise the orbit's ecliptic inclination with low-cost maneuvers. We expect the solution to slightly differ from ESA's optimized one, since the algorithm for now designs ballistic flybys only. In any case, reaching a solution similar to SolO's mission profile would give a strong proof of the proposed design approach.
The reference case is the solution proposed with launch in January 2017 in European Space Agency (ESA) (2011), even though the actual launch has eventually happened in February 2020. We report the just mentioned optimized trajectory in Fig. 9. We use our algorithm to reproduce the resonant flybys only, i.e., to find a solution similar to the trajectory between the close approaches marked as V 2 and V 6.
We report the main features of the four resonant orbits in Table 3, with the first flyby to happen on May 22, 2020 (European Space Agency (ESA) 2011). The resonances are also available in European Space Agency (ESA) (2011), (3/4, 3/4, 2/3, 3/5) in that order. We  use the same notation introduced in European Space Agency (ESA) (2011) to identify the different close approaches, with V i standing for the ith flyby of Venus and V i − V j referring to the heliocentric orbit between the flybys i and j. The orbits are uniquely defined, since the position of Venus at each specified time adds the three remaining parameters.

Unperturbed (interplanetary) level
For the sake of testing the fully ballistic design algorithm, we remove some of the mission constraints, especially keeping perihelion and aphelion altitudes within specified boundaries. Therefore, we ask the solver to perform the specified change in inclination in the same number of flybys, letting all the parameters adjust autonomously. Given the optimal resonances mentioned above, the algorithm selects (k, h) values only within the range 1 to 5. We report the initial and final orbits for a Solar Orbiter-like resonant phase with Venus in Table 4.
The algorithm converged to a ballistic optimal solution, obtaining a v residual after the last flyby whose norm is on the order of 10 −8 km/s. The four deflections are We report the orbits designed in Fig. 10 and their parameters in Table 5. The intermediate orbits feature resonances (3/4, 2/3, 3/5) and the last orbit complies with the final inclination, albeit the last orbit reported in Table 3 is still in a 3/5 resonance. As a visual check, the algorithm has produced orbits that are similar to the optimized ones ( Fig. 9) available in European Space Agency (ESA) (2011). Also, the orbital parameters obtained are similar, apart from the perihelion distance which appears to be the main objective of control maneuvers, given also the purposes of the mission itself (EADS-Astrium 2011).
Only the first one out of three resonances corresponds to the ones designed in European Space Agency (ESA) (2011). This may be due to the absence of maneuvers and mission constraints in the optimization algorithm, which, for instance, leads to perihelion distances too low compared to the constraints specified in European Space Agency (ESA) (2011). Also the Cartesian formulation of the objective function that features an automatic v split might affect the obtained solution, since the optimization does not include user-specified intermediate orbital parameters.  Figure 11 shows the perturbed b-plane points that are chosen for achieving a deflection that gets the closest to the unperturbed interplanetary target. Visually, their distance from the optimal unperturbed b-plane coordinates seems to be minimized, and this may be due to the higher difference between the b-plane circle shapes than the actual v deflections, as it is discussed in Sect. 3.4.1. For the sake of conciseness, all the resonant flybys are reported in the same plot, although the b-plane reference frame itself has a different orientation for each close approach. All the resonant deflections feature some residual, shown in Equation (41), which is explained by the constraint of lying on their b-plane perturbed circles. The residual on the last free flyby is in the order of 10 −6 km/s:

Perturbed b-plane level
Providing a further and final check to the deflection model, we obtain the maximum error on the third flyby, equal to 0.628% with respect to the simulated trajectory.

Performances
In the unperturbed design, the heaviest part from the computational point of view is the target selection. patternsearch.m implements an intelligent direct search method that needs to perform several function evaluations at each algorithm iteration. For large domains and small tolerances, this results in a rather high computational time, due to the optimal targets not known at all. Instead, we obtain an extremely small computational time when the algorithm is executed again for retrieving the complete solution, of only 1.5 seconds on a single core of a CPU Intel ® Core tm i5-5200U running at 2.20 GHz.
About the perturbed design, in order to underline the performances of the deflection model itself we focus on the performance of the local computations, whose time mentioned below includes the numerical simulation of the selected trajectory. The heaviest part is the map generation, with higher computational time required for higher resolution. For a given map, we obtain again an extremely fast convergence, for about 2.7 seconds in the last free flyby and 1.9 seconds for the resonant close approaches, with the same CPU and machine settings of the unperturbed case. Considering the pair precision level reached and software used, a full transition to better programming languages could easily tackle these design problems with higher precision and even lower computational times.

Conclusion
The precise identification of the resonant loci of point in the b-plane reference frame successfully brought a semi-analytic modification of the current model. We obtained an exact solution, in terms of simulation compliance, at the b-plane point where the perturbing coefficients have been computed, achieving a highly reliable approximation in the nearby b-plane regions. The immediately next developments already comply with the planetary protection application needs, i.e., showing a full Monte Carlo simulation on the same b-plane. The model should eventually detach from the patched conics approximation: The proposed semi-analytic solution should be generalized to account for the variations of the perturbing parameters among different b-plane regions. To this extent, some focus should be put toward both exploring optimal mapping strategies and new analytic solutions, preserving the efficiency features and ensuring an in-b-plane precise impact/resonance probability computation.
Finally, we have shown the deflection modeling capabilities of the b-plane formalism by implementing an efficient ballistic resonant trajectory design algorithm. The algorithm could already be used for preliminary mission analysis purposes, given that a good enough guess for the resonant deflection is known. The main outcome is perhaps having shown how embedding and connecting different orbital representations could enhance the algorithm performances, using the b-plane features for the deflection design, the traditional Keplerian elements for the objective definition and the simple Cartesian representation to perform simple computations. The results achieved by the algorithm, in terms of precision and required computational resources, suggest keeping developing both the model and the method. On the one hand, the restriction to resonant flybys drastically reduced the number of degrees of freedom left to the problem, which is the main reason for such low computational costs. This also highlights that proper analytic or semi-analytic developments, such as optimal perturbation mapping strategies and in-b-plane maneuver design, could lead to application extensions to the non-resonant case. On the other hand, optimal trajectory design techniques already exist; thus, our method could contribute to increasing the functionality of the current procedures.

Appendix B Perturbation and mapping strategy for flyby design
We map perturbing effects onto the b-plane before designing the flyby, with a finer mesh close to the planet, and then interpolate to obtain the perturbing values not exactly at the map nodes, exploiting the deflection model to avoid numerical propagations within the optimization process.
We perform a set of simulations at the beginning of the design itself: We describe the starting conditions with the b-plane polar coordinates (b min ≤ b ≤ b max , 0 ≤ α < 2π), thus with their correspondent (ξ, η, ζ ) positions at the entrance of the sphere of influence, and the incoming planetocentric velocity of the current flyby. All these points together build a mesh (Fig. 12), whose propagation brings the following set of parameters at each node: • The perturbation γ on the deflection magnitude.
• The perturbation ψ on the deflection direction.
• The perturbation θ on the outgoing angle between the planet's and the object's planetocentric velocities. • Mesh-related handle quantities for the interpolation.
The different α angles are linearly spaced, whereas b features a quartic distribution between b min and b max . As shown in Fig. 12, we arbitrarily chose the quartic distribution on b to reach a finer mapping nearby the lower limit, where the difference between two values of b has an higher impact on the deflection error.
Despite the reason of mapping, the whole domain is fully evident in the free flyby case; we would not necessarily need to propagate the whole mesh for the resonant close approaches. If we knew the b-plane region of a specified resonance, we could restrict the domain to the interest areas only. Nevertheless, in the current application, not only is the position of the resonant circles at a specified close approach not known, but also all the configurations of the intermediate encounters are let free, leading to even more uncertainties than the sole circle location. For the current implementation, the optimal two-body result could cut some computational cost, identifying the interest areas in the N-body analysis. However, aiming to detach from the patched conics approximation, a complete mapping is anyway needed to include the perturbing effects, since we would run the optimization already at the perturbed level.
The interpolation strategy is similar to a bi-linear Lagrange polynomials, interpolating over b and the arc bα. The four interpolation coefficients are with i, j the indexes identifying the nodes of the cell enclosing b and α.
Being in a design phase and having mapped the perturbing effects on the b-plane, we do not need to pick one node only as the reference to draw the perturbed resonant circles, whose approximation looses accuracy far from the reference. Since in general the position of the circle intersection with theζ axis is not known, some iteration method is required. We propose a fixed-point-iteration-like algorithm: 1. Initialize the search by computing the unperturbed circle parameters D and R. 2. Compute the highest magnitude intersection b with theζ axis (accounting for both positive and negative intersections withζ ). 3. Get the interpolated perturbing angles γ , ψ and θ on this point. 4. Compute the perturbed circle parameters D p and R p and get the updated intersection b new with theζ axis. 5. If the difference between the current b new and the previous iteration's one (equal to b when initializing the algorithm with the unperturbed case) is lower than some tolerance stop the algorithm, otherwise return to step 1.
We expect a fast convergence, given the map resolution and the deflection model validation results. In fact, despite having set a rather strict tolerance of 10 −6 km for the difference between two iterations and a maximum of 1000 steps, this number was never reached and convergence happened in all the test cases analyzed.
This model does not take into account possible different values of the perturbation θ and still assumes the resonant loci of points to be described by circles. A possible generalization performs the optimization over both the variables (b, α) of the free flyby and relaxes the circle approximation of the perturbed resonances, but still makes use of the perturbing coefficients to insert the following nonlinear constraint based on the perturbed resonant geometric deflection: cos θ * R − cos(γ + γ ) cos θ * − sin(γ + γ ) sin θ * cos(ψ + ψ) = 0.
We are removing the assumption that such a deflection must bring a circular locus of points in the b-plane, accounting for possible different angles γ , ψ and θ in the various b-plane positions to be obtained via map interpolation. We attempted to optimize the deflection with this generalized resonance, although the new setup has not converged to an optimal solution. The results reported in the main body involve only the circular resonance model, and a further analysis to find a suitable method for the generalized resonance optimization is required.