Common Envelope Wind Tunnel: Coefficients of Drag and Accretion in a Simplified Context for Studying Flows Around Objects Embedded Within Stellar Envelopes

TL;DR

This study uses hydrodynamic simulations in a wind tunnel setup to quantify drag and accretion coefficients for objects within stellar envelopes, revealing their dependence on flow parameters and implications for binary evolution.

Contribution

It introduces a 3D hydrodynamic wind tunnel model to systematically study flow coefficients around embedded objects in common envelopes, considering different density gradients and adiabatic exponents.

Findings

Coefficients of drag increase with steeper density gradients.

Drag and accretion are higher in more compressible flows with γ=4/3.

Rapid initial inspiral occurs near the stellar limb, affecting transient timescales.

Abstract

This paper examines the properties of flows around objects embedded within common envelopes in the simplified context of a "wind tunnel." We establish characteristic relationships between key common envelope flow parameters like the Mach number and density scale height. Our wind tunnel is a three-dimensional, cartesian geometry hydrodynamic simulation setup that includes the gravity of the primary and secondary stars and allows us to study the coefficients of drag and accretion experienced by the embedded object. Accretion and drag lead to a transformation of an embedded object and its orbit during a common envelope phase. We present two suites of simulations spanning a range of density gradients and Mach numbers -- relevant for flow near the limb of a stellar envelope to the deep interior. In one suite, we adopt an ideal gas adiabatic exponent of $\gamma=5/3$, in the other, $\gamma=4/3$. We find that coefficients of drag rise in flows with steeper density gradients and that coefficients of drag and accretion are consistently higher in the more compressible, $\gamma=4/3$ flow. We illustrate the impact of these newly derived coefficients by integrating the inspiral of a secondary object through the envelopes of $3M_\odot$ ($\gamma\approx5/3$) and $80M_\odot$ ($\gamma\approx4/3$) giants. In these examples, we find a relatively rapid initial inspiral because, near the stellar limb, dynamical friction drag is generated mainly from dense gas focussed from deeper within the primary-star's envelope. This rapid initial inspiral timescale carries potential implications for the timescale of transients from early common envelope interaction.

Figures

Equations

$$v_{\mathrm{k}}={\left(\frac{GM}{a}\right)}^{1/2}$$

$$\mathcal{M}=\frac{v_{\infty }}{c_{\mathrm{s},\infty }}.$$

$$R_{\mathrm{a}}=\frac{2G{M}_2}{v_{\infty }^2},$$

$$\frac{R_{\mathrm{a}}}{a}=\frac{2}{f_{\mathrm{k}}^2}\frac{M_2}{M}=\frac{2}{f_{\mathrm{k}}^2}\frac{1}{1+q^{-1}}.$$

$$H_{\rho }=-\rho {\left(\frac{d\rho }{dr}\right)}^{-1},$$

$${ϵ}_{\rho }=\frac{R_{\mathrm{a}}}{H_{\rho }}$$

$$\frac{d\rho }{dr}$$

$$\frac{dP}{dr}$$

$${\left(\frac{d\ln P}{d\ln \rho }\right)}_{\mathrm{envelope}}={\Gamma }_{\mathrm{s}},$$

$${\left(\frac{d\ln P}{d\ln \rho }\right)}_{\mathrm{ad}}={\gamma }_1,$$

$${\left(\frac{d\ln T}{d\ln \rho }\right)}_{\mathrm{ad}}={\gamma }_3-1,$$

$$\gamma ={\gamma }_1={\gamma }_3.$$

$$\frac{dP}{dr}=-g\rho .$$

$$\frac{dP}{d\rho }={\Gamma }_{\mathrm{s}}\frac{P}{\rho }={\gamma }_1\left(\frac{{\Gamma }_{\mathrm{s}}}{{\gamma }_1}\right)\frac{P}{\rho }=\left(\frac{{\Gamma }_{\mathrm{s}}}{{\gamma }_1}\right)c_s^2,$$

$$\frac{c_s^2}{\rho }\frac{d\rho }{dr}=-g{\left(\frac{{\Gamma }_{\mathrm{s}}}{{\gamma }_1}\right)}^{-1}.$$

$$\frac{c_s^2}{H_{\rho }}=\frac{G{M}_1}{a^2}{\left(\frac{{\Gamma }_{\mathrm{s}}}{{\gamma }_1}\right)}^{-1}.$$

$$c_s^2{H}_{\rho }^{-1}=\frac{v_{\mathrm{k}}^2}{a}\frac{M_1}{M}{\left(\frac{{\Gamma }_{\mathrm{s}}}{{\gamma }_1}\right)}^{-1}.$$

$$\frac{a}{H_{\rho }}=\frac{v_{\mathrm{k}}^2}{c_s^2}\frac{M_1}{M}{\left(\frac{{\Gamma }_{\mathrm{s}}}{{\gamma }_1}\right)}^{-1}=\frac{{\mathcal{M}}^2}{f_{\mathrm{k}}^2}\frac{M_1}{M}{\left(\frac{{\Gamma }_{\mathrm{s}}}{{\gamma }_1}\right)}^{-1}.$$

$${\mathcal{M}}^2={ϵ}_{\rho }\frac{(1+q{)}^2}{2q}{f}_{\mathrm{k}}^4\left(\frac{{\Gamma }_{\mathrm{s}}}{{\gamma }_1}\right)$$

$${ϵ}_{\rho }=\frac{2q}{(1+q{)}^2}{\mathcal{M}}^2{f}_{\mathrm{k}}^{-4}{\left(\frac{{\Gamma }_{\mathrm{s}}}{{\gamma }_1}\right)}^{-1}.$$

$$P=\left(\gamma -1\right)\rho e$$

$$a=f_{\mathrm{k}}^{2}GM=\frac{1}{2}{f}_{\mathrm{k}}^2(1+q^{-1}).$$

$${\vec{a}}_{\mathrm{grav},1}=-\frac{G{M}_1}{(y-y_1{)}^2}\hat{y}$$

$${\vec{a}}_{\mathrm{grav},2}=-\frac{G{M}_2}{|\vec{r}{|}^2}\hat{r}.$$

$$\dot{M}=\frac{1}{\Delta t}{\int }_{\mathrm{sink}}(\rho -{\rho }_{\mathrm{sink}})dV,$$

$${\dot{p}}_x=\frac{1}{\Delta t}{\int }_{\mathrm{sink}}(\rho -{\rho }_{\mathrm{sink}})v_{x}dV.$$

$$d{\vec{F}}_{\mathrm{grav}}=\frac{G{M}_2\rho dV}{r^2}\hat{r},$$

$$d{F}_{\mathrm{grav},x}=\frac{G{M}_2\rho dVx}{r^3}.$$

$$F_{\mathrm{df}}=\int d{F}_{\mathrm{grav},x}.$$

$${\dot{M}}_{\mathrm{HL}}=\pi R_{\mathrm{a}}^2{\rho }_{\infty }{v}_{\infty },$$

Full text

1 Introduction

Common envelope episodes occur in binary star systems when one star engulfs its companion (Paczynski, 1976). These episodes, during which drag on the surrounding gas transforms and tightens the binary orbit, are thought to be critical in the formation of compact binaries. In particular, binaries that are able to merge under the influence of gravitational radiation (like the recently discovered merging pairs in the LIGO experiment) must be brought to separations smaller than the size of the stars that created them (e.g. Kalogera et al., 2007; Abbott et al., 2016). This can only occur through a phase of orbital transformation like a common envelope episode. The details of this common envelope interaction phase therefore determine the exact nature of the resultant binary (Iben & Livio, 1993; Taam & Sandquist, 2000; Taam & Ricker, 2010; Ivanova et al., 2013b; Postnov & Yungelson, 2014).

