Orbital evolution of eccentric perturbers under dynamical friction: crossing the sound barrier

In a gaseous medium, dynamical friction (DF) reaches a maximum when the orbital speed of a (point-like) perturber moving on a circular orbit is close to the sound speed. Therefore, in a quasi-steady state, eccentric orbits of perturbers approaching the sound barrier (from below) should rapidly circularize as they experience the strongest drag at pericenter passage. To investigate this effect, we extend the solution of Desjacques et al. 2022 for circular DF in a uniform gaseous medium to eccentric Keplerian orbits. We derive an approximation to the steady-state DF force, which is valid for eccentricities as high as $e=0.9$ in a limited range of Mach number around the transition to supersonic regime. We validate our analytical result with 3-dimensional simulations of the gas density response. Although gaseous DF generally dissipates orbital energy, we find that it can be directed along the motion of the perturber near pericenter passage when the eccentricity is $e\gtrsim 0.9$. We apply our results to compute the long-time evolution of the orbital parameters. Most trajectories tend to circularize as the perturber moves into the supersonic regime. However, orbits with eccentricities $e\gtrsim 0.8$ below the sound barrier experience a slight increase in eccentricity as they loose orbital energy. Possible extensions to our analytical approach are also discussed.

⋆ E-mail: robinbuehler@campus.technion.ac.ilMost theoretical studies thus far have assumed that the perturber moves in linear motion.Exact solutions such as e.g.Ostriker (1999)'s are routinely applied to model the impact of DF on proto-planetary systems or on the dynamics of compact stellar binaries (see for instance Iben & Livio 1993;Grishin & Perets 2015;Staff et al. 2016;Grishin & Perets 2016;MacLeod et al. 2017;Antoni et al. 2019;Ginat et al. 2020;De et al. 2020;Everson et al. 2020;Rozner & Perets 2022).However, it would be very desirable to extend the scope and validity of the theoretical results to generic bound (eccentric) orbits.Several pieces of work have investigated the DF experienced by circularly-moving perturbers using a variety of analytical and numerical methods for both collisionless and collisional media (see for instance Tremaine & Weinberg 1984;Sánchez-Salcedo & Brandenburg 2001;Kim & Kim 2007;Kim et al. 2008;Kaur & Sridhar 2018;Sánchez-Salcedo 2019;Banik & van den Bosch 2021;Desjacques et al. 2022).Using linear response theory, Desjacques et al. (2022) developed an analytical approach to compute the DF for a circular motion in a gaseous medium.The salient differences with the corresponding linear motion formula are the absence of a far-field, logarithmic divergence and the appearance of a radial (i.e.perpendicular) component in the DF force.Like the linear-motion result however, the steady-state circular DF peaks for a Mach number M ≃ 1.Therefore, if the steadystate approximation to DF holds, the orbit of a perturber moving on a bound eccentric trajectory should rapidly circu-larize as the perturber looses orbital energy and increasingly moves at supersonic speed.
To investigate this issue further, we build on the approach of Desjacques et al. (2022) to explore Dynamical Friction when the orbital eccentricity is significant.The paper is organized as follows.Section §2 summarizes our computation of the friction coefficient for a generic elliptic orbit; Section §3 shows that our analytical approximation is valid for a range of Mach numbers M ∼ 1; In Section §4 we apply our results to eccentric orbits to study their evolution under the effect of DF; We summarize our results and conclude in Section §5.

General relations
Following Ostriker (1999); Desjacques et al. (2022), the DF force in Newtonian gravity can be generally expressed as where ρg is the density of the unperturbed (uniform) gaseous medium, u = r − rp(t) is the separation vector relative to the current position rp(t) of the perturber, and α(r, t) is the fractional gas density perturbation.In the linear response theory considered here, α(r, t) solves the driven, linearized sound wave equation Here, cs is the speed of sound, whereas h(t) is 1 if the perturber is active and zero otherwise.Transforming to Fourier space and applying Green's method, we can solve for the overdensity and, thereby, express the DF force as , after taking advantage of the Fourier transform d 3 u u u 3 e ik•u = 4π ik k 2 of the Coulomb potential.We have also defined in spherical coordinates for which k = (k, φ k , ϑ k ).Eq. ( 3) is still completely general as far as the orbital motion rp(t) is concerned.

