3.1. Uniform density lagging prolate spheroid model

The gravitational potential at a point (x, y, z) inside an ellipsoid of uniform density $\overline{{\unicode{x03C1}}}$ is given by

(4) \begin{equation} \Phi = {\unicode{x03C0}} G \overline{{\unicode{x03C1}}} \left({\unicode{x03B1}} x^2 + {\unicode{x03B2}} y^2 + {\unicode{x03B6}} z^2 - {\unicode{x03C7}} \right) \! ,\end{equation}

where the coefficients ${\unicode{x03B1}}$ , ${\unicode{x03B2}}$ , ${\unicode{x03B6}}$ , and ${\unicode{x03C7}}$ are given in Lamb (Reference Lamb1879), pp. 146–154, 700, Chandrasekhar (Reference Chandrasekhar1969), pp. 38–45, and Binney & Tremaine (Reference Binney and Tremaine2008), pp. 83–95. Take the origin to be the particles’ orbital centre of mass (CM) in the simulation and the x-axis to be the axis passing through both particles. For the special case of a constant-density prolate spheroid with major axis lagging this binary axis by the angle $\Delta\unicode{x03D5}$ , the z-component of torque at any point inside the spheroid is given by

(5) \begin{equation} \begin{aligned} {\unicode{x03C4}}_z = {\unicode{x03C0}} G \overline{{\unicode{x03C1}}} m &\big\{\sin \left(2\Delta \unicode{x03D5}\right) \left({\unicode{x03B2}} - {\unicode{x03B1}}\right)\left(y^2-x^2\right) \\ &+ 2xy\left[{\unicode{x03B1}}\cos\left(2\Delta\unicode{x03D5}\right) - {\unicode{x03B2}}\cos\left(2\Delta\unicode{x03D5}\right)\right] \big\}, \end{aligned}\end{equation}

where

(6) \begin{align} m = M_{\textrm{1,c}} + M_2 = 1.068\,{\textrm{M}}_{\odot},\end{align}

(7) \begin{align} {\unicode{x03B1}} = \frac{1-e^2}{e^3} \ln \frac{1+e}{1-e}-2 \frac{1-e^2}{e^2}, \end{align}

(8) \begin{align} {\unicode{x03B2}} = \frac{1}{e^2}-\frac{1-e^2}{2 e^3} \ln \frac{1+e}{1-e}, \end{align}

(9) \begin{align} e =\sqrt{1-\widetilde{B}^2 / \widetilde{A}^2}, \end{align}

$\widetilde{A}$ and $\widetilde{B}$ are the semi-major and semi-minor axes of the spheroid, respectively, and e is the spheroid eccentricity. As a consequence of Newton’s third theorem (Binney & Tremaine Reference Binney and Tremaine2008, p. 87), the ellipsoid’s mass outside the similar, concentric ellipsoid (of the same orientation) that passes through the particles exerts no net force on the binary. Therefore, the torque does not depend directly on the size of the ellipsoid. However, below we obtain $\overline{{\unicode{x03C1}}}$ by averaging the density inside an ellipsoidal contour, so the overall size of the contour does play a role. In our choice of coordinate system, the expression for the aggregate torque on both particles reduces to (Appendix A)

(10) \begin{equation} {\unicode{x03C4}}_z = - \frac{1}{2} {\unicode{x03C0}} G m \left( {\unicode{x03B2}} - {\unicode{x03B1}} \right) \cos (\Delta \unicode{x03D5}) \sin (\Delta \unicode{x03D5}) \overline{{\unicode{x03C1}}} a^2,\end{equation}

where a is the orbital separation.

3.2. Simplified double-perturber model

Another approach is to solve the governing equations after making various approximations. Reference Kim, Kim and Sánchez-SalcedoK08 solved numerically a linearized version of the hydrodynamics equations for two perturbers in a circular orbit, building on a study of a single perturber in a circular orbit, Kim & Kim (Reference Kim and Kim2007), which in turn extended the work of Ostriker (Reference Ostriker1999) for rectilinear motion. The Reference Kim, Kim and Sánchez-SalcedoK08 model is valid if the perturbers are weak and if various complicating effects – gas self-gravity, orbital eccentricity, the density gradient, and rotation of the gas in which the perturbers are embedded, the motion of the binary CM, and the time dependence due to orbital evolution – can all be safely neglected. Their model further assumes an ideal gas with adiabatic index $\gamma=5/3$ , which is also assumed in our CE simulation.

Below we restrict our discussion to the case of equal mass perturbers. While Reference Kim, Kim and Sánchez-SalcedoK08 also focuses on this case, we note that their general model includes the mass ratio as a parameter. The force on each perturber can be written as

(11) \begin{equation} F = -\mathcal{F}(\mathcal{I}_R{\hat{\boldsymbol{R}}} + \mathcal{I}_{\unicode{x03D5}}{ {\hat{\unicode{x1D6DF}}}}), \end{equation}

where $\mathcal{I}_R{\hat{\boldsymbol{R}}}$ and $\mathcal{I}_{\unicode{x03D5}}{\hat{\unicode{x1D6DF}}} $ characterize the radial and azimuthal components in the orbital plane,

(12) \begin{equation} \mathcal{F} \equiv \frac{4{\unicode{x03C0}}{\unicode{x03C1}}_0 (GM_{\textrm{p}})^2}{V_{\textrm{p}}^2},\end{equation}

$M_{\textrm{p}}$ and $V_{\textrm{p}}$ are the mass and speed of a perturber, and ${\unicode{x03C1}}_0$ is the unperturbed gas density. Further, the drag force is exerted both by the wake of the perturber and by the wake of its companion, so we can write

(13) \begin{equation} \mathcal{I}_R = \mathcal{I}'_{R} + \mathcal{I}''_{R}\end{equation}

and

(14) \begin{equation} \mathcal{I}_{\unicode{x03D5}} = \mathcal{I}'_{\unicode{x03D5}} + \mathcal{I}''_{\unicode{x03D5}},\end{equation}

where the first term, denoted by prime, is due to the perturber’s own wake and the second term, denoted by double-prime, is due to the wake of the other perturber.

