MNRAS 489, 3334–3350 (2019)

doi:10.1093/mnras/stz2312

Advance Access publication 2019 August 21

Roberto Iaconi,1,2,3‹ † Keiichi Maeda,1 Orsola De Marco,2,3 Takaya Nozawa4
and Thomas Reichardt2,3
1 Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan
2 Department of Physics & Astronomy, Macquarie University, Sydney, NSW 2109, Australia
3 Astronomy, Astrophysics and Astrophotonics Research Centre, Macquarie University, Sydney, NSW 2109, Australia
4 Division of Science, National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan

Accepted 2019 August 5. Received 2019 August 5; in original form 2019 July 8

ABSTRACT

We investigate the common envelope binary interaction, that leads to the formation of compact
binaries, such as the progenitors of Type Ia supernovae or of mergers that emit detectablegravitational waves. In this work, we diverge from the classic numerical approach that models
the dynamic inspiral. We focus instead on the asymptotic behaviour of the common envelope
expansion after the dynamic inspiral terminates. We use the SPH code PHANTOM to simulate
one of the set-ups from Passy et al., with a 0.88 M , 83 R RGB primary and a 0.6 Mcompanion, then we follow the ejecta expansion for 50 yr. Additionally, we utilize a tabulated
equation of state including the envelope recombination energy in the simulation (Reichardt
et al.), achieving a full unbinding. We show that, as time passes, the envelope’s radial velocities
dominate over the tangential ones, hence allowing us to apply an homologous expansion
kinematic model to the ejecta. The external layers of the envelope become homologous assoon as they are ejected, but it takes 5000 d (14 yr) for the bulk of the unbound gas to achieve thehomologously expanding regime. We observe that the complex distribution generated by the
dynamic inspiral evolves into a more ordered, shell-like shaped one in the asymptotic regime.
We show that the thermodynamics of the expanding envelope are in very good agreement withthose expected for an adiabatically expanding sphere under the homologous condition andgive a prediction for the location and temperature of the photosphere assuming dust to be themain source of opacity. This technique ploughs the way to determining the long-term lightbehaviour of common envelope transients.
Key words: hydrodynamics – methods: numerical – stars: AGB and post-AGB – binaries:
close – stars: evolution – dust, extinction.

1 I N T RO D U C T I O N
The common envelope interaction, hereafter CE (Paczynski 1976;
see also Ivanova et al. 2013b for a review), represents a shortevolutionary phase in the life of binary systems that takes place when
one of the components evolves to the giant stage whilst the other is amuch smaller object, such as a main sequence star or a white dwarf.
Compact binary white dwarfs, neutron stars, and black holes havelikely gone through one or more CE events during their evolution.
These systems can merge at a later time, potentially generating Type E-mail: roberto.iaconi@kusastro.kyoto-u.ac.jp† JSPS International Research Fellow (Graduate School of Science, Kyoto
University).

Ia supernovae (Webbink 1984), gamma-ray bursts (Fryer et al. 2007)
or emitting detectable gravitational waves (Abbott et al. 2016).
The CE interaction is characterized by the inspiral of the two stars
towards each other in a time-scale comparable to the dynamical
time of the giant star (1 yr), therefore it is usually named dynamic
inspiral. The dynamic inspiral is triggered by orbital instabilities,
such as the Darwin instability (Darwin 1879). This leads to a phaseof Roche lobe overflow that is unstable and eventually resultsin the dynamical inspiral, when energy and angular momentum
from the orbit are transferred to the primary’s envelope. The
physical mechanism regulating the energy exchange is primarily
the gravitational drag affecting the gas in the proximity of thecompanion star (for studies of the gravitational drag in simulationsof the CE interaction see Passy et al. 2012, Ricker & Taam
2012, Staff et al. 2016, Iaconi et al. 2017, MacLeod et al. 2017,
and Reichardt et al. 2019, while for an analytical study on theC 2019 The Author(s)Published by Oxford University Press on behalf of the Royal Astronomical Society

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

Properties of the post-inspiral common envelope ejecta – I. Dynamical
and thermal evolution

Post inspiral common envelope ejecta – I

reduces the post-dynamic inspiral expansion to a simple analyticalmodel, which can be easily evolved using initial conditions from
3D hydrodynamic simulations. This then allows us to focus on the
dynamic, thermodynamic and optical properties of the post-inspiral
CE ejecta over longer time-scales (up to 50 yr). Moreover, if thehomologous condition holds this method also allows the study ofthe ejecta on the time-scales for the formation of planetary nebulae(tens to hundreds thousands of years).
The calculation of the location and temperature of the photosphere as the envelope expands, accounting for radiation transport,
is also greatly simplified in such model and offers many advantages
compared to the approach of Galaviz et al. (2017).
This paper is structured as follows. In Section 2, we introduce
the kinematic law of the homologous expansion model and thetime dependence of the various thermodynamic quantities if thehomologous law is applied to an adiabatic expanding ideal gas. In
Section 3, we discuss our numerical choices and check whether thesimulation satisfies the basic conditions of homologous expansion.
In Section 4, we test how the simulation reproduces the analyticalmodel during its evolution, considering both the evolution of the
dynamic and thermodynamic quantities. In Section 5, we calculate
the location and temperature of the photosphere. Finally, wesummarize our findings and state our conclusions in Section 6.
2 H O M O L O G O U S E X PA N S I O N
Homologous expansion kinematic models are widely used to model
radiation transfer in supernova ejecta simulations (see e.g. Röpke
2005; Maeda, Kutsuna & Shigeyama 2014). The assumption ofhomologous expansion provides the codes with a way to determine
the kinematics of the ejecta so that most of the computational power
can be focused in solving the radiation transfer equations. Thismodel holds under the assumptions that the material is expanding
radially and that it has received an energy input that triggered the
expansion itself. After that, the amount of energy injected into andlost by the expanding material is assumed to be negligible (i.e.
adiabatic conditions). If these conditions are satisfied, then
r
r
= ,
(1)
vrad =t − t0thwhere v rad is the radial velocity, r is the radial location of a gas
parcel, t is the current time, and t0 is time at which homologous
expansion is assumed to start. We then define th = t − t0 as the‘homologous time’. The starting time of homologous expansion,t0 , can be somewhat arbitrary and we will discuss our choice inSection 4. By differentiating the previous expression it is possible
to obtain the radial location r at time th for the material located at
any initial radius ri and homologous time th, i as
r=

th
ri .
th,i

(2)

Therefore the radius at a certain time th is determined by the factor
th /th, i .
Assuming that our system is composed of an ideal polytropic gas,
for a sphere with a uniform temperature and energy, the thermalenergy (non-specific) is directly proportional to the temperature,according to the equation
3
nkB T
(3)
2
and the equation of state governing the thermodynamic variables is
Etherm =P =ρ

Etherm(γ − 1),M(4)
MNRAS 489, 3334–3350 (2019)

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

phenomenon see Ostriker 1999). The envelope is, as a result, liftedfrom the potential well of the binary and partly unbound, while
the orbital separation between the primary core and the companion
is dramatically reduced. Orbital energy is unlikely to be the only
source leading to the unbinding of the envelope. Additional sourcesof energy may contribute and probably the most discussed of them ishydrogen and helium recombination (Nandez, Ivanova & Lombardi
2015; Ivanova & Nandez 2016; Nandez & Ivanova 2016). The
efficiency of such energy source has been recently debated, with
Grichener, Sabach & Soker (2018), Soker, Grichener & Sabach(2018), and Wilson & Nordhaus (2019) arguing that the hydrogen
recombination energy might be partly radiated away and therefore
not available to do work on the envelope, while Ivanova (2018)
argue for a much higher efficiency. Other scenarios that could aid
envelope unbinding are envelope fall-back (Kuruwita, Staff & De
Marco 2016), stellar pulsations (Clayton et al. 2017), jets from thecompanion (Shiber, Kashi & Soker 2017; Shiber & Soker 2017;
Shiber et al. 2019), dust-driven winds (Glanz & Perets 2018) andconvection (Wilson & Nordhaus 2019).
Given the inherent complexity of the dynamic inspiral phase,
its treatment requires 3D hydrodynamic simulations. The firstsimulations of the common envelope interaction go back morethan three decades, though the number has increased dramaticallywithin the last decade. Notable works include Livio & Soker(1988), Terman, Taam & Hernquist (1994), Rasio & Livio (1996),
Sandquist et al. (1998), Sandquist, Taam & Burkert (2000), Ricker& Taam (2008), Ricker & Taam (2012), Passy et al. (2012), Nandez
et al. (2015), Nandez & Ivanova (2016), Ohlmann et al. (2016a),
Ohlmann et al. (2016b), Iaconi et al. (2017), Iaconi et al. (2018),
Reichardt et al. (2019), Chamandy et al. (2018), Chamandy et al.
(2019), and Prust & Chang (2019).
The numerical works cited above focus on the physics andthe immediate outcomes of the dynamical inspiral and the masstransfer phase immediately preceding it. Only those simulations
including the effects of recombination energy of the gas (Nandez
et al. 2015; Nandez & Ivanova 2016) succeed in unbinding theenvelope and only in lower mass stars (Nandez & Ivanova 2016;
Iaconi et al. 2018; Reichardt et al. in preparation). Whether or
not recombination energy can be an efficient way to remove theenvelope, alone it is not the solution to unbinding the common
envelope. Recently, a different picture is emerging: the envelope
is unbound during an extended phase that follows the dynamical
inspiral. This phase cannot be modelled by current 3D simulations,
unless some simplifying assumptions are made.
One approach to the modelling of time-scales much longer thanthe dynamical time-scale is with 1D codes. However, as shown
by Ivanova & Nandez (2016), reproducing 1D common envelopelike structures by using stellar evolution codes to deposit energy
into a giant star is challenging and produces results not particularly
similar to those obtained from 3D simulations. This caveat aside,
1D simulations of the common envelope interaction, including thelonger thermal time-scales have been carried out by, e.g. Clayton
et al. (2017), Fragos et al. (2019).
Here, we take another approach. We postulate that as theunbound envelope expands after all orbital and recombinationenergy is injected into it, its thermodynamic quantities graduallyredistribute, so that the expansion can be described by an adiabatic,
homologous kinematic model. Further motivation is provided
by the fact that a homologous expansion model can reproduce
relatively well the observed ejecta properties of red novae (e.g.
the object Vg4332; Kamiński et al. 2018), which are possible
observational counterparts of CE events (Ivanova et al. 2013a). This