A typical common envelope event is thought to pass through several phases (Podsiadlowski, 2001; Taam & Ricker, 2010). Over stellar evolution ($\gtrsim 10^{5}$ year) timescales, a star in a pair evolves from its main-sequence size onto the giant branch, growing in radius significantly. Along the way it may become so large that it starts to impinge on the orbit of its companion star, when the orbital separation is similar to the stellar radius, $a\sim R_{*}$. This growth initiates interaction between the pair, which exponentiates through exchange or loss of mass (and angular momentum) from the system (Pejcha et al., 2016b, a), or by tidal instability (the Darwin instability, which is relevant mainly in pairs of unequal mass) (e.g. MacLeod et al., 2016). Thus the onset of the interaction occurs over a timescale regulated either by mass loss or by tidal dissipation, perhaps lasting hundreds to thousands of orbital periods in either case (e.g. Tylenda et al., 2011; Nandez et al., 2014; MacLeod et al., 2016). Both of these processes desynchronize the orbit of the secondary star from that of the primary’s envelope.

Eventually, one object is engulfed within the envelope of the other and the common envelope phase begins. Supersonic relative motion between the engulfed object and the envelope gas leads to gravitational focussing and the buildup of a dense wake behind the embedded object, which exerts a gravitational drag on the orbital motion (Ivanova et al., 2013b). The result of this ‘dynamical friction’ drag is a rapid inspiral through the increasingly dense stellar envelope. This phase, sometimes called the dynamical plunge, has two possible conclusions. In some cases, the pair of stars merge. In others, the two stellar cores both heat a fraction of envelope material, resulting in subsonic relative motion between the embedded objects and the gas. With this transition to subsonic relative velocity, drag forces drop off dramatically (e.g. Ostriker, 1999) and the new binary’s orbit stabilizes (e.g. Ricker & Taam, 2012; Passy et al., 2012; Ohlmann et al., 2016b; Iaconi et al., 2017).

Despite significant recent effort and progress (e.g. Ricker & Taam, 2008, 2012; Passy et al., 2012; Nandez et al., 2014, 2015; Nandez & Ivanova, 2016; Ivanova & Nandez, 2016; Ohlmann et al., 2016b, c, a; Staff et al., 2016a, b; Iaconi et al., 2017), global simulations of common envelope remain challenging to perform with many potentially important processes and timescales at play. A particular concern is the resolution requirement of simulating the full spatial extent of the binary for many orbital timescales implies that the simulations are either very computationally expensive, allowing a small number of calculations to be performed (e.g. Iaconi et al., 2017), or they are performed with extremely low numerical resolution but can span some parameter space (e.g. Ivanova & Nandez, 2016). As a result of these numerical limitations and physical complexity, fully interpreting and learning from the results of these simulations has proven to be challenging.

In this paper, we adopt the complementary approach of studying a well-defined, but idealized scenario related to common envelope encounters in detail. We model flow past a gravitating object like one embedded in a common envelope phase in the context of a “wind tunnel” numerical setup. By restricting the scope of the problem from the global scenario, this idealized approach allows us to examine the importance of individual physical processes separately from the full, complex system. These results, in turn, can prove valuable in interpreting the findings of global calculations. Our work builds on a long history of study of supersonic flows past gravitating objects, starting with Hoyle & Lyttleton (1939) and Bondi & Hoyle (1944); Bondi (1952). The analytic scalings of these flows have also informed a large portion of our understanding of the hydrodynamics of common envelope interactions (see, for example Livio & Soker, 1988; Iben & Livio, 1993; Ivanova et al., 2013b).

Numerical studies have augmented this analytic understanding beginning with pioneering simulations by Hunt (1971). Later, work by Shima et al. (1985) was the first to solve for the properties of this flow with a finite volume computational method. This work also took the important step of calculating coefficients of dynamical friction drag due to the gravitational interaction of the object with its wake. Subsequent work studied flows in inhomogeneous media, making the results more directly applicable to understanding flow around objects embedded in the common envelope (Livio et al., 1986b; Soker et al., 1986; Livio et al., 1986a; Fryxell et al., 1987; Fryxell & Taam, 1988; Taam & Fryxell, 1989; Armitage & Livio, 2000).

Numerical advances facilitated what remains a benchmark series of simulations of Hoyle-Lyttleton flow by Ruffert (Ruffert, 1994a; Ruffert & Arnett, 1994; Ruffert, 1994b, 1995, 1996, 1997, 1999). This work is notable for its relatively high numerical resolution, and for being the first broadly successful attempt to span a wide parameter space of flow mach numbers, adiabatic exponents, and object sizes. Blondin & Pope (2009) and Blondin & Raymer (2012) show that with high resolution and modern numerics, three-dimensional flows in homogenous media are stable and reach a steady state with accretion rates on the order of the Hoyle & Lyttleton (1939) estimate.111See section 3.1 of MacLeod & Ramirez-Ruiz (2015a) for a more detailed discussion of the numerical assumptions and results of this recent work. While accretion and flow morphology are the focus of much of the above work, dynamical friction drag forces have also been a focus of recent numerical studies, some of which have conditions that are particularly relevant to common envelope flow (in particular, Sánchez-Salcedo & Brandenburg, 1999; Kim & Kim, 2009; Sánchez-Salcedo, 2012; Thun et al., 2016).

This paper adds to this history of numerical study of Hoyle–Lyttleton and related flows and their application to the common envelope phase of binary interaction. To do so, we expand on an idealized formalism for studying the dynamical inspiral phase of common envelope introduced in MacLeod & Ramirez-Ruiz (2015a) and model flow past an object embedded in a numerical “wind tunnel.” To determine the conditions of the wind, we consider stellar structures (and gas adiabatic exponents) relevant to two key regimes of common envelope encounters. In the convective envelopes of low-mass stars, gas pressure dominates and a $\gamma=5/3$ equation of state describes the gas well. In higher mass stars, radiation pressure is quite important and the gas response to compression is closer to $\gamma=4/3$. This paper examines both regimes (for encounters with a 1:10 mass ratio) and compares flow properties in each.

The remaining sections of this paper are organized as follows. In Section 2 we introduce key descriptive parameters for common envelope flows and derive relationships between them. These flow properties inform the setup of our numerical experiments – which we call the Common Envelope Wind Tunnel. Section 3 describes our numerical method. Section 4 examines the results of a set of numerical experiments comparing flows with $\gamma=5/3$ to those with $\gamma=4/3$. Section 5 illustrates the implications that these idealized results have for our understanding of the nature of typical common envelope inspirals. In Section 6 we conclude.

2 Flow Conditions During Common Envelope Inspiral

2.1 Characteristic Scales

Let us imagine the interaction between a giant-star primary with total mass $M_{1}$ and radius $R_{1}$ with a secondary object of mass $M_{2}$ and radius $R_{2}$, which will become embedded within the primary. The separation between these two objects is $a$, and during the interaction, $a<R_{1}$. We define the mass ratio of the system as $q=M_{2}/M_{1}$. The characteristic orbital velocity is

[TABLE]

where $M=M_{1}+M_{2}$. In general, the orbital motion of $M_{2}$ is desynchronized from the primary’s gaseous envelope and the relative velocity will be written as $v_{\infty}=f_{\rm k}v_{\rm k}$, where $f_{\rm k}$ is the fraction of keplerian velocity that describes the relative motion between the secondary object and the gas.

A long-standing conceptual framework for understanding flows during the dynamical plunge phase of the common envelope inspiral has been that of Hoyle & Lyttleton (1939) accretion flows (e.g. Taam et al., 1978; Meyer & Meyer-Hofmeister, 1979; Livio & Soker, 1988; Fryxell & Taam, 1988; Kato & Hachisu, 1991; Iben & Livio, 1993; Chevalier, 1993; Brown, 1995; Ivanova et al., 2013b). In these cases gas moves supersonically past a gravitating object. The Mach number is

[TABLE]

Gravitational focussing leads gas within an impact parameter,

[TABLE]

to be energetically bound to the accreting object of mass $M_{2}$. Note that the simple expression above of Hoyle & Lyttleton (1939), and later Bondi & Hoyle (1944), ignores the gas internal energy, which is added in the Bondi (1952) formalism (see Edgar, 2004, for a review).

The relationship between $R_{\rm a}$ and $a$ is dictated by the relative masses in the system and fraction of Keplerian rotation,

[TABLE]