Reference Kim, Kim and Sánchez-SalcedoK08 (see also Kim & Kim Reference Kim and Kim2007) solved numerically for the drag force exerted on the perturbers by their wakes and obtained fitting formulas that closely match the numerical solutions. The adiabatic sound speed is given by

(15) \begin{equation} c_{\textrm{s}} = \left(\gamma\frac{P}{{\unicode{x03C1}}}\right)^{1/2},\end{equation}

where $\gamma=5/3$ and P is the gas pressure. The sonic Mach number is given by

(16) \begin{equation} \mathcal{M}_{\textrm{p}} = \frac{V_{\textrm{p}}}{c_{\textrm{s,0}}},\end{equation}

where $c_{\textrm{s,0}}$ is the sound speed in the unperturbed medium. We shall see below that the subsonic regime, with $\mathcal{M}_{\textrm{p}}\lt1$ , is most relevant for this work. In this regime the Reference Kim, Kim and Sánchez-SalcedoK08 fitting formulas are

(17) \begin{equation} \mathcal{I}'_{\unicode{x03D5}} = 0.7706\ln\left(\frac{1+\mathcal{M}_{\textrm{p}}}{1.0004 - 0.9185\mathcal{M}_{\textrm{p}}}\right) - 1.4703\mathcal{M}_{\textrm{p}}\end{equation}

and

(18) \begin{equation} \mathcal{I}''_{\unicode{x03D5}} = -0.022\mathcal{M}_{\textrm{p}}^2(10-\mathcal{M}_{\textrm{p}})\tanh\left(\frac{3\mathcal{M}_{\textrm{p}}}{2}\right).\end{equation}

An analytical solution has been derived by Desjacques et al. (Reference Desjacques, Nusser and Bühler2022) that agrees with the numerical solutions of Reference Kim, Kim and Sánchez-SalcedoK08, but as the analytical solution is cumbersome we choose to compare our solutions with the fitting formulas (17) and (18).Footnote ^c

The z-component of the torque on a single perturber in the rest frame of the CM of the perturbers is given by $-(a/2)\mathcal{F}\mathcal{I}_{\unicode{x03D5}}$ , where $\mathcal{F}$ is given by equation (12) and $\mathcal{I}_{\unicode{x03D5}}$ by equation (14). Thus, for equal mass perturbers with equal Mach numbers, the torque on the binary is given by

(19) \begin{equation} {\unicode{x03C4}}_z = -a\mathcal{F}\mathcal{I}_{\unicode{x03D5}}.\end{equation}

However, in the simulation the azimuthal speeds of the equal-mass point particles are not precisely equal due to asymmetry in the gas distribution. One can generalize equations (17) and (18) by replacing $V_{\textrm{p}}$ with $V_{i,\unicode{x03D5}}$ in equation (12), and $\mathcal{M}_{\textrm{p}}$ with

(20) \begin{equation} \mathcal{M}_i = V_{i,\unicode{x03D5}}/c_{\textrm{s,0}}\end{equation}

in equations (17) and (18), where i represents the particle on which the force is being calculated, i.e. 1 for the AGB core particle and 2 for the companion particle. In this case, the total torque on the binary is given by

(21) \begin{equation} {\unicode{x03C4}}_z = -\frac{a}{2}\left(\mathcal{F}_1\mathcal{I}_{1,\unicode{x03D5}} + \mathcal{F}_2\mathcal{I}_{2,\unicode{x03D5}}\right).\end{equation}

Writing this general torque equation for equal mass perturbers explicitly, we have

(22) \begin{align} {\unicode{x03C4}}_z &= -2{\unicode{x03C0}} a{\unicode{x03C1}}_0(GM_{\textrm{p}})^2 \nonumber \\[2pt] &\quad\times\left[ \frac{I'_{\unicode{x03D5}}(\mathcal{M}_1)+I''_{\unicode{x03D5}}(\mathcal{M}_1)}{V_{1,\unicode{x03D5}}^2} +\frac{I'_{\unicode{x03D5}}(\mathcal{M}_2)+I''_{\unicode{x03D5}}(\mathcal{M}_2)}{V_{2,\unicode{x03D5}}^2} \right]. \end{align}

Equations (10) and (22) represent the approximate models for the torque on the binary that we will use to compare with the simulation.

Figure 3. Top left: Density contours at ${\unicode{x03C1}}=0.01{\unicode{x03C1}}_{\textrm{max}}(t)$ , $0.006{\unicode{x03C1}}_{\textrm{max}}(t)$ , and $0.005{\unicode{x03C1}}_{\textrm{max}}(t)$ in the orbital plane at the time $t=188.7\,\textrm{d}$ , with ellipse fitted to the contour ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}=0.006{\unicode{x03C1}}_{\textrm{max}}(t)$ , which was found to enclose effectively all of the gas producing significant torque (see also Figure 2). The ellipse is phase-shifted by an angle $\Delta\unicode{x03D5}$ with respect to the axis that passes through the particles. Top right: Similar to the top-left panel but now showing the plane perpendicular to the orbital plane and rotated clockwise by the angle $\Delta\unicode{x03D5}$ about the orbital axis, shown by the dashed line in the top left panel. The length of the ellipse major axis is set equal to that in the orbital plane, but the length of the minor axis is allowed to differ. Bottom left: Adapted from figure 1 of Reference Kim, Kim and Sánchez-SalcedoK08, showing the steady state for $\mathcal{M}_{\textrm{p}}=0.6$ sliced through the orbital plane in their idealized double-perturber model. The black circle shows the orbit of the point masses that perturb the background density. Colour shows the density contrast $\log\mathcal{D}$ , where $\mathcal{D}=(c_{\textrm{s,0}}^2a/Gm)\lambda$ with m the binary mass and $\lambda=({\unicode{x03C1}}-{\unicode{x03C1}}_0)/{\unicode{x03C1}}_0$ . Overplotted for comparison is the time-averaged best fit ellipse in the orbital plane from our simulation, with parameter values noted on the plot (see also Figure 5). Bottom right: Similar to the bottom left panel but now showing $\log\mathcal{D}$ for our simulation, at the same time as the top row, when ${\unicode{x03C1}}_0=0.44{\unicode{x03C1}}_{\textrm{c}}=3.16\times10^{-6}\,{\textrm{g}}\,{\textrm{cm}}^{-3}$ , and $c_{\textrm{s,0}}=\overline{c}_{\textrm{s}}=93.0\,{\textrm{km}}\,{\mathrm s}^{-1}$ . The region outside $0.44{\unicode{x03C1}}_{\textrm{c}}$ has negative values of $\mathcal{D}$ and is coloured grey. The contour ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ is plotted in yellow. The white contours near the softening radius (black circles, AGB core on the left and companion on the right) show contours of $\log\mathcal{D}=1.2\,\textrm{and}\,0.8$ , while the contours outside ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ show $\log\mathcal{D}=-0.4, -0.8\,\textrm{and}\,-1.2$ .

