Geometric analysis of a two-body problem with quick loss of mass

. We consider a two-body problem with quick loss of mass which was formulated by Verhulst in [14]. The corresponding dynamical system is singularly perturbed due to the presence of a small parameter in the governing equations which corresponds to the reciprocal of the initial rate of loss of mass, resulting in a boundary layer in the asymptotics. Here, we showcase a geometric approach which allows us to derive asymptotic expansions for the solutions of that problem via a combination of geometric singular perturbation theory [5] and the desingularisation technique known as “blow-up” [4]. In particular, we justify the unexpected dependence of those expansions on fractional powers of the singular perturbation parameter; moreover, we show that the occurrence of logarithmic (“switchback”) terms therein is due to a resonance phenomenon that arises in one of the coordinate charts after blow-up.

Here, the dot denotes differentiation with respect to time t P R`, while 0 ă ε ! 1 is a small parameter and n ą 1 is a real exponent. The function f is assumed to be continuously differentiable in pt, x, 9 x, εq with bounded derivatives for t P r0, Ls, where L is a positive, ε-independent constant; moreover, it may be assumed that 0 ă m ă f pt, x, 9 x, εq for all pt, x, 9 x, εq P R`ˆR 2ˆr 0, ε 0 s, with m and ε 0 constant.
Equation (1) arises in the study of modifications of the classical two-body problem from celestial mechanics, such as in the modelling of binary stars with rapid loss of mass. There, the two stars interact by Newtonian gravitational forces, whereby their motion is confined to a plane by conservation of momentum, with mass being ejected isotropically. The variables x and t are then associated with r´1 and θ, respectively, in that orbital plane, where pr, θq denotes standard polar coordinates. The variable u represents the total mass in the system; crucially, the large parameter ε´1 encodes the initial rate of loss of mass, while the exponent n in (1b) determines the rate at which mass decays in time. For a detailed derivation of Equation (1), the reader is referred to [13].
Equation (1) represents a classical singular perturbation problem that exhibits a boundary layer at u " 0 due to the presence of the small parameter ε in (1b). Correspondingly, in [14], Verhulst derives asymptotic solution expansions for (1) by applying matched asymptotics in the formulation of Vasil'eva [12], as well as the method of multiple scales that is due to O'Malley [9]. In the process, he observes that the case where n " 2 in (1b) is distinguished due to the presence of logarithmic ("switchback") terms in ε in the resulting expansions for px, 9 xqptq.
Here, we showcase how the classical asymptotic analysis of the boundary layer problem in Equation (1), as presented in [14], can be re-interpreted within a dynamical systems framework. Specifically, we combine the geometric singular perturbation theory due to Fenichel [5] with the desingularisation technique known as "blow-up" [4] to obtain rigorous asymptotic expansions for (1). In particular, we identify the reason for why the case of n " 2 in (1b) is distinguished, and we justify rigorously the structure of the corresponding expansions via a resonance phenomenon in one of the coordinate charts after blow-up. For definiteness, we focus on the case where f " 1 in Equation (1) which was also considered by Verhulst [14]. As will become apparent in the following, that choice of f allows for a relatively explicit analysis while being representative of the asymptotics of (1) under the more general assumption in Equation (3).
The article is structured as follows. In Section 2, we formulate the fast-slow system in standard form on which our analysis is based. In Section 3, we introduce the blow-up transformation which allows us to desingularise the flow of that system; in particular, we construct a singular orbit Γ, and we show the persistence thereof for ε sufficiently small. In Section 4, we derive asymptotic solution expansions for our system; we discuss the resonant case where n " 2, and we compare our results with those of Verhulst [14]. Finally, we conclude in Section 5 with a summary, and an outlook to potential future research.

