How Drag Force Evolves in Global Common Envelope Simulations

TL;DR

This study investigates how drag forces evolve during common envelope phases in binary star systems through global simulations, revealing when traditional models apply and explaining late-stage drag reduction.

Contribution

It provides detailed global simulation data on drag force evolution in common envelope events, highlighting discrepancies with analytic models and the impact of flow symmetry.

Findings

Drag force matches analytic predictions during initial in-spiral

Drag force significantly weaker after first periastron passage

Flow symmetry reduces net drag at late stages

Abstract

We compute the forces, torque and rate of work on the companion-core binary due to drag in global simulations of common envelope (CE) evolution for three different companion masses. Our simulations help to delineate regimes when conventional analytic drag force approximations are applicable. During and just prior to the first periastron passage of the in-spiral phase, the drag force is reasonably approximated by conventional analytic theory and peaks at values proportional to the companion mass. Good agreement between global and local 3D "wind tunnel" simulations, including similar net drag force and flow pattern, is obtained for comparable regions of parameter space. However, subsequent to the first periastron passage, the drag force is up to an order of magnitude smaller than theoretical predictions, quasi-steady, and depends only weakly on companion mass. The discrepancy is exacerbated for larger companion mass and when the inter-particle separation reduces to the Bondi-Hoyle-Lyttleton accretion radius, creating a turbulent thermalized region. Greater flow symmetry during this phase leads to near balance of opposing gravitational forces in front of and behind the companion, hence a small net drag. The reduced drag force at late times helps explain why companion-core separations necessary for envelope ejection are not reached by the end of limited duration CE simulations.

Equations

$${\mathbf{F}}_{2-\mathrm{gas},1}={\mathbf{F}}_{2-\mathrm{gas}}-(M_2/M_{1,\mathrm{c}}){\mathbf{F}}_{1-\mathrm{gas}},$$

$${\tau }_z=\frac{M_2}{M_{1,\mathrm{c}}+M_2}a{F}_{1-\mathrm{gas},\varphi }+\frac{M_{1,\mathrm{c}}}{M_{1,\mathrm{c}}+M_2}a{F}_{2-\mathrm{gas},\varphi }.$$

$$\begin{array}{rl}{J}_{1-2,\mathrm{z}}=& M_{1,\mathrm{c}}[({\mathbf{s}}_1-{\mathbf{s}}_{\mathrm{CM}})\times ({\mathbf{v}}_1-{\mathbf{v}}_{\mathrm{CM}}){]}_z\\ & +M_2[({\mathbf{s}}_2-{\mathbf{s}}_{\mathrm{CM}})\times ({\mathbf{v}}_2-{\mathbf{v}}_{\mathrm{CM}}){]}_z,\end{array}$$

$$\dot{W}={\mathbf{F}}_{1-\mathrm{gas}}\cdot {\mathbf{v}}_1+{\mathbf{F}}_{2-\mathrm{gas}}\cdot {\mathbf{v}}_2$$

$$E_{1-2}=\frac{1}{2}{M}_{1,\mathrm{c}}{v}_1^2+\frac{1}{2}{M}_2{v}_2^2-\frac{\mathrm{G}{M}_{1,\mathrm{c}}{M}_2}{a}$$

$$R_{\mathrm{a}}=\frac{2\mathrm{G}M}{c_{\infty }^2+v_{\infty }^2}$$

$$F\sim \dot{M}{v}_{\infty }\ln \left(\frac{r_{\max }}{r_{\min }}\right)\sim \frac{4\text{uppi}{\mathrm{G}}^2{M}^2{\rho }_{\infty }{v}_{\infty }}{(c_{\infty }^2+v_{\infty }^2{)}^{3/2}}\ln \left(\frac{r_{\max }}{r_{\min }}\right).$$

$$F_0\equiv \frac{4\text{uppi}{\mathrm{G}}^2{M}_2^2{\rho }_0{v}_0}{(c_0^2+v_0^2{)}^{3/2}},$$

$$R_{\mathrm{a},\mathrm{DM}}=\frac{R_{\mathrm{a},0}}{1+R_{\mathrm{a},0}^2/(4H_{\rho }^2)},$$

$$R_{\mathrm{a},0}=\frac{2\mathrm{G}{M}_2}{c_0^2+v_0^2}$$

$$F_{\mathrm{DM}}=F_0\ln \left(\frac{r_{\max }}{r_{\min }}\right){\left(\frac{R_{\mathrm{a},\mathrm{DM}}}{R_{\mathrm{a},0}}\right)}^2.$$

$$F_{O99}=-I\frac{4\text{uppi}{\mathrm{G}}^2{M}^2{\rho }_{\infty }}{v_{\infty }^2},$$

$$I=\frac{1}{2}\ln \left(\frac{1+{\mathcal{M}}_{\infty }}{1-{\mathcal{M}}_{\infty }}\right)-{\mathcal{M}}_{\infty }.$$

$${ϵ}_{\rho }=\frac{2\mathrm{G}{M}_2}{v_{\infty }^2{H}_{\rho }},$$

$$F_{M17}\simeq 0.6{\mathrm{e}}^{0.61{ϵ}_{\rho }}\times \frac{4\text{uppi}{\mathrm{G}}^2{M}_2^2{\rho }_0}{v_0^2}.$$

$$t_{\mathrm{d}}\sim \frac{Nl}{c}\sim \frac{R^2\rho {\sigma }_{\mathrm{T}}}{m_{\mathrm{H}}c},$$

Full text

1 Introduction

Common envelope evolution (CEE) is the most natural mechanism for rapidly tightening binary orbits and likely facilitates many phenomena, including gravitational wave-emitting mergers and type Ia supernovae. In CEE, the primary and secondary cores inspiral from drag, transferring orbital energy to the envelope until the latter ejects, or the cores merge.

Hydrodynamic simulations of this process generally do not eject the envelope. Although it is possible that the cores should merge for the parameter regime explored in some of the simulations (Iaconi et al., 2018), in other simulations the rate of decay of the inter-particle separation $a$ decreases dramatically at values of $a$ too large for a merger during the computation. There may also be missing physics in the simulations. For example, most employ an ideal gas equation of state (EOS), whereas a more sophisticated EOS should account for ionization and recombination. When recombination energy is injected locally, the envelope is found to eject or almost eject in at least some cases (Nandez et al., 2015; Nandez & Ivanova, 2016; Prust & Chang, 2019).