3335

3336

R. Iaconi et al.

Etherm ∝ th−2 ,(5)T ∝ th−2 ,(6)P ∝ th−5 ,(7)ρ ∝ th−3 ,

(8)

and
S = const.,(9)

where S is the entropy of the gas. Based on equations (5)–(9), one
can predict the behaviour of the main thermodynamic propertiesof an adiabatically expanding sphere of ideal gas following thehomologous kinematic.
3 NUMERICAL CHOICES AND MODELSANITY CHECKS
It is easy to see how the homologous kinematic model described in
Section 2 can be used to describe supernova ejecta, given the very
short energy injection time-scale and the high radial velocity of theejecta.
Here, we apply the homologous expansion model to the unboundportion of the CE ejecta obtained from numerical simulations
and assess the fidelity of this kinematic recipe on both shortand longer time-scales. However, CE interactions present some
substantial differences compared to supernova explosions. Namely,
the injection of energy cannot be considered instantaneous and at
least initially, the velocity field generated by the dynamic inspiral
is not dominated by the radial component of the velocities (i.e.
the tangential components are not negligible). However, as we willshow below, it is still possible to approximate the post inspiral CE
ejecta by a homologous model.
3.1 Grid versus SPH codes
To test the model proposed in this work, we have used a CE SPH
simulation along the lines of those performed by Iaconi et al. (2017)
and Reichardt et al. (2019) using the Smooth Particle Hydrodynamic(SPH) code PHANTOM (Price et al. 2018).
SPH has the ability to follow the ejecta out to large radii fromthe central binary, while grid simulations have a relatively small
computational domains. SPH therefore allows us to compare boththe initial and asymptotic behaviours of the expanding layers of theenvelope with our kinematic prescription. Increasing the domain
size of a grid CE simulation so as to follow the ejecta, is not advisable
because of the effect of the ‘hot vacuum’, a numerical expedientused in grid simulation to stabilise the initial star. As discussed by
Iaconi et al. (2018) this expedient can have a dynamical effect on the
MNRAS 489, 3334–3350 (2019)

simulation for large domain sizes. Last, the ‘hot vacuum’ heats theouter layers of the ejecta, altering the thermodynamics propertiesof the photosphere (Galaviz et al. 2017).
Reichardt et al. (in preparation) have upgraded PHANTOM toinclude recombination energy by using a tabulated equation of state,similarly to the implementation of Nandez et al. (2015) and Nandez& Ivanova (2016). In this implementation, the entire recombinationenergy payload is thermalized instantly and available to do work.
With the parameters of our simulation, the entire envelope isunbound.
3.2 Numerical set-up
Similarly to Passy et al. (2012) and Iaconi et al. (2017), we havesimulated a red giant branch primary with a total mass of 0.88 M
and a radius of 83 R . The stellar core is simulated by a point
mass particle with a mass of 0.39 M . The main sequence/WD
companion is simulated by a point mass particle of 0.6 M and
it is placed on the primary’s surface in circular Keplerian orbit at
the beginning of the simulation, something that triggers immediate
inspiral. The point masses interact with the gas only gravitationallyand their gravity is smoothed by a 3 R softening length (see Priceet al. 2018 for more information on the smoothing procedure).
We use  80 000 SPH particles. This is the same as the lowest
resolution simulations carried out with identical parameters by
Reichardt et al. (2019) and was deemed sufficient by their resolution
studies for the analysis of the inspiral and ejecta parameters.
Given the longer time-scales considered in this work, increasing
the resolution further would have been prohibitive. The simulationconserves energy at the 0.1 per cent level and angular momentum
at the 0.03 per cent level.
We carried out our CE simulation using a tabulated EoS, something that includes the effects of recombination energy. Full details
on the treatment of the recombination energy will be presented by
Reichardt at al. (in preparation), but the implementation is identicalto that of Nandez et al. (2015) and Nandez & Ivanova (2016).
The binary systems are simulated from their initial configuration
through the dynamic inspiral phase, and then they are further
evolved up to 18 250 d (50 yr) after the dynamic inspiral terminates.
3.3 Energy injection time-scales
In order to apply our homologous expansion approximation, we
need to assess the time at which the CE energy sources have been
fully injected. The first source of energy is the orbital energy.
This is fully injected by the end of the inspiral, estimated using
the criterion of Sandquist et al. (1998, see also Passy et al. 2012;
Iaconi et al. 2017). This happens at 300 d after the beginning of thesimulation. At the end of the dynamic inspiral, the core-companion
separation is 20 R .
The second source of energy is the recombination energy. The
entire envelope becomes unbound by 500 d after the start of thesimulation. However, gas is still recombining at that time. By
620 d, no more recombination takes place and the energy of theunbound gas remains constant. At the end of the recombinationenergy injection phase the core-companion separation is 17 R .
3.4 Bound versus unbound ejecta
Our simulation unbinds 94 per cent of the envelope. This value has
been calculated by deeming a parcel of gas to be unbound if thesum of its potential, kinetic, and thermal energies is positive. In

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

