Flow Morphology of a Supersonic Gravitating Sphere

Stars and planets move supersonically in a gaseous medium during planetary engulfment, stellar interactions and within protoplanetary disks. For a nearly uniform medium, the relevant parameters are the Mach number and the size of the body, $R$, relative to its accretion radius, $R_A$. Over many decades, numerical and analytical work has characterized the flow, the drag on the body and the possible suite of instabilities. Only a limited amount of work has treated the stellar boundary as it is in many of these astrophysical settings, a hard sphere at $R$. Thus we present new 3-D Athena++ hydrodynamic calculations for a large range of parameters. For $R_A\ll R$, the results are as expected for pure hydrodynamics with minimal impact from gravity, which we verify by comparing to experimental wind tunnel data in air. When $R_A\approx R$, a hydrostatically-supported separation bubble forms behind the gravitating body, exerting significant pressure on the sphere and driving a recompression shock which intersects with the bow shock. For $R_A\gg R$, the bubble transitions into an isentropic, spherically-symmetric halo, as seen in earlier works. These two distinct regimes of flow morphology may be treated separately in terms of their shock stand-off distance and drag coefficients. Most importantly for astrophysical applications, we propose a new formula for the dynamical friction which depends on the ratio of the shock stand-off distance to $R_A$. That exploration also reveals the minimum size of the simulation domain needed to accurately capture the deflection of incoming streamlines due to gravity.


INTRODUCTION
Supersonic motion of a gravitating body within a gaseous medium arises in many areas of astrophysics, and the shock waves which adorn such objects facilitate both their gravitational and hydrodynamic interactions with the gas.These problems vary in the physics responsible for the generation of the bow shock, such as wind-wind interactions (Meyer et al. 2014;Scherer et al. 2016;Tarango-Yong & Henney 2018) or MHD shocks from magnetospheres (Kotova et al. 2021).Bow shocks may also be generated through purely hydrodynamical processes, as in the supersonic motion of an accretor.Often, the gravitating body can be modeled as a point source (Hunt 1971), as an accretion boundary (Shima et al. 1985;Ruffert 1994;Ruffert & Arnett 1994;Blondin & Raymer 2012;Xu & Stone 2019) or with a Plummer potential (Sánchez-Salcedo & Brandenburg 2001;Kim & Kim 2009;Morton et al. 2021).In any case, the body feels a drag force due to dynamical friction, which has been analytically modeled for both subsonic and supersonic motion (Ostriker 1999).Accreting bodies feel an additional drag equal to the total momentum of the accreted gas.
Also of interest is the case in which the gravitating body has a rigid surface at a fixed radius (or cannot effectively accrete), which ★ LJP: ljprust@kitp.ucsb.eduhas received relatively little study.Rather than an accretion drag, such an object feels a net force due to the non-uniform pressure distribution over its surface.Shara & Shaviv (1980) computed the flow structure around a hypersonic rigid sphere, though their analysis focused on isodensity contours in the far field.The dearth of studies on rigid spheres was ameliorated by Thun et al. (2016, hereinafter T16), who carried out a thorough investigation of shock stand-off distance and dynamical friction on such objects in homogeneous media.T16 studied supersonic bodies with Mach numbers M ≥ 2 and proposed a new formula for dynamical friction which was demonstrated to perform well for objects with large accretion radii   relative to their physical radii .Here the accretion radius is defined as where  and  ∞ are the mass and velocity of the body.Though the gas accumulated at the surface of the sphere is too optically thick to radiatively cool, preventing actual accretion onto the surface,   is nevertheless useful as a metric of the strength of gravity relative to inertial forces.Flow around a rigid gravitating sphere is primarily relevant to problems in binary stellar interactions or planet-star interactions.One example is the interaction of a star with supernova ejecta from its companion, which was studied by Bauer et al. (2019) as a possible formation mechanism for runaway stars in unusual locations on the Hertzsprung-Russell diagram.Another is the migration and orbital evolution of planetesimals in a protoplanetary disk, which can take place due to aerodynamic drag and dynamical friction as well as other effects such as tides and Lindblad resonances (Grishin & Perets 2015).The atmosphere of a protoplanet is also affected by the presence of a bow shock (Mai et al. 2020), though protoplanets generally do not develop such shocks unless they are on eccentric orbits (Bailey et al. 2021).The use of a rigid surface to model the body in these scenarios is well-motivated by the large entropy contrast between the stellar or planetary material comprising the body and the shocked gas surrounding it.
A prime example is the engulfment of a planet by its host star.As a star ascends the red giant branch (RGB) or asymptotic giant branch (AGB), its increase in radius can result in the engulfment of its inner planets.The planets spiral into the stellar envelope due to drag and experience a wide range of Mach numbers and values of   /.This makes these events an ideal physical application for studies of a rigid gravitating sphere.As such, we focus our parameter space and analysis through the lens of planetary engulfment, though our results are applicable to other problems with similar parameters.
Analytical estimates suggest that engulfed planets of mass ≈10   (where   is one Jupiter mass) may eject a portion of the stellar envelope (Livio & Soker 1984;Nelemans & Tauris 1998;O'Connor et al. 2023).The numerical work of Yarza et al. (2022) confirmed this result, but found that most engulfed planets are instead destroyed through either tidal disruption (for planets which fall within tidal radius of the core) or ram-pressure stripping.Regardless of the mechanism of destruction, the energy that these planets deposit during their inspiral causes a temporary increase in stellar radius and luminosity which may be observable as a transient (MacLeod et al. 2018;Gurevich et al. 2022;De et al. 2023;O'Connor et al. 2023) at a galactic rate of ∼0.1 -1 per year (Metzger et al. 2012).
To determine the ultimate fate of the planet, as well as its energy deposition and eccentricity evolution, it is necessary to study the morphology of the flow near the planet and the resulting forces.Szölgyén et al. (2022) integrated the trajectory of an engulfed body to find that dynamical friction can cause the orbital eccentricity to either increase or decrease depending on the equation of state and the local density gradient, and that this could result in binaries exiting the common envelope phase with non-zero eccentricity.Their semi-analytic models used the dynamical friction prescribed by Ostriker (1999), which is a classic and well-tested result but leaves open the question of which length scales to choose for the Coulomb logarithm.Yarza et al. (2022) performed similar integrations assuming a circular orbit but determined the drag force using a suite of 3-D hydrodynamical simulations.They interpolated both the dynamical friction and hydrodynamical drag from the results of these simulations but remained agnostic as to their functional forms.
Despite the strides made in recent years, there are still several open questions about gravitating spheres such as the dynamical friction on low-mass objects and the form of the hydrodynamical drag.It is also unclear how the results of T16 generalize to lower masses or into the transonic regime.To this end, we perform a series of hydrodynamical simulations of engulfed bodies over a parameter space tailored to answer these questions.This paper is organized as follows.We describe the numerical setup and parameters of our wind-tunnel simulations in section 2. In section 3, we present our findings on flow morphology, including the shock structure as well as the presence of hydrostatic structures.The gravitational and hydrodynamic forces that result from this morphology are analyzed in section 4. We discuss these findings and conclude in section 5.