Chamandy et al. (2019) (hereafter Paper II) applied the CE energy formalism (van den Heuvel, 1976; Webbink, 1984; Livio & Soker, 1988) to show that for a reasonable energy parameter $\alpha_{\mathrm{CE}}\lesssim 0.3$, theory correctly predicts that the envelope will not eject in our simulation or in any other with very similar initial conditions (Ohlmann et al., 2016) because the simulations do not reach the predicted separation for ejection by the end of the runs. That simulations do not eject the envelope because they do not attain small enough separations is partly supported by observations, exhibiting small final separations (Iaconi et al., 2017; Iaconi & De Marco, 2019).

Although extra energy sources (e.g. recombination energy or energy released by accretion onto the companion) may help to eject the envelope, they should also result in larger final separations, since less transfer of orbital energy is then required for ejection. On the other hand, energy sinks, such as loss via radiation, may offset energy gain by the envelope gas (Sabach et al. 2017; Grichener et al. 2018, but see Ivanova 2018).

Separations at late times tend to be overestimated because of inadequate numerical resolution (Ohlmann et al. 2016; Iaconi et al. 2017, 2018; Paper II), but this is unlikely a dominant effect–the slow decrease of $a$ at late times needs to be explained physically. Ricker & Taam (2008, 2012); Staff et al. (2016) and Iaconi et al. (2017, 2018) include some explorations of the drag force in their global CE simulations. And recent work by Reichardt et al. (2019) showed that the decrease in the rate of orbital tightening at late times was consistent with reduction in the drag force on the companion measured in one of their simulations. The reduction was explained qualitatively by a reduction in the angular velocity of the cores relative to the gas in their vicinity.

The goal of this work is to analyze the drag force in three otherwise identical simulations, but each with a different companion mass, and to compare our results with results from analytic theory and local wind tunnel CE simulations of flow near the secondary. In Sec. 2 we summarize our numerical methods. Sec. 3 contains the results of our simulations for the net force. We compare these results to analytic theory in Sec. 4. The evolution of the flow around the secondary, with a focus on the simulation with largest companion mass, is explored in Sec. 5. The results for the net force and flow properties are then compared to wind tunnel simulations in Sec. 6. We summarize and conclude in Sec. 8.

2 Simulation parameters and methods

We employ the hydrodynamics code astrobear, which includes adaptive mesh refinement (AMR). The primary is an $M_{\mathrm{1}}=1.96\,\mathrm{M_{\odot}}$ red giant branch (RGB) star with radius $R_{\mathrm{1}}=48\,\mathrm{M_{\odot}}$ and core mass $M_{\mathrm{1,c}}=0.37\,\mathrm{M_{\odot}}$, and the secondary has mass $M_{\mathrm{2}}=0.98\,\mathrm{M_{\odot}}$ (Model A), $M_{\mathrm{2}}=0.49\,\mathrm{M_{\odot}}$ (Model B) or $M_{\mathrm{2}}=0.245\,\mathrm{M_{\odot}}$ (Model C). The primary and secondary are initialized in a circular orbit with separation $a_{\mathrm{i}}=49\,\mathrm{R_{\odot}}$. Aside from the companion mass and initial velocities, the three runs are identical. Model A is the same as Model A of Chamandy et al. (2018) (hereafter Paper I) and Paper II.

RGB core and companion are modeled as point particles (“particle 1” and “particle 2”, respectively) that interact with each other and gas via gravity only. The particle potential is smoothed according to a spline function (Springel, 2010) such that it is Newtonian for $r>r_{\mathrm{soft}}$ and shallower than Newtonian for $r<r_{\mathrm{soft}}$, where $r_{\mathrm{soft}}$ is the spline softening radius. The RGB model is adapted from a mesa (Paxton et al., 2015) 1D profile using a similar method to that of Ohlmann et al. (2017) to model the gas profile within the softening radius. The spline softening radius and smallest resolution element are respectively $r_{\mathrm{soft}}=2.4\,\mathrm{R_{\odot}}$ and $\delta=0.14\,\mathrm{R_{\odot}}$ from $t=0$ to $t=16.7\,{\rm d}$, and $r_{\mathrm{soft}}=1.2\,\mathrm{R_{\odot}}$ and $\delta=0.07\,\mathrm{R_{\odot}}$ thereafter. Refinement at the highest resolution is applied everywhere within a sphere of dynamically changing radius $r_{\mathrm{refine}}$ (see Fig. 1), centred on the primary core before $t=16.7\,{\rm d}$, and companion thereafter. The simulation domain size is $L_{\mathrm{box}}=1150\,\mathrm{R_{\odot}}$, with $512^{3}$ base cells of size $2.25\,\mathrm{R_{\odot}}$. (Four levels of AMR are used for $t<16.7\,{\rm d}$ and five levels thereafter, and going up one level halves the cell size.). Extrapolation boundary conditions are employed.

An ideal gas EOS with $\gamma=5/3$ is employed. The ambient density and pressure are $\rho_{\mathrm{amb}}=6.7\times 10^{-9}\,{\rm g\,cm^{-3}}$ and $P_{\mathrm{amb}}=1.0\times 10^{5}\,{\rm dyn\,cm^{-2}}$. The simulations are stopped after $t=40\,{\rm d}$. More details about the setup and methods can be found in Papers I and II.

In Model B of Paper I, the secondary was a sink particle that accreted mass at a rate which was an upper bound to the true accretion rate. Since the orbit, and hence the drag force, were not drastically affected by this accretion, we exclude accretion onto the companion in the present simulations.

3 Overall Evolution

3.1 Orbital separation

Fig. 1 shows orbital separation versus time for Models A, B and C in solid blue, dashed red and dash-dotted black, respectively. The quantities $r_{\mathrm{soft}}$ and $r_{\mathrm{refine}}$, which do not change between runs, are also shown for reference. As the companion mass is lowered, the initial orbital speed and separation decay rate are both reduced. At later times however, the separation decays more rapidly for lower mass, and the curves cross. This behaviour is consistent with other studies (e.g. Passy et al., 2012).

We have computed the tidal shredding radius $r_{\mathrm{shred}}$ for a main sequence secondary using the initial density profile of the primary along with the estimate of Nordhaus & Blackman (2006) and the mass-radius relation from Eker et al. (2018), and find $r_{\mathrm{shred}}<1\,\mathrm{R_{\odot}}$ for all three models. For a white dwarf secondary, $r_{\mathrm{shred}}$ would be smaller still. Thus, the secondary is not expected to tidally shred during any of our simulation runs.

3.2 Drag Force