Fast-Slow Analysis
In this section, we reformulate Equation (1) as a first-order system of ordinary differential equations that is amenable to our geometric approach: we introduce the new variable y " 9 x in Equation (1) to obtain 9 x " y, (4a) 9 y " u´x, (4b) ε 9 u "´f pt, x, y, εqu n , (4c) which will be the starting point for our analysis. Given our assumption in (3), Equation (4) has a unique steady state at the origin in px, y, uq-space. Henceforth, we will assume that for t P R`, the dynamics of (4) evolves on a compact set D :" r´X, Xsˆr´Y, Y sˆr´U, U s Ă R, with X, Y , and U positive and sufficiently large; in particular, p0, 0, 0q P D.
Equation (4) is a fast-slow system in standard form, written on the slow time scale t, where px, yq are slow variables and u is fast; the corresponding fast system on the fast time scale τ " t ε takes the form where the prime denotes differentiation with respect to τ .
Setting ε " 0 in Equations (4) and (5) yields two limiting systems, the "reduced problem" 9 x " y, (6a) 0 "´f pt, x, y, 0qu n (6c) and the "layer problem" The critical manifold for Equation (6), which we denote by S, is obtained from (6c): since f is bounded away from zero, by Equation (3), S is given by tu " 0u, i.e., by the px, yq-plane; linearisation of (7c) gives B Bu "´f p0, x, y, 0qu n ‰ˇˇˇt u"0u "´nf p0, x, y, 0qu n´1ˇt u"0u " 0 (8) due to n ą 1. Hence, the manifold S is not normally hyperbolic; in other words, the mass u of the system decays to zero algebraically in time, rather than exponentially. (We emphasise that S is still weakly attracting, as well as that the decrease in u is monotonic, by (7c).) It follows that standard geometric singular perturbation theory [5,6] is not applicable; to remedy the loss of normal hyperbolicity, we apply the blow-up technique in its formulation due to [4,7]. (For the reader's reference, the "regular" case where n " 1 in Equation (4) is discussed briefly in Appendix A.) Throughout the remainder of the article, we restrict to the case where f " 1 in (4): given that f is uniformly bounded away from zero in general, by (3), other choices of f can be considered in a similar fashion; details can be found in Section 5. Under that restriction, the general solution to Equations (6) and (7) then reads respectively. In particular, the solution of the reduced problem in (9) is usually known as the (leading-order) "outer" solution, while the solution to the layer problem in (10) is referred to as the "inner" solution. However, for ε positive and sufficiently small, neither solution perturbs to a uniformly valid solution for Equation (4); that deficiency is traditionally overcome by asymptotic matching in an "intermediate" region. Equation (4) is hence naturally studied in three regions I, II, and III; here, regions I and III correspond to the inner and outer regions, respectively, in the language of matched asymptotics, while region II represents the overlap between those two regions. We emphasise that, while the manifold S is non-hyperbolic, the existence of a slow manifold S ε " S follows from the invariance of S for any ε. (Correspondingly, Equation (9) then also gives the exact solution to (4) on S ε for arbitrary ε P r0, ε 0 s.) Therefore, blow-up does not yield the existence of S ε in our case; rather, it will allow us to desingularise the flow of Equation (4) in the vicinity of S ε , i.e., in regions II and III, and to justify rigorously the small-ε asymptotics derived in [13] .
In Figure 1, we visualise the fast-slow dynamics of Equation (4) for various values of n and ε and f " 1; the intermediate time scale which is observed is due to the non-hyperbolicity of S ε , and the resulting algebraic decay of u, for n ą 1. Figure 2 shows dynamics with different choices of f in (4); comparison with Figure 1 indicates that the dependence of f on x, y, and t, respectively, does not affect the qualitative dynamics of Equation (4).
Remark 1. In [14], Tikhonov's theorem [15,Theorem 8.1] is applied to show that, for ε Ñ 0, solutions of Equation (4) converge to those of the reduced problem, Equation (6), uniformly in time (t). In particular, the resulting leading-order outer approximation for the corresponding solution to (4) agrees with Equation (9) above. Here, we will argue that convergence to S ε under the fast flow of Equation (5) is, in fact, exponential in suitably chosen "blow-up" coordinates.
In the singular limit as ε Ñ 0, Equation (11) reduces to given the initial state Q´: px 0 , y 0 , u 0 , 0q, recall (2), we easily obtain the general solution to (12), as given in Section 2. Let us now introduce the following notation: we denote by Σ in the hyperplane tu " r 0 u in px, y, u, εq-space, where r 0 ą 0 is chosen sufficiently small, but fixed, i.e., we choose Σ in to define the boundary between regions I and II. In particular, since r 0 ă u 0 and as u is monotonically decreasing in τ , by (7c), we may determine the time τ in ą 0 such that upτ in q " r 0 : Hence, the segment Γ I of the singular orbit Γ in region I can be written as follows: Γ I :" px, y, u, 0qˇˇx " xpτ q, y " ypτ q, and u " upτ q for τ P r0, τ in s ( , where px, y, uqpτ q are as stated in Equation (10) and τ in is defined in (13). Moreover, we label the point of intersection of Γ I with Σ in as P in : px 0 , y 0 , r 0 , 0q.
Next, we study the persistence of the orbit Γ I , for ε positive and sufficiently small. Correspondingly, we denote by Qέ : px 0 , y 0 , u 0 , εq the initial point for the flow of Equation (11), with ε fixed; finally, we define the corresponding line of initial states in the extended px, y, u, εq-space as ´: " px 0 , y 0 , u 0 , εqˇˇε P r0, ε 0 s ( . By regular perturbation theory, Γ I will persist as an orbit Γ I ε that is initiated in Qέ : for n ‰ 2, we find xpτ, εq " x 0`ε y 0 τ`Opε 2 q, (15a) ypτ, εq " y 0´ε " x 0 τ´" pn´1qτ`u 1´n whereas for n " 2, we have which, together with (10c), implies that, for ε fixed, Γ I ε is given by Γ I ε :" px, y, u, εqˇˇx " xpτ, εq, y " ypτ, εq, and u " upτ q for τ P r0, τ in s and ε P p0, ε 0 s ( ; (17) here, xpτ, εq and ypτ, εq are defined in (15) and (16) for n ‰ 2 and n " 2, respectively, while upτ q is as in (10c). (We remark that we have stated explicitly the ε-dependence of x and y, for clarity.) Remark 2. We note that the union Ť εPr0,ε 0 s Γ I ε defines a manifold in px, y, u, εq-space, which is hence tracked from the inner region I into the outer region III via the intermediate region II. However, for ease of comparison with [13], we will take an orbit-focused approach in the following.