DF for eccentric orbits
Since we are interested in a perturber on a bound eccentric orbit, it is convenient to parameterize the latter with the eccentric anomaly η.Assuming that the motion takes place in the x − y plane, we have where a, e and ϑ are the semi-major axis, eccentricity and true eccentric anomaly respectively.For a perturber orbiting a (massive) companion (located at the origin of coordinates) counterclockwise, the position vector of its eccentric orbit is delineates a circular orbit (e = 0) with identical semi-major axis.A non-zero eccentricity thus perturbs the circular orbit in two ways: it changes i) the physical shape of the orbit (from a circle to an ellipse) and ii) the time lapse along the orbit.As we will see shortly, the second effect dominates across a range of Mach number for which it is possible to derive an accurate prediction for the DF force.
Parameterizing the orbit with the mean anomaly, the DF force can be expressed as The Rayleigh decomposition of e ik•(rp(η)−rp(η ′ ) is particularly powerful for the circular case (see Desjacques et al. 2022) since FDF(η) can then be conveniently expanded on the (spherical) helicity basis {ẑ, ê+, ê−} with ê± =1 √ 2 (iŷ ∓ x), This decomposition can also be used in the eccentric case, although the variation of the orbital radius r(η) makes the calculation tedious.On substituting into Eq.( 8) and performing the Gaunt integral, we arrive at Here, is a characteristic Mach number 1 and M• ≫ M is the mass of the companion.The (complex) friction coefficient I(M(a), e, η) encodes the dependence of the DF force on the nature of the medium and the value of the orbital elements.
Appendix §A outlines an approximation to the steady-state friction coefficient, which captures timing variation in the orbit (i.e.t(η)) relative to the circular case but neglect the change in the orbit radius (i.e.r(η)) The final expression of  16) for a perturber with characteristic Mach Number Ma = 0.9 and eccentricity e = 0.9.The semi-major axis factorizes out and is thus left unspecified.Snapshots of log(α) (represented by the color scale) are shown in the orbital plane at four successive times given by η = 3π 2 (top left panel), 2π, 5π 2 and 3π (bottom right panel).The orbit and the position of the perturber are indicated by a curve and a white circle, respectively.We omit the first half rotation because the wake has not fully developed by that time and is thus not very informative.In each panel, a zoomed-out inset shows the evolution of the far-field density wake as it moves outward the orbit.
I(M(a), e, η) is given by the multipole expansion (A5) and (A13).Appendix §A also demonstrates that this expansion has a short distance logarithmic divergence, which is regulated by truncating the series at some maximum multipole ℓmax.
Projecting the force onto the instantaneous radial and tangential directions êr(η) = cos ϑ(η) x + sin ϑ(η) ŷ and êϑ (η) = − sin ϑ(η) x + cos ϑ(η) ŷ, and using the relation we eventually obtain where are the radial and azimuthal components of the DF force along the trajectory of the perturber.Note that the instantaneous, radial unit vector êr(η) is directed outward, while the azimuthal unit vector êϑ (η) points in the direction of the (counterclockwise) motion.