The centre of mass of the particles accelerates during the simulation (Paper II) due to the gravitational interaction between the gas and each of the particles. To facilitate comparison with theory and local simulations that treat the primary as fixed and non-rotating, we compute the dynamical friction force on particle 2 in the non-inertial rest frame of particle 1, $\bm{F}_{\mathrm{2-gas,1}}$, where the subscript ‘1’ after the comma denotes this reference frame. As seen in the lab frame, this frame orbits with particle 1 but does not rotate. This introduces a fictitious force so that the force exerted on particle 2 by gas in this frame is given by

[TABLE]

where the terms on the right are computed in the lab frame (nearly the centre of mass frame of the entire system; see Paper II). We have not included the terms $(1+M_{\mathrm{2}}/M_{\mathrm{1,c}})\bm{F}_{2-1}$ because these terms involve forces between the particles, and we are interested in computing the forces between gas and particles. Note that a force in the $-\phi$ direction (a drag) on particle 1 in the lab frame contributes a drag on particle 2 in the frame of particle 1. To compute the terms on the right of equation (1), we simply integrate the force per unit volume on each particle, for example: $\bm{F}_{\mathrm{2-gas}}=\mathrm{G}M_{\mathrm{2}}\sum_{V}\rho(\bm{s})[(\bm{s}-\bm{s}_{\mathrm{2}})/|\bm{s}-\bm{s}_{\mathrm{2}}|^{3}]\mathrm{d}^{3}s$, where $\rho(\bm{s})$ is the gas density at position $\bm{s}$, $V$ is the volume of the simulation domain, and $\bm{s}_{\mathrm{2}}$ is the position of particle 2. We then compute the $\phi$-component $(\bm{s}_{\mathrm{2}}-\bm{s}_{\mathrm{1}})\times\bm{F}_{\mathrm{2-gas,1}}/a\cdot\hat{\bm{z}}$. Likewise, we compute the projection of the force along the velocity vector relative to particle 1: $\bm{F}_{\mathrm{2-gas,1}}\cdot(\bm{v}_{\mathrm{2}}-\bm{v}_{\mathrm{1}})/|\bm{v}_{\mathrm{2}}-\bm{v}_{\mathrm{1}}|$.

The force on particle 2, multiplied by $-1$, is presented in Fig. 2. We refer to positive values on the plot as ‘drag’ and negative values as ‘thrust’. The $\phi$-component is plotted as a solid black line, and the projection along the relative velocity is plotted as a dash-triple-dotted gold line, for Model A (top), Model B (middle) and Model C (bottom). The separation $a$ is plotted with respect to the right axis for reference.111The inter-particle separation never differs from its projection in the $xy$-plane by more than $0.2\%$.

For all models, the $-\phi$-component of the force steadily increases from $t=0$, attains a broad peak of a few days width (with a double humped morphology, at least for Models A and B), and then reduces before becoming quasi-constant until the end of the simulation. The broad peak roughly coincides with the first periastron passage, though for Models B and C it happens slightly earlier. The peak magnitude is roughly proportional to the companion mass: $\sim 32$, $\sim 16$ and $\sim 8$ in units of $10^{33}\,{\rm dyn}$ for Models A, B and C, respectively.

Periodicity emerges at later times, particularly in Models A and B, with the force magnitude greatest (smallest) when $a$ is smallest (greatest). The evolution is slower in Model C, so we expect such variations to become more regular only after $t=40\,{\rm d}$. At late times, the magnitude is only weakly dependent on the companion mass, being $\sim 7\times 10^{33}\,{\rm dyn}$ for Models A and B and closer to $\sim 4\times 10^{33}\,{\rm dyn}$ for Model C, but has not yet stabilized by the end of that run.

The dotted grey curve shows the contribution to the solid black curve from only the the first term on the right of equation (1). Ignoring the fictitious force exerted by gas on particle 1 would thus lead to the wrong conclusion that the $\phi$-component of $\bm{F}_{\mathrm{2-gas,1}}$ is sometimes positive.

The component of force along $\bm{v}_{\mathrm{2}}-\bm{v}_{\mathrm{1}}$ has a period-averaged magnitude similar to the $\phi$-component at late times, but varies strongly with orbital separation and oscillates between drag and thrust. The oscillations occur because the gas force exerted on each particle is dominated by gas in the direction of the other particle, so particle approach (recession) produces thrust (drag).

To explore the effect of changing the softening length and resolution at $t=16.7\,{\rm d}$, we compared the drag force in Model A to that of Model F from Paper II, for which both $r_{\mathrm{soft}}$ and $\delta$ retain their initial values for the full $t=37.3\,{\rm d}$ simulation. The separation $a(t)$ for Model F differs only slightly from that of Model A, as does the force (Appendix A).

3.3 Torque and Orbital Energy Dissipation

The torque on the particles about their centre of mass is plotted on the left side of Fig. 3. The solid black line shows the $z$-component of the torque computed using the forces obtained by integration over the simulation domain,

[TABLE]

Here the forces are in the lab frame because torques from fictitious forces cancel.

The dash-triple-dotted magenta line shows the $z$-component of the rate of change of the particle angular momentum. This is obtained by first computing the angular momentum of the particles about the particle centre of mass (denoted ‘CM’)

[TABLE]

and then numerically time-differentiating $J_{z}$ (with sampling interval $\approx 0.23\,{\rm d}$ per frame). The two methods of computation should in principle yield identical results, except for the sampling error on $\dot{J}_{z}$, and the level of agreement is indeed excellent. The $xy$-plane components of the torque are negligible in both methods.

Given that the particle orbital energy dissipation rate could be used to estimate the observed luminosity of potential transient CE events such as luminous red novae, we also compute this. However, such an estimate would need to consider radiative transfer and is left for future work. The right side of Fig. 3 shows strong agreement between the two different methods, first from computing the rate of work done by gas on particles,

[TABLE]

(black solid line), and second by numerical time-differentiating the total particle energy

[TABLE]

(dash-triple-dotted magenta line).

4 Comparison to Analytic Theory

4.1 Estimate for Uniform Density

The dynamical friction force can be estimated from Bondi-Hoyle-Lyttleton (BHL) theory (Hoyle & Lyttleton, 1939; Bondi & Hoyle, 1944; Bondi, 1952). Here, gas approaching with impact parameter less than the accretion radius

[TABLE]

accretes onto the star, where $c_{\infty}$ and $v_{\infty}$ are the sound speed of the unperturbed envelope, and its speed relative to the secondary. The accretion rate can be estimated as $\dot{M}\sim\uppi R_{\mathrm{a}}^{2}\rho_{\infty}(c_{\infty}^{2}+v_{\infty}^{2})^{1/2}$, where $\rho_{\infty}$ is the unperturbed density, resulting in a drag force