4.1. Orbit

The evolution of the distance a between the AGB core particle and companion is shown in Figure 1, and the orbit of the particles is shown in the inset. These results are similar to those in other CE simulations (e.g. Ohlmann et al. Reference Ohlmann, Röpke, Pakmor and Springel2016a; Chamandy et al. Reference Chamandy2019b; Chamandy et al. Reference Chamandy, Blackman, Frank, Carroll-Nellenback and Tu2020).

4.2. Torque exerted on the particles by the gas

In the top panel of Figure 2, we show the evolution of the z-component of the torque applied on the particles about their CM (left axis, solid black), from $t\approx115\,\textrm{d}$ onward. The evolution of the orbital separation is shown for comparison (right axis, dashed). The drastic reduction of the torque with time is consistent with previous studies involving different initial binary parameter values (Chamandy et al. Reference Chamandy2019a; Reference Chamandy, Blackman, Frank, Carroll-Nellenback and Tu2020).

In order to model this torque, we first need to determine which part of the gas contributes significantly to it. This is done by integrating the torque out to various trial values of ${\unicode{x03C1}}/{\unicode{x03C1}}_{\textrm{max}}$ , with ${\unicode{x03C1}}_{\textrm{max}}$ the maximum density in a slice through the orbital plane, for a given snapshot (throughout the simulation, the density maximum occurs near the position of the AGB core particle). It is found that using the trial value ${\unicode{x03C1}}/{\unicode{x03C1}}_{\textrm{max}}=0.006$ leads to a torque that deviates the least from the actual torque. In the top panel of Figure 2, it can be seen that the torque exerted by gas inside the contour ${\unicode{x03C1}}/{\unicode{x03C1}}_{\textrm{max}}=0.006$ (magenta) matches closely with that exerted by all the gas (black). Henceforth, we refer to this as the threshold or critical density ${\unicode{x03C1}}_{\textrm{c}}=0.006{\unicode{x03C1}}_{\textrm{max}}(t)$ . Further justification for the value $0.006$ is provided in Appendix B.

However, it is important to note that ${\unicode{x03C1}}_{\textrm{max}}(t)$ depends on the choices involved in modifying the initial mesa profile and softening the gravity around the AGB core particle (Section 2). Therefore, the value of ${\unicode{x03C1}}_{\textrm{max}}$ (and by extension, the ratio ${\unicode{x03C1}}_{\textrm{c}}/{\unicode{x03C1}}_{\textrm{max}}$ ) has limited physical significance. This is not a major concern because the role of ${\unicode{x03C1}}_{\textrm{max}}$ is simply to provide a convenient normalization for the density. We shall see below that the linear spatial scale of the contour ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}(t)$ is roughly proportional to the core-companion separation a(t) and that the morphological evolution is approximately self-similar, as also seen in Reference Escala, Larson, Coppi and MardonesE04. This self-similarity helps to explain why setting ${\unicode{x03C1}}_{\textrm{c}}$ to be a fixed fraction of ${\unicode{x03C1}}_{\textrm{max}}$ turns out to be a fruitful approach. A choice of reference density other than ${\unicode{x03C1}}_{\textrm{max}}$ , e.g., that at the particle CM or averaged over the circle with radius $a/2$ centred on the particle CM, could perhaps be used instead, and we plan to explore this possibility in future work.

4.3. Fitting an ellipsoid to the threshold density contour

To apply the uniform density ellipsoid model of Section 3.1, it seems natural to use the ellipsoid that best approximates the contour at the threshold density ${\unicode{x03C1}}=0.006{\unicode{x03C1}}_{\textrm{max}}$ . To simplify the fitting procedure, we choose to fit two ellipses in perpendicular planes (which uniquely define an ellipsoid) rather than directly fitting an ellipsoid. We first fit the ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ contour in the orbital plane with a co-planar ellipse, treating the semi-major axis A, aspect ratio $B/A$ , and phase shift $\Delta\unicode{x03D5}$ as fit parameters. We then fit the ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ contour in the plane perpendicular to the orbital plane that contains the ellipse major axis (i.e. the plane rotated by $\Delta\unicode{x03D5}$ from the line joining the particles), setting the length of the semi-major axis of this ellipse also equal to A but fitting for the semi-minor axis C. These two ellipses define an ellipsoid with major axis and one other axis in the orbital plane, and the third axis perpendicular to the orbital plane. The top two panels of Figure 3 show an example of such a fit, at the snapshot $t=188.7\,\textrm{d}$ , where the phase shift ( $\Delta \unicode{x03D5}$ ) is equal to $14.7^\circ$ (which also happens to be almost equal to the mean value of $\Delta\unicode{x03D5}$ over the time domain analyzed, as discussed below). The blue ellipses show the fits to the ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ blue contour in the orbital plane (left) and orthogonal plane containing the ellipse major axis (right).

Our interest is in modeling the slow spiral-in phase, so the models cannot be applied right from $t=0$ . To determine the appropriate time domain to apply the ellipsoid model, we qualitatively judged when the density contours in the orbital plane resemble ellipses, which led to the estimate $t\gt125\,\textrm{d}$ . This choice is supported by quantitative analysis shown in Figure 4, where we plot the root mean square error (RMSE) of our ellipsoidal fit, normalized to be between 0 and 1. It is seen that the RMSE decreases below $0.1$ after $t=125\,\textrm{d}$ and remains below this value thereafter.