VALIDATION WITH SIMULATIONS
To validate our approximation, we compute the DF force after solving the driven sound wave equation ( 2) on a 3dimensional grid.
Using the retarded Green's function, we calculate the overdensity α on a regular, 512 3 cubical mesh of length 16a centered on the massive companion, i.e.
where ri are discretized grid coordinates, rp(η) given by Eq. ( 6) is the position of the perturber and η plays the role of the clock.The second equality follows from approximating the Dirac-delta distribution with a Gaussian of width of σ = 0.01a.The simulations assume absorbing boundary conditions at the outer edge of the grid and no accretion on the perturber.They implement the finite time perturbation such that h(η) = 1 for η > 0 and zero otherwise.Fig. 1 displays the evolution of the gas fractional density fluctuation α(r, t) in the orbital plane for an elliptic orbit with (Ma, e) = (0.9, 0.9).Snapshots are shown at four different times corresponding to eccentric anomalies η = 3π/2, 2π, 5π/2 and 3π as indicated in the figure.The instantaneous Mach number is M(η) = 3.9 (resp.0.2) at pericenter (resp.apocenter).At η = 3π 2 the near-field density wake (in the vicinity of the perturber) is nearly circular, leading to a DF force which is close to zero.As the perturber passes through the pericenter, the near-field wake becomes asymmetrical and elongated while the supersonic motion of the perturber produces a Mach cone.All this causes the DF force to rise.The Mach cone lasts until apocenter passage, where the motion becomes subsonic again while the trailing density wake detaches from the perturber and propagates outwards as a spiral shock wave.The fairly symmetric distribution of the near-and far-field density wakes at apocenter minimizes the DF force.In the zoomedout insets of Fig. 1, the spiral shock wave which detached at the first apocenter passage (η = π) can be seen propagating outwards.Note also that the wake density always exceeds the average density, i.e. α(r, η) ≥ 0 everywhere.This arises from the fact that the Green's function is positive definite (∝ 1/r) and the perturber is an overdense perturbation.
Using Eq. 1, we calculate the DF force acting on the perturber for each 3-dimensional snapshot of the overdensity field α(r, t).First, as a consistency check, we tested our simulation setup for the circular case e = 0 to ensure that the size of the box and the resolution are sufficient enough to properly capture the DF force.Our simulation setup successfully recovers the analytical results of Desjacques et al. (2022) when the largest multipole ℓmax ∼ π/(∆/a) is matched to the mesh resolution ∆ = a/32.Next, we produced a suite of "simulations" for the parameter choices e ∈ [0.3, 0.6, 0.9] and Ma ∈ [0.8, 1.0].A comparison between the "simulated" DF force and the analytical approximation based on equations (A5) and (A13) is presented in Fig. 2. The latter assumes steady-state, and only takes into account the dependence of t(η) on eccentricity (i.e.rp(η) ≃ rc(η) as discussed in Appendix §A).
For the finite time perturbation implemented by the simulations, the asymmetry of the perturber's trajectory suggests that, unlike the circular case for which steady-state is achieved exactly after one sound-crossing time tsc = 2a/cs of the system (see Desjacques et al. 2022), convergence to steady-state may occur on a different timescale when e > 0. Notwithstanding, our theoretical predictions appear to reproduce the numerical results reasonably well for the parameter combinations considered here, although discrepancies can be seen at large eccentricities especially around pericenter passage.In general our solution tends to overestimate F ϑ while it underestimates Fr.For Mach numbers outside the range [0.75, 1.05], we have found that our analytical approximation to DF is a poor match to the numerical simulation regardless the eccentricity.
Fig. 2 also shows that Fr reaches a (positive) maximum (the radial force is thus directed outward) in the time interval 3π 2 ≲ η ≲ 2π, which coincides with the minimum of F ϑ .Fr changes then abruptly at pericenter passage and reaches a minimum for η ≃ π 2 , which is somewhat delayed relative to the maximum of F ϑ .The latter turns out to be positive for e = 0.9 so that the azimuthal component is directed along the direction of motion (and thus increases the kinetic energy of the perturber).Note that these extrema occur along the orbit approximately when the instantaneous Mach number of the perturber becomes larger or smaller than its orbit averaged value Ma (i.e.η = π/2 and 3π/2.)