[TABLE]

Typically, $r_{\mathrm{max}}$ is taken to be $R_{\mathrm{a}}$ and $r_{\mathrm{min}}$ as the radius of the star. Equation (7) was first derived by Dokuchaev (1964) and survives among different estimates (Edgar, 2004) subjected to refinements from numerical studies, e.g. Shima et al. (1985). We neglect turbulence (Krumholz et al., 2006) which may be important in general. We do consider the influence of a density gradient, as explained below.

To make contact with previous work, we plot the $\phi$-component of the drag force, as in Fig. 2, but now with additional lines representing theoretical predictions or results from local simulations, in Fig. 4. The dash-dotted red line shows the quantity

[TABLE]

where $v_{\mathrm{0}}\equiv|\bm{v}_{\mathrm{2}}-\bm{v}_{\mathrm{1}}|_{\mathrm{0}}$. In our notation, quantities with a ‘0’ subscript are computed from the initial envelope profile at radius $a(t)$, with velocity computed assuming a circular orbit with primary mass equal to the mass interior to the orbit, $m_{\mathrm{1}}(a)=M_{\mathrm{1,c}}+m_{\mathrm{env}}(a)$. We find that replacing $v_{\mathrm{0}}$ by its actual value measured in the simulation $|\bm{v}_{\mathrm{1}}-\bm{v}_{\mathrm{2}}|$ increases the amplitude of the oscillations in $F_{\mathrm{0}}$ but otherwise the results are similar, so we opt to use the relative velocity computed from the initial profile.222Defining the Mach numbers $\mathcal{M}=|\bm{v}_{\mathrm{2}}-\bm{v}_{\mathrm{1}}|/c_{\mathrm{0}}$ or $\mathcal{M}_{\mathrm{0}}=v_{\mathrm{0}}/c_{\mathrm{0}}$, these are found to be in the range $1.1<\mathcal{M},\mathcal{M}_{\mathrm{0}}<5.8$. Within this range, the values are larger for larger companion mass. Both peak at the same time, at about $6$, $7$ and $9\,{\rm d}$ for Models A, B and C, respectively, before leveling near to the minimum as the simulation progresses.

At early times, $F_{\mathrm{0}}$ is effectively zero due to the small $\rho_{\mathrm{0}}$. Subsequently, $F_{\mathrm{0}}$ rises to be comparable to $\phi$-component of $-\bm{F}_{\mathrm{2-gas,1}}$ just before the first periastron passage, before continuing to rise, in contrast to the $\phi$-component of $-\bm{F}_{\mathrm{2-gas,1}}$, which decreases and then levels off. Hence, equation (8) overestimates the magnitude of the drag force at late times. Qualitatively similar results were obtained by Staff et al. (2016). However, the phase and amplitude of the variability seen in the $\phi$-component of $-\bm{F}_{\mathrm{2-gas,1}}$ at late times in Models A and B is reproduced by this minimalist theoretical model.

4.2 Estimate Including Density Gradient

Using a more refined version of the theory that includes the logarithmic factor of equation (7) and accounts to some extent for gradients in the radial direction might be expected to produce better agreement. Hence, for a more general estimate to the drag force, we multiply $F_{\mathrm{0}}$ by $\ln(r_{\mathrm{max}}/r_{\mathrm{min}})$ and a correction factor Dodd & McCrea (1952) (hereafter DM) $R_{\mathrm{a,DM}}^{2}/R_{\mathrm{a,0}}^{2}$, that accounts for a linear or at most quadratic density gradient in the correction to the accretion radius, computed by

[TABLE]

where

[TABLE]

and $H_{\rho}=-\rho_{\mathrm{0}}/(\mathrm{d}\rho/\mathrm{d}r)_{\mathrm{0}}$ is the scale height. The modified force magnitude is then

[TABLE]

We adopt $r_{\mathrm{max}}=R_{\mathrm{a,DM}}$ and $r_{\mathrm{min}}=r_{\mathrm{soft}}|_{t=0}=2.4\,\mathrm{R_{\odot}}$. The dashed purple line shows the resulting corrected estimate.333Radial variations in the density gradient of the initial envelope profile cause the noise. This variation is present in the 1D mesa solution, which was retained for our initial condition outside of $r=2.4\,\mathrm{R_{\odot}}$.

We see from Fig. 4 that despite some differences, the level of agreement between theory which includes the density gradient (dashed purple) and the simulation results (solid black) is overall comparable to that obtained using $F_{\mathrm{0}}$ (dashed red). The DM correction marginally improves agreement for Models A and B, but marginally worsens agreement for Model C.444We also tried other variations, with only the $\ln(r_{\mathrm{max}}/r_{\mathrm{min}})$ factor or only the DM correction included, and found results that are generally similar to the cases plotted.

4.3 Theory works best at intermediate times

We do not expect good agreement between simulation and theory at early times because of three interrelated factors: (i) tidally drawn envelope material increases the density near particle 2 beyond ambient and outer envelope layer values, (ii) density scale heights are initially small compared to the accretion radius and (iii) the initial condition of a secondary placed just outside the spherically symmetric primary at $t=0$ is not fully realistic. Point (ii) is seen by comparing solid black and green lines in Fig. 5, where we plot the various length scales as a function of time, for each simulation. The DM correction in principle helps to account for (ii) but considers only the lowest order effect of the density gradient.

At late times, we also expect poor agreement. In Model A (q=1/2) $R_{\mathrm{a}}\sim a$ shortly after the first periastron passage, as seen in Fig. 5 (compare blue and red dashed lines with black solid line), so we do not expect good agreement. For Model B $(q=1/4)$, $R_{\mathrm{a}}$ remains marginally smaller than $a$, while for Model C (q=1/8), $R_{\mathrm{a}}\sim 0.5a$ by the end of the simulation. Theoretical predictions for late times improve slightly as $q$ decreases to $1/8$, but not dramatically.

At intermediate times, when $R_{\mathrm{a}}\ll a$ and $H_{\rho,0}\gtrsim R_{\mathrm{a}}$, we expect and find agreement to be much better. If we use $|\bm{v}_{\mathrm{2}}-\bm{v}_{\mathrm{1}}|$ measured directly from the simulation to compute $R_{\mathrm{a}}$ (dashed blue in Fig. 5), then the time range when $R_{\mathrm{a}}(t)<a(t)$ at the previous periastron passage becomes $12\,{\rm d}\lesssim t\lesssim 14\,{\rm d}$ for Model A, $13\,{\rm d}\lesssim t\lesssim 24\,{\rm d}$ for Model B, and $t\gtrsim 16\,{\rm d}$ for Model C. BHL/DM theory approximates the numerical results reasonably well in these time ranges. In particular, during the broad force peak, theory (equation 8 or 11) correctly predicts the force to within a factor of $\sim 2$ for all models.