NUMERICAL SETUP
We perform calculations using the ATHENA++ code, which is described in detail by Stone et al. (2020).ATHENA++ solves the Euler equations

𝜕 𝜌 𝜕𝑡 for a gas with density , fluid velocity , gravitational potential Φ, specific energy  and pressure .We model the gas as ideal with specific heat ratio  = 5/3.A second-order van Leer predictor-corrector scheme is used for time integration with the Courant-Friedrichs-Lewy number set to  CFL = 0.3.We use a Harten-Lax-van Leer contact (HLLC) approximate Riemann solver (Toro et al. 1994) and handle spatial reconstruction via a piecewise linear method.ATHENA++ utilizes a van Leer slope limiter, including the corrections described in Mignone (2014) to keep it total variation diminishing.We perform "wind-tunnel" simulations which model the immediate vicinity of the body.This approach is similar to that of MacLeod & Ramirez-Ruiz (2015), MacLeod et al. (2017) and De et al. (2020) which studied common envelope evolution by modeling the engulfed body as a sink particle within a wind tunnel.Yarza et al. (2022) used a similar setup with a reflective boundary rather than a sink.We use a spherical-polar geometry with the body at the centre of the spherical domain.The radius of the body sets a characteristic length scale and is held constant throughout this paper at  = 0.1R ⊙ , which is a value typical of Jovian planets.In contrast, the accretion radius   is allowed to vary according to the local flow conditions, as we discuss in section 2.3.The radius of the simulation domain is  out = 100, within which we assume that the gas is initially uniform.We use static mesh refinement, reducing the linear size of the mesh by a factor of 2 within a radius of 20 and by a factor of 4 within  = 5.Unless otherwise stated, our base grid before refinement is 640 × 40 × 10 cells in the radial, polar, and azimuthal directions, respectively.This results in 1.6 × 10 6 mesh cells after refinement.In Appendix A, we perform convergence tests for both the linear resolution and the size of our simulation domain.Because the mass of the body far exceeds that of the gas within the domain (by 1 -4 orders of magnitude), we neglect the self-gravity of the gas and treat the sphere as the only source of gravity.

Boundary Conditions
The surface of the body is treated as a reflective boundary.We do not impose the no-slip condition at the surface, as we take the Reynolds number to be effectively infinite and the boundary layer to be negligibly small.For the outer boundary, the upstream ( < /2) and downstream ( > /2) regions are treated separately: (i) At the upstream boundary, the fluid is set to some specified inflow conditions  ∞ ,  ∞ and  ∞ .We assume that the orbital separation is much larger than the domain, so that the curvature of the velocity field is negligible.The incoming streamlines are taken to be parallel so that  ∞ =  ∞ x, though  we investigate the domain size needed for this assumption to be valid in Appendix A.
(ii) The downstream boundary is set to "diode" conditions: gas is not allowed to enter the domain but may exit with zerogradient conditions.
We check that this methodology is valid by comparing against the results of laboratory tests.A key quantity in describing a bow shock is the stand-off distance of the nose of the shock from the body   .This is defined as the distance from the centre of the sphere to the bow shock in the upstream direction (along the  = 0 pole).Billig (1967) complied experimental data from various supersonic wind tunnels to obtain the fit We carry out several test simulations with M in the range 1.5 -6 and set  = 7/5 to be consistent with the laboratory tests in air.
Mach numbers above this range can lead to chemical changes in the air and are not of interest to us.We use a slightly different setup than that described above, using a smaller domain since gravity is absent: the radius of the simulation domain is  out = 50 and the base grid is 1800 × 80 × 40 cells in the radial, polar, and azimuthal directions.
We refine this grid by a factor of 2 within  = 10 and by 4 within  = 5 for a total of 4.27 × 10 7 cells.Fig. 1 shows density and M snapshots from the M = 4 test after it has reached a steady state condition.In Table 1, we compare the stand-off distances obtained from our simulation results against the fit of Billig (1967) through the fractional error We find good agreement between our results and the fit, with errors at or below a few percent.