In Figure 5 (left axis), we plot the evolution of key parameters of the fit: the ratio of the semi-major axis to orbital separation $A/a$ , the ratio of semi-minor to semi-major axis $B/A$ ( $C/A$ in the orthogonal plane), and the phase angle $\Delta\unicode{x03D5}$ (right axis). The best fit values oscillate on the orbital period; averaging over this variability we find that they are fairly constant, thought they also vary somewhat on longer timescales for reasons not yet understood. Furthermore, the values of B and C are comparable, which suggests that the ellipsoid can be approximated by a spheroid ( $B=C$ ). This is important because it allows us to make use of the simple expressions presented in Section 3.1. For the general case of ellipsoids, the quantities ${\unicode{x03B1}}$ , ${\unicode{x03B2}}$ , and ${\unicode{x03B6}}$ have an integral form, which would complicate modeling.

Figure 4. Relative root mean squared error (RMSE) of ellipse fitting in both planes for the ellipsoid model discussed in Section 3.1. (the highest value is set to unity and other values are calculated with respect to it). We observe a higher RMSE in the perpendicular plane as the major axis in this plane is forced to have the same value as that in the orbital plane. We start our analysis at $t=125\,\textrm{d}$ .

4.4. Mass evolution

In Figure 6 (left axis), we show the ratio of the mass of the gas enclosed by the ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}=0.006{\unicode{x03C1}}_{\textrm{max}}(t)$ surface to the mass $m=M_{\textrm{1,c}}+M_2$ of the particles. This ratio decreases with time as the particles inspiral and is always $\ll1$ . This is the regime for which Reference Escala, Larson, Coppi and MardonesE04 applied their uniform ellipsoid model. The mean density inside the ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ surface is shown in blue (right axis). The mean density $\overline{{\unicode{x03C1}}}$ is roughly constant between $t=125\,\textrm{d}$ and $t=270\,\textrm{d}$ , but decreases sharply thereafter. This sudden decrease may be an artefact caused by the ratio $r_{\textrm{soft}}/a$ becoming too large, which implies a larger fractional volume inside which the gravitational potential is underestimated (see also Section 5.2.1).

4.5. To what extent can the uniform spheroid model reproduce the torque?

To model the torque, we calculate the ${\unicode{x03B1}}$ and ${\unicode{x03B2}}$ parameters in equations (5)–(9) separately for each frame, taking

(23) \begin{equation} \widetilde{A} = A, \qquad \widetilde{B} = \frac{B+C}{2}.\end{equation}

For the mean density $\overline{{\unicode{x03C1}}}$ , we make the approximation that the torque applied by the heterogeneous ellipsoid defined by the surface ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}=0.006{\unicode{x03C1}}_{\textrm{max}}(t)$ can be replaced by that applied by a homogeneous spheroid with dimensions determined using equations (23), and with density equal to the mean density $\overline{{\unicode{x03C1}}}$ inside the ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ contour. The torque derived using equation (10), shown by the cyan line in the top panel of Figure 2, reproduces quite closely the actual torque (black line in the same figure).

Further, by approximating $\Delta\unicode{x03D5}$ and $B/A$ by their mean values over the whole time interval of analysis, $\langle\Delta\unicode{x03D5}\rangle=14.9^\circ$ and $\langle \widetilde{B}/\widetilde{A} \rangle=0.654$ , equation (10) can be written as

(24) \begin{equation} {\unicode{x03C4}}_z \approx - \frac{1}{2} {\unicode{x03C0}} G m \left(\widetilde{{\unicode{x03B2}}}-\widetilde{{\unicode{x03B1}}}\right) \cos \langle\Delta \unicode{x03D5}\rangle \sin \langle\Delta \unicode{x03D5}\rangle \overline{{\unicode{x03C1}}} a^2,\end{equation}

where $\widetilde{{\unicode{x03B1}}}=0.457$ and $\widetilde{{\unicode{x03B2}}}=0.771$ are the values of ${\unicode{x03B1}}$ and ${\unicode{x03B2}}$ when $\widetilde{B}/\widetilde{A}$ is replaced by $\langle \widetilde{B}/\widetilde{A}\rangle$ . Equation (24) was also used by Reference Escala, Larson, Coppi and MardonesE04 to model equal-mass inspiralling BSMBHs (see Section 1). We find that this approximation, shown as an orange line in the top panel of Figure 2, aligns quite well with the actual torque (black), though it introduces a phase shift of about a half period. The phase difference between the model and the actual torque can be explained by the direct dependence of equation (24) on the square of the separation a, which causes the modeled torque to be in phase with the separation.

Figure 5. Time evolution of key fit parameters for the lagging spheroid model: (i) ratio of semi-major axis A to separation a, (ii) ratios of semi-minor axis (B in the orbital plane and C in the perpendicular plane) to semi-major axis A, and (iii) lag angle $\Delta \unicode{x03D5}$ between the binary axis and the major axis of the fitted ellipsoid (right axis). The simulation output (thin lines) oscillates rapidly with time. Thick lines show $10\,\textrm{d}$ -moving averages and the dotted lines show mean values over the time domain of the analysis (125 days onward).

4.6. To what extent can solutions of the idealized binary perturbers problem reproduce the torque?

In this section, we compare the torque evolution in the simulation with that predicted by the Reference Kim, Kim and Sánchez-SalcedoK08 model, i.e. equation (19). The Reference Kim, Kim and Sánchez-SalcedoK08 model assumes that the unperturbed medium is uniform, and in order to apply this model, one needs to specify the unperturbed density ${\unicode{x03C1}}_0$ and unperturbed sound speed $c_{\textrm{s,0}}$ . It is not clear how this should be done, so we tried different choices. One possible choice is to set ${\unicode{x03C1}}_0(t)$ to the value of the density on the surface that bounds the region that contributes significantly to the torque, ${\unicode{x03C1}}_{\textrm{c}}=0.006{\unicode{x03C1}}_{\textrm{max}}(t)$ , and $c_{\textrm{s,0}}(t)$ to the mean sound speed interior to this surface. In practice, we introduce an adjustable scale factor multiplying ${\unicode{x03C1}}_0(t)$ , which is expected to be of order unity, so that we finally take