4.4 Improved theory is needed for late times

Reichardt et al. (2019) suggested that a reduction in the relative velocity between the particles and gas at late times in their simulation with $q=0.68$ might help explain the reduction of the drag force. Might replacing our initial values of the gas density, velocity with respect to particle 2, and sound speed with those measured directly from the simulation help reconcile theory and simulation at late times?

From the right column of Fig. 6 (see Sec. 5 for details), we can estimate the parameters $\rho_{\infty}$, $v_{\infty}$ and $c_{\infty}$ at $t=22.0\,{\rm d}$. Since the orbital separation is small, one option is to choose a location in the vicinity of particle 2, namely in the region below and to the left of particle 2 and within a few $\!\,\mathrm{R_{\odot}}$ of particle 2 in the plots. The velocity vectors show that some of this gas is being sling-shotted around particle 2 and this produces a drag force. Here, $\rho_{\infty}\sim(5$–$10)\rho_{\mathrm{0}}$, $v_{\infty}\sim\tfrac{1}{2}v_{\mathrm{0}}$, and $c_{\infty}\sim 2c_{\mathrm{0}}$. This leads to a force $F_{\infty}\sim(1.9$–$3.8)F_{\mathrm{0}}$, using the analogue of equation (8), so the predicted force is even larger than before, instead of smaller, as needed. Perhaps more reasonable is to choose a location farther away from particle 2 on the circle centred on particle 1 with radius equal to $a$, near coordinates $(-10,-5)$ in the right column of Fig. 6. In this vicinity, we find $\rho_{\infty}\sim 3\rho_{\mathrm{0}}$, $v_{\infty}\sim 0.3v_{\mathrm{0}}$, and $c_{\infty}\sim c_{\mathrm{0}}$, giving $F_{\infty}\sim 4.9F_{\mathrm{0}}$. Pushing the values as far as seems reasonable to obtain a smaller estimate for the force, we could instead choose $\rho_{\infty}\sim\rho_{\mathrm{0}}$, $v_{\infty}\sim\tfrac{1}{2}v_{\mathrm{0}}$ and $c_{\infty}\sim 2c_{\mathrm{0}}$, which leads to $F_{\infty}\sim 0.38F_{\mathrm{0}}$. This value is still not small enough to explain the factor of $\sim 10$ between the solid black and dash-dotted red lines in the top panel of Fig. 4 at $t=22.0\,{\rm d}$.

A second possibility is to note that the conventional theoretical drag force formula changes as $\mathcal{M}_{\infty}$ becomes small (Ostriker, 1999). Could this account for the small values seen at late times, as suggested by Staff et al. (2016) for their simulations? In the subsonic regime Ostriker (1999) derives the solution

[TABLE]

where

[TABLE]

This leads to the factor $[(1+\mathcal{M}_{\infty}^{2})^{3/2}/\mathcal{M}_{\infty}^{3}]I\approx 0.48$ multiplying our above estimate of the force if $v_{\infty}=\tfrac{1}{2}v_{\mathrm{0}}$ and $c_{\infty}=2c_{\mathrm{0}}$, or $\approx 0.56$ if $V_{\infty}=0.3v_{\mathrm{0}}$ and $c_{\infty}=c_{\mathrm{0}}$ are used instead. Thus, we obtain the new estimates $F_{\infty}\sim(1$–$2)F_{\mathrm{0}}$, $3F_{\mathrm{0}}$ and $0.18F_{\mathrm{0}}$ for the three cases described above. Even the smallest of these is about a factor of two too large to explain the force measured in the simulation. We obtain better agreement at $t=40\,{\rm d}$. However, it is not clear whether equation (12) is even applicable in the present context.

Thus, such theoretical estimates are not well-motivated for late times owing to the small inter-particle separation, are sensitive to arbitrary choices, and cannot reproduce the force measured from the simulation.

Another possibility is to note that at late times the particles accrete their own quasi-static, quasi-spherical “bulges” of gas (Paper I). The bulge around particle 2 is mainly pressure supported but partially rotation supported, and has size $r_{\mathrm{b}}\sim 2r_{\mathrm{soft}}=2.4\,\mathrm{R_{\odot}}$. Although particles interact with gas only through gravity, and thus can only experience dynamical drag, the composite particle-bulge “system” in addition experiences a hydrodynamic drag as it moves through the surrounding envelope gas. The hydrodynamic drag force can be estimated as $F_{\mathrm{h}}\sim\rho v^{2}\uppi r_{\mathrm{b}}^{2}$, where $\rho v^{2}$ is the ram pressure exerted by the gas encountered by the bulge around particle 2 as it orbits. For Model A at $t=22.0\,{\rm d}$, using the above estimate $\rho\sim(5$–$10)\rho_{\mathrm{0}}$ and $v\sim 0.5v_{\mathrm{0}}$. gives $F_{\mathrm{h}}\sim(0.5$–$1)\times 10^{34}\,{\rm dyn}$, which is just the order of magnitude needed. However, this formula underestimates the drag force at $t=40\,{\rm d}$, and hence is inadequate.

We leave further exploration of the force at late times for future study.

5 Evolution of flow properties

Fig. 6, shows snapshots of various quantities, sliced through $z=0$, for Model A. The snapshots are taken at $t=6.9$, $11.1$, $16.7$ and $22.0\,{\rm d}$. Defining the end of the dynamical plunge-in phase as the time of first periastron passage (Paper II) the four times represent, roughly speaking, the flow near the beginning of plunge-in, the end of plunge-in, the transition to slow spiral-in, and at the beginning of slow spiral-in. The orbital motion is counter-clockwise. In each panel, particle 2 is at the centre and the view is rotated so that particle 1 is on the negative $y$-axis. Circles show the softening spheres around the particles. We focus on Model A because it displays the strongest deviation from theory, and because other aspects of the run were extensively studied in Papers I and II.

5.1 Force density