This ratio describes the relative size of an accretion structure to the size of the binary orbit. Note that in the simplifying case of $f_{\rm k}=1$ and $M_{1}\gg M_{2}$, then $R_{\rm a}/a\approx 2q$.

A final length scale that plays a role in defining the common envelope interaction is the density scale height,

[TABLE]

where $d\rho/dr$ describes the density profile within the primary star’s envelope (e.g. MacLeod & Ramirez-Ruiz, 2015a). Another important ratio describes the density gradient across the accretion radius. We define the ratio

[TABLE]

to describe the number of density scale heights subtended by the accretion radius ($\epsilon_{\rho}\rightarrow 0$ describes a homogeneous density structure, where $\epsilon_{\rho}\rightarrow\infty$ describes a very steep density gradient).

Taken together, the accretion radius, Mach number, and density gradient describe the hydrodynamic properties of common envelope flows, which we will focus on in this paper.

2.2 Polytropic Stellar Envelopes

To obtain approximate profiles of the stellar envelope structure into which the secondary star plunges in the common envelope event, we will consider polytropic envelope profiles in hydrostatic equilibrium. Here we will also make the approximation of a mass-less envelope (which implies that the bulk of the mass is concentrated in the giant-star core). This is a crude approximation of stellar structure, but one that still yields useful results (as we will show in the following sections).

In this case, coupled differential equations of pressure and density profile describe the envelope structure,

[TABLE]

where $g={GM_{1}/r^{2}}$. The parameter $\Gamma_{\rm s}=1+1/n$ is the polytropic index of the stellar profile such that

[TABLE]

where the subscript denotes that this expression is evaluated along the envelope profile – a change in density in the stellar profile implies a change in pressure, $P_{1}\propto{\rho_{1}}^{\Gamma_{\rm s}}$.

2.3 Gas Equation of State

The envelope gas may have a different response to compression than its arrangement in the hydrostatic profile. General equations of state have four adiabatic indices, which describe their thermodynamic behavior (see, for example, chapter 3 of Hansen et al., 2004), in these cases there may be departures between $\gamma_{1}$, defined by

[TABLE]

and $\gamma_{3}$, which is defined by

[TABLE]

where the subscript indicates partial derivatives along an adiabat (at constant entropy). The gas’s adiabatic behavior is particularly relevant because any compression induced by a companion will happen on a timescale much shorter than the stellar envelope thermal timescale. The first exponent, $\gamma_{1}$, is relevant in computing the gas sound speed, $c_{\rm s}^{2}=\gamma_{1}P/\rho$, and the third, $\gamma_{3}$, enters into the equation of state relationship between pressure, density, and internal energy as $P=(\gamma_{3}-1)\rho e$. Constant entropy stellar envelope structures (for example, a convective envelope) have $\Gamma_{\rm s}\approx\gamma_{1}$, while other structures (like a radiative envelope) may have $\Gamma_{\rm s}<\gamma_{1}$.

For an ideal gas, all of the adiabatic exponents are identical. In this case, there is a single adiabatic exponent, $\gamma$, which is

[TABLE]

This implies that when ideal gas is compressed (or allowed to expand) adiabatically, the pressure follows $P\propto\rho^{\gamma}$.

2.4 Relationships between Flow Parameters

In a common envelope encounter, the secondary star plunges into the envelope of the primary. We use the (simplified) polytropic description above to show that there are relationships between the characteristic flow parameters described in Section 2.1.

These relationships directly result from the fact that the envelope is in hydrostatic equilibrium in opposition to the same gravitational forces that determine the secondary object’s orbit. Therefore, these relationships will hold for any hydrostatic envelope structure, which needs not be the pre-encounter stellar envelope structure.

To derive the relationships between the characteristic scales of common envelope flows, we start with the pressure gradient of the envelope,

[TABLE]

We re-write the left-hand side as $dP/dr=(dP/d\rho)\times(d\rho/dr)$. We can then use the derivative of pressure with respect to density within the envelope to find,

[TABLE]

because $c_{s}^{2}=\gamma_{1}{P/\rho}$. Substituting this in, our expression becomes,

[TABLE]

We substitute in the definition of the density scale height (equation 5), the definition of $g$, and set $r=a$ (the separation of the pair) to find

[TABLE]

Further substitutions into this expression are useful. We use the definition of $v_{\rm k}$, equation (1), and total system mass $M=M_{1}+M_{2}$ to write

[TABLE]

Then, rearranging and introducing the Mach number, equation (2), based on a flow relative velocity $v_{\infty}=f_{\rm k}v_{\rm k}$,

[TABLE]

We can then substitute in for the accretion radius $R_{\rm a}$ and $\epsilon_{\rho}$, equations (3) and (6), to express the relationships between the flow parameters,

[TABLE]

or

[TABLE]

We note that where the enclosed primary-star mass, $m_{1}(a)$, is substantially less than $M_{1}$, the enclosed mass may be replaced into the above equations by using $q_{\rm enc}=M_{2}/m_{1}(a)$ in place of $q$. These expressions are extremely useful because they reduce the multi-parameter space of common envelope flows down to a plane of allowed combinations on the basis of the hydrostatic equilibrium nature of the stellar envelopes.

2.5 Example Profiles for Two Primary Stars

We illustrate these relationships between flow parameters using a secondary object $(q=0.1)$ embedded within the unperturbed envelopes of two giant stars in Figure 1.

To compute these envelope profiles, we use the MESA stellar evolution code, version 8845 (Paxton et al., 2011, 2013, 2015). We show a $3M_{\odot}$ red giant with $Z=0.01$ that has evolved to have a $31R_{\odot}$ radius and an $0.43M_{\odot}$ helium core. The input list for this model is based on the 7M_prems_to_AGB test suite input list, with a change to $3M_{\odot}$. This input includes a mixing length $\alpha$ of 1.73, and a Reimer’s red giant branch wind prescription with coefficient 0.5. We also show a lower metallicity massive star, which is $80M_{\odot}$ with $Z=0.001$ and has evolved to $720R_{\odot}$, with a $41.2M_{\odot}$ helium core. This model was run using the 150M_z1m4_pre_ms_to_collapse test suite example inlist, modified to $80M_{\odot}$ initial mass and $Z=0.001$, and no other modifications to the inlist_massive_defaults parameters, which include a mixing length $\alpha$ of 1.5, and a semiconvection $\alpha$ of 0.01.222input lists available upon request to the corresponding author.

The panels of Figure 1 map out profiles of gas compressibility, density, and pressure within the stellar envelope along with profiles of Mach number and density gradient. The panels are normalized to the location of a hypothetical secondary, embedded within the stellar envelope, and the x-axes show distance in units of the accretion radius of this object, $R_{\rm a}$. Within $\pm R_{\rm a}$, we also show a polytropic reconstruction of the local profile. These panels adopt $q=0.1$ (which implies secondary masses of $0.3M_{\odot}$ and $8M_{\odot}$, respectively) and $f_{\rm k}=1$.

The flow in these and other common envelope encounters is described by these profiles of pressure and density across the accretion radius – but this description can be compactly represented in the parameters of the Mach number, density gradient, gas compressibility ($\gamma_{1}$, $\gamma_{3}$), and structural gamma ($\Gamma_{\rm s}$). The panels of Figure 1 show that equations (18) and (19) reproduce the Mach number and density gradient at a given position, and that a polytropic profile reproduces the approximate slope of these parameters around the central value.

The two examples in Figure 1 show overall similarity despite originating in relatively different stars: the highest Mach numbers and density gradients are found near the stellar limb where radiative losses contribute to a reduction of the scale height. Typical Mach numbers are $\mathcal{M}\sim 1-5$ and density gradients, $\epsilon_{\rho}$, are of order unity. Gas adiabatic exponents, $\gamma_{1}$ and $\gamma_{3}$, and the structural parameter, $\Gamma_{\rm s}$, are both $\approx 5/3$ in the interior of the $3M_{\odot}$ star’s convective envelope, but drop to lower values in zones of partial ionization nearer to the surface. The values of $\gamma_{1}$ and $\gamma_{3}$ diverge from $\Gamma_{\rm s}$ in the radiative interior of the $3M_{\odot}$ star. The more massive, $80M_{\odot}$, star has an extended convective envelope (so $\Gamma_{\rm s}\approx\gamma_{1}$) and a more compressible equation of state, with $\gamma_{1}\sim 1.4$ and $\gamma_{3}\sim 1.35$ due to a partial contribution to the pressure from radiation (see, e.g. Sanyal et al., 2017, for more details).