(25) \begin{equation} {\unicode{x03C1}}_0(t) = \xi{\unicode{x03C1}}_{\textrm{c}}(t), \qquad c_{\textrm{s,0}}(t) = \overline{c}_{\textrm{s}}(t),\end{equation}

where $\xi=0.44$ and bar represents average inside the contour ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}(t)$ . Of the methods tried, this is the one that best reproduces the torque measured in the simulation.

The Reference Kim, Kim and Sánchez-SalcedoK08 torque on the perturbers is then calculated using equation (22) and is shown by magenta lines in the bottom panel of Figure 2, dashed for $\xi=1$ and solid for $\xi=0.44$ . The modeled torque using $\xi=1$ is seen to agree quite well with the simulation from $t=270\,\textrm{d}$ onward, but this may be a coincidence since numerical effects may cause the torque to be overestimated in the simulation at late times (see Section 5.2.1). For $t\lt250\,\textrm{d}$ , the Reference Kim, Kim and Sánchez-SalcedoK08 torque obtained using $\xi=1$ is too large by a factor of order unity; after multiplying the modeled torque by $\xi=0.44$ , the model reproduces the simulation torque remarkably well, including (to varying degrees) the magnitude, amplitude and temporal phase evolution.

We also plot the particle Mach number $\mathcal{M}_{\textrm{p}}$ , equal to the particle tangential speed in the particle CM frame divided by the volume-averaged mean sound speed $\overline{c}_{\textrm{s}}$ averaged inside the contour ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ , in green using the right axis of the bottom panel of Figure 2. The value of $\mathcal{M}_{\textrm{p}}$ can be seen to vary on the orbital period and also on longer timescales, but is restricted to the range $0.5$ – $0.7$ , so the particles’ azimuthal motion is subsonic.

For completeness we now mention alternative choices for specifying ${\unicode{x03C1}}_0(t)$ and $c_{\textrm{s,0}}(t)$ . As a second approach, we tried setting $c_{\textrm{s,0}}(t)$ to the value of $c_{\textrm{s}}$ on the ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ surface, which leads to a larger Mach number and values of the torque that are about a factor of two larger than above. A third approach we tried was to take both ${\unicode{x03C1}}_0(t)$ and $c_{\textrm{s,0}}(t)$ to be equal to the average values inside the ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ surface. This produces results more similar to the chosen method (first approach) but the torque is similar in magnitude to what is obtained in the second approach. From the bottom left panel of Figure 3, one can see that the density at the contour ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ exceeds the background density. Thus, to the extent that our simulation lines up with the Reference Kim, Kim and Sánchez-SalcedoK08 model, it makes sense that the best-fit value of ${\unicode{x03C1}}_0$ for reproducing the torque is smaller than ${\unicode{x03C1}}_{\textrm{c}}$ .

Fourth, we considered using the original profile of the AGB star by setting ${\unicode{x03C1}}(t) = {\unicode{x03C1}}(r)$ with $r=a(t)$ , and likewise for the sound speed. This would be convenient if it worked, because then one would only need to know the initial conditions. However, this choice produces torque values that are higher than those measured in the simulation by about an order of magnitude, and a torque evolution that does not resemble the simulation even after dividing by a constant numerical factor. The failure of this fourth approach suggests that it is not straightforward to model accurately the torque at late times using the initial stellar profile and Reference Kim, Kim and Sánchez-SalcedoK08 model alone, at least not for companions whose masses are appreciable compared to that of the giant’s core. This is perhaps not surprising because transfer of orbital energy and angular momentum to the envelope causes the profiles of density and sound speed to be drastically affected, as shown in Figures C1 and C2 of Appendix C, and, as we argue in Section 5.2.2, the system loses memory of its previous state after about one sound-crossing time.

To summarize, the Reference Kim, Kim and Sánchez-SalcedoK08 model reproduces the simulation remarkably well between $t\approx125\,\textrm{d}$ and $t\approx250\,\textrm{d}$ when one adopts ${\unicode{x03C1}}_0=0.44{\unicode{x03C1}}_{\textrm{c}}(t)$ and $c_{\textrm{s,0}}=\overline{c}_{\textrm{s}}(t)$ , whereas at late times choosing ${\unicode{x03C1}}_0={\unicode{x03C1}}_{\textrm{c}}(t)$ produces closer agreement. This is discussed further in Section 5.2.1.

Figure 6. Time evolution of the ratio of the gas mass enclosed by the equipotential surface ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}=0.006{\unicode{x03C1}}_{\textrm{max}}(t)$ to the binary particle mass $m= M_{\textrm{1,c}} + M_2$ (left axis) and the mean density inside that surface (right axis).

4.7. Comparison of the two idealized models

As we have shown, the lagging spheroid model of Reference Escala, Larson, Coppi and MardonesE04 and the double perturber model of Reference Kim, Kim and Sánchez-SalcedoK08 both reproduce the torque in the simulation reasonably well. Thus, they would be expected to be consistent with one another. Indeed, Reference Kim, Kim and Sánchez-SalcedoK08 applied their model to the Reference Escala, Larson, Coppi and MardonesE04 simulation, arguing that their model can be used to explain the Reference Escala, Larson, Coppi and MardonesE04 results if $\mathcal{M}_{\textrm{p}}\sim0.6$ . The bottom-left panel of Figure 3 is adapted from Reference Kim, Kim and Sánchez-SalcedoK08 (their figure 1a). The quantity plotted in colour is $\log(\mathcal{D})=\log[c_{\textrm{s,0}}^2a/Gm)\lambda]$ , which is equal to a constant plus the logarithm of the fractional density perturbation, $\lambda=({\unicode{x03C1}}-{\unicode{x03C1}}_0)/{\unicode{x03C1}}_0$ (with $\lambda\ll1$ assumed in Reference Kim, Kim and Sánchez-SalcedoK08), for $\mathcal{M}_{\textrm{p}}=0.6$ . From the bottom panel of Figure 2 the relevant part of our simulation typically has $\mathcal{M}_{\textrm{p}}\sim0.6$ as well. To compare the ellipsoid model with that of Reference Kim, Kim and Sánchez-SalcedoK08, we annotate their figure by drawing an ellipse with parameters set equal to the time-averaged values from our simulation: $\langle A/a\rangle=1.47$ , $\langle B/A\rangle=0.649$ , and $\langle\Delta\unicode{x03D5}=14.9^\circ\rangle$ (see also Figure 3). The fairly close correspondence helps to explain why the Reference Escala, Larson, Coppi and MardonesE04 and Reference Kim, Kim and Sánchez-SalcedoK08 models produce similar values of the torque. The time-averaged best fit ellipse in the orbital plane from our simulation is seen to align quite well with the structure obtained for $\log(\mathcal{D})$ by Reference Kim, Kim and Sánchez-SalcedoK08, though the latter is not precisely elliptical as the wakes form elongated tails.