The top row shows the magnitude of the $\phi$-component of the force per unit volume exerted by gas on particle 2 in the co-orbiting but non-rotating rest frame of particle 1. Positive (negative) drag contributions are indicated by solid (dashed) contours, spaced by the values on the color bar. As there are positive and negative contributions from both terms in equation (1), the plot contains four sets of contours. Forces between gas and particle 2 dominate the contours in the upper part of the plot, while forces between gas and particle 1 (fictitious forces on particle 2) dominate the lower sets of contours. A drag (thrust) on particle 1 in the lab frame produces a fictitious drag (thrust) on particle 2 in the reference frame orbiting with particle 1. The black arrow shows the relative magnitude and direction of $\bm{F}_{\mathrm{2-gas,1}}$, while the blue arrow shows the same for the velocity of particle 2 in the same reference frame.

At $t=6.9\,{\rm d}$ (column 1), a low density tidal tail from the primary wraps around the secondary from behind it (second row), providing a dynamical friction force. The density of the gas in front of particle 2 is much smaller so the force pulling particle 2 backward dominates. The contribution from the force on particle 1 is small because this contribution is dominated by the gas near particle 1, which is distributed symmetrically in the trailing and leading directions.

At $t=11.1\,{\rm d}$ (column 2), the gas near particle 2 now has a higher density, leading to a greater drag force. If logarithmic intervals in radius with respect to particle 2 contributed equally to the net force, as suggested by equation (7), the force per unit volume contours would be separated by a factor $\sim 10^{1/3}\approx 2.15$ in radius.

This is roughly valid out to about the third contour: in each of the snapshots, the ratio in radius between the innermost contours is $\sim 2$, but the ratio between adjacent contours decreases as one moves outward from particle 2. The ratio of radii of the fourth to third contour is generally $\sim 1.5$, implying a decrease in contribution by a factor $\sim 1.5^{3}/10\approx 0.3$.

Calculating $R_{\mathrm{a}}$ we obtain $28$–$38\,\mathrm{R_{\odot}}$ for $t=6.9\,{\rm d}$ and $15$–$16\,\mathrm{R_{\odot}}$ for $t=11.1\,{\rm d}$; the first value uses equation (6) with $c_{\infty}$ and $v_{\infty}$ replaced by $c_{\mathrm{0}}$ and $v_{\mathrm{0}}$, while the second value uses $|\bm{v}_{\mathrm{2}}-\bm{v}_{\mathrm{1}}|$ directly from the simulation. These values correspond to roughly the third contour, so adopting $r_{\mathrm{max}}=R_{\mathrm{a}}$ is reasonable; the contribution to the force beyond this radius is small.

By the next snapshot, at $t=16.7\,{\rm d}$ (column 3), the force density about particle 2 has become much more symmetric, and the drag force is approximately zero (consistent with the top-left panel of Fig. 3). Correspondingly, the contours show a high degree of left-right symmetry. The force $\bm{F}_{\mathrm{2-gas,1}}$ now has comparable contributions from $\bm{F}_{\mathrm{2-gas}}$ and $-(M_{\mathrm{2}}/M_{\mathrm{1,c}})\bm{F}_{\mathrm{1-gas}}$.

For the final snapshot at $t=22.0\,{\rm d}$ (column 4), the $\phi$-component of $\bm{F}_{\mathrm{2-gas,1}}$ is quasi-steady. The force remains small, owing to the high degree of symmetry in the force density, so that the thrust and drag contributions balance except for a small net drag. The force density pattern is quasi-steady thereafter, consistent with the net force being quasi-steady.

5.2 Gas density evolves toward symmetry

The second row of Fig. 6 shows the gas density normalized to $\rho_{\mathrm{0}}(a)$ of the initial density profile of the RGB star. Vectors show the gas velocity in the frame co-orbiting and co-rotating with particle 2, as this is the appropriate reference frame for comparison with theory and local simulations. The second to fifth rows are zoomed in by a factor of four compared to the top row. At $t=6.9\,{\rm d}$, the density is $\sim 10^{3}$ times larger than $\rho_{\mathrm{0}}(a)$ as deeper gas is pulled tidally around particle 2. The value of $\rho(a)/\rho_{\mathrm{0}}(a)$ then reduces to $\sim 10$ at $t=22.0\,{\rm d}$ and $\sim 1$ by $t=40\,{\rm d}$. The flow around particle 2 becomes more axisymmetric, rotating $\sim 20\%$ of the Keplerian value at late times (Paper I).

5.3 Mach number and turbulence

In the third row, we plot the gas Mach number computed in the frame co-orbiting and co-rotating with particle 2. In this frame, the gas around particle 2 is mostly supersonic during plunge-in, except in the bow shock. For these snapshots we obtain, from earliest to latest, $\mathcal{M}_{\mathrm{0}}=v_{\mathrm{0}}/c_{\mathrm{0}}=5.3$, $2.1$, $1.7$, and $1.7$. At $t=16.7\,{\rm d}$, a shock is still seen, now above particle 2 in the plot. By $t=22.0\,{\rm d}$, the gas about the particles and within the orbit is not only subsonic, but turbulent. The turbulence is well-developed by $t=18\,{\rm d}$.

To verify that the turbulence is not produced by the sudden change in softening length and resolution at $t=16.7\,{\rm d}$, we compared the 2D snapshots to a run (Model F of Paper II) for which the softening length and maximum AMR level are not changed in this way. Turbulent eddies are conspicuous by $t=18\,{\rm d}$ in that run as well, even though the smallest scales of the turbulence are larger.

The onset of turbulence roughly coincides with the transition from plunge-in to slow spiral-in once $R_{\mathrm{a}}\sim a$, as shown in Fig. 5, and the particles have completed a full orbit since the first periastron passage. The gas they encounter no longer moves out supersonically due to the deeper potential and confinement by overlying layers (Paper II) so is continually being ‘reprocessed’.

5.4 Azimuthal velocity and sound speed

Finally, we present plots of the velocity and sound speed, which were used above to obtain modified parameter values for theoretical estimates (Sec. 4.4). The vectors in the fourth row of Fig. 6 show the gas velocity in the corotating frame of particle 2 as in the third row. The color shows the $-\phi$-component of this velocity with respect to particle 1, normalized to $v_{\mathrm{0}}$. A value of unity (orange) corresponds to the relative tangential velocity estimated from the initial stationary envelope profile, while negative values (yellow to purple) denote oppositely moving gas. In all snapshots, the streamlines curl counter-clockwise around particle 2 so that the $\phi$-component of the velocity reverses sign. The speed of gas approaching particle 2 is of order $v_{\mathrm{0}}$ in the first three snapshots, but only about $\tfrac{1}{2}v_{\mathrm{0}}$ by the fourth snapshot (Sec. 4.4).

Finally, in the fifth row we plot the sound speed normalized to $c_{\mathrm{0}}$, along with the same velocity vectors plotted in rows 2 and 4. The gas flowing toward and deflected by the secondary has sound speed $\sim 2c_{\mathrm{0}}$ (blue).