LONG-TERM ORBITAL EVOLUTION
In spite of its limited range of validity, our approximation to the eccentric DF force can be used to calculate the evolution of orbital eccentricity as the perturber crosses the sound barrier.
It is convenient to use dimensionless units in order to calculate the evolution of the orbital parameters.For this purpose, we introduce a characteristic semi-major axis a0 and frequency Ω0, which are related through Kepler's third law Ω0 = (GM•) 1/2 a −3/2 0 .They define the dimensionless variables which we use in the numerical implementation below.
Here, q = M/M• ≪ 1 is ratio of the perturber's to the massive companion's mass and Fr,ϑ = ) for all combination of e = {0.3,0.6, 0.9} (rows top to bottom) and M = {0.8,1.0} (columns left to right).The simulation results implement the finite time perturbation and are shown for the first two rotations of the perturber.The theoretical, steady-state prediction matches well the overall behaviour of the numerical data for all parameter combination, although it is not able to always reproduce the exact values.These discrepancies grow with eccentricity and are most pronounced around pericenter passage.
Eq. ( 15), averaging the rate of change of the orbital elements over one period gives where ρg = ρg Figure 3. Evolution of the semi-major axis and eccentricity across 1000 orbits assuming an initial eccentricity e i = 0.3 and Mach number M i = 0.9, a binary mass ratio q = 10 −3 and a uniform gas density ρg = 10 −3 .The solid line was obtaining by evolving the orbit with the N-body integrator REBOUND with the instantaneous components Fr(η) and F ϑ (η) of the DF force given by Eq. ( 14).The dashed line represents the solution to the coupled ODEs (20) obtained from the orbit-averaged frictions Ia and Ie; the dashed-dotted line shows the effect of ignoring Fr in the calculation of Ia and Ie; the dotted line shows the effect of using the linear motion solution of Ostriker (1999) for the computation of Ia and Ie.The zoomed-in inset focuses on the first 20 rotations. and • ℜ e i(η−ϑ) I Mã, e, η sin ϑ + ℑ e i(η−ϑ) I Mã, e, η (cos ϑ + cos η) , with Mã = M a(ã) = Ωã(Ω0a0/cs) (12).Since these orbit averaged quantities must be evaluated numerically, we found prudent to check our results with the high-precision N-body integrator REBOUND (Rein & Liu 2012).
For this purpose, we set it up with one central mass and a perturber with q = 10 −3 .The initial (t = 0) position and velocity match an unperturbed, elliptic orbit with eccentricity e = 0.3 and orbit averaged Mach number Mã = 0.9.For t > 0, we apply, in addition to the gravitational pull of the central mass, the DF force the perturber would experience if it were moving in an uniform gaseous medium of density ρg = 10 −3 .The smallness of the product q ρg = 10 −6 ensures that the orbital parameters e and ã vary on a timescale significantly longer than the dynamical time, so that steadystate Dynamical Friction holds.Therefore, we shall assume the steady-state approximation to the DF force given in Appendix §A throughout.The component of the DF force are calculated according to Eq. ( 15) using the instantaneous eccentricity and semi-major axis provided by REBOUND.
The results of this simulation are displayed in Fig. 3 as the solid curves.These are compared to the solution to the coupled ODEs Eq. ( 20) with Ia and Ie calculated i) following equations ( 21) and (22) (dashed line) (ii) ignoring the imaginary part (dotted line) and (iii) using the (purely complex) friction coefficient I = I(Mã) derived by Ostriker (1999) in the linear motion case (dashed-dotted line).Unsurprisingly, case (i) matches best the instantaneous evolution given by REBOUND: the eccentric evolution is accurately reproduced, while the evolution of the semi-major axis deviates only by ≈ 1.5% after 1000 orbits.Case (ii) demonstrates that discarding only the real part or, equivalently, the radial component Fr already leads to a noticeable deviation in the evolution of the orbital parameters.The discrepancy is even larger for case (iii), for which the real part is zero while the imaginary part is computed from the linear-motion solution of Ostriker (1999).