The bottom-right panel of Figure 3 shows $\log(\mathcal{D})$ from our simulation at $t=188.7$ d, which is the same snapshot plotted in the top two panels. To calculate $\mathcal{D}$ , we assume ${\unicode{x03C1}}_0=0.44{\unicode{x03C1}}_{\textrm{c}}=0.44\times0.006{\unicode{x03C1}}_{\textrm{max}}(t)$ , as in the best-fit torque model discussed above. The morphology is similar to that obtained by Reference Kim, Kim and Sánchez-SalcedoK08, but the wakes do not have tapered, elongated tails in the simulation. We can also see that $\lambda$ is higher near the particles in the simulation than it is for Reference Kim, Kim and Sánchez-SalcedoK08. This makes sense because Reference Kim, Kim and Sánchez-SalcedoK08 assume the perturbers to be weak ( $\lambda\ll 1$ ), but in reality they are not (see Section 5.2.2 below).

4.8. Spatial variation of key quantities and its evolution

In Figure 7, We show the time evolution of four different quantities, one per row. Columns show, from left to right, snapshots of slices through the orbital plane at $t= 138.9$ , $208.3$ , and $277.8\,\textrm{d}$ . At these times, the particle separations are respectively $a=19.5$ , $13.1$ , and $9.4\,{\textrm{R}}_{\odot}$ .

Figure 7. From left to right, columns show snapshots in the orbital plane at times $t=138.9$ , $208.3$ , and $277.8$ days. Primary core (left) and companion (right) softening spheres are shown in black or yellow in the bottom row. Rows from top to bottom are: (i) Mass density ( ${\unicode{x03C1}}$ , thick white contour for ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ , and thin grey contours showing ${\unicode{x03C1}}=4{\unicode{x03C1}}_{\textrm{c}}, 2{\unicode{x03C1}}_{\textrm{c}}, 0.5{\unicode{x03C1}}_{\textrm{c}}\,\text{and}\, 0.25{\unicode{x03C1}}_{\textrm{c}}$ ); (ii) The local Mach number of the particles $V_{\unicode{x03D5}}/c_{\textrm{s}}$ , where $V_{\unicode{x03D5}}$ is the $\unicode{x03D5}$ -component of the particle velocity in the particle centre of mass frame ( $V_{1,\unicode{x03D5}} = V_{2,\unicode{x03D5}}$ ) and $c_{\textrm{s}}$ is the sound speed, which depends on position (blue denotes the subsonic region and red the supersonic region); (iii) The difference over the sum of the angular speeds of the perturbers and gas $(\Omega-\Omega_{\textrm{gas}})/(\Omega+\Omega_{\textrm{gas}})$ (red being the region dominated by $\Omega$ , blue being the region dominated by $\Omega_{\textrm{gas}}$ ); (iv) The torque density on both particles. Black and white contours respectively show the negatively and positively contributing regions to the $\unicode{x03D5}$ -component of the drag. The highest contour levels near the particles are $\pm10^{12}\,{\textrm{dyn}}\,{\textrm{cm}}^{-2}$ , and contours are also plotted for $\pm10^{11}\,{\textrm{dyn}}\,{\textrm{cm}}^{-2}$ , $\pm10^{10}\,{\textrm{dyn}}\,{\textrm{cm}}^{-2}$ , and so on. Movies starting from $t=115.7$ days to the end of the simulation are available at https://doi.org/10.5281/zenodo.17575148.

The top row shows the mass density in the orbital plane. The white contour is the equidensity surface ( ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}=0.006{\unicode{x03C1}}_{\textrm{max}}$ ) that we fit with an ellipse, whereas the grey contours show a few other equidensity surfaces. The softening spheres of the AGB core (left) and companion (right) are represented by black circles. Note that the axis units are normalized to the current particle separation. Hence, the particle separation is always equal to 1 in these units, but the softening radius, which is constant ( $r_{\textrm{soft}}=2.41\,{\textrm{R}}_{\odot}$ ), appears to grow with time in these normalized units. The left and right snapshots are very similar despite 29 orbits elapsing between them, which implies a high degree of self-similarity in the problem as the separation decreases with time. Note that the AGB core particle has a somewhat higher concentration of gas around it than the companion particle, which is reasonable given that the AGB core particle starts the simulation surrounded by dense gas, whereas the companion particle starts the simulation outside the AGB star.

The second row shows the spatially dependent Mach number of the particles, equal to the particle tangential speed in the rest frame of the particle CM divided by the local sound speed $c_{\textrm{s}}$ . The gas in the torque-dominating region, delineated by the ${\unicode{x03C1}}={\unicode{x03C1}}_{\textrm{c}}$ contour, is subsonic. Note the spiral shock, particularly evident in the middle panel, where a sudden transition from blue (high sound speed) to red (low sound speed) is visible across the shock. A jump in density across the shock can be seen in the middle panel of the top row where density contours almost overlap.

The third row shows the normalized difference between the angular velocity of the perturber and gas $(\Omega-\Omega_{\textrm{gas}})/(\Omega+\Omega_{\textrm{gas}})$ , where $\Omega=(\Omega_1+\Omega_2)/2$ is the mean angular speed of the particles and $\Omega_{\textrm{gas}}$ is the local angular speed of the gas. Gas coloured red (blue) is rotating slower (faster) than the perturbers and gas coloured white is in corotation. Note that there is a fairly large region around each particle where the gas is corotating with the particles in their orbit. This corotating gas exerts a torque on the particles because it preferentially lags them. This assertion is supported by the fourth row, which shows the torque density (colour and contours). The black (white) contours show regions that contribute positively (negatively) to the drag in the azimuthal direction. We can see that while there is a large degree of symmetry, the positive contributions are slightly stronger (as evidenced by the slightly larger contours), leading to a net drag, rather than a net thrust.