where n is the number of particles in the system, kB is the Boltzman
constant, M is the mass of the system, T is the temperature, P is
the pressure, ρ is the density, and γ is the adiabatic index. If wecouple this with the assumption that the expansion is spherical(trivially correct in the homologous expansion case; V ∝ r3 , where
V is the volume of the sphere) and that the expansion is adiabatic(i.e. dEtherm + PdV = 0) with γ = 53 , for an homologous expanding
system it is possible to derive the proportionality laws

Post inspiral common envelope ejecta – I

3.5 Radial and tangential velocities
Homologous expansion assumes gas velocities to be only radial.
In Fig. 1, we compare the radial velocity, v rad , and the tangential
velocity v tan for all unbound SPH particles at various sample times
between 300 and 18 250 d to show that this assumption is accurate
in our case. The values are reported in Table 1.
At 300 d, when the orbital separation between the primary’s
core and the companion has stopped decreasing, only half of theenvelope has v rad > v tan . This fraction quickly increases to 86 percent by 620 d, when almost all recombination energy has been
injected in the envelope gas. By 1800 d, 98 per cent of the envelope
has a radial velocity larger than the tangential component. A more
stringent condition, namely that v rad > 10v tan , shows that the amountof envelope with a radial velocity at least 10 times greater than thetangential velocity is 11 per cent at 300 d, when the separation
has just stopped decreasing. The amount then becomes 78 percent at 1800 d. In both cases, the SPH particles with the highest
tangential velocities reside in the inner envelope, near the central
binary formed by the primary’s core and the companion, as can beseen in the second and third columns of Fig. 1.
We conclude that at 300 d homologous expansion does not
reproduce closely a substantial portion of the envelope. By 1800 d
about 20 per cent of the envelope still does not satisfy our stronger
criterion. However, after 1800 d the number of particles satisfying
our stronger criterion starts to decrease rapidly and by 5000 d
90 per cent of the particles have v rad > 10v tan . Finally, at 9125 d(equivalent to half of our total simulation time), all SPH particles
satisfy the criterion.

4 H O M O L O G O U S E X PA N S I O N I N C E E J E C TA
In this section, we apply the homologous expansion kinematic
recipe to the unbound ejecta and analyse the similarities anddifferences between the analytic model and our simulations.

4.1 Analytic versus numerical radial velocities
To verify how well the homologous expansion reproduces the
kinematics of the ejecta determined by the SPH simulation, we
compare the actual SPH particles’ radial velocity distribution withthat derived from the analytical expression of the homologous
expansion model (equation 1), at fixed homologous times.
In supernovae, the injection of energy happens over a negligible
time compared to the time-scale of the expansion of the ejecta,
therefore the initial time, t0 , can be set at any time during the energy
injection phase, without affecting how the analytical model fits theejecta kinematics. However, in the case of the CE interaction, theinjection of energy happens over a longer time and the choice of t0
must be made more carefully.

We investigate three values for t0 : the beginning of the SPH
simulation at t0 = 0 d, the end of the rapid inspiral at t0  300 d(i.e. when the orbital energy has been completely injected into theenvelope) and the moment when all recombination energy has been
injected into the envelope gas at t0  620 d. In Fig. 2, the black dots
are the unbound SPH particles, the red, cyan, and gray lines are thehomologous expansion analytic distributions with t0 = 0, 300 and
620 d, respectively.
In the case with t0 = 0 d, we observe that as time passes the
SPH particles’ distribution tends to flatten closer and closer to theanalytical distribution, marked as a red line. This is particularly
evident for the external portion of the gas distribution (i.e. the partwith radius greater than 5 au in the first panel of Fig. 2), that
gradually moves out in radius and converges with the upper portionof the analytical solution as the simulation advances. The SPH
particles closer to the core (i.e. with radius smaller than 5 au in
the first panel of Fig. 2) also move out in radius and slowly flatten
against the analytical line.
Next let us consider t0 = 300 d and t0 = 620 d. In the second andthird panels of Fig. 2, we see that the SPH particles do not tend to theanalytical distribution. However, the analytical lines pass throughthe bulk of the particle distribution and possibly could also represent
the homologous behaviour at those times. The asymptotic behaviourof the SPH particles towards the cyan and gray lines in the fourthand fifth panels of Fig. 2 shows that given enough time all the threelines tend to converge. We therefore conclude that, although theenergy injection in the case of CE is not instantaneous, on a longenough time-scale the choice of the initial time for homologous
expansion does not have a great impact on the results. We will
therefore adopt t0 = 0 d as the starting point for homologous
expansion for the remainder of this paper. As a result of this choice,the homologous time, th , coincides with the elapsed simulation time,
t. From now on, we will use the variable t to identify both simulationand homologous times.
A more detailed study of the correlation between the homologous
expansion zero-point, t0 , and the beginning of the unbinding
process would require a set of simulations starting with different
configurations tailored not to trigger the dynamic inspiral as soon as
the simulations starts. Examples of such possible configurations are
simulations starting at or immediately beyond Roche lobe overflow(e.g. Reichardt et al. 2019; Iaconi et al. 2017). However, suchconfigurations require a large amount of computational time, hence
we defer this analysis to future work.
4.2 Total mass and energy distributions
In Fig. 2, first panel, we notice a concentration of SPH particles atradii 2 AU, with radial velocities larger than the homologous
model. The concentration can be seen flattening and moving
outwards in the remaining panels of Fig. 2 with the outer edge of theconcentration moving out to radii 8, 50, 350, and 700 AU
for progressive snapshots. We marked this locations with verticaldashed lines in the panels of Fig. 2. These particles’ distribution
tends to the homologous distribution as time passes. We find thatmost of the envelope mass resides in this concentration for theentire duration of the simulation. However, we notice that slowly
the external portion of the concentration starts to spread out. Thispart of the distribution is identifiable as the portion of particles
between radii of 5 and 10 au in the second panel of Fig. 2 and25 and 60 au in the third panel. These particles eventually fully
spread out, in terms of mass and energy, mixing with the rest of
the envelope at those locations. In Fig. 3, we plot the radial density
MNRAS 489, 3334–3350 (2019)

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

this simulation, all the SPH particles tend to the same asymptotic
behaviour. In this work, we will not argue whether all the recombination energy should be allowed to do work, or whether there should
be other sources of energy at play. By using the gas distribution of
our recombination energy simulation, we are implicitly assumingthat the envelope is (almost) fully ejected because of the orbitalenergy injection and a second source of energy that operates duringor shortly after the dynamical inspiral. We therefore make the tacit
assumption that the kinematic and thermal parameters of the gas
are those that result from such assumptions.

3337

3338

R. Iaconi et al.