6 Comparison to Wind Tunnel Simulations

6.1 Drag force comparison

Here we compare our results with the local CE wind tunnel simulations of MacLeod et al. (2017) (hereafter M17) for which a particle representing the secondary was fixed at the centre of the grid and a wind was launched from the $-x$ boundary with a prescribed $x$-velocity and density gradient in the $-y$ direction. By approximating the gas as a polytrope and assuming that the upstream wind velocity equals the local Keplerian orbital speed, the upstream Mach number is determined once the dimensionless density gradient parameter

[TABLE]

and the mass ratio $q_{\mathrm{enc}}=M_{\mathrm{2}}/m_{\mathrm{1}}(a)$ are specified. M17 chose $q_{\mathrm{enc}}=0.1$ and explored the dependence on $\epsilon_{\rho}$. In our simulations, $q_{\mathrm{enc}}=q=1/2$, $1/4$ or $1/8$ at $t=0$, but then increases with time as $m_{\mathrm{1}}(a)$ decreases: computing $m_{1,0}(a)$ from the initial envelope profile, we obtain $q_{\mathrm{enc,0}}=2.0$, $1.2$ and $0.6$ at $t=40\,{\rm d}$ for Models A, B and C, respectively. M17 results are likely to be sensitive to their fixed choice of $q_{\mathrm{enc}}$; nevertheless we proceed with the comparison for Model C, which proves to be fruitful.

We first fit the drag force in M17 for $\gamma=5/3$ in their Fig. 10. Replacing $\rho_{\infty}$ and $v_{\infty}$ by $\rho_{\mathrm{0}}$ and $v_{\mathrm{0}}$ we obtain

[TABLE]

We apply this for $0.4<\epsilon_{\rho}<2.7$ (consistent with the range of parameter space explored by M17) and plot the resulting force only for times in our simulation when $\epsilon_{\rho,0}=2\mathrm{G}M_{\mathrm{2}}/(v_{\mathrm{0}}^{2}H_{\rho,0})$ is within this range. The result is the dashed green line in the bottom panel of Fig. 4.

At intermediate times, between $t=12$–$26\,{\rm d}$, the agreement is excellent. At $t=26\,{\rm d}$, when $-\dot{a}$ attains its second maximum (the first having occurred between $t=13$–$14\,{\rm d}$), the $\phi$-component of $-\bm{F}_{\mathrm{2-gas,1}}$ decreases from its peak value, but equation (15) predicts the force to continue rising. At this time, $q_{\mathrm{enc,0}}=0.28$, or almost three times larger than that assumed by M17, which likely contributes to this discrepancy.

6.2 Flow structure comparison

Flow structure of global and local simulations can also be compared. For Model C, we choose the time $t=20.8\,{\rm d}$, at which $\epsilon_{\rho}=0.80$ (and $q_{\mathrm{enc}}=0.15$), to compare with the lower left panels of Fig. 2 of M17. We plot the mass density normalized to $\rho_{\mathrm{0}}(a)$ and velocity vectors in the corotating frame of particle 2 in the top panel of Fig. 7, and the Mach number in the corotating frame of particle 2 in the bottom panel. The unit of the M17 axes is $2\mathrm{G}M_{\mathrm{2}}/v_{\mathrm{0}}^{2}=7.4\,\mathrm{R_{\odot}}$, so our plotting region is slightly larger than theirs. The level of agreement is remarkable. All of this suggests that the local simulations approximate global simulations for this window of parameter space.

Curiously, there is also some correspondence between the flow in Model A and that of the local simulations, even though the mass ratio of the former is much larger. Whereas M17 adopted $q_{\mathrm{enc}}=0.1$ and $\epsilon_{\rho}=0.2$–$2$, Model A has $q_{\mathrm{enc,0}}=0.50$, $0.75$, $1.23$ and $1.53$ and $\epsilon_{\rho,0}=20.3$, $2.1$, $1.7$ and $1.4$ at $t=6.9$, $11.1$, $16.7$ and $22.0\,{\rm d}$, respectively. Despite differences in $q_{\mathrm{enc}}$ and the force seen in Fig. 4, the panels of Fig. 6 showing $\rho/\rho_{\mathrm{0}}$ and $\mathcal{M}$ in Model A at $11.1\,{\rm d}$ show similarities to those of $\epsilon_{\rho}=2.00$ in Fig. 2 of M17, as seen by comparing the second column, third and fourth rows of Fig. 6 with Fig. 2 of M17. The flow pattern is similar and both methods exhibit a thin spiral shock. However, the normalized gas density in Fig. 6 is $1$–$1000$ (whereas in Fig. 7 we used $0.1$–$100$, as in M17). Thus in Model A our normalized densities are almost an order of magnitude larger than those of M17, likely because our $q_{\mathrm{enc}}$ is $7.5$ times larger than theirs at that time. The size of the region plotted in units of $2\mathrm{G}M_{\mathrm{2}}/v_{\infty}^{2}$ differs from M17: for Model A we obtain $2\mathrm{G}M_{\mathrm{2}}/v_{\mathrm{0}}^{2}=29$, $19$, $16$ and $14\,\mathrm{R_{\odot}}$, respectively, for the four snapshots.555 When $\epsilon_{\rho}$ and $2\mathrm{G}M_{\mathrm{2}}/v_{\infty}^{2}$ are estimated using the actual velocity $\bm{v}_{\mathrm{2}}-\bm{v}_{\mathrm{1}}$, rather than $v_{\mathrm{0}}$, the values are larger by $37\%$ for $t=6.9\,{\rm d}$ but hardly differ for the other snapshots of Model A.

As expected, snapshots at other times hardly resemble those of M17. At $t=6.9\,{\rm d}$, $\epsilon_{\rho}$ is an order of magnitude larger than that explored by M17. We do see increasing density contrast and larger rotation angle of the bow shock with increasing $\epsilon_{\rho}$, as in M17, but our shock is thick and morphologically complex. By $t=16.7\,{\rm d}$, $R_{\mathrm{a}}$ has already become comparable to $a$, as shown in the top panel of Fig. 5. The assumptions of M17, namely that (i) the envelope gas encountered by the secondary had not been previously affected; (ii) their $\rho_{\infty}$ smoothly and monotonically decreases with distance from the RGB core, and (iii) the gravity force from the RGB core can be approximated as everywhere downward, are no longer valid. Moreover, in our final snapshot, turbulence likely affects the dynamics.