Geometric desingularisation
In this section, we desingularise the dynamics in a neighbourhood of the critical manifold S for Equation (11). Our analysis will proceed in two steps: first, we will blow up the px, yqplane in the extended px, y, u, εq-space, which will allow us to give a rigorous description of the dynamics for ε " 0; in particular, we will define a singular orbit Γ in that limit. The dynamics that is obtained in the various coordinate charts will then be combined into a global description of the flow of Equation (11) near the degenerate px, yq-plane which is uniformly valid in ε, and which hence shows the persistence of Γ for ε P p0, ε 0 s, with ε 0 positive and sufficiently small.
We introduce the quasi-homogeneous blow-up transformation Φ : B Ñ R 4 in Equation (11): x "x, y "ȳ, u "rū, and ε "r n´1ε ; (20) here, B " R 2ˆS1ˆr 0, r 0 s, where pū,εq P S 1 , with S 1 the unit circle, and r 0 is positive and fixed. As the preimage of the hyperplane tu " 0u under Φ equals R 2ˆS1ˆt 0u, S is hence blown up to a cylinder in R 4 ; see Figure 3 for a visualisation.
We will require two coordinate charts in our analysis, which we denote by K 1 and K 2 ; these charts are obtained forū " 1 andε " 1 in (20), respectively, which implies for the blow-up transformation Φ therein. Loosely speaking, the phase-directional chart K 1 covers a neighbourhood of the equator of the blow-up cylinder R 2ˆS1ˆt 0u, while the top of that cylinder is covered by the rescaling chart K 2 ; cf. again Figure 3. The change-of-coordinates transformation κ 12 between charts K 1 and K 2 is given by (23) its inverse κ 21 " κ´1 12 reads κ 21 : px 2 , y 2 , u 2 , r 2 q Þ Ñ´x 1 , y 1 , ε 1 1´n 1¯.