Eccentric evolution for Mach numbers M ∼ 1
Fig. 4 shows the integral curves defined by the flow equations ( 20) assuming an initial eccentricity in the range 0 < ei < 1 but a unique, initial semi-major axis ãi ≈ 1.8 corresponding to a Mach number Mã = 0.75.Furthermore, since the vector flow is independent of the product ρgq (which can be absorbed into a redefinition of the time coordinate), we have set ρgq = 1 without loss of generality.Since the perturber loses energy regardless of the choice of ei and ãi (DF transfers orbital energy to the density wake), the orbit always shrinks to smaller semi-major axes.As a result, the characteristic Mach number eventually exceeds the upper bound above which our approximation ceases to be accurate.This occurs when ã(t) ≈ 0.9, at which point we stop the computation of the integral curves.
The eccentric evolution is sensitive to the choice of ei.For ei ≲ 0.7, the orbit tends to circularize by the time Ma exceeds unity, with an effect strongest in the range 0.2 ≲ ei ≲ 0.4.For ei ≳ 0.7, the orbit becomes more eccentric as can be seen from the solid (black) curve, which marks the locus for which de/dt = 0. Fig. 2 suggests a simple, intuitive explanation: near pericenter passage, the azimuthal component F ϑ can be positive at high eccentricities.This increases the kinetic energy (i.e. the orbital energy) of the perturber and, thereby, the distance of the apocenter.As a result, the orbit becomes more elliptic.The converse is true at low eccentricities: F ϑ is negative and thus slows down the perturber near pericenter passage, which tends to circularize the orbit.
In order to quantify this further, we follow the analytical argument of Szölgyén et al. (2022) and introduce the specific angular momentum h = GM•a(1 − e 2 ) and orbital energy ε = − GM• 2a .This allows us to express the eccentricity as A body subject to dynamical friction experiences a change of energy where the velocity is given by Furthermore, DF generates a torque which changes the angular momentum by Integral curves in the (e, ã) plane obtained by solving the system of ODEs 20.We have assumed cs/a 0 Ω 0 = q ρg = 1.0, which implies M ã = ã−1/2 .The range of ã is chosen such that our analytical approximation is viable.Colors represent the rate ⟨de/d t⟩ of eccentricity change.The overall magnitude is arbitrary since ⟨de/d t⟩ ∝ q ρg.However, the locus where ⟨de/d t⟩ = 0 shown as the black curve is robust to the choice of q ρg.On the left of it, orbits circularize while, on the right of it, their eccentricity increases.The dotted curve is an analytical estimate of this boundary based on the change of orbital energy and angular momentum (see text for details).
Therefore, DF changes u = e 2 − 1 (which is a proxy for the eccentricity) by Rather than integrating over a whole orbit, the change of ∆u u can be estimated from the empirical observation that ∆u/u reaches a positive maximum at η = π/2 and negative minimum at η = π.In other words, the loss of eccentricity is maximum at η = π/2, while the gain of eccentricity is largest at η = π.We thus write ∆u Approximating the time interval during which the DF force acts on the body as and using our analytical solution to the DF force provides an estimate for ∆ϵ and ∆h when the gain/loss of eccentricity is maximum and, thereby, an estimate for ∆u/u|tot as given by Eq. 28.Setting ∆u/u|tot = 0 gives the locus shown as the dotted line in Fig. 4, for which the gain and loss of eccentricity balance each other, i.e. de/dt = 0.This prediction is in good agreement with that inferred from the computation of the integral curves (solid black curve).Summarizing, most trajectories will tend to circularize as the sound barrier is crossed.However, orbits with e ≳ 0.8 at characteristic Mach number Ma ∼ 0.8 experience a (slight) increase in eccentricity while the perturber looses orbital energy and moves into the supersonic regime.

DISCUSSION AND CONCLUSIONS
We have investigated the effect of dynamical friction (DF) for a perturber moving on a bound eccentric orbit in a gaseous medium.We have extended the multipole approach of Desjacques et al. (2022) to capture timing variations relative to the circular case through a perturbative expansion in the orbital eccentricity (the "small" parameter) e.However, we have not succeeded in capturing the physical deformation of the orbit (which breaks the planar symmetry) and have thus neglected it.
We have validated our analytical (steady-state) approximation based on timing variations with measurements of the DF force extracted from 3-dimensional simulations of the gas density response.We have found good agreement for characteristic Mach numbers Ma = Ωa/cs (a is the ellipse semimajor axis) in the range 0.75 ≲ Ma ≲ 1.05, even for eccentricities as large as e = 0.9.The reason why the timing variation dominates in this range of characteristic Mach number has remained elusive.The observed agreement indicates also that the finite time perturbation implemented by the 3-dimensional simulations approaches steady-state on a dynamical timescale, that is, the sound-crossing time of the system as in the circular case (see Kim & Kim 2007;Desjacques et al. 2022).Furthermore, snapshots of the gas density response show that, at high eccentricities, the trailing density wake induced by the perturber is a series of concentric, incomplete ring-like patterns produced in "bursts" around pericenter passage.
Like the linear and circular motion case, the DF force with e > 0 exhibits a short-distance, logarithmic divergence when the instantaneous Mach Number is supersonic, regardless of the choice of orbital parameters.This Coulomb (logarithmic) divergence is encoded in our perturbative approach and regularized with the introduction of a maximum multipole (set to match the resolution of the simulations).By contrast, the DF force always converges when the orbital velocity is locally subsonic.
We have also investigated the impact of DF on the longtime evolution of the eccentricity in the range 0.75 ≲ Ma ≲ 1.05 where our theoretical approximation is a reasonable description of the true DF force.The latter leads to orbital decay and the inspiraling of the perturber, such that the characteristic Mach number grows with time.Therefore, initial conditions are laid down at Ma = 0.75 and the system is evolved until Ma = 1.05.We have checked that the time evolution of the orbit-averaged orbital parameters closely matches that obtained from a numerical integration of the instantaneous DF across 1000 orbits.The eccentric evolution depends on the initial eccentricity ei (set when Ma = 0.75): for ei ≲ 0.8, the orbit tends to circularize by the time Ma = 1.05 is achieved while, for ei ≳ 0.8, it becomes more eccentric.At a qualitative level, this behaviour reflects the fact that the tangential component of the DF force can be directed along the motion near pericenter passage when the eccentricity is high.At a quantitative level, the limit between orbit circularization and eccentricity growth is reasonably predicted by comparing the relative loss of specific orbital energy and angular momentum at those orbital positions where the gain and loss of eccentricity are largest.
Our approach, which has focused on a single perturber in an eccentric orbit, can be readily extended to a binary system along the lines of Desjacques et al. (2022).It can also include the self-gravity of the medium, be it gaseous or not.However, extending the scope of this perturbative expansion to any (characteristic) Mach number requires that we can take into account the deformation of the orbit (from a circle to an ellipse).At a technical level, this looks challenging since this contribution implies both a time variation in the separation r(η) between the perturber and its companion as well as a preferred direction in the orbital plane, which make the plane wave expansion (in spherical harmonics) less appealing.Alternatively, for moderate eccentricities e ≲ 0.5 and outside the range 0.8 ≲ Ma ≲ 1 explored here, substituting the instantaneous Mach number of the eccentric orbit into the circular solution of Desjacques et al. (2022) yields a better match to the simulation results (see Fig. 5), but it performs worse than the perturbative approach for Ma ∼ 1. l max = 1/2 q max (M a , e, η) (0.9, 0.3, 0) (0.9, 0.3, π) (0.9, 0.6, 0) (0.9, 0.6, π) (1.0, 0.6, 0) .For all the parameter combinations corresponding to an instantaneous Mach number larger than unity, ℑ(I(Ma, e, η)) exhibits a logarithmic divergence, whereas convergence is achieved when the instantaneous Mach number is subsonic Right : The upper panel shows DF force extracted from two simulations of the linear response density α(r, t) with resolution corresponding to ℓmax = 32 (square symbol) and 64 (cross symbol).Our theoretical predictions with ℓmax = 32 (solid curve) and 64 (dotted curve) are overlaid for comparison.The lower panel displays the fractional difference between the two simulation results.Results are shown for (Ma, e) = (0.9, 0.6).