MNRAS 489, 3334–3350 (2019)

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

Figure 1. First column: Radial velocity versus radius for the unbound particles. Second column: Tangential velocity versus radius for the unbound particles.
Third column: Radial over tangential velocity ratio versus radius for the unbound particles; the horizontal blue lines mark v rad > v tan and v rad > 10v tan , also
note that the vertical scale is logarithmic in this column. Rows are in order of increasing time from top to bottom and the snapshots are taken at 300, 620, 1800,
9125, and 18 250 d, respectively.

Post inspiral common envelope ejecta – I
Table 1. Number, mass, and mass/Menv of the unbound SPH particles withv rad > v tan (2nd, 3rd, and 4th columns) and v rad > 10v tan (5th, 6th, and7th columns). Here, Menv is the envelope mass. The values are given at the
sample times used for Fig. 1.

(#)

300
620
1 800
9 125
18 250

34012
65247
73971
75780
75778

v rad > v tan(M )(Menv %)
0.22
0.41
0.47
0.48
0.48

45
86
98
100
100

(#)
8584
22 57759 50175 64575 729v rad > 10v tan(M ) (Menv %)
0.05
0.14
0.38
0.48
0.48

11
30
78
100
100

distributions. The concentration of SPH particles visible in Fig. 2
corresponds to the central density peak.
As the radial velocity plots in Fig. 2 show, initially the SPH
particles in the central concentration and those farther out aredecoupled, in the sense that the dynamics of the outskirts is alreadythat of a free expanding, homologous system, while the innerparticles in the concentration are still interacting with each other andtheir dynamics is moving towards that of an homologous expanding
medium. As time passes, the entire central concentration spreadsout (note that the range of the abscissa of the panels in Fig. 2
increases with time) and moves towards a more organized radial
velocity distribution, even though its density is higher than that therest of the distribution.
In the last two panels of Fig. 3, we observe that at a fixedradial location within the SPH particle concentration there are SPH
particles with a wide range of densities. This is an indication that,as the velocities redistribute, the envelope gas becomes mixed.
We see that this mixing process does not happen on the sametime-scale as the dynamical inspiral and that at t = 18 250 d it isstill partially ongoing. Therefore, the process that brings the CEejecta from an initial velocity distribution, dictated by the complex
dynamics of the inspiral to an ordered distribution such the one
described by the homologous expansion approximation, takes placeat least over tens of dynamical times. As introduced in Section 1,
the dynamical inspiral is powered by a gravitational drag force
acting primarily on a local scale approximately the size of theBondi accretion radius of the companion. This force dissipatesenergy and angular momentum from the orbit and distributes themto the envelope, also evidenced by the spiral structure imprinted
on the ejecta, with gaseous layers expanding away from thecentral binary primarily on the orbital plane. Therefore, the local
mechanism that generates a turbulent distribution in both densityand velocity, slowly moves towards a symmetric, global distribution.
In addition, recombination takes place during the first 620 d of thesimulation and adds an additional expansion force that is reasonably
symmetric (possibly more with a symmetry axis than sphericalsymmetry). This force may contribute to circularize the distributionof the envelope and promote the evolution towards a homologous
expansion.
4.3 Shell formation
Let us now consider the evolution of the mass and energy containedin the overdensity of SPH particles discussed in Section 4.2. By
9125 d, we are already in the asymptotic regime where all theenergy, both from the orbit and from the gas recombination, hasbeen delivered and the radial velocity distribution is approximately
homologous (Fig. 2). In the next 9125 d, the central particles’
overdensity stretches out to 350 au, therefore the velocity at which

this particle distribution stretches outwards is 66 km s−1 . However,
the mass contained within it, i.e. the mass enclosed within 300 au
at 9125 d or, 650 au at 18 250 d, remains approximately the same,
0.46 M , and the shape of the distribution itself only shows smallchanges. An analogous behaviour can be observed for the totalenergy, which we do not plot.
This is a result of the fact that velocities are almost entirely
redistributed at the time of the second-to-last snapshot (9125 d),following the homologous expansion kinematics. Hence, we areleft with a distribution of material containing a given mass andenergy, travelling homologously outwards. In other words, there isan overdensity of expanding material shaped like a shell, whosemean density decreases ∝t−3 (see equation 8, and discussion in
Section 4.4).
Figs 4 and 5 show a set of density slices on the orbital plane at z
= 0 and on the x–z plane (perpendicular to the orbital plane) at y =0, respectively. The slices are taken at 300 d (panel a), 620 d (panelb), 9125 d (panel c), and at 18 250 d (panel d). Panels (c) and (d) of
Fig. 4 reveal the disappearance of the spiral pattern typical of the
dynamic inspiral (clear in panels a and b), the cicularization of theenvelope ejecta and the appearance of expanding shell-like featuresthat decrease in density with time. Therefore, the various outburstsof ejecta generated by each orbit during the dynamic inspiral tend
to separate from one another, as expected from a distribution ofgas expanding homologously. In panels (d) of Figs 3 and 5, we canalso observe that not all the shells manage to perfectly circularize
before reaching the homologous regime. This is clear by looking at
the inner shell, which maintains an identical shape between 9125 d(panel c) and at 18 250 d (panel d). Note that while in panel (c) it
is possible to see three shells, the most external, low-density one
falls below the density threshold of the figure in panel (d) and istherefore no longer visible.
The overall shape of the expanding shells is elongated in the z
direction. This can be explained by the mechanics of the common
envelope ejection. The envelope gas ejected during the dynamic
inspiral is mainly concentrated along the orbital plane (Fig. 5, panela), where most of the interaction between sequentially ejected layersof gas happens. The envelope layers ejected at later times impact
with the previously ejected ones forming mild shocks (see e.g.
Ricker & Taam 2012; Iaconi et al. 2018). This interaction contributes
to slow down the overall expansion velocity on the orbital plane
and to diffuse gas along the polar direction. This gas leaves thebinary system at different angles outside the orbital plane, allowing
the formation of the shell we observe. Additionally, when moving
towards the polar direction the gas can expand more freely andmaintain a higher velocity, resulting in a shape elongated towards thez- axis.
Kamiński et al. (2018) observed the three red novae V4332 Sgr,V1309 Sco, and V838 Mon with ALMA after 22, 8, and 14 yr
from their eruptions, time-scales similar to that simulated in thiswork. The shape of the ejecta of two of them, V4332 Sgr and
V1309 Sco, appears bipolar, suggesting therefore that some levelof asymmetry such as the one produced in our simulation might
be present in systems which are possible post-CE candidates. A
more detailed comparison is unfeasible due to many uncertainties
and free parameters but it is still interesting to observe such
similarity.
The formation of bipolar structures can easily be explained bysubsequent phases of ejections where a later ejection impacts theequatorially enhanced CE outflow. This is well known to give rise
to bipolar planetary nebulae around post-CE binaries (Frank et al.
2018).
MNRAS 489, 3334–3350 (2019)

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

Time
(d)

3339

3340

R. Iaconi et al.

4.4 Adiabatic cooling of the ejecta
Let us now examine how the adiabatic condition is satisfied whenthe ejecta cool down as they approach the homologous expansion

MNRAS 489, 3334–3350 (2019)

regime. The results are shown in Fig. 6, where we plot thermalenergy (panel a), average temperature (panel b), average pressure(panel c), average density (panel d), and specific entropy (panel e)for the unbound ejecta. All the averages are arithmetical, i.e. we

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

Figure 2. Radial velocity versus radius for the unbound SPH particles at the same sample times of Fig. 1 (black dots). Shown as solid lines are the analyticalcurves representing the radial velocities predicted by the homologous expansion model, with a range of initial homologous times. The red line corresponds to
t0 = 0 d, the cyan line corresponds to t0  300 d (i.e. when the orbital energy has been completely injected into the envelope), the gray line corresponds to t0 
620 d (i.e. when both orbital and recombination energies have been completely injected into the envelope). Note that the initial homologous time corresponding
to the blue model is larger than the time of the first panel, while that corresponding to the gray model is larger than the times of the first two panels. The verticaldashed lines mark the radius containing the central overdensity of SPH particles described in Section 4.2.

Post inspiral common envelope ejecta – I

3341

summed the quantities over all the particles then divided by theparticle number (this has the effect of mass-weighing the averages
because SPH particles have equal mass). We overplotted (dashedcurves) the analytical behaviour of the various physical quantities
according to the adiabatic expansion law under homologous conditions, as shown in equations (5)–(9).

The total thermal energy of the unbound ejecta (Fig. 6, panel a)increases more and more as mass becomes progressively unbound in
the first 620 d. This increase is followed by a decrease, steeper at first
and then asymptotically flattening. We find a very good agreement
between the computed energy and the analytic expression (dashedblack line) after the end of the steep decrease at 5000 d. ThisMNRAS 489, 3334–3350 (2019)

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

Figure 3. Density versus radius for the unbound SPH particles at the same sample times of Fig. 1.

3342

R. Iaconi et al.

(b)

(c)

(d)

Figure 4. Density slices along the orbital plane at z = 0 for the recombination simulation at 300 d (a), 620 d (b), 9125 d (c), and 18 250 d (d).

indicates that the expansion becomes adiabatic after the initial partof the interaction.
The average temperature shows a behaviour in line with theenergies and a good agreement with the analytic expectation, in
accordance with the fact that temperature and energy are directly
proportional for an ideal polytropic gas (equation 3). Finally, alsopressure and density have a similar pattern as those of the previous
thermodynamic variables, since their analytic evolution can bederived from the equation of state once temperature is known.
The numerical values tend to converge with the dashed lines after
5000 d.
A special mention goes to the entropy evolution. Entropy depends
on the behaviour of all the other thermodynamic quantities discussed
above, in fact, in our model it can be derived from the fundamental
thermodynamic relation
dEtherm = T dS − P dV ,(10)where V ∝ M/ρ. Therefore entropy is very sensitive to deviations from an adiabatic behaviour. We observe that also entropy
converges towards the analytic solution after 5000 d, confirming
the previous results. We are hence witnessing an expansion thatMNRAS 489, 3334–3350 (2019)

tends to the behaviour of an ideal, polytropic gas under adiabatic
conditions.
Our simulation, started with the companion on the stellar surface,ignores the possible effects that circumstellar material ejected prior
to the dynamic inspiral might have on the entropy generation. The
envelope lost through the Lagrangian point L2 before the dynamic
inspiral (Pejcha et al. 2017, Reichardt et al. 2019) may be shocked
by the higher velocity ejecta generated by the dynamic inspiraland result in an additional entropy increase with respect to the one
observed in Fig. 6 (panel e). Such entropy increase would in any casehappen on time-scales of at least one order of magnitude smaller
than those required to reach homologous expansion and would notaffect the asymptotic evolution of the system.
The entropy evolution is also in line with what was observedby Ivanova & Nandez (2016) and Ivanova (2018), with an initial
entropy increase at the hand of the drag force energy deposition andshocks, followed by a decrease of the gas entropy while recombination takes place during the expanding phase after the dynamic
inspiral. Eventually, the gas enters the asymptotic regime, whereentropy becomes constant, which takes place when the alreadyrecombined layers of expanding gas gradually stop interacting witheach other and tend towards homologous expansion.

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019(a)

Post inspiral common envelope ejecta – I

(b)

(c)

(d)

Figure 5. Density slices along the x–z plane at y = 0 for the recombination simulation at 300 d (a), 620 d (b), 9125 days (c), and 18 250 d (d).

5 C H A R AC T E R I Z AT I O N O F T H E
P H OT O S P H E R E
In this section, we provide a simple estimate of the location of thephotosphere along the orbital plane. The calculation we perform
below could be carried out on the computational data without assuming that the envelope expands following the homologous expansion
law. However, here we couple the procedure for determining thephotosphere’s location with homologous expansion, which, as wehave shown in Sections 3 and 4, approximates the ejecta evolution
at large times after the dynamic inspiral has completed. This is toshow how using the analytical homologous pattern to model theejecta evolution opens the possibility of locating and characterizing
the photosphere and the ejecta emission properties for much longer
times than can be calculated in 3D.
A more detailed determination of the photosphere location would
require a dedicated radiation transfer simulation assuming the same
homologous pattern. We defer such study to future work.

5.1 Method
Assuming now that after 5000 d the envelope kinematics follows
the homologous expansion, we can utilize the physical quantities

from the code at 5000 d, coupled with the analytical recipe, tocreate a simple model to estimate location and temperature of thephotosphere as a function of time. We assume that the light emitted
by the central close binary system cannot be seen directly as it is wellwithin the ejecta, so that the light is reprocessed by the ejecta andemitted as a blackbody at the photosphere. Under such conditions
we can apply the Stefan–Boltzmann law to our photosphere to
determine its temperature: L  14
,
(11)
Tph =
2
4π Rphσwhere L is the central object’s luminosity, Rph is the radius of thephotosphere and σ is the Stefan–Boltzmann constant.
For the luminosity of the central object, we assumed a minimumof Lcore and a maximum of LEdd . Lcore is the luminosity of the RGB
core at the moment when the CE terminates. The luminosity of thecompanion can be neglected if the object is a main sequence staror a WD. LEdd is the Eddington luminosity of the RGB core at the
moment when the CE is initiated and, assuming accretion processeson the core by the companion, is the maximum possible luminositythat can be achieved. Both the values are taken from the MESA
model that has been mapped into PHANTOM for the hydrodynamic
simulation and correspond to Lcore  102.81 L and LEdd  104.32 L .
MNRAS 489, 3334–3350 (2019)

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

(a)

3343

3344

R. Iaconi et al.

(b)

(c)

(d)

(e)
Figure 6. Thermal energy (a), average temperature (b), average pressure (c), average density (d), and average entropy (e) versus time for the unbound particles.
Additionally, we report as dashed lines the curves representing the temporal evolution of the quantities according to the analytic derivation of Section 3. We
determined the proportionality constants of the analytic curves so as to match the numerical values at the end of the simulations. Namely, tf = 18 250 d, Etherm, f
 2.8 × 1041 erg, Tf  89.7 K, Pf  2.4 × 10−7 dyne × cm−2 , ρ f  5.8 × 10−16 g × cm−3 and Sf  1.1 × 1017 erg × K−1 g−1 , where the subscript ‘f’denotes the values from the last output of the simulations.

MNRAS 489, 3334–3350 (2019)

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019(a)

Post inspiral common envelope ejecta – I

where ti , ρ i , and ri are time, density and pressure at the moment
when homologous expansion kicks in. In our case at 5000 d. In thisway, the optical depth is a function of both the optical propertiesof the gas and of the initial homologous quantities, as a result Rph
depends on the same quantities.
For the low-temperature situation considered here, we assume
that the mass absorption coefficient κ is dominated by newly formed
dust within the CE ejecta. In the accompanying paper, we will
justify this assumption by showing that the unbound CE ejecta
satisfy the conditions for dust formation after 5000 d (Iaconi et al.,in preparation; Reichardt et al., in preparation). Similarly to thecase of dust formation in supernova ejecta, the main contributions
to the mass absorption coefficient may well come from MgSiO3(silicates) or carbon dust (Nozawa et al. 2008; Maeda et al. 2013).
We will therefore ignore any other source of opacity and focus
on the MgSiO3 and carbon dust grains. We will analyse two
possible regimes: dust component dominated by MgSiO3 and dust
component dominated by carbon. RGB star chemistry is dominatedby oxygen-based dust because the gas-phase C/O ratio is lower thanunity. However, we experiment with the case of carbon dust because
AGB stars, in many ways similar to RGB stars, can have C/O ratios
above unity and be dominated by carbon type dust. Using both types
will therefore be informative for future work.
We first verify that dust could actually form and remain stablein the expanding envelope by using the procedure of Nozawa &Kozasa (2013) and dust parameters from Nozawa et al. (2003).
We find that, at 300 d, when the dynamic inspiral terminates, thesome SPH particles cool down enough that dust can form at theoutskirts of the expanding envelope. After that, dust keeps forming
until 5000 d. The dust that forms is stable according to the criterion
from Nozawa & Kozasa (2013) we used here. At the time whenhomologous expansion sets in, dust has therefore formed over theentire envelope. This ensures that the entire amount of dust thatcan possibly form, actually forms, both in the case of MgSiO3 andcarbon. Details of the dust formation processes and calculations
will be presented elsewhere (Iaconi et al., in preparation). Assumingthat all the dust-forming metals are confined either to MgSiO3 orcarbon grains, and using a Solar composition, the amount of dust in
the whole ejecta is 2.9 × 10−4 (0.06 per cent of the envelope mass)and 2.5 × 10−3 M (0.5 per cent) in the former and latter cases,
respectively.
All the quantities at the photosphere location for the range of dust
grain sizes and as a function of time are determined by sampling
the density distribution at 5000 d into a series of concentric radial
shells (a total of 106 shells), numerically performing the integral inequation (12) on the initial distribution and by evolving over timethe value of τ according to its dependence on the homologous
factor t/ti up to 18 250 d (for a total of 104 steps). The mass
absorption coefficients necessary for the calculations are takenfrom the tables used by Nozawa et al. (2008), which report mass

absorption coefficient as a function of wavelength for four different
dust grain sizes (0.001, 0.01, 0.1, 1 μm). For each time-step, we
compute all the possible Tph from all the wavelengths available in
the tables through the process described above. The (wavelengthaveraged) photospheric temperature is then determined by using
the Wien’s displacement law, which provides a relation between
Tph and the peak wavelength of the blackbody.

5.2 Location and temperature of the photosphere
We show the estimated location of the photosphere in Fig. 7, while
in Fig. 8 we plot the envelope mass enclosed in the photosphereas a function of time. In both cases, the zig-zag appearance ofthe curves is a numerical artefact due to the simplified treatmentof the temperature determination coupled with the wavelength
discretization in the mass absorption coefficient tables.
The difference between the behaviour of MgSiO3 and carbon
is easily noticeable and depends on how the dust absorption
coefficients of the two dust types depend on wavelength, the former
showing a more complex evolution as a function of time. The
location of the photosphere is at larger radii for the carbon dust withrespect to MgSiO3 . This depends in part on the optical propertiesof the carbon dust and in part on the abundance of carbon in theenvelope gas.
If we compare the radii achieved by the particles’ distribution(Fig. 1, last two panels) with the photospheric radii for the two dust
types, we can see that in both cases the location of the photosphereis inside the region where the bulk of the SPH particles are. The gas
residing in the external tail of the radial distribution has an extremely
low density and generates a minimally contributes to the total opticaldepth. It is interesting to see that the photosphere is not at a fixed
homologous coordinate, i.e. even though its location increases overtime, it falls behind the homologous expansion velocity and moves
gradually inwards. This is especially true in the case of MgSiO3
dust, for which the increase in the photospheric radius becomes lesssteep as time passes. As a result of this the mass enclosed in thephotosphere decreases over time (Fig. 8). In the case of MgSiO3 , theinitial mass enclosed in the photospheric radius is smaller than for
carbon and the percentage decrease over time is larger for silicatesthan carbon.
In Fig. 9, we show the temperature at the photosphere location
for the 0.1 and 1 μm dust grains. We do not plot the results for the0.01 and 0.001 μm dust grains because their behaviour is identicalto the 0.1 μm case. This is because for the range of temperatures
we consider (100–400 K), the mass absorption coefficient behaves
similarly for the 0.001, 0.01, and 0.1 μm grains and as a result also
the optical depth is independent of the grain radius.
The receding photosphere discussed in the context of Fig. 8 is also
apparent in Fig. 9 by observing that the temperature evolves with ashallower slope with respect to that expected from equation (11) and1
homologous expansion, namely Tph ∝ t − 2 (black line in Fig. 9). The
1
Tph ∝ t − 2 evolution assumes that the photosphere does not movein homologous coordinates. This would be true if τ  2/3 withthe photosphere residing at the outer edge of the mass distribution.
However, the location of τ = 2/3 resides inside the envelope and,as the density decreases during the expansion, the photosphere
automatically starts to move inwards. This causes the trend we
observe in Figs 7–9.
With this in mind, we observe that for all the grain sizes the
combination of carbon dust and a central luminosity of LEdd yields
the temperature evolution closest to the homologous evolution.
MNRAS 489, 3334–3350 (2019)

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

The dependence of Tph on the the homologous expansion kinematics and on the local properties of the medium comes into playwhen determining the radius of the photosphere. Rph can be obtained
by numerically integrating the optical depth inwards, τ = κρdR,
where κ is the opacity (more accurately, κ is the mass absorption
coefficient of dust, nomenclature we will adopt instead of using
the word ‘opacity’). Rph is then obtained as the radius where thegas becomes optically thick (τ = 2/3). Assuming homologous
expansion, the optical depth equation can be rewritten as t −2 (12)κρi dri ,
τ=
ti

3345

3346

R. Iaconi et al.

(b)

Figure 7. Location of the photosphere as a function of time for two different dust types: MgSiO3 (a) and carbon (b). The calculations have been performed
for Lcore and LEdd and for four dust grain sizes, as described in the legend.

(a)

(b)

Figure 8. Envelope mass enclosed in the photosphere as a function of time for two different dust types: MgSiO3 (a) and carbon (b). The calculations havebeen performed for Lcore and LEdd and for four dust grain sizes, as described in the legend.

This also corresponds to the case where the mass enclosed inside
the photosphere decreases by the smallest amount. In the othercases, since the photosphere moves closer to the central binary, thetemperature shows a flatter decrease.
If we observe the behaviour of the density (Fig. 3), we noticethat at 9125 d, the density at the location of the photosphere tends
to be flatter inside the silicates dust photosphere with respect to thephotosphere obtained for carbon dust, this supports the difference

MNRAS 489, 3334–3350 (2019)

in recession we observe between silicates and carbon (i.e. thephotosphere in the case of silicates dust moves inwards faster thanin the case of carbon dust). The photosphere is more inside forsilicates than carbon dust and if the density distribution is flatthen the speed of the photosphere recession is higher at smaller
radii.
Typical photospheric temperatures correspond to blackbody
emission in the IR band. However, the dust mass absorption

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019(a)

Post inspiral common envelope ejecta – I(b)Figure 9. Temperature of the photosphere as a function of time plotted for two different dust grain sizes [0.1μm (a) and 1 μm (b)]. As specified in the text,
the remaining grain sizes show a behaviour exactly identical to that of panel (a), therefore have not been plotted here. For each grain size MgSiO3 dust and C
dust at both Lcore and LEdd are represented. To highlight the range of possible temperatures we shaded the respective areas between Lcore and LEdd . The black1
think line is a t − 2 curve draw for comparison.

coefficient is wavelength dependent and so is the optical depth.
As a sanity check, we performed the same integration at a fixedwavelength of 0.4 μm, in the optical band. At this wavelength,
MgSiO3 has a very low-mass absorption coefficient, while thecarbon one is very high. Indeed this results in an optically thinenvelope in the former case and in a optically thick one in thelatter (Figs 10 and 11). In the carbon case, the location of thephotosphere is in line with the various mass absorption coefficientsat different grain sizes reported in the tables. This confirms thatour procedure of integration behaves correctly. On the other hand,in case of MgSiO3 , only a fraction of the luminosity from thecentral system may indeed be absorbed within the ejecta, which
may reduce the expected temperature of the emission component
associated with the ejecta.
Our estimate of the photospheric temperature is not very precise,and a more detailed study would require a full radiation transfer
calculations. However, based on the geometry of the ejecta we
can predict features that may appear in post-dynamic inspiralobservations. In the first few years after the dynamic inspiral theejected material has a larger density on the orbital plane than inother directions. Therefore, we expect dust obscuration primarily
along the orbital plane, while we would see radiation from a more
central location or even from the binary itself when observing thesystem face-on. In this situation, the IR emission from the dusty
torus would be in addition to that emitted by the central object andresult in an IR excess. This is essentially an IR echo by the opticallythick material added to the unobscured optical/NIR light from thecentral system. As time passes the envelope distribution becomes
gradually more symmetric and the shells discussed in Section 4.3
form. We can therefore expect that the IR excess would gradually
die out and that the system would look similar when observed fromdifferent angles, once it has achieved symmetry.
There are already several observations of transients, thought toderive from common envelope mergers, that developed the signatureof newly formed dust. Examples are V 1309 Sco that displayed a

1 mag dip just before the outburst that marked the moment of themerger (Nicholls et al. 2013), V 838 Mon (Bond et al. 2003), andmore recently the luminous red nova AT 2017jfs in NGC 4470,
which is a massive counterpart of the first two (Pastorello et al.
2019). Ultimately, surveys such as the Spitzer survey SPIRITS(Kasliwal et al. 2017) and telescopes such as the dedicated Palomar
Gattini-IR (Moore et al. 2016) will allow us to discover the IR
signature of common envelope transients due to the dust that forms
in the expanding envelope.

6 S U M M A RY A N D C O N C L U S I O N S
In this work, we analysed the behaviour of the ejecta from a
CE interaction simulated with the 3D SPH hydrodynamic code
PHANTOM over a post dynamic inspiral time-scale much longer thanusually achieved. We carry out a CE simulation with the same setup used by Passy et al. (2012), namely, a primary RGB star withan initial mass of 0.88 M and an initial radius of 83 R , plus apoint-mass companion of 0.6 M . The companion is placed on the
primary surface in circular orbit at the beginning of the simulation.
Additionally, we utilize a tabulated equation of state that allows us
to include the recombination energy into the gas and achieve a fullunbinding of the envelope.
During the long post-inspiral time-scale the ejecta have time toreach a regime where they evolve without any injection or lossof energy. We propose that, since this gas does not receive any
additional energy and is unbound, the evolution of the ejecta can bemodelled as a homologously expanding ideal gas under adiabatic
conditions. Our main results are the following:(i) A key condition for the homologous expansion model torepresent the data is that the velocity field must be dominated by
the radial component of the velocities. We observe that during the
dynamic inspiral (which lasts 300 d) and until 1800 d from the
beginning of the simulation, the tangential velocities remain nonMNRAS 489, 3334–3350 (2019)

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

(a)

3347

3348

R. Iaconi et al.

(b)
Figure 10. Same as Fig. 7, but for a fixed wavelength of 0.4 μm.

(a)

(b)
Figure 11. Same as Fig. 8, but for a fixed wavelength of 0.4 μm.

negligible. However, after 1800 d the fraction of gas with high radial
velocities rapidly increases and dominates the gas distribution.
Because of this we can apply a homologous model to the numericaldata.
(ii) Since the energy injection time-scale in CE interactions is
not instantaneous as is the case in supernovae, envelope layers
ejected at different times tend to be best approximated by different
choices for the initial homologous time, t0 . This is mostly evident
for the envelope layers ejected early in the dynamic inspiral, whose
distribution rapidly converges towards the analytical relation in the

MNRAS 489, 3334–3350 (2019)radial velocity versus radius plane. We find that if t0 is chosen
during the first 620 d from the beginning of the simulation, time
during which all the possible energy is injected into the envelope,
any choice of t0 becomes a viable approximation to describe thehomologous expansion of the envelope. This becomes evident once
the average thermodynamic quantities of the envelope start evolving
in accordance to the homologous regime [5000 d, more on this in
point (iv)].
(iii) When the system achieves homologous expansion, we observe that morphology and density profiles do not change shape

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019(a)

Post inspiral common envelope ejecta – I

Concluding, we have studied a part of the CE evolution that had
not been considered previously in detail. Interestingly, a complex
phenomenon like the CE dynamic inspiral can be described with arather simple analytic model in its asymptotic regime. We here show
how this simplifies the determination of the photosphere locationand temperature. A similar approach can be applied alongsidea full radiation transfer calculations, as is done for supernovae.
This would allow us to calculate a synthetic light curve for thepost dynamic inspiral ejecta to be compared with observations oftransients, bypassing many of the obstacles encountered by theapproach of Galaviz et al. (2017).

AC K N OW L E D G E M E N T S
RI is grateful for the financial support provided by the Postodoctoral
Research Fellowship of the Japan Society for the Promotion of Science (JSPS P18753) and by the International Macquarie University
Research Excellence Scholarship. KM acknowledges the support
provided by the JSPS KAKENHI grant 18H04585, 18H05223, and
17H02864. TN acknowledges the support provided by the JSPS
KAKENHI grant 18K03707. We are grateful to the referee, Noam
Soker, for the useful comments and for helping to improve the paper.
We also thank Brian Metzger and Tomasz Kaminski for their feedback after posting this work on astro-ph. We acknowledge the computational facilities of the Department of Physics and Astronomy ofMacquarie University, on which the simulations have been carriedout. The computations described in this work were performed using
the PHANTOM code (https://phantomsph.bitbucket.io/), which isthe product of a collaborative effort of scientists under the leadof A/Prof. Daniel Price (Monash University, Melbourne, VIC,
Australia). All simulation outputs are available upon request by
e-mailing roberto.iaconi@kusastro.kyoto-u.ac.jp.
REFERENCES
Abbott B. P. et al., 2016, Phys. Rev. Lett., 116, 061102
Bond H. E. et al., 2003, Nature, 422, 405
Chamandy L. et al., 2018, MNRAS, 480, 1898
Chamandy L., Tu Y., Blackman E. G., Carroll-Nellenback J., Frank A., Liu
B., Nordhaus J., 2019, MNRAS, 486, 1070
Clayton M., Podsiadlowski P., Ivanova N., Justham S., 2017, MNRAS, 470,
1788
Darwin G. H., 1879, Proc. R. Soc. I, 30, 1
Fragos T., Andrews J. F., Ramirez-Ruiz E., Meynet G., Kalogera V., Taam
R. E., Zezas A., 2019, preprint (arXiv:1907.12573)
Frank A. et al., 2018, Galaxies, 6, 113
Fryer C. L. et al., 2007, PASP, 119, 1211
Galaviz P., De Marco O., Passy J.-C., Staff J. E., Iaconi R., 2017, ApJS, 229,
36
Glanz H., Perets H. B., 2018, MNRAS, 478, L12
Grichener A., Sabach E., Soker N., 2018, MNRAS, 478, 1818
Iaconi R. et al., 2017, MNRAS, 464, 4028
Iaconi R., De Marco O., Passy J.-C., Staff J., 2018, MNRAS, 477, 2349
Ivanova N., 2018, ApJ, 858, L24
Ivanova N., Nandez J. L. A., 2016, MNRAS, 462, 362
Ivanova N., Justham S., Avendano Nandez J. L., Lombardi J. C., 2013a,
Science, 339, 433
Ivanova N. et al., 2013b, A&AR, 21, 59
Kamiński T. et al., 2018, A&A, 617, A129
Kasliwal M. M. et al., 2017, ApJ, 839, 88
Kuruwita R. L., Staff J., De Marco O., 2016, MNRAS, 461, 486
Livio M., Soker N., 1988, ApJ, 329, 764
MacLeod M., Antoni A., Murguia-Berthier A., Macias P., Ramirez-Ruiz E.,
2017, ApJ, 838, 56
Maeda K. et al., 2013, ApJ, 776, 5
Maeda K., Kutsuna M., Shigeyama T., 2014, ApJ, 794, 37
Moore A. M. et al., 2016, in Hall H. J., Gilmozzi R., Marshall H. K., eds,
Proc. Conf. Ser. Vol. 9906, Ground-based and Airborne Telescopes VI.
SPIE, Bellingham
Nandez J. L. A., Ivanova N., 2016, MNRAS, 460, 3992
Nandez J. L. A., Ivanova N., Lombardi J. C., 2015, MNRAS, 450, L39
Nicholls C. P. et al., 2013, MNRAS, 431, L33
Nozawa T., Kozasa T., 2013, ApJ, 776, 24
Nozawa T., Kozasa T., Umeda H., Maeda K., Nomoto K., 2003, ApJ, 598,
785
Nozawa T. et al., 2008, ApJ, 684, 1343
Ohlmann S. T., Röpke F. K., Pakmor R., Springel V., 2016a, ApJ, 816,
L9
Ohlmann S. T., Röpke F. K., Pakmor R., Springel V., Müller E., 2016b,
MNRAS, 462, L121
Ostriker E. C., 1999, ApJ, 513, 252
Paczynski B., 1976, in Eggleton P., Mitton S., Whelan J., eds, Proc. IAU
Symp. 73, Structure and Evolution of Close Binary Systems. Dublin,
Ireland, p. 75
Passy J.-C. et al., 2012, ApJ, 744, 52
Pastorello A. et al., 2019, A&A, 625, L8
Pejcha O., Metzger B. D., Tyles J. G., Tomida K., 2017, ApJ, 850, 59
Price D. J. et al., 2018, PASA, 35, e031
Prust L. J., Chang P., 2019, MNRAS, 486, 5809
Rasio F. A., Livio M., 1996, ApJ, 471, 366
Reichardt T. A., De Marco O., Iaconi R., Tout C. A., Price D. J., 2019,
MNRAS, 484, 631

MNRAS 489, 3334–3350 (2019)

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019as time passes, in line with the behaviour expected from a homologously expanding system. The result of this behaviour is theformation of a series of expanding shells with decreasing averagedensities from the central binary (Fig. 4, bottom panels). Since thegas is expanding homologously, the shells tend to move apart fromeach other as time passes, with the external shells moving fasterthan the internal ones.
(iv) Thermal energy, temperature, pressure, density, and entropy
evolution initially deviate from the homologous evolution as the
dynamic inspiral and energy injection by recombination take place(in the first 620 d of the simulation). The difference persists until
5000 d, during which the bulk of the ejecta are still carrying theresiduals of the previous, turbulent interaction. After this the chaoticgas distribution becomes more ordered. After 5000 d the main
thermodynamic quantities move towards their homologous values.
Eventually all quantities evolve homologously until the end of thesimulation, showing therefore a decrease as a power law of time
for the main thermodynamic quantities (Etherm ∝ t−2 , T ∝ t−2 , P ∝t−5 , ρ ∝ t−3 , and S ∝ const.). This shows that the time-scales forthe envelope to become homologous are a factor of 10 longer thanthose of the dynamic inspiral.
(v) Utilizing the homologous model proposed here and the dust
formation model by Nozawa & Kozasa (2013), we calculate thelocation of the photosphere on the orbital plane of the ejecta startingat 5000 d after the beginning of the inspiral. We find that, irrespectiveof the dust type used (silicates and carbon dust) and of the dust grainsize, the photosphere is always located within the outer edge of thedense part of the ejecta, where the bulk of the envelope mass resides.
The photosphere does not expand homologously, but slowly moves
inwards. This happens faster for silicates than for carbon dust. From
5000 d onwards, the photosphere has temperatures between 100 and
450 K, compatible with emission in the IR band. We also highlightthat the behaviour of the photosphere does not really rely on thehomologoues expansion, and qualitatively a similar behaviour is
expected even in the phases where the homologous expansion does
not apply.

3349

3350

R. Iaconi et al.

MNRAS 489, 3334–3350 (2019)

Soker N., Grichener A., Sabach E., 2018, ApJ, 863, L14
Staff J. E., De Marco O., Wood P., Galaviz P., Passy J.-C., 2016, MNRAS,
458, 832
Terman J. L., Taam R. E., Hernquist L., 1994, ApJ, 422, 729
Webbink R. F., 1984, ApJ, 277, 355
Wilson E. C., Nordhaus J., 2019, MNRAS, 485, 4492

This paper has been typeset from a TEX/LATEX file prepared by the author.

Downloaded from https://academic.oup.com/mnras/article-abstract/489/3/3334/5552668 by National Astronomical Observatory Japan user on 26 September 2019

Ricker P. M., Taam R. E., 2008, ApJ, 672, L41
Ricker P. M., Taam R. E., 2012, ApJ, 746, 74
Röpke F. K., 2005, A&A, 432, 969
Sandquist E. L., Taam R. E., Chen X., Bodenheimer P., Burkert A., 1998,
ApJ, 500, 909
Sandquist E. L., Taam R. E., Burkert A., 2000, ApJ, 533, 984
Shiber S., Soker N., 2017, MNRAS, 477, 2584
Shiber S., Kashi A., Soker N., 2017, MNRAS, 465, L54
Shiber S., Iaconi R., De Marco O., Soker N., 2019, MNRAS, 488, 5615