3 Numerical Approach: Common Envelope Wind Tunnel

We study flows under typical common envelope conditions using an idealized three-dimensional hydrodynamic setup, which we call the Common Envelope Wind Tunnel. This section outlines the details of our numerical method.

3.1 Hydrodynamic Implementation

We solve the equations of inviscid hydrodynamics using the FLASH code (Fryxell et al., 2000). FLASH is a grid based code with adaptive mesh refinement. We use the directionally split Piecewise Parabolic Method Riemann solver in the calculations presented here (Colella & Woodward, 1984). We use an ideal gas, gamma-law equation of state,

[TABLE]

and take different values of the compressibility, $\gamma$ in different simulations. As noted in equation (11), the ideal gas approximation assumes that $\gamma=\gamma_{1}=\gamma_{3}$, which is generally a reasonable (but inexact) approximation for thermodynamic conditions of interest for stellar envelopes – see Figure 1 to note the small departures between $\gamma_{1}$ and $\gamma_{3}$ in the MESA stellar models due to their more sophisticated treatment of the equation of state.

Like the simulations of MacLeod & Ramirez-Ruiz (2015a), our Common Envelope Wind Tunnel calculations have a 3D cartesian computational domain with a point mass representing the embedded object at the coordinate origin. The simulations are performed in dimensionless units, in which $R_{\rm a}=v_{\infty}=\rho_{\infty}=1$, where $\rho_{\infty}$ is the density of the primary-star envelope at $r=a$. These units imply a time unit of $R_{a}/v_{\infty}=1$, or one flow crossing time of the accretion radius. The mass of the embedded object is $M_{2}=(2G)^{-1}$ in these simulation units, and the primary has mass $M_{1}=q^{-1}M_{2}$. Because $v_{\infty}=f_{\rm k}v_{\rm k}=1$, we can solve for the binary separation, $a$ in code units,

[TABLE]

The orbital plane defined in the simulation is the $x-y$ plane. We locate the primary at $y_{1}=-a$. The gravitational force from the primary therefore acts in the $-y$ direction. As a concession to the cartesian geometry of our domain, the primary-star gravity only depends on the $y$-coordinate,

[TABLE]

where $\vec{a}_{\rm grav,1}$ is the gravitational acceleration from the primary star’s gravity. The acceleration from the point mass at the coordinate origin (representing $M_{2}$) is

[TABLE]

where $\vec{r}=(x,y,z)$ is the distance from the coordinate origin to the cell, where we are calculating the force within the computational domain.

The simulation $-x$ boundary feeds a wind into the wind tunnel and past the point mass. The wind has a gradient of pressure and density in the $y$ direction, and is uniform in the $z$ direction. The conditions of this wind are specified by an upstream Mach number, $\mathcal{M}$, a density gradient, $\epsilon_{\rho}$, and the pressure and density at $y=0$ (along the $x$ axis). We begin by specifying a corresponding pair of $\epsilon_{\rho}$ and $\mathcal{M}$ given $q$, $f_{\rm k}$, $\Gamma_{\rm s}$, and $\gamma$ using equations (18) and (19). We then assign $P_{\infty}$ to generate the sound speed that satisfies the selected flow Mach number, based on $\rho_{\infty}$ and $c_{s}^{2}=\gamma P/\rho$. We therefore have $P_{\infty}=\mathcal{M}^{-2}\gamma^{-1}\rho_{\infty}v_{\infty}^{2}$ (note that $\rho_{\infty}=v_{\infty}=1$ in our code units).

Once the values at $y=0$ are set, we integrate to both positive and negative $y$ using the expressions of hydrostatic equilibrium for a mass-less atmosphere, see equation (2.2). Here the relevant differentials in our code units become $dP/dy$ and $d\rho/dy$. We extend the hydrostatic profile to the ghost zones in the $-y$ boundary such that the hydrostatic pressure gradient is preserved. The wind fed into the box is therefore in hydrostatic equilibrium with the primary-star gravitational force, and since it is supported on its lower boundary, it does not rise or fall unless an additional force is applied. On the $+x$, $+y$, and $\pm z$ boundaries, we apply ‘diode’ boundary conditions, which allow material to freely leave but not enter the grid.

The initial condition is uniform in the $x$ and $z$ directions with flow properties based on the integrated profile of pressure and density in the $y$ direction. The velocity everywhere is set to $\vec{v}=v_{\infty}\hat{x}$. We turn the central point mass on progressively over the first code time unit, so $\vec{a}_{\rm grav,2}$ is fully active for $t>1R_{\rm a}/v_{\infty}$.

We create an absorbing central “sink” surrounding the point mass, with radius $R_{\rm s}$. The calculations in this paper use $R_{\rm s}=0.05R_{\rm a}$. Each timestep, the average pressure and density of a spherical shell, which extends from $R_{\rm s}$ to $2R_{\rm s}$ are computed. The conditions inside the sink are reset to a fraction, usually $10^{-3}$ of these values, creating an effective vacuum – and deleting (accreting) mass and energy from the grid every timestep. This prescription represents accretion without feedback on the surrounding flow (MacLeod & Ramirez-Ruiz, 2015a).

Our computational domain extends from $\pm 4R_{\rm a}$ in the $x$ direction and $\pm 3.5R_{\rm a}$ in the $y$ and $z$ directions. This domain is covered by 8 blocks in $x$ and 7 each in $y$ and $z$ of $8^{3}$ cells in each direction. We employ the PARAMESH adaptive mesh refinement package (v4), and base refinement choices on the second derivative of gas internal energy (erg g*-1*) (MacNeice et al., 2000). We set the minimum refinement level to 2 (so all blocks are refined at least once) and the maximum refinement to 6. The maximum cell size is therefore $R_{\rm a}/16$ and the minimum is $R_{\rm a}/256$. As in MacLeod & Ramirez-Ruiz (2015b), to focus the highest resolution cells in the center of the computational domain (near $M_{2}$), we drop the maximum refinement with distance. Blocks with size less than $\alpha r$ (we adopt $\alpha=0.3$) are not allowed to refine further. The first drop in refinement occurs at $\gtrsim 1.5R_{\rm s}$ and drops one level each time the distance from the point mass doubles.

3.2 Diagnostics of Flow Properties

Several key diagnostics and integral quantities of the flow are computed at runtime and recorded every timestep in our setup.

3.2.1 Accreted Quantities

We record the properties of material that falls into the central sink just prior to deleting it. These quantities represent the accreted mass, angular momentum, and linear momentum. Each timestep, we perform a volume integral over the sink cells and sum the total accreted quantity. We convert this sum to an accretion rate by dividing by the timestep. For example, the accretion rate of mass is given by

[TABLE]

where $\rho_{\rm sink}$ is the density the sink cells were set to on the previous timestep. Similarly,

[TABLE]

is the accretion rate of linear momentum along the direction of motion. This accretion of momentum also represents a force that we will call $F_{\dot{p}_{x}}$ in what follows.

3.2.2 Gravitational (Dynamical Friction) Drag Forces

If the mass distribution around the embedded object is not spherically symmetric, it experiences a net gravitational force. The component of this force directed along the direction of motion of the object constitutes a “gravitational drag” or gas dynamical friction that modifies the motion of the gravitating object (Chandrasekhar, 1943; Ostriker, 1999). In the case of supersonic flows past a gravitating object, an overdense wake is generated that exerts a stronger gravitational force than the upstream gas (Ostriker, 1999). The net force decelerates motion, and is thus termed a drag.

The gravitational force on the secondary, $M_{2}$ by a volume of gas $dV$ is

[TABLE]

(note the inversion of the sign in this expression as compared to equation 23). The component of this force along the direction of motion is,

[TABLE]