(24)
We define several sections for the flow of Equation (11) -or, rather, of the corresponding blown-up systems in charts K 1 and K 2 : ! px 1 , y 1 , r 0 , ε 1 qˇˇx 1 P r´X, Xs, y 1 P r´Y, Y s, and ε 1 P r0, 1s ! px 2 , y 2 , 1, r 2 qˇˇx 2 P r´X, Xs, y 2 P r´Y, Y s, and r 2 P Here, X and Y , as well as r 0 , are positive constants that are defined as before. The sections Σ in 1 and Σ out 1 clearly correspond to the respective boundaries between the inner region I and the intermediate region II and between the intermediate region II and the outer region III, respectively, when expressed in the original px, y, u, εq-coordinates; recall Section 1. For future reference, we note that the section Σ in 1 is equivalent to Σ in , after blow-up and transformation to chart K 1 , while Σ out 1 corresponds to the section Σ in 2 under the change of coordinates κ 12 defined in (23): Σ out 1 " κ 12 pΣ in 2 q. We emphasise that K 1 is the natural entry chart for the flow of Equation (4) after blow-up due to up" r 1 q decaying monotonically to zero; recall Section 2.1, where "inner" expansions for the corresponding solutions were derived. Hence, any orbit of the associated blown-up vector field will intersect the section Σ in 1 in finite time. As will become apparent below, orbits will enter chart K 2 through Σ in 2 after passage through K 1 ; that step is the direct geometric analogue of asymptotic matching between the "outer" and "inner" solution expansions for (4). In particular, the asymptotics in the rescaling chart K 2 will be regular in r 2 , and will correspond to the "outer" solution expansion. Remark 4. For any object˝in px, y, u, εq-space, we denote the corresponding blown-up object by˝. Moreover, in chart K i pi " 1, 2q, that object will be denoted by˝i.
Remark 5. The blow-up transformation defined in (20) is homogeneous inr. In general, one may make a quasi-homogeneous ansatz of the form x "r α 1x , y "r α 2ȳ , u "r α 3ū , and ε "r α 4ε , where α i pi " 1, . . . , 4q are positive integers; see e.g. [4]. One typically determines α i by finding a distinguished limit in the resulting rescaling chart, i.e., by balancing powers ofr there. An alternative, more systematic approach is produced by the method of Newton polygons [1].
3.1. Dynamics in region III. The dynamics in the "inner" region III is naturally described in the rescaling chart K 2 , where the blow-up transformation defined in (20) is expressed as in (22), corresponding to a rescaling of px, y, uq in powers of ε; substituting into Equation (11), we find in the coordinates of that chart.
By dividing out a factor of r n´1 2 from the right-hand sides in Equation (26), we obtain the desingularised equations where the prime now denotes differentiation with respect to a new rescaled time τ 2 . Next, we define the line 2 " ; for r 2 " ε 1 n´1 positive and fixed, the associated point on 2 corresponds to the point Qὲ before blow-up.
A direct calculation yields Lemma 1. Any point on 2 is a non-hyperbolic steady state for Equation (27), with eigenvalues i and 0.
We remark that the degeneracy of 2 is immaterial to us, as the flow in chart K 2 will not be affected thereby; rather, generic orbits of Equation (27) will bypass 2 entirely, as follows.
For r 2 " 0p" εq, the solution to (27) is given as in (9); hence, orbits reduce to concentric circles in the px 2 , y 2 q-plane, i.e., on S 2 , in the singular limit. For u 2 positive and r 2 " 0, we observe a drift of that family with decreasing u 2 as τ 2 Ñ 8. In fact, we can easily solve Equation (27c) to find u 2 pτ 2 q " " pn´1qpc`τ 2 q ‰ 1 1´n , where c is a constant that is determined by matching -or, rather, "patching" -with the incoming orbit from chart K 1 , as Σ in 2 " κ 12 pΣ out 1 q. Let Γ 2 denote the segment of the singular orbit Γ that is located in chart K 2 , and label the point of intersection of that orbit with Σ in 2 as P in 2 : px in 2 , y in 2 , 1, 0q, which implies u in 2 " 1 and In sum, Γ 2 can be written as Γ 2 :" px 2 , y 2 , u 2 , 0qˇˇx 2 " x 2 pτ 2 q, y 2 " y 2 pτ 2 q, and u 2 pτ 2 q for τ 2 P r0, 8q where x 2 pτ 2 q " x in 2 cospτ 2 q`y in 2 sinpτ 2 q and y 2 pτ 2 q " y in 2 cospτ 2 q´x in 2 sinpτ 2 q, by (9), and where u 2 pτ 2 q is as given in (28).
The smooth persistence of Γ 2 for r 2 positive and sufficiently small then follows from regular perturbation theory: substituting u 2 pτ 2 q into (27b) and taking r 2 pτ 2 q " r 2 constant, we can solve for x 2 pτ 2 , r 2 q and y 2 pτ 2 , r 2 q as follows: x 2 pτ 2 , r 2 q " x in 2 cospτ 2 q`y in 2 sinpτ 2 q`r 2 The orbit Γ 2 hence perturbs in a smooth fashion, for r 2 P r0, r 0 s and τ 2 P r0, T 2 s, with T 2 ą 0 arbitrary, to an orbit Γ 2ε for Equation (27) Remark 6. Equation (30) can alternatively be obtained from Green's function theory [3], as Gpτ 2 , σq " sinpτ 2´σ q is a Green's function for the differential operator d 2 . Alternatively still, one may make use of the fact that Equation (27) is near-integrable to derive the asymptotics in (30); see, e.g., [8] for details and references.
The geometry in chart K 2 is illustrated in Figure 4.  (21); substituting into Equation (11) and dividing out a factor of r n´1 1 from the resulting equations, we obtain where the prime denotes differentiation with respect to a new rescaled time τ 1 .
The steady states of Equation (31) are located in the plane π 1 :" px 1 , y 1 , 0, 0qˇˇx 1 P r´X, Xs and y 1 P r´Y, Y s ( ; here, X and Y are suitably chosen, positive constants, as before. Linearisation of (31) about π 1 shows Lemma 2. Any point in π 1 is a partially hyperbolic steady state for Equation (31), with eigen-values´1, 0 (double), and n´1.
Blow-up has hence resulted in a partial desingularisation of the flow of (11) in a neighbourhood of S ε , as we have gained two hyperbolic directions, by Lemma 2.
Here, we note that Γ1 corresponds to the leading-order inner solution, initiated in the point Q´, after transformation to chart K 1 : in particular, after blow-up, we thus observe exponential decay towards tr 1 " 0u, i.e., to the equivalent of the critical manifold S in K 1 .
Similarly, in the invariant hyperplane tr 1 " 0u, we identify the segment of the singular orbit Γ which corresponds to the leading-order outer solution through the point Q 1 " px 0 , y 0 , 0, 0q; we label that segment Γ1 . To derive a representation for Γ1 , we rewrite Equation (32) by introducing ε 1 as the independent variable: which can be solved explicitly with x 1 p0q " x 0 and y 1 p0q " y 0 to give Hence, the orbit Γ1 can be expressed as Γ1 :" px 1 , y 1 , 0, ε 1 qˇˇx 1 " x 1 pε 1 q and y 1 " y 1 pε 1 q for ε 1 P r0, 1s ( . In summary, the portion of Γ that lies in the intermediate region II consists of the union of the two orbits Γ1 and Γ1 and the steady state at Q 1 . The geometry in chart K 1 is illustrated in Figure 5. 3.2.2. Transition through region II. Next, we establish the persistence of the singular heteroclinic orbit Γ for ε positive and sufficiently small. To that end, we need to connect the two persistent orbits Γ I ε and Γ 2ε , as constructed in Sections 3.1 and 3.2, respectively, via the orbit segment Γ 1ε that is located in region II. Thus, we need to describe the transition between the sections Σ in 1 and Σ out 1 under the flow of (31); recall (25). Specifically, we will approximate the transition map Π 1 , with in other words, Π 1 maps the point P in 1ε :" Γ 1ε X Σ in 1 to the point P out 1ε :" Γ 1ε X Σ out 1 , where x in 1 " x 0`O pεq and y in 1 " y 0`O pεq are defined as in Equations (18) and (19) for n ‰ 2 and n " 2, respectively.
Proof. In a first step, we introduce ε 1 as the independent variable in Equation (31), which gives clearly, Equation (40c) can be solved explicitly for r 1 pε 1 q, as it decouples from the px 1 , y 1 qsubsystem in (40): imposing r 1 pr 1´n 0 εq " r 0 in Σ in 1 , we find which immediately implies (39c), as ε 1 " 1 in Σ out 1 , by (25). In particular, the flow of (40) hence enters chart K 1 through the section Σ in 1 , evolves with increasing ε 1 , and leaves via the section Σ out 1 . Substituting (41) into (40), we can solve the resulting equations for x 1 pε 1 q and y 1 pε 1 q, with initial condition px 1 , y 1 qˇˇt ε 1 "r 1´n 0 εu " px in 1 , y in 1 q; thus, we find x 1 pε 1 q " x 1h pε 1 q`x 1p pε 1 q and y 1 pε 1 q " y 1h pε 1 q`y 1p pε 1 q for the solution to (40a) and (40b), where the homogeneous solution is given by We first assume that n ‰ k`1 k , with k " 2, 3, . . . integer; then, the particular solution reads where 1 F 2 denotes the generalised hypergeometric function [2, Chapter 16]. Given the definition of Π 1 in Equation (38), we need to evaluate (42) and (43) in Σ out 1 , i.e., for ε 1 " 1, to obtain the sought-after asymptotics of x out 1 and y out 1 , as stated in (39). We only outline the argument for x out 1 here, leaving the approximation of y out 1 to the reader: expanding the contribution to x 1 p1q from the homogeneous solution x 1h in (42a) and recalling the expansions for x in 1 and y in 1 from (18), we find It remains to estimate the particular solution in (43a); to that end, we expand Since, however, 1 ă n ă 2 implies that Opε 1 n´1 q " Opεq, only the first term in (45) contributes to the order considered here. Combining Equations (44) and (45), we hence obtain Equation (39a), as claimed; we note that the error is O`ε mint2, 1 n´1 u˘d ue to the fact that 2 ă 1 n´1 for 1 ă n ă 3 2 , whereas the reverse is true for 3 2 ă n ă 2. Moreover, we remark that (39a) is r 0 -independent, as is to be expected.
Finally, the case where n " k`1 k in (40) has to be studied separately, as was also done in [14], since the hypergeometric functions occurring in (43) are not defined then. It can be shown that the particular solutions x 1p pε 1 q and y 1p pε 1 q can be expressed in terms of the Sine Integral and the Cosine Integral [2, Chapter 6] in that case, with Since the corresponding expressions are algebraically involved, we do not reproduce them; however, a straightforward adaptation of the above argument show that the resulting asymptotics of x out 1 and y out 1 is again given by Equation (39) to the order considered here, which completes the proof.
Remark 7. Equation (40) can equally be solved by combining (40a) and (40b) into the secondorder equation the Green's function for the above differential operator is given by Gpsq " pn´1q sin`s n´1˘, which implies the particular solution Expanding that integral, as well as the corresponding homogeneous solution, in ε 1 and evaluating the result in Σ out 1 , one may then again approximate x out 1 and y out 1 as above.
Finally, we consider the case of n ą 2.
Proof. The proof is similar to that of Proposition 1; in particular, the closed-form solution in (42) is still valid when n ą 2. As was done there, we evaluate that solution at ε 1 " 1, i.e., in Σ out 1 . Again, we only outline the argument for x out 1 here, leaving the approximation of y out 1 to the reader: the contribution to x out 1 from the homogeneous and particular solutions is given by (44) and (45), respectively, as before; however, since n ą 2 now, it follows that the Opεq-correction therein is Opε to the order considered here, which implies Equation (48a), as claimed. (A similar observation was made in Section 5 of [14].) The transition through chart K 1 is the geometric analogue of the matching procedure performed by Verhulst in [14]; correspondingly, the occurrence of logarithmic ("switchback") terms in the resulting asymptotic expansions for n " 2 was observed there already. From a geometric point of view, these terms can be shown to be due to resonances between the eigenvalues of the linearisation of Equation (31) about π 1 ; recall Lemma 2. To motivate that observation, we note that the term r 1 ε 1 p" εq in (31b) is resonant of order 2. We can then perform a sequence of near-identity normal form transformations in Equation (31) in order to remove any non-resonant terms therefrom, which yields an approximation for the transition map Π 1 ; see [16] for details. Lemma 3. Let n " 2; then, there exists a sequence of near-identity transformations, with px 1 , y 1 q Þ Ñ pξ, ηq, such that Equation (31) can be written as Solving Equation (49) with a transformed initial condition pξ in , η in , r 0 , r´1 0 εq, we find pξ, η, r 1 , ε 1 qpτ 1 q "`ξ in , η in`ε τ 1 , r 0 e´τ 1 , r´1 0 εe τ 1( 50) to the order considered here. Now, the transition time T 1 between the (transformed) sections Σ in 1 and Σ out 1 can be determined by solving 1 " ε 1 pT 1 q " ε r 0 e T 1 , cf. Equation (38), which gives T 1 "´ln ε`Op1q. Thus, ηpT 1 q " η 0´ε ln ε`Opεq contains a logarithmic switchback term of the order Opε ln εq, in agreement with Proposition 2 above.
Remark 8. Incidentally, logarithmic terms of the order Opε k ln εq are also found in the asymptotics resulting for n " k`1 k , with k " 2, 3, 4, . . . , in the statement of Proposition 1; for k " 2, for instance, we have Evaluating the above at ε 1 " 1, i.e., in Σ out 1 , and expanding the result in ε, we find that the asymptotics of x out 1 and y out 1 contain the switchback terms 4 cosp2qε 2 ln ε and´4 sinp2qε 2 ln ε, respectively. These terms can be shown to be due to resonances at higher order in the governing equations in chart K 1 via a higher-order normal transformation that specifies the Op3q-terms in Equation (49).