5.1. Ellipsoid model

The simplest version of the ellipsoid model is the constant aspect ratio, constant lag angle prolate spheroid model, akin to the model used by Reference Escala, Larson, Coppi and MardonesE04. This minimalistic model depends on three parameters, viz. $\widetilde{A}/a$ , which is time-dependent, and $\widetilde{B}/\widetilde{A}$ and $\Delta\unicode{x03D5}$ , which are both treated as constants. Recall that $\widetilde{A}/a$ is relevant only because it affects the value of $\overline{{\unicode{x03C1}}}$ (see the discussion of Newton’s third theorem in Section 3.1). This model does a reasonable job of producing the simulation results, particularly the period-averaged magnitude of the torque, which deviates from the torque measured in the simulation by $\lt16\%$ in the time domain studied (compare orange and black lines in the top panel of Figure 2). The model also does a good job of reproducing the amplitude and period of the torque variation, though it does a poor job of reproducing the temporal phase. This minimalistic version of the model is successful because the values of the parameters are fairly constant (Figure 5). However, it is not clear whether this is generally true, given the large diversity of common envelope parameter values in nature. Moreover, we have not tried to motivate the values of the parameters from first principles and knowledge of the initial conditions (but a place to start would be the Reference Kim, Kim and Sánchez-SalcedoK08 model).

In any case, it is rather remarkable that the values of the parameters reported by Reference Escala, Larson, Coppi and MardonesE04, viz. $\widetilde{B}/\widetilde{A}=1/2$ and $\Delta\unicode{x03D5}=22.5^\circ$ , are close to the time-averaged values we obtained, viz. $\langle\widetilde{B}/\widetilde{A}\rangle=0.654$ , and $\langle\Delta\unicode{x03D5}\rangle=14.9^\circ$ . In addition, we obtain $\langle \widetilde{A}/a\rangle = 1.47$ ; Reference Escala, Larson, Coppi and MardonesE04 do not report this value, though they do find the behaviour to be self-similar, i.e. $\widetilde{A}/a\approx\textrm{const}$ . In fact, when we measure the properties of the ellipse in their figure 8, we find $A/a\approx1.43$ , $B/A=0.658$ , and $\Delta\unicode{x03D5}=16.1^\circ$ , values that are remarkably similar (almost identical) to the ones we obtain. This near-perfect correspondence is likely a coincidence, but it could suggest that the parameter values tend to converge to values that are fairly universal, in the sense that they only depend on a few basic quantities, such as $\mathcal{M}_{\textrm{p}}$ . This possibility should be investigated in future work.

A correspondence of merger phenomena on vastly different scales is further suggested by the morphology in gas density seen in the merging galaxy simulations of Bortolas et al. (Reference Bortolas, Gualandris, Dotti and Read2018); see the pre-merger snapshots shown in their figure 2, which are strikingly similar to those seen in our simulation.

For real CE events the ratio $q=M_2/M_{\textrm{1,c}}$ generally differs from unity, but the ellipsoid model cannot be directly applied to the case of unequal perturber masses. For $q\ne1$ , the density distribution is roughly pear-shaped (e.g. Gagnier & Pejcha Reference Gagnier and Pejcha2025), making it more difficult to model the resulting gravitational potential. However, the $q=1$ case provides the opportunity to gain physical insight and interpret the simulation results with the simplest model. Our result – that dynamical friction torque can be modeled as that due to a roughly constant-shape, $\sim a$ -sized mass distribution that lags the particle orbit by a $\sim$ constant phase angle – is a principle that may carry over to unequal mass cases, albeit with a more complicated shape. A generalized approach analogous to our ellipsoid modeling would be interesting to investigate for $q\ne 1$ cases. But the Reference Kim, Kim and Sánchez-SalcedoK08 model and extensions thereof (Kim Reference Kim2010) are already more promising as a template for modeling general cases since they already allow for $q\ne1$ .

5.2. Double perturber model

Unlike the ellipsoid model which is partly phenomenological, the Reference Kim, Kim and Sánchez-SalcedoK08 model is motivated from first principles but relies on various assumptions, discussed below in Section 5.2.2. Aside from these assumptions, an important limitation is that the parameters of the model, viz. the density and sound speed of the unperturbed medium, are difficult to constrain.

5.2.2. Consequences of neglecting various effects, as in the Reference Kim, Kim and Sánchez-SalcedoK08 model

The good agreement between the Reference Kim, Kim and Sánchez-SalcedoK08 model and the simulation suggests that several approximations of the former are valid in the region of parameter space simulated. For one, the Reference Kim, Kim and Sánchez-SalcedoK08 model neglects gas self-gravity, which is justified in the regime considered since the point masses dominate the potential. It is reasonable to neglect self-gravity if the ratio of the orbital period to the gas free-fall time, $P_{\textrm{orb}}/t_{\textrm{ff}}\sim\sqrt{G\overline{{\unicode{x03C1}}}}/\Omega\sim\sqrt{\overline{{\unicode{x03C1}}} a^3/m}$ is small. From Figure 6, $\overline{{\unicode{x03C1}}}\lesssim 1.8\times10^{-5}\,{\textrm{g}}\,{\textrm{cm}}^{-3}$ , and from Figure 1, $a\lesssim23\,{\textrm{R}}_{\odot}$ . Thus, we obtain $P_{\textrm{orb}}/t_{\textrm{ff}}\lesssim0.2$ at $t\approx125\,\textrm{d}$ and smaller thereafter, so self-gravity likely plays only a minor role.