The net dynamical friction drag, $F_{\rm df}$, is the volume integral of the contributions to the gravitational force along the direction of motion $dF_{{\rm grav},x}$,

[TABLE]

The sign of $F_{\rm df}$ in our coordinate setup is such that a positive value represents a drag force (deceleration).333Alternatively, one could also calculate the drag force by measuring the momentum and pressure change of the gas (rather than the gravitational net force on the particle). We note here that Ricker & Taam (2008, equation 3) and MacLeod & Ramirez-Ruiz (2015a, equation 13) used this approach and measured the net momentum transport by gas passing through a spherical surface to evaluate drag forces generated within the enclosed volume. However, in their analysis of the momentum equation in steady state, Thun et al. (2016) show that these expressions are incomplete because they do not include the difference in pressure across the surface. In their section 3.5, Thun et al. (2016) show that the sum of a surface integral of momentum transport and a surface integral of (net) pressure balance dynamical friction (see Thun et al. (2016), section 3.5 for a complete derivation). For this application, that of supersonic Hoyle-Lyttleton flow, the pressure term is opposite in sign and smaller in magnitude than the momentum transport term (see, for example, Figure 3 of Thun et al., 2016). This suggests that the drag forces derived from Ricker & Taam (2008, equation 3) and MacLeod & Ramirez-Ruiz (2015a, equation 13) are of the correct magnitude, but are likely moderate overestimates of the drag force generated within a particular volume.

4 Numerical Results

In this section, we present and analyze two suites of calculations, each spanning a range of density gradients and corresponding Mach numbers. Each calculation adopts $R_{\rm s}=0.05R_{\rm a}$. In one suite, we take $\gamma=\Gamma_{\rm s}=5/3$, relevant for the convective envelopes of low-mass stars, as shown in Figure 1. In a second suite, we take $\gamma=\Gamma_{\rm s}=4/3$, as exemplifying the high-compressibility limit of massive star envelopes in which radiation pressure becomes increasingly important. Here we compare some key flow properties, including rates of accretion and the generation of drag forces, realized in these simulations.

4.1 Flow Morphology

The morphology of flows around objects embedded in the common envelope is distorted and asymmetric in response to the gradient of density in the upstream, stellar envelope material (MacLeod & Ramirez-Ruiz, 2015a). Our new calculations are consistent with this result, but the inclusion of the primary-star properties – in the form of gravitational force, hydrostatic equilibrium pressure gradients, and the relationships between flow parameters discussed in Section 2.4 – impacts the expected nature of common envelope flows.

In Figure 2 through Figure 5, we show slices of density and Mach number through the orbital ($x$-$y$) plane, and perpendicular to the orbital $(x-z)$ plane for the simulation suites with $\gamma=\Gamma_{\rm s}=5/3$ and $\gamma=\Gamma_{\rm s}=4/3$, respectively. In the each case, we run simulations for six log-spaced values of the density gradient, $\epsilon_{\rho}$. For $\gamma=\Gamma_{\rm s}=5/3$, this corresponds to $\epsilon_{\rho}=0.2,0.32,0.5,0.8,1.26,2.0$. For $\gamma=\Gamma_{\rm s}=4/3$, we have $\epsilon_{\rho}=0.3,0.48,0.75,1.19,1.89,3.0$. The ranges of density gradient parameters were chosen such that in the steepest gradient cases, the object was embedded to a depth of approximately one accretion radius within the envelope of the primary.

The corresponding upstream Mach numbers also vary across these simulations, following equation (18) with $f_{\rm k}=1$, from $\mathcal{M}=1.1$ for $\epsilon_{\rho}=0.2$ to $\mathcal{M}\approx 3.48$ for $\epsilon_{\rho}=2$ and $\mathcal{M}\approx 4.26$ for $\epsilon_{\rho}=3$. By comparison to Figure 1 and the associated discussion, we can see that the steepest density gradients and the highest Mach numbers correspond to flow near the limb of the primary-star envelope, while shallower gradients (and lower Mach numbers) are found deep within the common envelope. Therefore, in our panels of Figures 2 through 5, the upper left panels correspond to flow around a deeply embedded object, while the lower right panels correspond to flow around an object nearer to the envelope limb.

The flow Mach number is plotted in the lower panel sets of Figures 2 through 5. A dramatic transition occurs here with upstream density gradient. In the case of shallow density gradients, the symmetry of the bow shock is nearly preserved. Upstream flow is supersonic, while downstream flow is subsonic after crossing the shock and meeting a pressure gradient imposed by the convergence of flow into the post-shock region. As the density gradient steepens the portion of material in the post-shock region that remains supersonic increases dramatically. In flows with the steepest density gradients, the material moves nearly ballistically with only a small fraction having $\mathcal{M}\ll 1$.

One implication of the changing Mach number can be seen in the bow shock morphology. Bow shocks in homogenous flow exhibit an opening angle, $\Theta$, proportional to the Mach number, where $\sin\Theta\sim\mathcal{M}^{-1}=c_{\rm s,\infty}/v_{\infty}$, because the disturbance from the shock wave moves laterally at approximately the sound speed while the stream motion is supersonic with $v_{\infty}$. As the Mach number (and density gradient) increases in these simulations we see a narrowing of shock opening angles – particularly along the well-defined edge facing the stellar center (and the flux of the densest material). While the shock for $\epsilon_{\rho}=0.2$ in Figures 2 and 3 is nearly planar, by the time $\epsilon_{\rho}=2.0$, the trailing shock opens in a much narrower cone. This consequence of choosing corresponding combinations of density gradient and Mach number can be contrasted to the simulations of MacLeod & Ramirez-Ruiz (2015a), which adopted $\mathcal{M}=2$ for all density gradients.

The shock morphology shows an interesting secondary effect perpendicular to the orbital plane in cases of the steepest density gradient, particularly for $\epsilon_{\rho}=1.19,1.89,3.0$ in Figure 5. In these snapshots, we see that the shock opening angle is not constant, but, in fact, widens with increasing displacement into the wake. What we observe from the streamlines in Figure 4 is that material focused onto the wake at larger $+x$ displacements comes from a larger impact parameter in the $-y$-direction. Recalling the profiles of Figure 1, this material, originating from deeper within the stellar envelope, has higher sound speed. As a result, there is a gradient of upstream Mach number in the $y$-direction (which can be observed in the lower panels of Figure 2 and 4). The shock opening angle, which depends inversely on this upstream Mach number, thus broadens as the focussed material is drawn from deeper in the stellar envelope potential well. This effect is observable primarily in cases of steep gradient (near the envelope limb), where the derivatives of $\mathcal{M}$ and $\epsilon_{\rho}$ become large.

The equation of state of the stellar envelope gas also plays a role in determining flow structure. The flow in the $\gamma=\Gamma_{\rm s}=4/3$ shown in Figures 4 and 5 is more compressible than the flow in the $\gamma=\Gamma_{\rm s}=5/3$ shown in Figures 2 and 3. This results in higher densities in the immediate wake of the embedded object because the pressure does not build up as rapidly upon compression in the focused material. In the steeper-gradient cases of $\gamma=\Gamma_{\rm s}=4/3$, we see a nested shock outside of an accretion line, which differs from the much broader fan of material in the $\epsilon_{\rho}=2$, $\gamma=\Gamma_{\rm s}=5/3$ simulation.

In all cases, the secondary’s gravitational focus lifts some dense material from the stellar interior against the primary star’s gravity. This gravitational force leads some material (with impact parameter $\gg R_{\rm a}$) to rise and fall in a “tidal bulge” trailing the embedded object. In material with impact parameter $\lesssim R_{\rm a}$, as shown in the streamlines overplotted on the upper panels of Figures 2 and 4, this gravitational force leads to a slingshot around the embedded object. Some of this gas leaves the simulation box after being deviated through a large angle then expelled toward the lower-density of the primary star’s limb ($+y$-direction in our simulation setup).

4.2 Rates of Accretion