Summary.
To summarise, the singular heteroclinic orbit Γ that connects the initial point Q´to the critical manifold S -or, rather, its analogueΓ in the blow-up space -is defined as the union of the orbits Γ I , Γ1 , Γ1 , and Γ 2 obtained in Sections 2.1, 3.1, and 3.2, respectively, as well as of the steady state at Q 1 . Similarly, for ε positive and sufficiently small, the orbit Γε is obtained as the union of the corresponding persistent orbits Γ I ε , Γ 1ε , and Γ 2ε . The global geometry in the blow-up space is illustrated in Figure 6.

Asymptotic expansions
Finally, we derive asymptotic solution expansions for Equation (4) in regions I, II, and III; to that end, we recall that any orbit of (4) will terminate in the rescaling chart K 2 (region III) after passage through chart K 1 (region II) and the inner region I. It follows, in particular, that an expansion for px, y, uqpt, εq which is valid for sufficiently large times t can be derived by tracking a given orbit Γ ε that is initiated in region I through region II and into region III. (Correspondingly, the solution to the governing equations in K 2 , as given in Equation (30), is bounded for τ 2 Ñ 8.) We denote by Γ II ε and Γ III ε the orbits that are obtained from Γ 1ε and Γ 2ε , respectively, after blow-down, i.e., in the original px, y, u, εq-space.
In a first step, we determine the value of the time variable τ in the section Σ out , which is obtained from Σ out Lemma 4. Let ε P p0, ε 0 s, with ε 0 ą 0 sufficiently small, and let τ out be defined by upτ out q " ε 1 n´1 in Σ out ; then, Proof. Recalling the exact solution for upτ q from Equation (10c), we solve upτ out q " ε 1 n´1 for τ out to obtain (51), as claimed.
Proof. The proof is analogous to that of Proposition 4: substituting the expansions for x out 1 and y out 1 from (46) into Equation (53), applying the same trigonometric identities as there, and noting that the particular solution in (53) contributes an Opεq-correction, we obtain (54).
Proof. As in the proof of Proposition 4, we consider the expressions for px, yqpt, εq in Equation (53); however, the particular solution now cannot be neglected, as Opεq " Opε 1 n´1 q due to n ą 2. To account for the contribution therefrom to (53a), for instance, we expand the integrand and then evaluate the resulting leading-order integral to find Substituting (56), as well as the expansions for x out 1 and y out 1 from (48), into (53a) and expanding the result in ε, we obtain the expansion for xpt, εq in Equation (55a), as claimed. The expansion for ypt, εq in (55b) is derived in an analogous fashion. Finally, Equation (55c) is obtained as in the proof of Proposition 4; recall (52c).
We emphasise that the expansions derived in Propositions 1 through 3 do not depend on the arbitrary definition of the intermediate sections Σ in 1 p" Σ in q and Σ out 1 p" Σ in 2 q, as is to be expected; in particular, there is no r 0 -dependence in these expansions.
We summarise our asymptotic expansions for the solutions to Equations (4) and (5) in regions I, II, and III in the following Table 1; here, we recall that τ in " 1 n´1`r 1´n 0´u 1´n 0˘a nd t out " ετ out " 1 n´1 p1´εu 1´n 0 q. Moreover, we note that, in order to be able to write px 1 , y 1 , r 1 qpε 1 q in terms of the original variables pτ, εq from Equations (42) and (47), we need to express ε 1 as a function of those variables: Lemma 5. Let ε P p0, ε 0 s, with ε 0 ą 0 sufficiently small, and let ε 1 be defined as in Equation (21); then, Proof. We recall that ε 1 " ε r n´1 1 " ε u n´1 , by (21); substituting for upτ q from Equation (10c), we obtain (57), as claimed.
To reiterate, the "inner" expansion obtained in region I is a regular perturbation expansion for the solution to the fast system, Equation (5); it is valid on an Op1q-interval in the fast time τ or, equivalently, for t " Opεq, capturing the rapid initial decay towards tu " 0u. The long-time behaviour of solutions is captured by the "outer" expansion in region III, which is uniformly valid in the slow time t, for t " Op1q as ε Ñ 0. These expansions are matched in the "intermediate" region II; in particular, the expressions for x out 1 and y out 1 that are obtained in the transition through that region correspond to the matching conditions given in [14]. III t P rt out , 8q px, y, uqpt, εq: (52) (1 ă n ă 2), (54) (n " 2), (55) (n ą 2) . Absolute error between the asymptotics in Equations (52) (1 ă n ă 2), (54) (n " 2), and (55) (n ą 2) and a numerical solution of Equation (4) for ε P r0.0023, 0.012s and t˚" 50.

Conclusions
In this article, we have presented a geometric analysis of a singularly perturbed two-body problem with quick loss of mass that was first proposed by Verhulst in [14]. Our aim was two-fold: first, we wanted to show that geometric singular perturbation theory [5,6], in combination with geometric desingularisation ("blow-up") [4,7], is well-suited to the analysis of such problems. Second, we wanted to explain geometrically the matching conditions imposed in [14], in the derivation of Verhulst's matched asymptotic expansions. Our conclusions can be summarised as follows.
(1) In contrast to the asymptotic matching performed in [14], whereby "inner" and "outer" expansions are assumed without much justification, the unexpected-seeming scaling in powers of ε resonances between the eigenvalues of the linearisation of the governing equations in the phase-directional chart K 1 in the blow-up space [11]. The matching procedure applied in [14] is thus justified rigorously via our description of the transition through region II in K 1 , which establishes the connection between regions I and III. That region corresponds to an "intermediate" scaling in (5) which is equivalent to the domain of overlap that is conventionally introduced in matching. Our analysis complements that in [16], in the sense that no explicit small-ε approximation for x out 1 and y out 1 was derived there, which precluded a detailed comparison with [14]. Correspondingly, our asymptotic expansions for 1 ă n ă 2 and n " 2 agree with those of Verhulst after blow-down, i.e., after transformation to the original state variables x and 9 x; for n ą 2, we obtain more explicit expansions than are given in [14]. (2) Numerical simulation confirms the expected error of our "outer" asymptotics in ε, as stated in Equations (52), (54), and (55), respectively; cf. Figure 7, where ε up to Op10´2q is considered. Here, ∆xpt˚, εq and ∆ypt˚, εq denote the absolute error in x and y, respectively, between our asymptotics for xpt, εq and ypt, εq, respectively, and a numerical solution of (4) that was obtained in Mathematica with 8-digit relative and absolute precision, where we have chosen t˚" 50. Throughout, we observe agreement with the expected order of the error in ε. We emphasise that, when 1 ă n ă 2, the error behaves like Opε 2 q as ε Ñ 0 for 1 ă n ă 3 2 and like Opε 1 n´1 q for 3 2 ă n ă 2, which is in agreement with Equation (52); see again Figure 7. (The magnitude of the absolute error is due to the numerical value of the corresponding coefficient, which can well be of the order 10 for some values of n.) For illustrative purposes, we also visualise representative time series of ypt, εq for various values of n and ε in Figure 8; here, our asymptotics is displayed in magenta, whereas the numerical solution, which was obtained in Matlab using the standard solver ode45 with default parameters, is plotted in red. We note good agreement on the scale of the figure for moderate values of ε, as well as an apparent increase in the accuracy of our asymptotics with increasing n.
(3) On a related note, we emphasise that we do not address the question of the optimal truncation of the asymptotic expansions in the various regions, as summarised in Table 1.
Since the derivation of these expansions is based on singular perturbation techniques, with ε the (small) perturbation parameter, they are only guaranteed to approximate the solution to (4) well for sufficiently small ε; being asymptotic series, they are expected to diverge. Unlike in the theory of convergent power series, one hence does not automatically obtain a better approximation by including more terms in the truncation. The optimal truncation point, which will typically depend on ε, can be determined by investigating the so-called Gevrey properties of these series expansions; see, e.g. [10] for details and references. (4) While the singularly perturbed system in Equation (5) is atypical, in that it admits an explicit solution in terms of (generalised) hypergeometric functions in the phasedirectional chart, the geometric approach showcased here can equally be applied to systems where no such solutions are available; as is common practice in the application of blow-up, the corresponding transition maps through those charts can be approximated using well-adapted normal forms and invariant manifold techniques. (5) Finally, it may appear that our simplifying assumption of f " 1 in Equation (1) is restrictive. However, it can be shown that our analysis will remain valid for more general f provided (3) holds, in agreement also with [14]. If f is time-independent, a regular expansion as in Equation (16) can still be found in region I; due to the fact that no explicit expression for u will be available in general, that expansion will depend on the Taylor series of f about px, y, εq " px 0 , y 0 , 0q. In region III, the f -dependence of u 2 -and, correspondingly, of the integral terms in Equation (30) -will merely introduce higher-order corrections (in r 2 ) in the resulting asymptotics. Similarly, in region II, the dynamics of r 1 will be f -dependent, affecting our description of the transition through that region at higher order in ε 1 ; in particular, the corresponding Equation (40) may no longer be solvable exactly in terms of special functions. If, additionally, f is permitted to depend on time in (3), further complications arise due to the vector field in (4) being non-autonomous then; however, those can be remedied by including time as an artificial dependent variable both in Equation (4) and in the blow-up transformation in (20), and by expanding the function f about pt, x, y, εq " p0, x 0 , y 0 , 0q. In sum, the qualitative dynamics of (1) will remain unaffected by the specific choice of f given Equation (3), which is corroborated by numerical simulation; recall Figures 1 and 2.