Non-Linearity Parameter
We perform a suite of simulations in which we vary the mass and Mach number of the body.Here an important quantity is the "nonlinearity parameter" as defined by Kim & Kim (2009), who found that  can be used to describe the shock stand-off distance   via a broken power law: though they did not perform any runs with  < 1.This result was partially confirmed by T16, who found that T16 also introduced a new formula for dynamical friction when 1 <  < 50 and demonstrated that it is a good approximation to their numerical results within this range, but that it greatly overestimates the dynamical friction for  < 1.We discuss the dynamical friction further in section 4.1.
It is clear that further investigations of the  < 1 regime are necessary.To this end, we perform a suite of simulations with 0 ≤  ≤ 1.We also perform several runs with  > 1 to probe the transition into the high- regime, which is well-understood due to the works mentioned above.We consider only Mach numbers greater than one due to the fact that for subsonic flow, there is no bow shock and no local solution for the hydrostatic halo.Mach numbers greater than 4 are also not considered since they are of little importance to planetary engulfment, as we discuss below.

Simulation Realm Applicable to Planetary Engulfment
Exoplanet surveys have identified a population of planets orbiting subgiants at separations less than 1 AU which are expected to be engulfed (Schneider et al. 2011;Wright et al. 2011;Villaver et al. 2014).Lau et al. (2022) carried out global simulations of a hot Jupiter engulfed by an early red giant, which is a numerically tractable problem as the inspiral proceeds on the dynamical timescale of the star.However, many planets are expected to be engulfed at a later evolutionary stage where the envelope has expanded to hundreds of solar radii, requiring many hundreds of orbits to inspiral.This motivates our assumption that the gas within the simulation domain is uniform, as for near-sonic M the local density scale height is comparable to the orbital separation.
To set the stage, we choose conditions typical of the outer envelope of an AGB star:  = 5/3,  ∞ = 6.76 × 10 −9 g/cm 3 and  ∞ = 9109 dyn/cm 2 .The velocity  ∞ is then determined from the specified Mach number M ∞ .The choices of  ∞ and  ∞ do not affect the flow morphology except for their role in determining the sound speed  ,∞ = √︁  ∞ / ∞ .The Mach number of an engulfed planet in a convective envelope is generally expected to be supersonic, though subsonic velocities are possible.Consider the specific energy of (static) stellar material where  enc is the stellar mass enclosed within radius  and the specific internal energy  can be expressed as For a circular orbit appropriate to a slow inspiral,  enc / =  2 .Then the specific energy of the stellar material is If the stellar material is locally bound ( < 0), we arrive at the condition O 'Connor et al. (2023) find that in practice the Mach number always exceeds unity for a range of AGB stellar models generated in the stellar evolution code MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2023).Notably, they also find that  falls in the range 0.05 ≲  ≲ 2 for a Jupiter-mass planet.However, Szölgyén et al. (2022) find subsonic behavior over a portion of the inner envelope for a range of  values including  = 5/3, using solutions of the Lane-Emden equation to obtain the stellar profile.This discrepancy in M could be attributable to differences in stellar structure between these two works, though they also use different drag prescriptions.Regardless, it is clear that the Mach number is low yet supersonic over a large portion of the envelope.
Based on the considerations presented in this section, we choose a parameter space of 1.25 ≤ M ≤ 4 and 0 ≤  ≤ 2, as well as several runs with higher  and M = 4.It is necessary to carry out several runs with varying M for each value of , as there is no guarantee that the flow morphology and drag forces depend solely on .The alert reader will note from (7) that for fixed , an increase in M corresponds to an increase in   .A summary of the initial conditions and derived quantities for our calculations can be found in Table 2.

FLOW MORPHOLOGY
Each simulation is run until it reaches steady state, as determined by the time variation of its derived quantities.This takes approximately one fluid crossing time of the domain   = 2 out / ∞ .Typically a vortex ring briefly forms behind the sphere before being shed downstream.The exception is our non-gravitating M = 1.5 test case with  = 7/5, which retains its vortex ring.Density snapshots from our  = 5/3 runs are shown in Fig. 2, with  increasing to the right and M increasing down each column.The colormap bounds are unique to each column for the purpose of illustrating flow morphology.These snapshots show the  = 0 plane, though we do not find significant azimuthal dependence in any of our 3-D results.The qualitative features of these results are illustrated in Fig. 3 as a primer for the following discussion of these features.

Detached Shock Stand-Off
In Fig. 4 we show the shock stand-off distance as a function of  for our M = 4 runs.We find that   / =  above  = 1 and is constant below, in agreement with T16.However, when we decrease the Mach number toward unity we find that   depends strongly on M, as shown in Fig. 5. Specifically, it is proportional to the cosine of the Mach angle cos  = √︁ M 2 − 1/M, though its scaling with  is ambiguous.On the other hand,  = 2 exhibits altogether different behavior, generally increasing   with M. Because an increase in M corresponds to an increase in   at fixed , it is not guaranteed that increasing M will decrease   in contrast to the non-gravitating case.
The variation of   for a single value of  calls into question the usefulness of the  parameter as a metric, as opposed to the accretion radius.However, we find that cases with  ≤ 1 exhibit a distinctly different flow morphology than those with  > 1 despite overlapping in   , as we now discuss.