Our numerical approach replaces the embedded object with a sink on the grid, which absorbs convergent flow. The sink has a radius of $R_{\rm s}=0.05R_{\rm a}$. We note that this sink could be of a similar scale to that of a main-sequence star embedded in a typical common envelope (see Table 1 of MacLeod & Ramirez-Ruiz, 2015a), but is certainly much larger than the size of an embedded compact object like a white dwarf, neutron star, or black hole. Here we study rates and properties of material accreting through this inner boundary of our computational domain, but note that the accretion rate is dependent on the size of the sink boundary compared to the accretion radius (for example, MacLeod & Ramirez-Ruiz, 2015a, found lower accretion rates for $R_{\rm s}=0.01R_{\rm a}$ than for $R_{\rm s}=0.05R_{\rm a}$).

We begin by examining the mass accretion rate into the sink boundary as a function of density gradient in our $\gamma=\Gamma_{\rm s}=4/3$ and $\gamma=\Gamma_{\rm s}=5/3$ simulation suites, shown in Figure 6. Accretion rates in Figure 6 are normalized to the Hoyle-Lyttleton accretion rate,

[TABLE]

which is the flux of material passing through a cross-section of area $\pi R_{\rm a}^{2}$ assuming a uniform density background. One feature of the accretion rate is that when density gradients are introduced into the flow, the flow morphology becomes less laminar and variability is introduced into the mass accretion rate. Therefore, we plot the median values (pink and blue lines) along with the 5-th and 95-th percentile ranges (shaded regions) for mass accretion rate, $\dot{M}$, as a function of density gradient, $\epsilon_{\rho}$.

As density gradients steepen, the accretion rate into the sink drops dramatically and becomes more variable. We see accretion coefficients ($\dot{M}/\dot{M}_{\rm HL}$) spanning more than an order of magnitude as density gradient changes across typical values. In all regions, the accretion efficiency is substantially lower than accretion from a uniform medium. The imposition of a density gradient breaks the symmetry of the inflowing material (as seen in Figures 2 and 4). As opposed to the uniform medium case, where momenta of opposing streamlines cancel, there is net angular momentum in the flow, which forms a barrier to efficient accretion when the circularization radius is significantly outside the sink radius (MacLeod & Ramirez-Ruiz, 2015a). The increased variability in cases of steep density gradient can be attributed to the increased turbulence of the post-shock regions, as seen in Figures 2 through 5. The more compressible $\gamma=4/3$ flow accretes at higher rates, particularly in cases of mild density gradient, $\epsilon_{\rho}\lesssim 1$, where radial pressure gradients oppose flow convergence less strongly than in the $\gamma=5/3$ case.

Figure 6 also compares our new accretion rates to a fitting formula to the results of MacLeod & Ramirez-Ruiz (2015a) for $R_{\rm s}=0.05$, labeled M2015. The M2015 simulations all used $\mathcal{M}=2$, $\gamma=5/3$, and the density gradient in the background material, $\epsilon_{\rho}$, was uniform rather than polytropic. Finally, there was no corresponding pressure gradient (a uniform pressure background was assumed). We find that despite these differences, accretion rates of similar order of magnitude are found. However, differences appear to lie in the functional form of $\dot{M}(\epsilon_{\rho})$ and in the accretion rate for mild values of $\epsilon_{\rho}\lesssim 1$. The two simulation suites presented here show higher mass accretion rates for $0.2\lesssim\epsilon_{\rho}\lesssim 1.5$ by a factor of a few than M2015. One likely contribution to this difference is the lower Mach number in our current simulations for these density gradients.

Turning now to the accretion of angular momentum, Figure 7 evaluates the distributions of the magnitude of specific angular momentum, $|l|$, of material absorbed by the sink boundary. These are normalized to the Keplerian specific angular momentum at the sink surface, $l_{\rm kep}$ (for details see MacLeod & Ramirez-Ruiz, 2015a, section 4.3). In mild density gradient cases, accreted material has a relatively narrow distribution of specific angular momenta, with typical values much less than Keplerian. In these cases, the mass accretion rate is high (compare to Figure 6) because the net angular momentum of the flow does not substantially oppose accretion when $|l|\ll l_{\rm kep}$. At higher values of the density gradient, the distributions of specific angular momenta of accreted material are much broader, with typical values of $|l|/l_{\rm kep}\sim 0.5$. None of our simulations show signs of accreting material with nearly complete rotational support $|l|\sim l_{\rm kep}$, which makes sense because in a given timestep, any material that is fully rotationally supported will be unable to accrete. However, the $\gamma=4/3$ simulations show systematically higher specific angular momenta than those of the $\gamma=5/3$ simulations.

4.3 Drag Forces

With the inclusion of the gravitational restoring force of the primary star (and hydrostatic equilibrium pressure gradient) our simulations are much better suited to the evaluation of drag forces than those of MacLeod & Ramirez-Ruiz (2015a). This section studies contributions to the force on an embedded object from gaseous dynamical friction and from the accretion of linear momentum by the sink boundary. The details of both terms are outlined in Section 3.2.

We integrate the gravitational drag force (dynamical friction), given by equation (28), over ten different volumes each timestep. The integration volumes are spherical shells with inner radius equal to $R_{\rm s}$ and outer radius evenly spaced in $\ln(r_{\rm out})$ between $0.1R_{\rm a}$ and $3.5R_{\rm a}$. The outer integration radii therefore are, $r_{\rm out}/R_{\rm a}\approx 0.1$, 0.15, 0.22,0.33, 0.49, 0.72, 1.06, 1.59, 2.36, 3.5. We find that in the steepest gradient cases, our results for the outer two radii are somewhat sensitive to the boundary position (whether the $y,z$ domain extends $\pm 3.5R_{\rm a}$ or $\pm 4R_{\rm a}$) because of the diode (no inflow) boundary conditions imposed. We therefore show results for integration radii only out to $1.59R_{\rm a}$ in what follows, which are converged with respect to boundary location in even the steepest gradient cases.

The motivation for logarithmically spaced integration bins is that dynamical friction forces (in both gaseous and collisionless systems) in uniform media grow as $\propto\ln(r_{\rm out})$ (e.g. Chandrasekhar, 1943; Ostriker, 1999). By spacing our integration bins in this manner, each bin contains a similar contribution to the total force. The ‘advective’ force due to the rate of accretion of linear momentum by the sink boundary is given by equation (25) and is labeled $F_{\dot{p}_{x}}$ here. In our coordinate system the dynamical friction force acts in the positive direction (and is therefore a drag), while net accretion generally acts in the negative direction (and is therefore a thrust) because most material is accreted from behind the embedded object.

We begin by examining the net drag forces as a function of time in some example simulations. Figure 8 plots the net force on the embedded object for each value of the outer integration radius of the dynamical friction force, $r_{\rm out}$, as denoted by line color. We show results for two simulations from the $\gamma=\Gamma_{\rm s}=5/3$ simulation suite. These simulations have differing density gradients, $\epsilon_{\rho}\approx 0.5$ and $\epsilon_{\rho}\approx 0.8$. In both cases, we see an initial transient behavior while the flow sets up (the box crossing time is $\approx 8R_{\rm a}/v_{\infty}$), followed by a settling to a steady state.

In the $\epsilon_{\rho}\approx 0.5$ simulation, the flow visualizations in Figures 2 and 3 show relatively smooth post-shock structure without substantial small-scale vorticity. As discussed in relation to Figure 7, this is likely because the sink is large enough to swallow much of the circularizing material. As a result of the smooth wake, the dynamical friction force experienced by the embedded object is also relatively smooth. The net force is negative (a thrust) when the dynamical friction is only integrated to very small radii $\lesssim 0.3R_{\rm a}$. When we include the contributions from progressively larger radii, $F_{\rm df}$ outweighs $F_{\dot{p}_{x}}$ and the net force is positive (a drag).