Secondly, the Reference Kim, Kim and Sánchez-SalcedoK08 model solves a linearized set of equations which depend on the assumption that the density perturbations are small, i.e. the perturbers are weak. This constraint is relaxed and the full nonlinear equations are solved by Kim & Kim (Reference Kim and Kim2009), who study a single strong perturber in rectilinear motion, and by Kim (Reference Kim2010), who studies a single strong perturber in circular motion. However, the double perturber case has not yet been studied in the nonlinear regime. According to those papers, to be in the linear regime one requires $\mathcal{A}=GM_{\textrm{p}}/c_{\textrm{s,0}}^2r_{\textrm{s}}\ll1$ , where $r_{\textrm{s}}$ is the Plummer softening radius. In our CE simulation, $\mathcal{A}\approx12$ (taking $r_{\textrm{soft}}=2.8r_{\textrm{s}}$ ), so we are in the nonlinear regime. However, for $\mathcal{M}_{\textrm{p}}\lt1$ , which is the regime in which we find ourselves (bottom panel of Figure 2), solutions of the drag force are found not to differ very much from the linear case (Kim & Kim Reference Kim and Kim2009).Footnote ^d This likely helps to explain why the Reference Kim, Kim and Sánchez-SalcedoK08 model works so well to explain our CE simulation, up to a scaling factor of order unity. Nonlinear effects can also become important as $\mathcal{B}=Gm/c_{\textrm{s,0}}^2a$ increases, but for Mach numbers $\gt1$ (Kim Reference Kim2010). Even if the drag torque in the subsonic regime is well approximated by the linear equations, the morphology of the wakes produced by solving the full nonlinear equations, including weak shocks outside the orbit, matches better with what is obtained in our simulations (c.f. the top two panels of figure 5 of Kim Reference Kim2010).

A third idealization of Reference Kim, Kim and Sánchez-SalcedoK08 is that the unperturbed density and sound speed are spatially constant, which at first glance seems inconsistent with the CE case, where the density scale height of the original AGB profile can be comparable to the inter-particle separation during the simulation (Chamandy et al. Reference Chamandy2019a). However, the initial density profile gets drastically modified during the CE phase (e.g. Iaconi et al. Reference Iaconi2017; Chamandy et al. Reference Chamandy, Blackman, Frank, Carroll-Nellenback and Tu2020), and hence seems fairly irrelevant for the late times explored in this work. There seems to be little memory of the original profiles and flow pattern, aside from their influence on how much gas is present near the perturbers at later times.

Another idealization of Reference Kim, Kim and Sánchez-SalcedoK08 is that the CM of the binary is not in motion with respect to the envelope. In our CE simulation, the particle CM moves with a speed of a few ${\textrm{km}}\,{\textrm{s}}^{-1}$ (as can be estimated from the inset of Figure 1; see also Chamandy et al. Reference Chamandy2019b). As this speed is small compared to the orbital speed of ${\sim}60\,{\textrm{km}}\,{\textrm{s}}^{-1}$ (Figure C2), the CM motion is not expected to affect significantly the torque on the binary (Sánchez-Salcedo & Chametla Reference Sánchez-Salcedo and Chametla2014).

In the Reference Kim, Kim and Sánchez-SalcedoK08 model the semi-major axis of the binary orbit is fixed, whereas in our CE simulation it is decreasing with time. The timescale $\overline{a/\dot{a}}$ , where bar here denotes average over an orbital period, is much larger than the orbital period $P_{\textrm{orb}}$ . The flow is expected to reach a new steady state within a sound-crossing time (Kim & Kim Reference Kim and Kim2007; Desjacques et al. Reference Desjacques, Nusser and Bühler2022),

(26) \begin{equation} t_{\textrm{s}} \sim \frac{a}{c_{\textrm{s}}} \sim \frac{\mathcal{M}_{\textrm{p}} a}{V_{\textrm{p}}} \sim \frac{2\mathcal{M}_{\textrm{p}}}{\Omega_{\textrm{p}}} = \frac{\mathcal{M}_{\textrm{p}} P_{\textrm{orb}}}{{\unicode{x03C0}}} \lesssim 0.2P_{\textrm{orb}},\end{equation}

where the numerical estimate is made by taking $\mathcal{M}_{\textrm{p}}\sim0.6$ (bottom panel of 2). Thus, to a good approximation, the flow is able to continuously adjust to the local conditions and can be treated as quasi-steady.

In fact, the parameters ( ${\unicode{x03C1}}_0$ and $\mathcal{M}_{\textrm{p}}$ ) also vary on the shorter timescale $P_{\textrm{orb}}$ due to the orbital eccentricity (e.g. Szölgyén, MacLeod, & Loeb Reference Szölgyén, MacLeod and Loeb2022), but relation (26) shows us that even this timescale is large compare to $t_{\textrm{s}}$ . Hence, to a good approximation, one can treat the flow as quasi-steady and apply the Reference Kim, Kim and Sánchez-SalcedoK08 model at any point in time during the time interval studied. This is only possible because the perturbers move subsonically at these times. For recent studies of gaseous dynamical friction on perturbers in eccentric orbits see Buehler et al. (Reference Buehler, Kolyada and Desjacques2024) and O’Neill et al. (Reference O’Neill, D’Orazio, Samsing and Pessah2024).Footnote ^e

5.3. Non-applicability to very early or very late times

The duration of our simulation is about 40 orbits and the models are applied for the final 36 orbits. The symmetry that makes the Reference Escala, Larson, Coppi and MardonesE04 and Reference Kim, Kim and Sánchez-SalcedoK08 models applicable over the last roughly 36 orbits is not present during the initial plunge-in. At those times models based on wind-tunnel simulations (e.g MacLeod & Ramirez-Ruiz Reference MacLeod and Ramirez-Ruiz2015; MacLeod et al. Reference MacLeod, Antoni, Murguia-Berthier, Macias and Ramirez-Ruiz2017; De et al. Reference De2020), or single perturbers (e.g. Kim & Kim Reference Kim and Kim2007; Kim Reference Kim2010) may be effective. At very late times (perhaps $\gtrsim10^3$ orbits) the envelope will be largely ejected, and self-similar evolution is unlikely to persist. Some CE simulations reveal an asymmetric, partly evacuated region in the orbital plane just outside the binary at these late times (Gagnier & Pejcha Reference Gagnier and Pejcha2025; Vetter et al. Reference Vetter2025), so the radial density profile can be non-monotonic. Different torque models may be required to explain these late stages of evolution.