Thus, we would not expect wind tunnel simulations to approximate the results of Model A at late times even if several wind tunnel simulations of different $q_{\mathrm{enc,0}}$ and $\epsilon_{\rho,0}$ were patched together to accommodate dynamically changing values of these parameters in the global simulation. However, given the excellent agreement at intermediate times for Model C, it would be interesting to compare local and global simulations using such dynamical patching of the local simulations to refine the temporal range over which this approach could be useful and computationally efficient.

7 How Important is Radiative Transfer?

Radiative transfer is neglected in our simulations (and in virtually all global CE simulations to date). A large diffusive flux of radiation out of the central region might lead to a different flow structure at late times. Averting the buildup of thermal energy might allow the flow there to retain a structure closer to that which it had originally, and thus closer to that assumed in wind tunnel experiments, where the flow is assumed to be unaffected by previous orbital passages of the particles.

Thus, we estimate the diffusion time to determine whether cooling would be significant at late times in our simulations. We consider the flow properties at $t=22.0\,{\rm d}$ in Model A (right column of Fig. 6), focussing on the region around particle 2 with $\mathcal{M}<1$. The distance from particle 2 to the boundary of this region is estimated as $R\sim 7\,\mathrm{R_{\odot}}$. Within this region, a typical gas density is $\rho\sim 2\times 10^{-4}\,{\rm g\,cm^{-3}}$ and the temperature $T>10^{6}\,{\rm K}$. At this temperature, hydrogen gas is ionized and the opacity is dominated by electron scattering, with cross section $\sigma_{\mathrm{T}}=6.65\times 10^{-25}\,{\rm cm}^{2}$.

The diffusion time is estimated by multiplying the number of scatterings by the mean free path and dividing by the speed of light in vacuum $c$. The mean free path is given by $l=(n_{\mathrm{e}}\sigma_{\mathrm{T}})^{-1}$, the electron number density by $n_{\mathrm{e}}=\rho/m_{\mathrm{H}}$ where $m_{\mathrm{H}}$ is the mass of the hydrogen atom, and the number of scatterings is given by $N\simeq R^{2}/l^{2}$. Thus, we have

[TABLE]

or $t_{\mathrm{d}}\sim 20\,{\rm yr}$ for Model A (with corresponding optical depth $\tau\sim\sigma_{\mathrm{T}}n_{\mathrm{e}}R\sim 40$). For Models B and C, at a comparable time in the evolution, $R$ is somewhat smaller while $\rho$ is slightly larger than in Model A, and the value of $t_{\mathrm{d}}$ is of the same order of magnitude. Since $t_{\mathrm{d}}\gg t$, the neglect of radiative transfer in the $\mathcal{M}<1$ region around the particles is justified.

8 Conclusions

We computed the drag force in three global CE simulation runs of $40\,{\rm d}$ in which a companion point particle is placed in circular orbit around a 2$\,\mathrm{M_{\odot}}$ RGB star. The runs are identical except for the value of the companion mass, $\,\mathrm{M_{\odot}}$, $\tfrac{1}{2}\,\mathrm{M_{\odot}}$ or $\tfrac{1}{4}\,\mathrm{M_{\odot}}$. We found that:

• •

The drag force on the particles at late times, during the slow spiral-in phase, has mean magnitude $\sim 7\times 10^{33}\,{\rm dyn}$, depending only weakly on companion mass, and varies periodically with the orbit (Figs. 2, 8).

• •

BHL/DM theory overestimates the drag force at late times by at least an order of magnitude for the run with initial mass ratio $q=1/2$ (Fig. 4 top panel), and cannot reproduce the late time force for any of the three runs.

• •

BHL/DM theory and local wind tunnel simulations are particularly inapplicable at late times for large $q_{\mathrm{enc}}=M_{\mathrm{2}}/m_{\mathrm{1}}(a)$ because the accretion radius becomes comparable to the inter-particle separation. The gas encountered by the particles forms a turbulent, thermalized, highly symmetric region around the particles (Fig. 6 rightmost column). Hydrodynamic drag may even dominate over dynamical friction during this phase, but further work is needed.

• •

At earlier times, the drag force peaks at or just before the first periastron passage with value approximately proportional to the companion mass (Fig. 2). Near this peak, the drag force is reasonably well matched by BHL/DM theory and particularly well matched by local wind tunnel simulations (Fig. 4 bottom panel), which also reproduce various features of the 2D slices at that time (c.f. Fig. 7 of this work and Fig. 2 of M17).

Thus, for low $q_{\mathrm{enc}}$, BHL/DM theory and local wind tunnel simulations approximate the drag in global simulations during the intermediate plunge-in phase, but not before or after. Since $q_{\mathrm{enc}}$ evolves temporally in global simulations, different fixed $q_{\mathrm{enc}}$ wind tunnel simulations must be patched together to increase the fidelity of comparison with global simulations over a larger temporal range. This has not yet been done.

Finally, more general theoretical approaches are needed to account for the high degree of symmetry and turbulence in the flow once $R_{\mathrm{a}}\sim a$, and the associated reduced drag at late times. This reduced drag dramatically slows the inward evolution and explains why numerous CE simulations do not reach tight enough orbits by the end of runs to eject the CE envelope.

Acknowledgements

We thank Orsola De Marco, Rosa Everson, Robert Fisher, Morgan MacLeod, Jason Nordhaus, Enrico Ramirez-Ruiz, and Jan Staff for useful discussions. We are grateful to the referee for a helpful report. EB acknowledges the Aspen Center for Physics, supported by NSF grant PHY-160761. This work used the computational and visualization resources in the Center for Integrated Research Computing (CIRC) at the University of Rochester and the computational resources of the Texas Advanced Computing Center (TACC) at The University of Texas at Austin, provided through allocation TG-AST120060 from the Extreme Science and Engineering Discovery Environment (XSEDE) (Towns et al., 2014), which is supported by National Science Foundation grant number ACI-1548562. Financial support for this project was provided by the Department of Energy grant DE-SC0001063, the National Science Foundation grants AST-1515648 and AST-181329, and the Space Telescope Science Institute grant HST-AR-12832.01-A.

Appendix A Effect of changing softening length and resolution

Here we compare the $-\phi$-component of the force exerted on particle 2 by the gas in the frame of particle 1 in Model A and Model F of Paper II. Model F restarts from Model A at $t=16.7\,{\rm d}$ but the softening radius and smallest resolution element are not halved as in Model A. The evolution of the force is very similar, confirming that the halving of $r_{\mathrm{soft}}$ and $\delta$ does not importantly affect the overall evolution of the force.