The $\epsilon_{\rho}\approx 0.8$ panel of Figure 8 shows many similarities to the $\epsilon_{\rho}\approx 0.5$ panel, but exhibits substantially greater time variability. It is interesting to note that features in the variability overlap at many scales in the cumulative drag force plotted. Some variability is imprinted at the smallest scales, in particular, the short timescale (but relatively small amplitude) variation. The majority of the variability in the net drag is imposed at larger scales, of order the shock standoff distance. This occurs as vorticies shed in the wake cause some breathing and instability of the position of the bow shock. Since this large-scale flow instability is not present in the shallower-gradient cases, the drag force is much steadier. An interesting caveat to this point is that the size of the sink boundary likely plays a role in both the time variability of the drag force and which flow parameters generate highly variable post-shock regions. As shown in MacLeod & Ramirez-Ruiz (2015a, Figure 7), smaller sink boundaries result in more vorticity in the post-shock region and more variable accretion. We can speculate that the boundary condition might imprint itself on the time variability of the drag force in a similar manner.

The relative displacement of the lines in Figure 8 is also informative. To first order, the spacing in these simulations is relatively uniform, indicating that the contribution to the dynamical friction drag is growing approximately $\propto\ln(r_{\rm out})$. A transition to slightly larger spacings among the largest integration radii comes after passing the approximate standoff distance of the bow shock – in this simulation, $\sim 0.5R_{\rm a}$. Within the shock standoff radius, the density field is more symmetric (though still not entirely so) than integration radii that include the shock (e.g. Thun et al., 2016).

Figure 9 shows the cumulative drag within different (logarithmically spaced) integration radii. As seen in the time series of Figure 8, the net drag is negative (a thrust) when the dynamical friction is only integrated out to a small radius. In the $\gamma=\Gamma_{\rm s}=5/3$ case, we see that logarithmically spaced bins contribute roughly equally to the cumulative drag force. This implies that the net drag is growing according to $\sim\ln\left(r_{\rm out}/r_{\rm in}\right)$, as is the case for gas dynamical friction in homogenous media (Ostriker, 1999; Thun et al., 2016). The contribution of each increasing radius bin tells us something about the characteristic scale $r_{\rm in}$ in the dynamical friction force. A numerical comparison quickly reveals that the appropriate $r_{\rm in}$ is not $R_{\rm s}$, the radius of the inner boundary of our computational domain, but is instead something of order $\sim 0.5R_{\rm a}$. Thun et al. (2016) find something qualitatively similar in their analysis, and they point out that this minimum scale is the standoff distance of the bow shock, because the density field becomes markedly more asymmetric outside this distance. This implies that our drag force results could be scaled to different maximum radii using a factor of $\sim\ln\left(r_{\rm out}/0.5R_{\rm a}\right)$.

In the $\gamma=\Gamma_{\rm s}=4/3$ simulation suite, the first characteristic we note is that the overall coefficients of drag are substantially larger than those in the $\gamma=\Gamma_{\rm s}=5/3$ simulations. This is likely because the more compressible equation of state results in a higher density wake behind the embedded object, which then exerts a stronger gravitational deceleration on the embedded object’s motion. The shallow gradient cases appear to grow roughly logarithmically, with $r_{\rm min}\sim 0.5R_{\rm a}$ as in the $\gamma=\Gamma_{\rm s}=5/3$ simulations. The steepest gradient cases of the $\gamma=\Gamma_{\rm s}=4/3$ simulation suite show somewhat different behavior: the growth of $F_{\rm df}$ is superlogarithmic. Ostriker (1999)’s equation 13 shows that the logarithmic behavior comes, in part, from the constant opening angle of the Mach cone. The cases that grow more rapidly than $\ln(r_{\rm out})$ in our wind-tunnel calculation are those that show a flared wake due to the Mach number gradient discussed in Section 4.1. With a wider wake opening angle with increasing distance in these cases, the integrated dynamical friction drag grows faster than logarithmically in $r_{\rm out}$.

Finally, we summarize our results for coefficients of dynamical friction as a function of density gradient and Mach number in Figure 10 (not including the advected momentum term). We see immediately that the coefficient of drag rises with increasing density gradient by a factor of several across the span of values we have simulated. Qualitatively, what we see is that in cases of steep density gradient, the drag force depends not only on the value of the density at the embedded object’s position ($\rho_{\infty}$), but also the sweep of higher densities within the accretion radius. The flow focusses this dense material into the wake of the embedded object, causing it to contribute to the net dynamical friction. We also can observe a sharp downturn in the drag coefficient in the shallowest gradient $\gamma=\Gamma_{\rm s}=5/3$ case. This is readily explained by the low Mach number of this simulation, $\mathcal{M}=1.1$, which allows pressure to partially resist the density asymmetry of the wake. Ostriker (1999) discusses this effect extensively, and shows that the drag should behave $\propto\ln\left(1+\mathcal{M}^{-2}\right)$ in the supersonic limit, thus decreasing steeply as $\mathcal{M}\rightarrow 1$ (equation 15 in Ostriker, 1999). The lowest Mach number simulation in our $\gamma=\Gamma_{\rm s}=4/3$ suite has $\mathcal{M}=1.35$, so we do not expect (or see) as dramatic of a correction due to low flow Mach number.

We note here that the coefficients of drag and their dependence on $\epsilon_{\rho}$ derived here are different (though similar order of magnitude) from those in MacLeod & Ramirez-Ruiz (2015a, Figure 13), both because of our updated formalism and different flow parameters (section 2), and because of our corrected dynamical friction diagnostics described in section 3.2. These updates represent a significant improvement in our ability to correctly asses the dynamical friction acting on the embedded object, and the difference of our new results reflects these changes.

5 Implications for Common Envelope Inspiral

We have used idealized numerical simulations to study flow morphologies, as well as coefficients of drag and accretion for objects embedded in the common envelope. These quantities describe the transformation of an object and its orbit through the common envelope episode. Drag forces drive the orbital tightening, while flow convergence and mass accretion might transform the object itself.

There are, of course, caveats associated with the simplifications we have made here. We have isolated particular flow conditions and measured steady-state rates of drag and accretion, but it is worth considering that steady state might not be realized during the complex and violent flow of a common envelope interaction. Among the potential concerns with extrapolating the results of these simulations is that the geometry of our simulations does not match that of the large-scale common envelope: we have adopted a cartesian geometry, where stars are spherical. We similarly disregard the effects of the rotating frame that co-moves with the embedded object. These simplifications almost certainly affect the exact numerical values derived for our coefficents of drag, particularly on scales $>R_{\rm a}$, which become similar to the binary separation, $a$, the scale where curvature becomes very important. Similarly, by fixing the gas compressibility, $\gamma$, and studying two representative values of $4/3$ and $5/3$, we ignore thermodynamic transitions that might result from the gas’s passage through shocks and compression as it passes near the embedded object.

Our coefficients of accretion have dependence on the size of the sink boundary, as documented in MacLeod & Ramirez-Ruiz (2015a). These rates should thus be treated as rates of flow convergence through a boundary of a particular size: if we are considering an embedded compact object, which might be orders of magnitude smaller, it is not obvious that all of the converging material will reach the embedded object’s surface. Secondly, not all objects are thermodynamically able to accrete from the common envelope gas. Accretion onto white dwarfs or main-sequence stars has no obvious cooling channel (photons will be trapped in the very dense flow) and therefore we probably should not expect mass accumulation on these objects despite flow convergence. On the other hand, for high enough accretion rates neutrinos can likely mediate the accretion luminosity of accretion onto neutron stars (Houck & Chevalier, 1991; Chevalier, 1993; Fryer et al., 1996; Brown et al., 2000; MacLeod & Ramirez-Ruiz, 2015b), and, lacking a surface, black holes will certainly accrete material passing through their horizons.

Despite the remaining uncertainties, the coefficients of drag and accretion derived here carry lessons for the dynamics of common envelope episodes. In MacLeod & Ramirez-Ruiz (2015a, b) we argued that the ratio of drag coefficient to accretion coefficient that arises from asymmetric flows in thecommon envelope implies that objects grow by at most a few percent during their inspiral. This qualitative conclusion remains unchanged despite our improved derivations of drag and accretion coefficients.