Shock Structure
It was shown above that the stand-off distance can be bifurcated into two distinct regions separated by  = 1.We find that the reason for this is that the flow morphology is also distinctly different in these regions.
At low , the bow shock is near the surface of the body, with a shape that is largely determined by the upstream Mach number.Trailing the sphere is a wake comprised of a low-density bubble sheathed in a second shock, visible in many of the panels of Fig. 2. The bubble is quite high in entropy and is fairly isentropic, as shown in Fig. 6.There is a large contrast in velocity between the bubble and the rest of the wake, separated by a slip line.This phenomenon also occurs in the absence of gravity and is well-understood: fluid which crosses the bow shock at near its strongest point at its nose is slowed to subsonic speeds, though this region of subsonic flow in front of the body is fairly small.Indeed, our runs with  = 0 (in which gravity is absent) strongly resemble the non-gravitating test cases of section 2.1, differing only in the value of .Streamlines crossing the outer oblique regions of the bow shock remain supersonic, and it is these streamlines which crest the top of the sphere ( = /2) and undergo supersonic expansion until they separate from its surface.Increasing the Mach number delays boundary layer separation and thus decreases the size of the separation bubble (Loth 2008).The turning of the streamlines as they encounter the separation bubble necessitates a recompression shock which slows the fluid to nearsonic (or possibly subsonic) speed.
Illustration of the qualitative features of our results for low, intermediate and high .For low  (left), the morphology resembles the non-gravitating case, with a small separation bubble on the trailing edge.As  is increased (centre), the separation bubble increases in size and in its extent along the surface of the sphere, driving the recompression shock outward to intersect with the bow shock.Further increase in  (right) causes the separation bubble to envelop the body and transition into an isentropic, hydrostatic halo.As a case study, consider our run with  = 0.4 and M = 2, which we show in Fig. 7.The top left panel shows the Mach number, which exhibits the features discussed above.To quantify the isolation of the bubble from the rest of the flow, we use the Bernoulli constant where is the enthalpy. is conserved along each streamline in steady flow as a consequence of Crocco's theorem (Chapman 2000).Furthermore, since the Bernoulli constant is uniform over the inflow boundary, it is also uniform for all flow originating from that boundary.We see in the top right panel that  is constant except for the separation bubble, indicating that the bubble is isolated from the rest of the flow.It is also significantly hotter than the surrounding gas (bottom left panel) while the pressure is continuous across the slip line (bottom right panel).This pressure can exert significant force on the body, which is discussed further in section 4.2.
As  is increased from zero to unity, the gravitational influence of the body causes the wake to increase in size dramatically.This widening of the wake can lead to an intersection between the two shocks, which is particularly pronounced for  = 1 at high M (Figs. 6 and 8).In Fig. 9 we show the boundary of the separation bubble along the surface of the sphere   , recalling that  = 0 corresponds to the upstream direction.We see that the non-gravitating result - Table 3. Density and pressure at the base of the hydrostatic halo as well as the mass of the halo for our high- runs.The power-law scalings are also listed along with their expected values.that increasing M reduces the size of the bubble -reverses as  is increased.In fact, for  = 1, M = 4 the bubble nearly crests the top of the sphere ( = /2).If  is further increased beyond unity, this double-shock structure disappears entirely (shown in the rightmost panels of Fig. 2) as the bubble is replaced by a hydrostatic halo, as we show below.We see that the transition is characterized by  rather than   , as  = 1 and  = 2 exhibit different morphologies despite overlapping in their coverage of   .