Figure 1 .
Figure1.Evolution of the gas overdensity α(r, η) computed from Eq. (16) for a perturber with characteristic Mach Number Ma = 0.9 and eccentricity e = 0.9.The semi-major axis factorizes out and is thus left unspecified.Snapshots of log(α) (represented by the color scale) are shown in the orbital plane at four successive times given by η = 3π 2 (top left panel), 2π, 5π 2 and 3π (bottom right panel).The orbit and the position of the perturber are indicated by a curve and a white circle, respectively.We omit the first half rotation because the wake has not fully developed by that time and is thus not very informative.In each panel, a zoomed-out inset shows the evolution of the far-field density wake as it moves outward the orbit.

Figure 5 .
Figure5.Radial and azimuthal components Fr and F ϑ in units of 4π ρgM 2 a (GM/Ωa) 2 for two different parameter choices (Ma, e) = (4.0,0.3) (top panel) and (0.6, 0.3) (bottom panel).Data points represent the DF force extracted from 3-dimensional simulations of the gas density response.The solid curve is the prediction obtained by the method outlined here, whereas the dotted curve follows from substituting the instantaneous Mach number (on the eccentric orbit) in the circular DF solution.

Figure A1 .
Figure A1.Left : Imaginary part ℑ(I(Ma, e, η)) of the friction coefficient (Eq.A5) for several combinations of (Ma, e, η).Results are shown as a function of lmax and qmax such that qmax = 1 2 lmax (see text for details).For all the parameter combinations corresponding to an instantaneous Mach number larger than unity, ℑ(I(Ma, e, η)) exhibits a logarithmic divergence, whereas convergence is achieved when the instantaneous Mach number is subsonic Right : The upper panel shows DF force extracted from two simulations of the linear response density α(r, t) with resolution corresponding to ℓmax = 32 (square symbol) and 64 (cross symbol).Our theoretical predictions with ℓmax = 32 (solid curve) and 64 (dotted curve) are overlaid for comparison.The lower panel displays the fractional difference between the two simulation results.Results are shown for (Ma, e) = (0.9, 0.6).