In Figures 11 and 12, we illustrate the effect of including a coefficient of drag that varies with the flow parameters of the material that it is passing through. We use the primary-star profiles of Figure 1, and (as elsewhere in this paper), take a mass ratio $q=0.1$. We assume the primary star is initially non-rotating. We integrate the equation of motion of the secondary star relative to the enclosed mass of the primary, and add a drag force, $F_{\rm d}=C_{\rm d}\pi R_{\rm a}^{2}\rho v^{2}$, where $C_{\rm d}=F_{\rm df}/\pi R_{\rm a}^{2}\rho_{\infty}v_{\infty}^{2}$ is the coefficient of drag. We illustrate the influence of two choices: $C_{\rm d}=1$, a Hoyle-Lyttleton drag force, and a numerically-derived $C_{\rm d}$ from our simulations (which comes from Figure 10; here we take the force generated from $r<1.6R_{\rm a}$, and use the $\gamma=\Gamma_{\rm s}=5/3$ case for the $3M_{\odot}$ primary and the $\gamma=\Gamma_{\rm s}=4/3$ case for the $80M_{\odot}$ primary). Our example inspirals are initialized at $a=0.95R_{1}$ and are integrated until $a=0.1R_{1}$.

A priori, we might imagine that the initial inspiral of common envelope episodes is slow, taking many orbits while the secondary passes through the low-density atmosphere of the primary’s envelope. Instead, with the realistic coefficients, the initial common envelope inspiral is substantially more rapid than with the Hoyle-Lyttleton force alone. In the late inspiral, the drag force drops, and the orbits wrap tighter. This result can be qualitatively understood in the context of our simulations: when an embedded object lies along a steep density gradient (where the scale height is small compared to $R_{\rm a}$), the object gravitationally focusses dense material from deeper in the stellar interior into its wake. This denser material (compared to the density at the secondary’s position within the primary star) leads to a more massive wake, and a higher dynamical friction drag force. In terms of the flow streamlines shown in Figures 2 and 4, the envelope gas contributing to the wake comes largely from dense material with impact parameters in the $-y$-direction in simulation coordinates – toward the primary-star interior.

One potential impact of the increased rapidity of early inspiral is on transients from the onset of a common envelope. The emergent class of luminous red novae transients has been associated with mass ejection in stellar merger and common envelope encounters (see, e.g. Tylenda et al., 2011; Ivanova et al., 2013a; Williams et al., 2015; Kurtenkov et al., 2015; Blagorodnova et al., 2017; Smith et al., 2016; MacLeod et al., 2016, for recent examples). A rapid early inspiral would match the rapid lightcurve rise of some of these transients. For example, the M31 LRN 2015 outburst rose from detection to peak brightness in a timescale of order one binary orbital period. With a $\sim 3-5M_{\odot}$, $\sim 35R_{\odot}$ progenitor giant, this system had a primary star broadly similar to that shown in Figure 11 (Williams et al., 2015; MacLeod et al., 2016). This is a surprisingly rapid timescale when we compare to the slow early inspiral predicted by Hoyle-Lyttleton drag coefficients ($C_{\rm d}=1$), but it is perhaps more consistent with our numerically derived coefficients, which show substantial inspiral in a single orbit. There remains much work to be done, though, to establish the mappings between orbit evolution, mass ejection, and light-curve generation in these events.

The $q=0.1$ inspirals of Figure 11 and 12 differ qualitatively between the $3.0M_{\odot}$ primary and the $80M_{\odot}$ primary in the number of orbits elapsed during the inspiral. For the $3.0M_{\odot}$ primary, the secondary spirals to $a=0.1R_{1}$ in $\sim 4$ orbits, while in the $80M_{\odot}$ case, the plunge takes $\sim 13$ orbits (with the interpolated drag coefficients). This difference reflects the difference in primary-star envelope structure. The density of the $80M_{\odot}$ red supergiant envelope is very low, because radiation pressure (and the fact that the star is nearly at the Eddington limit) inflates the envelope (e.g. Sanyal et al., 2017). One consequence of this difference might be that the embedded star orbits through material that it has disturbed (or shock heated) in previous passages if the change in separation between orbits is not greater than the typical bow shock scale. In other words, when $\dot{a}P_{\rm orb}\lesssim R_{\rm a}$, we can expect that the envelope is disturbed from its initial state on subsequent orbits.444see, Livio & Soker (1988) and Iben & Livio (1993) for similar considerations based on stellar envelope profiles. In these cases, we might expect some departure from our common envelope wind tunnel flow relations. The exact extent to which this effect is important will depend on the spherical geometry of the flow (e.g. Kim & Kim, 2007b, a; Kim et al., 2008; Kim, 2010, 2011) and is difficult to asses within the context of the simulations presented here.

6 Conclusions

In this paper, we have shown characteristic relationships between the density scale height and Mach number in the common envelope based on the primary star’s structure, and we have studied three-dimensional realizations of these gas flows in an idealized “wind tunnel” setup. We draw several key conclusions from this work:

We have derived relationships for dimensionless flow scales that generically characterize common envelope flows. In particular, equations (18) and (19) relate flow Mach numbers and density gradients in terms of binary mass ratio, envelope structure, and relative velocity. 2. 2.

These relationships between flow parameters affect common envelope flow morphologies in a correlated way. Low Mach number flows tend to have mild density gradients, while high Mach number flows are always accompanied by steep gradients. 3. 3.

Density gradients in common envelope flows limit mass accretion toward the embedded objects to a fraction of $\dot{M}_{\rm HL}=\pi R_{\rm a}^{2}\rho_{\infty}v_{\infty}$, where density and velocity are defined at the location of the embedded object within the envelope. 4. 4.

Dynamical friction (gravitational) drag forces are enhanced by steep density gradients compared to the Hoyle-Lyttleton drag force, $\dot{M}_{\rm HL}v_{\infty}=\pi R_{\rm a}^{2}\rho_{\infty}v_{\infty}^{2}$, because of the contribution from dense material $(\rho\gg\rho_{\infty})$ that is focussed into the wake of the embedded object from deeper within the stellar interior (the $-y$-direction in our wind tunnel setup). These conditions are particularly relevant near stellar envelope limbs, implying more rapid orbital evolution at the onset of common envelope interactions than predicted from the Hoyle-Lyttleton force alone, with potential implications for the timescale of associated transients.

There remain many future questions to address, even in the context of simplified studies of flow within a common envelope “wind tunnel.” In future simulations, we imagine it will be particularly worthwhile to consider flow properties in cases of partial synchronization between the primary-star envelope and the secondary’s orbital motion ($f_{\rm k}<1$), the role of equation of state, and of non-accreting secondary stars. To understand open questions about the transition from dynamical plunge to subsonic motion and stabilized inspiral (Podsiadlowski, 2001; Ivanova et al., 2013b; Ivanova & Nandez, 2016; Kuruwita et al., 2016; Iaconi et al., 2017), it is likely critical to move beyond the wind tunnel formalism established here to capture the details of the passage of objects through envelope gas, which they have already perturbed. However, even shock-heated material will retain the relationships between density gradient and flow Mach number described in Section 2 if it is in (approximate) hydrostatic equilibrium.

We acknowledge helpful discussions with G. Laughlin, D. Lee, D. Lin, E. Ostriker, M. Rees, J. Samsing, and J. Stone, which helped shape and improve this work. We also acknowledge helpful feedback from anonymous referees. The software used in this work was in part developed by the DOE-supported ASCI/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. Simulation visualizations and analysis were made possible using the yt toolkit (Turk et al., 2011). This research made use of astropy, a community developed core Python package for Astronomy (Astropy Collaboration et al., 2013). The calculations for this research were carried out in part on the UCSC supercomputer Hyades, which is supported by the National Science Foundation (award number AST-1229745) and UCSC. MM is grateful for support for this work provided by NASA through Einstein Postdoctoral Fellowship grant number PF6-170155 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. AA gratefully acknowledges support from the NSF REU program LAMAT at UCSC, a UCSC Undergraduate Research in the Sciences Award, and the California Space Grant Consortium (CaSGC) Undergraduate Research Opportunity Program. AMB acknowledges UCMEXUS-CONACYT Doctoral Fellowship. PM is supported by an NSF Graduate Research Fellowship and a Eugene Cota-Robles Graduate Fellowship. ER-R acknowledges financial support from the Packard Foundation, Radcliffe Institute for Advanced Study and NASA ATP grant NNX14AH37G.