Hydrostatic Structures
Density snapshots from our calculations with  > 2 are shown in Fig. 10.Here, the gravity of the body is strong enough to cloak it in a halo of gas adorned with a single bow shock.The entropy and velocity gradients created by the bow shock have been shown to cause Rayleigh-Taylor and Kelvin-Helmholtz instabilities in the shocked material.This so-called "flip-flop" instability leads to oscillations in the surface of the shock cone (Matsuda et al. 1987;Foglizzo & Ruffert 1999).The instability is more violent for small and absorbing bodies and high Mach numbers (Matsuda et al. 1989(Matsuda et al. , 1991(Matsuda et al. , 1992;;Benensohn et al. 1997;Shima et al. 1998) and is suppressed in 3-D simulations compared to similar work in 2-D (Ruffert & Arnett 1994).A thorough survey of these instabilities and their relationship with dimensionality and numerical artefacts can be found in Foglizzo et al. (2005).Indeed, eddies and vortices are visible in the low-speed flow within the halo in Fig. 10, causing modest time variation in the flow even after a steady state has been reached.Because of this, we time-average our derived quantities for these runs.Displayed in Fig. 11 are radial slices of the density, pressure, and pseudoentropy for  = 10.The nose of the bow shock is sharply defined in the  = 0 slice.At a sufficient distance behind the shock, we see that the hydrodynamic variables in each of the three radial slices are identical, demonstrating spherical symmetry of the halo.T16 derived an analytical model for the radial profile of the halo in hydrostatic equilibrium (HSE) given the stand-off distance and post-shock fluid state (subscript 2): We recast this using   = , the definition of  and thermodynamic relations as This prediction and its associated pressure slightly underestimate our pressure and density profiles (Fig. 11).The pressure and density at the surface of the sphere  = , which may determined from ( 17), far exceed those of the ambient medium.For  ≫ 1,  ≫ ℎ 2 , the scaling relation   ∝  1/(−1) =  3/2 holds at constant M. For pressure, we similarly obtain   ∝    ∝  5/2 .The total mass contained within this halo can be obtained by integration of (17) over the halo, taking   as its radial extent.We perform this integration in Appendix B and show that  halo ∝  ∞  3 .In Table 3, we list the values of   ,   and  halo measured from our simulation output at high .We also compute the power-law scaling for each variable and compare to the values predicted above.The halo mass scaling is in good agreement with our prediction, whereas the pressure and density scalings fall slightly short.This indicates that the assumptions  ≫ 1,  ≫ ℎ 2  are not valid over this range of parameter space, which is unsurprising as the largest  we have considered is 10.We expect these scaling relations to be obeyed at much larger values of .
As mentioned in section 3, the time required for these simulations to reach a steady state is approximately the fluid crossing time of the simulation domain   = 2 out / ∞ .Another relevant timescale is that required to form a halo of mass  halo .Approximating the halo cross-section as  2  , the mass flux feeding the halo is  2   ∞ .Then the time needed to form the halo is Since  halo scales as  3  , this timescale goes as  halo ∝   / ∞ .This is the fluid crossing time of the halo, which is clearly less than the crossing time of the simulation domain, so the limiting timescale in our simulations is   .
In addition to the halo, HSE can also occur within the separation bubbles discussed in section 3.2.Such structures cannot be described by ( 17), as the assumption underlying (17) that the bow shock provides a well-defined outer boundary condition is no longer valid.Instead, we turn to the fundamental condition for HSE The left and right sides of ( 19) are computed from our simulation data and compared in Fig. 12 for a separation bubble and a halo.We see that for  = 0.4 the separation bubble is in hydrostatic equilibrium within several  of the body, whereas for  = 10 the equilibrium is maintained out to at least the stand-off distance   = 10.Fig. 11 also shows that the entropy within the halo is constant and uniform, varying at the percent level with the exception of some fluctuations near the shock front.This constant value is set by the entropy generation at the shock, located at  1 (), which depends on the incoming Mach number at the bow shock M 1 .From conservation of energy, assuming a constant sound speed prior to the shock.If one assumes that the shock is located at  1 ≈   to obtain a lower bound, it follows that Values of  1 greater than   are also possible in the outer regions of the shock far from the body, but here the shock is weak and does not generate significant entropy.Furthermore, since the shock is oblique, the Mach number normal to the shock M 1, will be no greater than that predicted by (21).It is not possible to analytically determine M 1, , but we show our numerical results in Fig. 13.We determine M  = M sin  by fitting the shape of the shock to determine its angle relative to the local velocity field .Without gravity,  is equal to the local shock angle and varies between  and the angle of the Mach cone  = arccos( √︁ M 2 − 1/M).The presence of gravity alleviates this variation somewhat by turning the streamlines normal to the shock.We find that M 1, falls off more slowly with  as gravity is increased.We also see that the position of the sonic line   increases with mass.Thus, for increasing   the strength of the bow shock and the entropy generated by that shock become more uniform.
Once M 1 has been determined, the entropy ratio can be found using ( 21) and the Rankine-Hugoniot shock jump conditions where the subscript 2 denotes post-shock conditions.The pseudoentropy ratio is then We can use our lower bound on M 1 (21) in conjunction with the jump conditions to set a lower bound on the entropy.For M ∞ = 4, we predict that From Fig. 11, we see that the numerical value is ≈6, which falls above our bound.This bound can be applied to the rest of our runs as well, as shown in Fig. 14.For our low- runs, which do not have a hydrostatic halo, we take the entropy at the base of the separation bubble.We see that the entropy converges toward our lower bound as  is increased.

DRAG FORCES
Both gravitational and hydrodynamic forces are at work on the body.The overdense wake and Mach cone trailing the sphere lead to a net gravitational force, known as dynamical friction.The non-uniform pressure distribution over the surface of the body owing to the ram pressure, shock structure and possible vortices also provides a net force.

Calculating Dynamical Friction
The dynamical friction on a supersonic body can be expressed as (Ruderman & Spiegel 1971), where the Coulomb logarithm depends on the effective linear size of the body  min and the maximum extent of the wake  max .For an infinite uniform medium, the extent of the wake is  max =  ∞ , so (26) predicts that  max → ∞ and  DF → ∞ as  → ∞ (Ostriker 1999).This is unphysical, as the increasing dynamical friction will slow the object into the subsonic regime where ( 26) is no longer valid.Thus, there is no steady-state solution for dynamical friction in an infinite medium.In a planetary engulfment, the size of the giant provides an upper bound for the extent of the wake, but  max remains largely unconstrained and is often approximated as the orbital separation.The choice of  min is ambiguous, with authors finding a variety of fits (Sánchez-Salcedo & Brandenburg 1999;Kim & Kim 2009;Cantó et al. 2011;Bernal & Sánchez-Salcedo 2013;Morton et al. 2021).However, these works use either a point mass or a Plummer potential and take the softening length   as a measure of the size of of the engulfed body.It is not clear how   relates to the physical radius  of a solid sphere and whether the scaling relations found for a Plummer potential also hold for a solid body, particularly when the    shock is near the surface of the body.T16 showed that the stand-off distance   is an appropriate choice for  min in the range 1 <  < 50, but that this choice overestimates the dynamical friction for  ≤ 1 (see their Fig.16).Due to our use of a uniform medium, we are ill-equipped to inform the choice of  max , but we can determine the effective linear size  min .We measure the dynamical friction from our simulation output as where   is the volume of the -th cell.This can also be expressed as a drag coefficient via so that We show   in the top panel of Fig. 15 for  <  max <  out (solid lines) as well as the prediction of T16 (dotted lines), which asserts that  min =   .We find that   ∝  max , as expected, but that for low  and M the curves are offset from their predicted values.Note from (30) that ln( min ) sets the "intercept" of the curve.Therefore, an offset between the two curves indicates that  min differs from   .Near the outer boundary, the slope of   vs ln( max ) drops below 4 ∞ (/ ∞ ) 2 , the value expected from (30).We find that this discrepancy is attributable to the finite size of the domain, as variation of the domain size  out reveals that this feature always occurs near  =  out .This and other effects of the finite domain are discussed in Appendix A.
For a given  max , we can determine  min by inverting (26): We see that  min is fairly uniform throughout the domain except for  ≈  and in the outer regions of the domain (bottom panel of Fig. 15).This excess in  min at large  max is unphysical and is due to the and  min (bottom) as a function of wake extent  max for  = 1 with M = 4 (blue lines) and M = 1.5 (orange lines).For intermediate radii, the slope of   agrees well with theoretical predictions so that  min plateaus.We choose to measure  min at  = 20 (black vertical line), as this typically falls within the plateau.
finite domain size, as mentioned above.Thus, to obtain an accurate  min we need only choose  max significantly greater than   and well within  out .For consistency, we choose  max = 20, which lies within the plateau evident in the bottom panel of Fig. 15.We find that  min can be well-approximated by the fit  min /  = 2/  , which can be rewritten in the alternative forms where  = arccos( √︁ M 2 − 1/M) is the Mach angle.We check this result by determining  DF using ( 26), where we plug in (32) and  max = 20.Here the shock stand-off distance is the only quantity we need enter "by hand."This prediction is shown in the top panel of Fig. 16 (dotted lines) along with the values measured from our simulations (solid lines).We also show the ratio between our predictions and measurements in the bottom panel.We find agreement between our predictions and measurements, particularly for  = 1.To obtain an estimate of the dynamical friction across all , we combine our results with those of T16 for  > 1: At  = 1, these predictions differ by a factor of M 2 /(M 2 − 1), which is near unity for sufficiently large M.   is often negative within a few  of the sphere, but far from the body is always positive and obeys the expected   ∝ ln( max ) behaviour.Indeed, the point  = 0.4, M = 1.25 is an outlier in Fig. 16 because the dynamical friction is negative even at  max = 20.Thus, the choice of  max determines not only the magnitude but also the sign of  DF .Rephaeli & Salpeter (1980) showed that the dynamical friction on a supersonic body contains another term in addition to the Coulomb logarithm which becomes large near M = 1: If this term is included, we need only modify our formula (33) via the substitution The pressure is uniform for high , while the separation bubble is clearly defined for low .

Pressure Drag
An engulfed body also experiences hydrodynamic drag as a result of a non-uniform pressure distribution over its surface.The pressure distribution for our M = 4 runs is shown in Fig. 17.The pressure is uniform for high , while for low  the boundary of the separation bubble is sharply defined.The force provided by the bubble can rival or even (for  = 1) exceed that on the front of the sphere.The pressure drag can be expressed as for some pressure drag coefficient   .The dependence of   on Mach number and Reynolds number at very high Reynolds number is not precisely known, but is usually quoted as anywhere from 0.5 to 1 (Bailey & Hiatt 1972).We measure   from our simulation output via the surface integral of the pressure distribution Then ( 36) can be used to calculate   , which we plot in the top panel of Fig. 18.We find that   is of order unity for  < 1, while for  = 1 it transitions from 1 all the way to −1.For  = 2 the drag coefficient goes to zero as   is increased, and is near zero for larger values of .This is unsurprising due to the sphericallysymmetric halo at high , though persistent fluctuations in the flow due to instabilities cause   to oscillate around zero.
For the purposes of comparison with the dynamical friction, Yarza et al. (2022) plotted the quantity   (/  ) 2 against /  , as the prefactors in ( 36) and ( 26) differ by a factor of (/  ) 2 .We follow suit in the bottom panel of Fig. 18, including   for comparison, though we emphasize that the behaviour of   depends on the choice of  max which we have chosen arbitrarily.We find slightly negative values of   over a range of /  , though dynamical friction dominates the total drag in this regime regardless.These results indicate that the standard form for the hydrodynamical drag (36) with a drag coefficient near unity is sufficient for the regime in which   is significant.For  < 1,   remains near unity, for  = 1 it varies from 1 down to −1.At higher ,   tends toward zero and may be slightly positive or negative.(Bottom) Gravitational drag coefficient (+ signs) and pressure drag coefficient multiplied by (/  ) 2 (circles) for the purpose of comparison.
For moderate   , we find slightly negative values, though they are much smaller than   regardless.

Domain Size
In Appendix A, we perform a convergence test in which the domain size  out is varied.This is necessary in the presence of a gravitating body, as the upstream flow is deflected by gravity in the far field.Thus, it is necessary to use a domain large enough that our assumption that incoming streamlines are parallel along the x-direction is valid.The results of this test are listed in Table 2, where the drag coefficients are shown to vary by ∼50 per cent over our parameter space, converging at  out ≈ 50  .We conclude that this is the minimum domain size needed to obtain accurate measurements of the drag coefficients in the presence of gravity.

DISCUSSION AND CONCLUSIONS
In this paper we investigated the morphology of the flow around a rigid supersonic gravitating sphere engulfed in a homogeneous gaseous medium using 3-D hydrodynamical simulations in ATHENA++.We focused our parameter space and analysis on the case of a planet engulfed by a giant star, though our results are applicable to any problem with a range of Mach numbers and   / comparable to what we consider here.In particular, we focus on the  ≤ 1 regime, which has received relatively little study.We also extend the supersonic results of T16 to lower M, as engulfed planets are expected to have Mach numbers near unity over a large portion of the stellar envelope.The shock stand-off distance is found to agree with T16 at high Mach number, but departs significantly as M approaches unity.We find that there are two distinct regimes of flow morphology.At low , the morphology resembles that of a non-gravitating sphere: a bow shock forms near the surface of the body, while downstream a recompression shock is provided by flow separation from the surface.However, the flow departs from the non-gravitating result due to gravitational acceleration and deflection of the flow upsteam of the bow shock.Furthermore, the base of the separation bubble is hydrostatically-supported and the bubble grows in size as the mass of the gravitating body is increased, though it is not entirely clear if this growth of the bubble is due to its own hydrostatic support or to changes in morphology upstream.The transition between these two morphologies occurs at  ≈ 1.Here the wake grows large enough that the recompression shock intersects with the bow shock, and a further increase in  causes these two shocks to combine into a single bow shock.Based on the range of  over the inspirals computed by O'Connor et al. (2023), this transition is expected to occur within the stellar envelope during a planetary engulfment event.
For high , the sphere is cloaked in a spherically-symmetric hydrostatic halo with a stand-off distance described by (9) and density profile (17), as shown by T16.The entropy of the halo is set at the shock front, and is fairly uniform due to the deflection of the streamlines by gravity prior to reaching the shock.The entropy gradients which remain are near the bow shock and produce instabilities which cause oscillations in the shock cone, but are suppressed in 3-D and for the relatively small values (< 20) of   / we consider here (Ruffert & Arnett 1994;Foglizzo et al. 2005).The hydrostatic profile of the halo (17) can be determined by setting an outer boundary condition at the stand-off distance using the jump conditions.The fluid state at the centre of the halo (i.e. at the surface of the body) as well as the halo mass can be determined from this profile, and scale as powers of .
These morphologies have consequences for the drag forces experienced by the body.For the dynamical friction, the two key unknowns are the maximum extent of the wake and the effective linear size of the perturber.We find that the effective linear size can be approximated by the ratio of the stand-off distance to half the accretion radius for  ≤ 1.The theoretical underpinnings of this result are not clear, though it performs well for  ≤ 1.When  > 1, we defer to the result of T16 that the effective linear size is equal to the stand-off distance.Together, these results offer a complete model of dynamical friction for a supersonic gravitating sphere (33).
The pressure drag is similarly bifurcated: in the case of a hydrostatic halo, the spherical symmetry of the fluid near the body results in a small pressure drag.This compounded by the fact that such objects have high dynamical friction, so we conclude that the pressure drag is negligible for  > 1.For low , the pressure drag is well-approximated by the standard non-gravitating case (36) with a drag coefficient near unity.This is in agreement with studies on the migration of planetesimals by Grishin & Perets (2015).
Finally, we investigate the effects of the simulation domain size on our results.Because we assume that the incoming streamlines are parallel, the use of a small domain risks underestimating the deflection of the streamlines due to gravity which can cause errors of ∼50 per cent in the drag coefficients.We find that a domain size comparable to ∼50  is necessary for convergence.We have assumed that the accretion radius is small compared to the length scale of any gradients in the medium (e.g. the local scale height of a stellar profile) so that we can assume a uniform medium, though this may not be valid for massive planets or small orbital separations in a planetary engulfment event.
We have not considered subsonic velocities, as they are not well-suited to our methods.There is no a priori reason why subsonic velocities cannot occur in any of the physical scenarios we have described, so future study on this regime is needed for a complete understanding of flow around a gravitating sphere.vary the number of cells in the radial direction   .The cells in the polar and azimuthal directions   and   are scaled proportionally.We show in Fig. A1 the change in our principal derived quantities   ,   and   as a function of   .We see that for   ≥ 640 the shock stand-off distance and gravitational drag coefficient vary on the percent level.The coefficient of pressure drag exhibits significant fluctuations, though these do not disappear even at high resolution.
We conclude that our results for   and   are converged, while the error in   may be of order ≈10 per cent.
We have assumed in this work that the streamlines of the fluid entering the domain from the boundary are initially parallel.However, in the presence of gravity the streamlines should have deflected by some amount even at arbitrarily large .One solution is to approximate the deflection angle of the streamlines at some radius from the gravitating body, though it is not possible to write an exact solution for this deflection.We choose to tackle this issue by ensuring that the radius of our domain is large enough that the deflection at the outer boundary is negligible.To this end, we perform a series of runs with  = 1 and M = 4 with varying values of  out , as shown in Table 2.As the domain is enlarged, we increase the number of cells in the radial direction such that the linear resolution is preserved.Thus, any discrepancies between these simulations can be attributed purely to the domain size.Our results in Fig. A2 show that the drag coefficients vary by ∼50 per cent before convergence, illustrating the large errors that can result from a domain of insufficient size.The flow morphology is also affected by the streamline deflection.From Table 2 we see that the extent of the separation bubble   changes by over 12 degrees as the domain size is varied, even cresting the top of the sphere ( = 90 deg).We do not find any change in the shock stand-off distance.The majority of the calculations presented in this paper use  out = 100, which appears to be near convergence for both   and   .For high  values, we extend the domain size to accommodate the large accretion radii, scaling by  such that  out = 100 ≈ 50  .The refinement boundaries are also scaled by .

APPENDIX B: HALO MASS
The density profile in the hydrostatic halo, as shown in T16 and confirmed in section 3.3, is .

(B1)
Here the subscript 2 denotes post-shock conditions.Conservation of energy gives the pre-shock Mach number and because   =  at high , this simplifies to This can then be used in conjunction with the jump conditions and ambient medium parameters to find the post-shock fluid state.Using the definitions of   and ℎ, we rewrite (B1) as From the discussion above, we conclude that the quantity in brackets in (B8) does not scale significantly with .Thus, the scaling is determined by the prefactor: Furthermore, since the jump conditions give we can plug  back in to obtain Using the definitions of   and , this becomes which scales as or Therefore, for given upstream flow conditions the mass of the hydrostatic halo within   scales as the ambient density times the cube of the accretion radius, or equivalently as  ∞  3 .
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Density and Mach number slices from our  = 7/5, M = 4 test case in its steady state.

Figure 2 .
Figure 2. Density snapshots from our  = 5/3 runs with  ≤ 2, omitting M = 3.5 for space.The colormap bounds for each column are unique for the purpose of illustrating the morphology.

Figure 4 .
Figure 4. Variation of the shock stand-off distance   with the non-linearity parameter  for M = 4.Our results (blue dots) are consistent with the result of T16 that   =  for  > 1 (orange line) and constant for  ≤ 1.

Figure 5 .
Figure 5. Shock stand-off distance   as a function of M for several values of .Significant deviations from the high-M value are seen for M approaching unity, which scales with the cosine of the Mach angle.

Figure 6 .
Figure 6.Entropy slices for  = 0 (top),  = 0.4 (middle) and  = 1 (bottom) with M = 4, normalized to the entropy of the inflowing gas.The entropy of the separation bubble is much higher than the surrounding gas.All plots have identical colormap limits for the purpose of comparison.

Figure 7 .
Figure7.Slices of various quantities for  = 0.4, M = 2.The Mach number (top left) contrasts regions of the flow which are subsonic (blue), supersonic but less than the inflow Mach number (orange) and exceeding the inflow Mach number (green).Several features of the morphology are visible such as the sonic line, supersonic expansion just upstream of the recompression shock and near-zero velocities within the separation bubble.The Bernoulli constant (top right) is uniform over the domain except for the separation bubble, indicating that it is isolated from the rest of the flow.The enthalpy (bottom left) is higher within the bubble than outside the bubble despite the pressure (bottom right) being continuous across the contact discontinuity.

Figure 9 .
Figure 9.The boundary of the separation bubble along the surface of the sphere ( = 0 is along the upstream direction).For  = 0 we reproduce the expected result that the separation bubble shrinks with increasing M. As  is increased, this trend reverses.

Figure 11 .Figure 12 .
Figure 11.Density (top), pressure (middle) and pseudoentropy (bottom) profiles along radial rays from the centre of the sphere for  = 10.Pressure and density predictions for the hydrostatic halo by T16 are also shown, which slightly underestimate the values we find here.Within the halo, the three radial slices are nearly identical, demonstrating spherical symmetry.The entropy fluctuates near the bow shock but is constant over much of the halo.

Figure 13 .Figure 14 .
Figure 13.Mach number M and its component normal to the bow shock M  immediately pre-and post-shock (subscripts 1 and 2, respectively) for  = 0 (top),  = 1 (middle) and  = 2 (bottom) with M ∞ = 4.The solid black line denotes the sonic line   , which moves to higher  as  is increased.

Figure 16 .
Figure 16.(Top) Comparison of the dynamical friction measured at  = 20 (solid lines) against the theoretical value using our definition of  min (dotted lines).(Bottom) The ratio between the numerical and theoretical values.

Figure 17 .
Figure17.Pressure distribution over the surface of the engulfed body at M = 4 for all .The pressure is uniform for high , while the separation bubble is clearly defined for low .

Figure 18 .
Figure 18.(Top) Coefficient of pressure drag as a function of Mach number.For  < 1,   remains near unity, for  = 1 it varies from 1 down to −1.At higher ,   tends toward zero and may be slightly positive or negative.(Bottom) Gravitational drag coefficient (+ signs) and pressure drag coefficient multiplied by (/  ) 2 (circles) for the purpose of comparison.Simulations with  = 0 are excluded from the bottom panel since they have /  = ∞.For large   ,   (/  ) 2 → 0 as   → ∞, as expected.For moderate   , we find slightly negative values, though they are much smaller than   regardless.

Figure A1 .Figure A2 .
Figure A1.Convergence of several of our key derived quantities with linear resolution.The stand-off distance and gravitational drag coefficient converge nicely, while the pressure drag exhibits some fluctuations even at high resolution.

)
Then the total mass contained within the halo out to  =   is

Table 2 .
Input parameters and derived quantities for our calculations performed in this work.From left to right, the quantities are: non-linearity parameter , inflow Mach number M ∞ , accretion radius   , mass of the gravitating body  in units of Jupiter mass   , specific heat ratio , radius of simulation domain  out , shock stand-off distance   , boundary of the separation bubble   , pseudoentropy in the separation bubble or halo   , pressure drag coefficient   , gravitational drag coefficient   measured at  = 20 and effective linear size of the perturber  min .The sound speed  ,∞ = √︁   ∞ / ∞ depends only on , as  ∞ and  ∞ are constant.