11institutetext: Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str.
1, D-85748 Garching, Germany 22institutetext: Center for Computational
Relativity and Gravitation, Rochester Institute of Technology, 85 Lomb
Memorial Drive, Rochester, NY 14623, USA 33institutetext: Steward Observatory,
University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA
44institutetext: Department of Physics, Technion, Haifa, 3200003, Israel
55institutetext: Department of Astronomy & Astrophysics, University of
California, Santa Cruz, CA 95064, USA 66institutetext: Department of Physics
and Astronomy, University of North Carolina at Chapel Hill, 120 E. Cameron
Ave, Chapel Hill, NC 27599, USA

# Population synthesis of Thorne-Żytkow objects

Rejuvenated donors and unexplored progenitors in the common envelope formation
channel

K. Nathaniel 1122<EMAIL_ADDRESS>A. Vigna-Gómez 11<EMAIL_ADDRESS>A. Grichener 113344 R. Farmer 11 M. Renzo 33 R. W. Everson 5566

###### Abstract

Aims.

Methods.

Results.

###### Abstract

Context.  Common envelope evolution of a massive star and a neutron star
companion has two possible outcomes: formation of a short-period binary (a
potential gravitational wave source progenitor) or a merger of the massive
star with the neutron star. If the binary merges, a structure with a neutron
star core surrounded by a large diffuse envelope, a so-called Thorne-Żytkow
object (TŻO), may form. The predicted appearance of this hypothetical class of
star is very similar to red supergiants, making observational identification
difficult.

Aims. Our objective is to understand the properties of systems that are
potential TŻO progenitors, e.g., binary systems that enter a common envelope
phase with a neutron star companion. We also aim to distinguish those that
have been through a previous stable mass transfer phase, which can rejuvenate
the accretor.

Methods. We use the rapid population synthesis code COMPAS at Solar
metallicity and with common envelope efficiency parameter set to unity to
determine the population demographics of TŻOs. We use one-dimensional
evolutionary TŻO models from the literature to determine a fit for TŻO
lifetime in order to estimate the current number of TŻOs in the Galaxy as well
as to assess core disruption during the merger.

Results.  We explore the progenitors in the Hertzsprung-Russell diagram,
calculate formation rates, and investigate kinematics of the progenitor stars.
We find that the vast majority ($\approx$92) of TŻO progenitors in our
population have experienced mass transfer and become rejuvenated before their
formation event. Using a constant star formation rate we estimate $\approx
2\text{\times}{10}^{-4}$ TŻOs per in our Galaxy, corresponding to $\approx
5(1)$ TŻOs in the Milky Way at present.

###### Key Words.:

binaries: general – stars: neutron – stars: massive

##  1 Introduction

In the late 1930s, Landau (1938) and Oppenheimer & Volkoff (1939) first
theorized that stars could have dense neutron cores. Four decades later,
Thorne & Zytkow (1975, 1977) predicted and calculated the first structural
models of a stable configuration of a star with a neutron star (NS) core
surrounded by a diffuse envelope, later designated as a Thorne-Żytkow object
(TŻO).

TŻOs are a theoretical class of star that appear like an ordinary M supergiant
on the surface, despite the NS core hiding underneath (e.g., Thorne & Zytkow
1977; Farmer et al. 2023). There are multiple proposed formation channels: the
merger of a non-degenerate star and NS via common envelope (CE) evolution
(e.g., Taam et al., 1978; Ablimit et al., 2022), direct collision of the NS
with the stellar core via a kick (e.g., Leonard et al., 1994; Hirai &
Podsiadlowski, 2022), or dynamical mergers within a cluster or triple system
(e.g., Ray et al., 1987; Eisner et al., 2022). The collision and dynamical
formation channels are thought to have very low rates comparatively
(Podsiadlowski et al., 1995; Renzo et al., 2019; Grichener, 2023), thus the
focus of this paper is the CE channel.

TŻOs forming via CE evolution share a common channel with most NS binaries
(NS-NS) in the Galaxy (Tauris et al., 2017). The channels are the same up to
the formation of the CE; if the CE fails to eject a TŻO may form, otherwise a
NS-NS forms with successful envelope ejection. This same channel was suggested
to account for GW170817 (e.g., Vigna-Gómez et al. 2018; Kruckow et al. 2018),
the first gravitational wave signal observed by LIGO and Virgo from the
inspiral of a NS-NS (Abbott et al., 2017). An outline of the CE channel is
shown in Fig. 1, starting with the binary at the zero-age main sequence (ZAMS)
in step 1. After undergoing mass transfer (MT) through Roche lobe overflow
(RLOF) and stripping (steps 2 and 3), the first supernova (SN) occurs (step
4). At this point, when there is just one NS and one non-degenerate star (step
5), the binary has several possible outcomes. If the explosion ejects enough
mass, the natal kick received by the NS is high enough, or both, the binary
may become unbound (e.g., Blaauw, 1961; Boersma, 1961) –such is the fate of
the vast majority of massive binaries (e.g., Eldridge et al., 2011; Renzo et
al., 2019). If the binary remains bound and the orbital evolution imposes it,
a CE may follow (step 6; e.g., Paczynski 1976; Taam et al. 1978; Webbink 1984;
Iben & Livio 1993; Ivanova et al. 2013, 2020; Röpke & De Marco 2023), which is
thought to be crucial to the formation of NS-NSs (e.g., Tauris et al., 2017;
Vigna-Gómez et al., 2018; Gallegos-Garcia et al., 2023).

Figure 1: Schematic showing the CE formation channel for TŻOs, based on Fig. 1
of Vigna-Gómez et al. (2020); Vigna-Gómez & Rebagliato (2020). The chronology
and the nomenclature is the following. 1. Zero-age main sequence (ZAMS); 2.
Roche lobe overflow (RLOF); 3. Post-RLOF, or envelope detachment; 4. Supernova
(SN) of the primary; 5. Neutron star (NS) forms from the SN explosion; 6.
Common envelope evolution (CE evolution); 7. Post-CE evolution, which can lead
to either the ejection of the envelope, forming a stripped star + NS, or the
formation of a TŻO from the merger of the NS and the stellar core; 8. RLOF of
the stripped binary, further stripping the companion star; 9. Ultra-stripped
SN of the companion star; 10. NS binary. The solid arrow from step 6 indicates
the CE formation channel for TŻOs and the dashed arrow from step 8 shows the
alternative channel to a TŻO with a HeMS donor star.

During the CE phase, a non-degenerate star overflows its Roche lobe and begins
dynamically unstable MT with its NS companion. A CE forms around the NS and
the core of the companion. Drag and frictional forces in this configuration
lead to a rapid inspiral phase that decreases the orbital separation of the
binary. If the envelope is ejected before the NS reaches the core, then the
binary becomes a short-period NS + stripped-star binary (e.g. Götberg et al.,
2020), the ideal candidate for a NS-NS progenitor (steps 7 through 10 in Fig.
1). Alternatively, the NS and core can merge, resulting either in a luminous
transient (e.g., Chevalier 2012; Soker & Gilkis 2018; Soker et al. 2019; Dong
et al. 2021; Grichener & Soker 2023) or a TŻO. In this study, we use the terms
“TŻO” and “TŻO candidate” interchangeably to describe these merger products,
while remaining agnostic to the details of establishing a stable post-merger
structure.

Following the pioneering work on structural models by Thorne & Zytkow (1975,
1977), evolutionary models of TŻOs were first created by Biehle (1991), Cannon
et al. (1992), and Cannon (1993) and more recently revisited by Farmer et al.
(2023). High mass TŻOs ($\gtrsim\qty{12}{}$) have a completely convective
envelope and extremely high temperatures close to the NS (Eich et al., 1989;
Cannon, 1993). Consequently, they undergo the interrupted rapid-proton (irp)
process, which is expected to result in a unique chemical abundance profile
(Biehle, 1994; Cannon, 1993) that may possibly be used as an observational
signature. Lower mass TŻOs are supported by accretion onto the NS and are not
expected to have a discernible observational signature from red giants; low
mass TŻOs have a radiative envelope and therefore no irp-process. For a deeper
discussion of TŻO structure and chemical profile analysis, see Cannon et al.
(1992), Cannon (1993), Levesque et al. (2014), and O’Grady et al. (2020,
2023).

TŻOs are an exciting prospect in stellar physics, as they are part of an
exotic class of theoretical star-like objects. Their interiors and unique
nucleosynthetic processes are an excellent test of the capabilities of state-
of-the-art stellar evolution codes (Levesque et al., 2014; Farmer et al.,
2023). They are excellent candidates for multimessenger astronomy, as the
asymmetries and rapid rotation of the NS core are expected to emit continuous
gravitational waves (DeMarchi et al., 2021). In practice, TŻOs are very
difficult to distinguish from red (super)giants through electromagnetic only
observations (Biehle, 1991, 1994; O’Grady et al., 2020, 2023). The current
candidates, HV 2112 (Levesque et al., 2014), VX Sgr (Tabernero et al., 2021),
and HV 11417 (Beasor et al., 2018), have all been identified through unusual
chemical abundance profiles, although there is still much debate about their
classification (e.g., Coe & Pightling 1998; Vanture et al. 1999; Tout et al.
2014; Maccarone & de Mink 2016; Tabernero et al. 2021; Farmer et al. 2023).

In this study, we investigate TŻO progenitor systems in detail. Though there
are many processes in binary evolution that have been studied previously or
merit further investigation, for simplicity we focus on one important aspect:
the role of rejuvenation, which can change the stellar structure (Hellings,
1984; Renzo et al., 2023). Rejuvenation occurs when a star in a binary
experiences stable RLOF and transfers mass onto its main sequence (MS)
companion (e.g., Neo et al. 1977, Hellings 1983, 1984, see step 2 of Fig. 1).
Accretion causes the MS accretor to increase its mass, and consequently leads
to an increase in the size of the convective core. As the core mass increases,
unburnt material from outside the core is mixed through the convective
boundary, refueling the core with fresh hydrogen and rejuvenating the accretor
star. This process is not driven by rotation, as even non-rotating models can
become rejuvenated (Renzo et al., 2023; Wagg et al., 2024).

The change in structure in the core-envelope boundary region due to
rejuvenation affects the binding energy of the accretor after the MT phase.
The core-envelope boundary region contains most of the mass of the envelope,
therefore structure changes in that layer have a long-lasting impact on the
envelope, which is amplified later on by the evolutionary contraction of the
core. Renzo et al. (2023) show that the binding energy at a fixed total mass
and radius is still lower after rejuvenation than that of an unrejuvenated
star of the same mass and radius, implying that if a binary system consisting
of a NS and a rejuvenated star engages in a CE phase, it is possible that the
less bound envelope of the rejuvenated star could be ejected and result in a
short-period stripped binary instead of a TŻO. Most massive binaries leading
to NS-NSs experience stellar rejuvenation (Vigna-Gómez et al., 2018, 2020).
However, the impact of rejuvenation is not yet fully understood, and it is
conceivable that it might have contradictory effects on the formation rate of
TŻOs.

In this paper, we use rapid population synthesis in order to quantify and
characterize binary systems that enter a CE phase with a NS companion. We then
distinguish between those that merge (i.e., TŻO candidates) and those that do
not (i.e., NS-NS candidates). We determine the potential impact of
rejuvenation of the companion at the population level for NSs interacting with
non-degenerate stars. We calculate the formation rate of TŻOs per unit solar
mass () and per SN and estimate the predicted number of TŻOs in the Galaxy, as
well as explore the potential effect of rejuvenation on the formation and
merger rate of double compact objects (DCOs). In Sect. 2 we discuss our
synthetic population and how we identify and quantify systems of interest. In
Sect. 3 we present our results. We place them in context and discuss their
broader meaning in Sect. 4. We finalize with a brief summary in Sect. 5.

##  2 Methods

###  2.1 Data and initial distributions

We use the synthetic stellar populations from Grichener (2023)111Data publicly
available via zenodo (DOI 10.5281/zenodo.11237180, Grichener 2024)., generated
with the rapid binary population synthesis code
COMPAS222https://compas.science/, version 02.31.06 (Stevenson et al., 2017;
Vigna-Gómez et al., 2018; Team-COMPAS et al., 2022). For emulating single-star
evolution, COMPAS uses the analytical fits of Hurley et al. (2000), based on
the stellar models of Pols et al. (1998); for binary evolution, COMPAS follows
Hurley et al. (2002). The stellar models extend up to 50, therefore COMPAS
extrapolates for stars of higher mass. During binary evolution, COMPAS
accounts for mass loss through stellar winds, MT via RLOF, CE episodes, SNe,
and gravitational wave radiation. COMPAS version 02.31.06 does not account for
rotation, magnetic braking, nor tidal interactions. For more details, we refer
to the COMPAS methods paper (Team-COMPAS et al., 2022).

Grichener (2023) studies NSs and black holes (BHs) that merged with giant
secondaries during CE evolution using synthetic populations of $10$ million
massive binaries at various metallicities and CE efficiency parameters (see
equation 1). In this study, we use the population created with solar
metallicity ($Z=Z_{\odot}=0.0142$, Asplund et al. 2009) and CE efficiency
parameter set to unity. The setup is based on the Fiducial model from Vigna-
Gómez et al. (2018). We define the star with the larger ZAMS mass as the
primary and denote its mass by $M_{1}$ throughout the evolution. In this
synthetic population, the primary mass at the ZAMS is drawn using a Kroupa
(2001) initial mass function (IMF) in the form
$\mathrm{d}N/\mathrm{d}M_{1,\mathrm{ZAMS}}\propto M^{-2.3}_{1,\mathrm{ZAMS}}$
from 5100. The mass $M_{2}$ of the secondary star is then calculated from a
flat distribution in the initial mass ratio ($q=M_{2}/M_{1}$) between 0.1 and
1, largely consistent with observational constraints (Sana et al., 2012;
Kobulnicky et al., 2012; Offner et al., 2023). The initial separation,
$a_{\mathrm{i}}$ is drawn from a flat-in-log distribution from 0.11000au (Abt,
1983; Sana et al., 2012). The eccentricity is set to zero at ZAMS. We
elaborate on our assumption on the CE parameters in Sect. 2.1.1 and compact
object remnants in Sect. 2.1.2.

####  2.1.1 Mass transfer and common envelope evolution

MT is one of the most common ways through which stars in a binary can
interact. MT occurs when the (donor) star fills its Roche lobe and begins to
deposit mass onto its companion, the accretor star. The point at which a donor
star has filled its Roche lobe is determined using the fitting formulae for
the Roche radii from Eggleton (1983). The dynamical stability of MT, i.e.,
whether it results in stable RLOF or a CE event, is determined solely by the
stellar properties at the onset of MT. Although the analytic prescriptions
determining the stability are meant to capture the long term stability based
on the instantaneous initial conditions, it is possible for stable MT to
become unstable over time (e.g., Ivanova et al. 2020). However, this is not
currently implemented in COMPAS.

In order to assess the stability of the MT episode, we use the criterion of
Soberman et al. (1997), which compares how the radius of the donor star
changes with mass loss ($\zeta_{\mathrm{star}}\equiv\mathrm{d}\ln
R_{\mathrm{star}}/\mathrm{d}\ln M$) versus how the radius of the Roche lobe
changes with MT ($\zeta_{\mathrm{RL}}\equiv\mathrm{d}\ln
R_{\mathrm{RL}}/\mathrm{d}\ln M$). This criteria is sensitive to the amount of
mass accreted by the companion (accretor) and how much specific angular
momentum is removed from the system by the non-accreted mass. If
$\zeta_{\mathrm{star}}$ is greater than $\zeta_{\mathrm{RL}}$, then the MT is
assumed to be stable. The amount of mass removed from the donor depends on its
stellar type. For donors on the MS, mass is removed until they are able to fit
within their Roche lobe, at which point the MT episode ends. For donors with a
defined core/envelope structure (i.e., post-MS), the whole envelope is removed
on a thermal timescale. The amount of accreted mass depends on the stellar
type of the accretor. In the synthetic population, we assume no accretion onto
compact objects during CE evolution, as it does not impact the results
(Grichener, 2023). Material that is not accreted is removed from the system
with the specific angular momentum of the accretor (for additional details,
see Sect. 4.2 of Team-COMPAS et al. 2022).

Binary systems with NSs tend to experience dynamically unstable MT, which
leads to CE evolution. During CE evolution, the orbital separation can
drastically decrease within several orbits; the CE evolution can result in the
ejection of the envelope or merger. In order to estimate the outcome of a CE
episode, COMPAS uses the “$\alpha_{\mathrm{CE}}$-$\lambda$” energy formalism
to calculate the orbital separation following the CE phase (Webbink, 1984;
Livio & Soker, 1988; Iben & Livio, 1993; Ivanova et al., 2013). The energy
formalism equates the change in orbital energy ($\Delta E_{\mathrm{orb}}$)
before and after CE evolution with the binding energy of the CE
($E_{\mathrm{bind}}$) as

$E_{\mathrm{bind}}=\alpha_{\mathrm{CE}}\Delta E_{\mathrm{orb}},$ (1)

where $\alpha_{\mathrm{CE}}$ is the CE efficiency factor that parameterizes
the fraction of orbital energy that can be used to eject the CE, which is set
to unity in our simulations (i.e., $\alpha_{\mathrm{CE}}=1$). COMPAS estimates
the binding energy of the envelope by calculating the structure parameter
$\lambda$ (on the order of unity) using the “Nanjing lambda” prescriptions of
Xu & Li (2010b, a), so that

$E_{\mathrm{bind}}=-\frac{GMM_{\mathrm{env}}}{\lambda R},$ (2)

where $G$ is the gravitational constant, $M$ is the mass of the (secondary)
star, $M_{\mathrm{env}}$ is the mass of the envelope, and $R$ is the stellar
radius. For our case of a NS engulfed by a star, we write Equation 1 as

$E_{\mathrm{bind}}=\alpha_{\mathrm{CE}}\left(\frac{GMM_{\mathrm{NS}}}{2a_{\mathrm{pre-
CE}}}-\frac{GM_{\mathrm{core}}M_{\mathrm{NS}}}{2a_{\mathrm{post-CE}}}\right),$
(3)

where $M_{\mathrm{core}}$ is the core mass of the non-degenerate star at the
onset of CE evolution, $M_{\mathrm{NS}}$ is the mass of the NS, and
$a_{\mathrm{pre-CE}}$ and $a_{\mathrm{post-CE}}$ are the orbital separations
before and after CE evolution, respectively. We can then solve for the post-CE
separation by substituting Equation 2 for $E_{\mathrm{bind}}$ in Equation 3 in
order to determine if a system merges. The binary is assumed to be
circularised throughout the CE episode, and therefore the post-CE eccentricity
is zero. In our synthetic population, we define a binary system as merged if
it fits the condition $a_{\mathrm{post-CE}}\leq R_{\mathrm{core}}$. Contrary
to Grichener (2023), we take $R_{\mathrm{core}}$ as the radius of the stripped
star as computed by COMPAS. If a star enters a CE with a NS companion and
merges, we consider it to be a TŻO candidate.

####  2.1.2 Compact object remnants

A SN must occur before the CE to form the NS that makes the core of the TŻO.
COMPAS considers four SN scenarios: regular core-collapse SNe (CCSNe),
electron-capture SNe (ECSNe), ultra-stripped SNe (USSNe), and pair-instability
SNe (PISNe). It uses the pre-SN core and envelope masses to determine SN type.
Our synthetic population, at solar metallicity, does not contain any PISNe
(Langer et al., 2007; Farmer et al., 2019).

When a SN occurs, some material can fall back onto the core remnant. The
fraction of material that falls back is known as the “fallback fraction” and
determines how much material is not ejected during the SN. A fallback fraction
of unity means that all material falls back onto the remnant, there is no
ejecta and the stellar collapse is not observable. To determine the SNe rate
($N_{\mathrm{SNe}}$), we find all the CCSNe and ECSNe events with a fallback
fraction less than unity. Out of the $10$ million simulated binaries, we find
$4\,222\,994$ core-collapse events predicted to be associated with an
observable SN, and thus define
$\mathcal{R}=N/N_{\mathrm{SNe}}=N/$4\,222\,994$$, which normalizes our event
rates to the number of SNe (we do not consider Type Ia nor Type Ibc supernovae
in our estimates) in our synthetic population. We only follow systems that
remain bound after the first SN and we do not follow the evolution of unbound
stars and DCOs. Therefore we do not account for the SN explosions that occur
in these systems in our SN rate estimate. However, we estimate this to change
the total SN rate by only a factor of $\lesssim 2$.

In our population, the maximum NS mass is set to 2 (Özel et al., 2010), which
reproduces the hypothetical mass gap between NSs and BHs (although this is not
confirmed, see, e.g., Siegel et al. 2023 and Zhu et al. 2024). The NS natal-
kick distribution is bimodal Maxwellian, with CCSNe at the higher mode of 265
and ECSNe/USSNe at the lower mode of 30 (Vigna-Gómez et al., 2018). The BH
natal kicks are drawn from the same distributions and then the kick magnitude
is reduced using the fallback model from Fryer et al. (2012).

We assume 100 primordial binarity, although the previously mentioned 1000au
upper limit on the semi-major axis means that approximately half of all
binaries in the synthetic population do not interact at any point in their
evolution. If the binary has an eccentricity equal to or greater than unity
after the natal kick, the binary becomes gravitationally unbound.

###  2.2 Identifying subpopulations of interest

We filter our synthetic population in order to identify which systems have an
evolutionary channel that leads to TŻO candidate formation. Fig. 1 shows an
overview of the TŻO evolutionary channel studied in this paper, as explained
in detail in Sect. 1. Our goal is to characterize TŻO demographics, which
involves also quantifying systems with similar evolutionary pathways. To this
effect, we first identify systems consisting of a star with a NS companion. We
then distinguish between systems that later undergo CE evolution and those
that become disrupted following the SN (the unbound stars). Systems that
undergo CE evolution can result in either the ejection of their envelope (the
stripped star and NS candidates, hereafter referred to as star + NS) or
mergers of the NS with the stellar core (TŻO candidates).

As star + NS systems continue to evolve, the naked core can have its own core-
collapse event (the second one in the system) and become a NS or a BH. As
previously shown, TŻOs and DCOs have similar evolutionary channels (see, e.g.,
Belczynski et al. 2016; Vigna-Gómez et al. 2018; Broekgaarden et al. 2021;
Grichener 2023 and Fig. 1), so we also quantify the DCOs in the population,
filtering as NS-NS, BH-BH, and NS-BH; and sub-categorizing by whether the
binary will merge within the age of the Universe. The coalescence time of DCOs
via gravitational waves is based on Peters (1964). If it is less than the
Hubble time the system is a merger candidate.

Rejuvenated stars or their progenitors can be present in any of the
aforementioned subpopulations, as rejuvenation only requires stable MT onto a
MS star. When rejuvenation occurs, the internal structure of the accretor is
reorganized on a timescale comparable to a few convective turnover timescales
of the accretor’s core or the MT timescale, whichever is shorter. In most
cases, this means a timescale comparable to or shorter than the global thermal
timescale of the accretor (Renzo et al., 2023). As the thermal timescale of
the accretor is much shorter than the typical remaining lifetime of the donor
post-MT, it is likely that the accretor is in thermal equilibrium by the time
the primary collapses and undergoes a SN. We identify stellar companions in
all of our subpopulations that had accreted via stable MT while on the MS and
flag them as rejuvenated.

###  2.3 Statistical uncertainties

In order to assess the statistical uncertainty in the outcomes of our
population synthesis simulations, we use the Clopper-Pearson method (Clopper &
Pearson, 1934). This approach, also known as the exact binomial confidence
interval, is based on the cumulative distribution function of the binomial
distribution. The statistical error arises from the finite sample size of 10
million stellar binary systems, within which we identify our different
subpopulations. If we were to draw different samples of 10 million systems,
the numbers would vary due to random sampling variability. The Clopper-Pearson
intervals account for this variability by providing a range within which the
true results are likely to fall, with a specified confidence level. The
intervals were computed using the beta.ppf function from the SciPy library,
providing a robust measure of the statistical uncertainty associated with our
estimated outcome rates. Within our synthetic population, we find the
likelihood of having a statistically different number of TŻO progenitors, star
+ NS progenitors, or unbound stars under a 1 confidence level.

###  2.4 Galactic rate estimates

In order to estimate the number of TŻO candidates presently in the Milky Way,
we require the time of TŻO formation, as well as the time the TŻO remains
observable in that configuration, and the Galactic star formation rate (SFR).
For the duration in which TŻOs remain observable, i.e., the TŻO lifetime, we
use the models from Farmer et al. (2023) to approximate the TŻO lifetime as a
function of the total TŻO mass. Their models have a mass range of 520, so we
extrapolate the fit to be flat outside of this range. The resulting linear fit
is

$\tau_{\mathrm{T\dot{Z}O}}=\left(9.29m+7.15\right)\times\qty{e3}{yr},\quad
5\leq m/\unit{}\leq 20,$ (4)

where $m$ is the mass of the TŻO in solar masses. Farmer et al. (2023)
determined a TŻO lifetime of $\tau_{\mathrm{T\dot{Z}O}}\approx 10^{4-5}$
years, which is comparable or shorter to models from Cannon et al. (1992);
Cannon (1993). We calculated our fit as accurate within $\lesssim$ 6.

For the formation rate we follow the procedure outlined in Appendix B of
Chruslinska et al. (2018). In order to estimate when a TŻO progenitor is
formed, we assume a constant SFR in the Galaxy,
$\mathrm{SFR}_{\mathrm{Gal}}=2.0\pm 0.7\,\unit{\mathrm{/}}$ (Elia et al.,
2022). In practice, this means that each TŻO progenitor system is assigned a
random birth time ($\tau_{\mathrm{bin}}$) between 0Gyr and 10Gyr
($T_{\mathrm{Gal}}$). We then calculate the delay time for each TŻO
($\tau_{\mathrm{delay}}$), which is the sum of their binary progenitor birth
time and their formation time ($\tau_{\mathrm{form}}$),
$\tau_{\mathrm{delay}}=\tau_{\mathrm{bin}}+\tau_{\mathrm{form}}$. We add the
TŻO’s lifetime from Equation 4 to the delay time,
$T=\tau_{\mathrm{delay}}+\tau_{\mathrm{T\dot{Z}O}}$. If this sum is greater
than or equal to the age of the Galaxy, $T_{\mathrm{Gal}}$, then the TŻO
currently exists and we increment $n_{\mathrm{T\dot{Z}O}}$, the unnormalized
number of TŻOs in the Galaxy.

The final step to obtain estimations for the number of TŻOs in the Galaxy as
well as TŻOs per solar mass is to normalize our values to the total stellar
mass formed in the Galaxy over its lifetime. We account for the full IMF,
which includes stars less massive than the 5 lower limit in our synthetic
population. We calculate a total stellar mass of
$M_{\mathrm{sim}}=\qty{6e8}{}$, assuming that all stars are created in
binaries111As the upper limit for the initial period is 1000au, the population
includes many non-interacting binaries that mimic stars that are born and
evolve as single stars, see Sect. 2.1.2. and use a total mass range of
0.08100. The synthetic population uses a mass range of 5100, as a TŻO requires
a NS to form and therefore progenitor systems must have a massive primary that
experiences core collapse. While stars with an initial mass of $\lesssim$5
will not experience core-collapse, NS formation is still possible due to MT.
We then normalize $n_{\mathrm{T\dot{Z}O}}$ with
$n_{\mathrm{norm}}=T_{\mathrm{Gal}}\times\mathrm{SFR}_{\mathrm{Gal}}/M_{\mathrm{sim}}$.
The final number of TŻOs expected in the Galaxy is
$N_{\mathrm{T\dot{Z}O}}=n_{\mathrm{norm}}n_{\mathrm{T\dot{Z}O}}$. TŻOs per
solar mass is calculated as the raw number of TŻOs in our synthetic population
divided by $M_{\mathrm{sim}}$. We present the results in Sect. 3.1.

##  3 Results

###  3.1 Formation rates

Here we present the formation rates of systems in the synthetic population
(Fig. 2). Understanding the event rates in our population is key for
constraining formation channels and predicting the likely end product of a
binary. We find that the rate of systems in the synthetic population that
undergo MT is $\mathcal{R}\approx 1.2$, from which $\mathcal{R}\approx 0.85$
for systems that have a CE evolution phase. There are more MT events than SNe,
as the synthetic population has a lower mass limit of 5, thus not all stars in
the population will explode (see Sect. 2.1.2). While not all binaries in our
synthetic population interact with their companions, those that do are likely
to experience CE evolution.

Figure 2: Summary of the event rates in our synthetic population at Solar
metallicity ($Z=0.0142$, Asplund et al., 2009). We quantify the rate of
binaries that experience the following (from left to right): mass transfer
(MT), common envelope (CE) events (a subset of previous, can only lead to
ejections or mergers), CE ejections, and mergers. For each column, the outline
(black) indicates all binaries in the population that go through the event,
the bar (light blue) indicates binaries with a NS, and hatching indicates
companions to the NS rejuvenated by previous stable MT. The $y$-axis is in
logarithmic scale and it shows the number of systems of interest ($N$)
normalised by the number of SNe ($N_{\mathrm{SNe}}$, see Sect. 2.1.2 for more
details).

Systems with a NS that undergo CE evolution are possible progenitors for
either TŻOs or star + NS systems; in our synthetic population systems with a
NS that undergo CE evolution have a formation rate of $\mathcal{R}\approx
0.040$ (second column of Fig. 2). If we limit ourselves to systems with
rejuvenated companions, the rate only drops by 14, showing the importance of
understanding the impact of rejuvenation on massive interacting binaries.

Systems with a NS that merge with the core are defined as potential TŻOs and
we find a rate of $\mathcal{R}_{\mathrm{T\dot{Z}O}}\approx 0.029$. If the
system does not merge, it is defined as a star + NS system and we find a rate
of $\mathcal{R}_{\mathrm{star+NS}}\approx 0.011$. We also estimate the number
of TŻOs present in our Galaxy at the current moment (see Sect. 2.4). We use
our value for $\mathcal{R}_{\mathrm{T\dot{Z}O}}$ to find $\approx
2\text{\times}{10}^{-4}$ TŻOs per at solar metallicity. From this we estimate
$$\approx 5$\pm 1$ TŻOs in the Milky Way at the present.

In our synthetic population, we find that a NS-NS binary is an alternative to
TŻO formation (Fig. 1). However, due to the large uncertainties in MT prior to
the second SN explosion and mass accretion onto the NS, we believe that NS-BH
and BH-BH might also be an alternative outcome of the TŻO CE formation
channel. Fig. 3 presents the DCO rates in our synthetic population, at solar
metallicity. These systems are categorized by whether they will merge within a
Hubble time or remain a binary, with each bar in the figure representing the
total number of systems of the stated type. With the exception of BH-BH
systems that do not merge, the DCOs in our population (mergers and not) are
evenly split between staying a binary and merging due to gravitational wave
emission. NS-NS systems are the most rare, with $\mathcal{R}_{\mathrm{NS-
NS}}\approx 0.002$, and the vast majority (91) of NS-NS systems had a
rejuvenated companion before both SNe, regardless of whether or not they
merge. We find a BH-BH rate of $\mathcal{R}_{\mathrm{BH-BH}}\approx 0.004$,
with 96 of merging systems having had a rejuvenated companion before both SNe,
and only 23 of non-merging systems with a rejuvenated companion before both
SNe. This is because merging systems must interact, increasing the chances of
previous rejuvenation. Wider systems do not have this requirement and can
evolve from MS to BH without ever interacting. NS-BH systems have a rate of
$\mathcal{R}_{\mathrm{NSBH}}\approx 0.003$ and 96 had a rejuvenated companion
before the second SN.

Figure 3:  Summary of the double compact object (DCO) rates in our synthetic
population at Solar metallicity ($Z=0.0142$). We quantify the rates of DCOs
(from left to right): NS-NS binaries, BH-BH binaries, and NS-BH binaries. For
each column, the upper portion (light pink) indicates non-merging DCOs, the
lower column (dark pink) indicates DCOs that merge in a Hubble time, and the
hatching indicates DCOs that had a rejuvenated companion before both SNe. The
$y$-axis is scaled as in Fig. 2.

We find that over 90 of DCO binaries have rejuvenated stars, indicating that
including the effects of rejuvenation in COMPAS simulations and post-
processing is pertinent.

###  3.2 Progenitors in the Hertzsprung-Russell diagram

Now that we have quantified the formation rates of our subpopulations of
interest (Sect. 3.1) we focus our attention on the location of TŻO and star +
NS progenitor systems in the Hertzsprung-Russell diagram (HRD). Fig. 4
presents systems with a star and a NS at the onset of RLOF leading to CE. They
are further divided by whether they merge (TŻO candidate, light blue) or eject
the CE (star + NS candidate, dark blue). The histograms above and to the right
of the main panel show that the majority of donors in both of these
populations have experienced a previous stable MT event and become rejuvenated
(denoted by the pink lines in the histograms).

Figure 4:  HRD showing the progenitors of TŻO candidates (light blue) and star
+ NS systems (dark blue) at Solar metallicity ($Z=0.0142$). These progenitors
are shown at the onset of RLOF leading to a CE event. Two representative
single star tracks (Hurley et al., 2000) for 520 are shown as solid lines
(gray). ZAMS and the terminal-age main sequence (TAMS) are shown as dotted
lines (black). The Humphreys-Davidson limit (Humphreys & Davidson, 1994) is
shown as a dotted line (gray). Above and to the right: histograms showing the
rate density of TŻO and star + NS progenitors in their respective colors. The
outline (pink) indicates the subset of systems with a rejuvenated companion.
The $y$-axes of the histograms are scaled by the number of observable SNe in
the synthetic population (see Sect. 2.1.2 for more details).

At the onset of RLOF, TŻO donor progenitors have a mass range of ≈2≈110 and
star + NS donor progenitors have a mass range of ≈2≈68. Previous TŻO models
have had an upper mass limit of ≈60 (e.g., Cannon et al. 1992; Farmer et al.
2023; Hutchinson-Smith et al. 2023), meaning this population could potentially
indicate TŻOs more massive than have ever been modeled. We find that hotter
stellar donors always result in mergers, while cooler donors might result in
envelope ejection (e.g., Klencki et al., 2021). Examining the features of the
HRD, we see stellar donors of TŻOs can be on the MS, even close to the ZAMS,
when the CE begins. Stellar donors for star + NS systems, however, are limited
to the end of the Hertzsprung Gap (HG) and within the Asymptotic Giant Branch
(AGB). There is a horizontal gap at $10^{4.5}$ L_⊙, a numerical feature from
the fitting formulae which enhances the gap between stars that begin core He
burning while in the HG and stars that do not reach the required core
temperature until after crossing the HG into the AGB (for more details, see
Vigna-Gómez et al. 2020). There is also a lack of TŻO progenitors past the MS
for donors below 5. These donor stars have white dwarf companions rather than
NS companions. These systems may result in interesting transients, such as
late SNe (Zapartas et al., 2017).

###  3.3 Kinematics

SN explosions lead to NS (and BH) formation, and the violent explosion is
usually accompanied by a natal kick (e.g., Fryer et al., 2012; Vigna-Gómez et
al., 2018). The sudden change in the velocity of the newly born compact object
alters the orbit and changes the center-of-mass velocity ($v_{\mathrm{CoM}}$)
of the system. We investigate these center-of-mass ”spatial velocities” as
well as their relationship to the mass of the secondary (the donor) in Fig. 5.
Here we explore if there are distinct, unequivocal signatures to differentiate
between unbound stars (which become unbound after the first SN explosion and
do not undergo CE evolution) and TŻOs.

Figure 5: The mass of the secondary plotted against the center-of-mass speed
for TŻO progenitors (light blue) and unbound stars (green) on the MS (top),
post-MS (middle), and HeMS (bottom). The contours enclose 90 of the systems.
Above and to the right of the main panel are histograms showing the density,
normalized so that each integrates to unity.

Overall, the unbound stars (see step 4 of Fig. 1) dominate the parameter
space, as the TŻO rate is $\mathcal{R}_{\mathrm{T\dot{Z}O}}\approx 0.029$ and
the unbound star rate is $\mathcal{R}_{\mathrm{unbound}}\approx 0.83$, i.e.
the formation rate of unbound stars in our synthetic population is
approximately 27 times greater than that of the TŻOs. It is useful to try and
differentiate between TŻOs and unbound stars. We do so by comparing TŻO
progenitors on the MS, past the MS, and at a stage where they become stripped
stars depleted of hydrogen and rich in helium (HeMS), to unbound stars in the
same respective stages. In general, unbound stars are more massive and their
velocity has negligible dependence on the natal kick of the NS or BH formed
when the companion explodes, thus they move more slowly (Eldridge et al.,
2011; Renzo et al., 2019). TŻOs are formed from binaries that do not become
unbound when the primary undergoes a SN, meaning the magnitude of the natal
kick of the NS is likely smaller (or it is in a tight binary), and
$v_{\mathrm{CoM}}$ is strongly dependent on it (Eldridge et al., 2011; Renzo
et al., 2019; van der Meij et al., 2021).

For the MS subpopulation (top panel of Fig. 5), we find
$\mathcal{R}_{\mathrm{MS\>T\dot{Z}O}}\approx 0.0071$ for the MS TŻO donor
progenitors and $\mathcal{R}_{\mathrm{MS\>unbound}}\approx 0.71$ for the MS
unbound stars. The top panel of Fig. 5 shows nearly complete degeneracy in
donor mass and $v_{\mathrm{CoM}}$.

For the subpopulation with stars beyond the MS (middle panel of Fig. 5), we
find $\mathcal{R}_{\mathrm{post-MS\>T\dot{Z}O}}\approx 0.022$ for the TŻO
progenitors and $\mathcal{R}_{\mathrm{post-MS\>unbound}}\approx 0.045$ for the
unbound stars. The TŻO donor and unbound star progenitors show less overlap in
both mass and $v_{\mathrm{CoM}}$. The TŻO progenitors have a higher mass
range, reaching up to ≈65, while the unbound stars have a maximum mass of ≈30.
The $v_{\mathrm{CoM}}$ of the TŻO progenitors suggest they move slower than
the unbound stars, as higher velocities imply a higher chance of being
disrupted. Mass, which is strongly correlated to luminosity, is a prominent
distinguishing feature in these distributions; i.e., massive (greater than
≈30), luminous objects might serve as an observable indicator for TŻOs.

For the subpopulation with HeMS stars (bottom panel of Fig. 5), we find
$\mathcal{R}_{\mathrm{HeT\dot{Z}O}}\approx 4.0\times 10^{-5}$ for the HeMS TŻO
donors (dashed arrow in Fig. 1) and
$\mathcal{R}_{\mathrm{HeMS\>unbound}}\approx 0.008$ for the HeMS unbound
stars. This subpopulation shows strong differentiation between the TŻO and
unbound star progenitors. While the number of HeMS TŻO donor progenitors in
our synthetic population is two orders of magnitude smaller than that of the
unbound stars, they have markedly higher donor masses, lower
$v_{\mathrm{CoM}}$, and are depleted in hydrogen; these properties make them
an extremely interesting future study prospect (see Sect. 4.4 for a
discussion).

In total, we find that $v_{\mathrm{CoM}}$ and stellar mass are largely
degenerate in the subpopulation of MS TŻO donor and unbound star progenitors,
but are more clearly distinguished in the post-MS and HeMS subpopulations. The
post-MS TŻO progenitors have a much higher maximum mass, which might serve as
an observable indicator for TŻOs. The HeMS TŻO donor progenitors are small in
number, but show clear differences to the unbound stars and are an interesting
prospect for future studies. As the unbound stars have a much higher rate than
the TŻO progenitors, finding TŻOs amongst observations of massive stars is
challenging. It is likely that to indisputably identify a TŻO would require a
combination of an unusual abundance profile, kinematics, and perhaps stellar
pulsations (e.g., Levesque et al. 2014; Farmer et al. 2023).

##  4 Discussion

We have presented population demographics of binary systems that enter a CE
phase with a NS companion, using a solar metallicity-synthetic population
generated with COMPAS and publicly available at Zenodo (see Sect. 2.1). In the
discussion, we focus on the subset of the population that is predicted to
result in a merger during CE evolution, which are the potential progenitors of
TŻOs.

###  4.1 Rates

Throughout this work, we have focused on the formation rates of TŻO
candidates, particularly CE evolution in field binaries under the premise that
TŻOs can be successfully assembled through this channel (see Sect. 4.6 for a
discussion on formation from collisions following a SN). Therefore, our quoted
rates likely represent an optimistic estimate. Under this premise, we find a
TŻO formation rate per core-collapse SN of
$\mathcal{R}_{\mathrm{T\dot{Z}O}}\approx 0.029$, which is equivalent to
$\approx 2\text{\times}{10}^{-4}$ TŻOs formed per in the Galaxy. From these
values, we estimate $\approx 5(1)$ TŻOs presently in the Milky Way (Sect.
3.1). We now compare these rates with those in the literature.

Grichener (2023) also calculates rates of NS and stellar core mergers using
the same population, although they exclude progenitor systems with a MS donor
and calculate $R_{\mathrm{core}}$ using the fitting formulae of Hall & Tout
(2014) for the $a<R_{\mathrm{core}}$ merger condition. Moreover, they
normalize by number of core-collapse SNe in the synthetic population, not
filtering by fallback fraction, and find a rate of $\approx 0.011$ relative to
core-collapse SNe, approximately one third of our rate. As TŻO progenitors
with MS donors make up ≈25 of our TŻO subpopulation and the SN filtering
discussed in Sect. 2.1.2 takes out ≈7 of the SN population, our differing
rates do not impact our conclusions.

Cannon (1993) estimates $10^{7}$ TŻOs in the history of the Galaxy. Using our
value for the current number of TŻOs and the mean TŻO lifetime in our
population (8.8e4yr), we estimate $\approx 3\text{\times}{10}^{4}$ TŻOs in the
history of the Galaxy. The large discrepancy is due to the fact that Cannon
(1993) assumes that a TŻO is a later evolutionary stage for massive X-ray
binaries, with a significantly larger formation rate of ∼0.001yr (Meurs & van
den Heuvel, 1989). It is common for massive X-ray binaries to have a CE
evolution phase, but not all CE evolution phases lead to mergers, which is
illustrated by our subpopulation of star + NS systems (see Fig. 2). However,
the dichotomy between envelope ejection and stellar merger resulting in a TŻOs
is not clearly understood, and will be dependent on the evolutionary history
of the progenitor star (see Sect. 4.2 and 4.3).

Using our rate of ≈2e-4/ TŻOs and multiplying by the previously mentioned SFR
gives ≈4e-4/, which matches the formation rates estimated by Podsiadlowski et
al. (1995) and Hutilukejiang et al. (2018). They both consider the TŻO
lifetime is e5e6yr, based on the calculations by Cannon (1993) and Biehle
(1994). Under these assumptions, they both predict similar rates:
Podsiadlowski et al. (1995) predicts between $20$ and $200$ TŻOs in our Galaxy
at any given moment, while Hutilukejiang et al. (2018) determines there are
$15-150$ TŻOs. In contrast, we use the e4e5yr estimated TŻO lifetime based on
the models from Farmer et al. (2023), which results in only $5(1)$ TŻOs
presently in the Milky Way. The major discrepancy between our number and the
number found in Podsiadlowski et al. (1995) comes from the difference in
assumed TŻO lifetime. Farmer et al. (2023) analyzed the pulsational stability
of TŻOs and found modes with growing amplitudes, which lead the models to shed
their envelopes and stop being TŻOs, impacting the lifetime (see also
Romagnolo et al., in prep.).

###  4.2 Progenitors

From our synthetic population, we determine: i) in which binary star systems
there is MT from a non-degenerate star onto a NS companion, ii) if these MT
episodes are dynamically unstable and lead to CE evolution, and iii) if CE
evolution will ultimately lead to a merger and thus possibly the formation of
a TŻO. There are several uncertainties in the stellar evolution leading to the
onset of MT (RLOF) and arguably even more in the merger process and the
potential assembly of a TŻO (Langer, 2012; Pavlovskii & Ivanova, 2015; Everson
et al., 2023; Ivanova et al., 2024). In this paper, we focus on providing
insights in the properties of the TŻO stellar progenitor rather than on the
TŻOs themselves. Understanding the evolutionary pathways of TŻOs, and
providing their more common configurations will be crucial in order to assess
whether or not a TŻO can form (see Sect. 4.3).

Generally, the non-stripped stellar progenitors of TŻOs can be anywhere
between the ZAMS and the (super)giant phase, with masses between 2 and 110,
and effectively anywhere in the HRD (Fig. 4). There is, however, a slight
preference for less evolved and more massive hot donors with envelopes that
are significantly more gravitationally bound than those of their less massive
cool counterparts. In our models, there is a non-negligible fraction
($\approx$15) of progenitors with masses above 20 (middle panel of Fig. 5),
which we consider as potential candidates for TŻOs (but see Sect. 3.3 for a
discussion about unbound stars).

We find that the vast majority (92) of TŻO progenitors have donor stars which
have been rejuvenated. This means that the donor star has previously stably
accreted mass during a MT episode (Fig. 2). Here, we analyzed the synthetic
solar metallicity population from Grichener (2023). Metallicity influences the
stellar radii across the evolution through different physical effects (e.g.,
Xin et al., 2022, and references therein). The possible variation with
metallicity of the fraction of rejuvenated donors for CE with a NS remains to
be studied.

Rejuvenation alters the structure of the accretor star, especially at the
core-envelope boundary (see Sect. 1 for more details). The structural changes
are highly uncertain and are likely to vary depending on the evolutionary
stage of the star when rejuvenation occurs. Renzo et al. (2023) propose that
the altered structure could lower the binding energy of the envelope. Under
the energy formalism used throughout this work (Sect. 2.1.1), the binding
energy of the envelope must be lesser than $\alpha_{\mathrm{CE}}\Delta
E_{\mathrm{orb}}$ for envelope ejection to occur. In our TŻO progenitor
population, $E_{\mathrm{bind}}$, while large enough to prevent envelope
ejection, was less than $1.5$ times $\alpha_{\mathrm{CE}}\Delta
E_{\mathrm{orb}}$. Assessing the impact of structural changes in rejuvenated
stars will require 3D simulations (Landri et al. in prep.).

A decrease in the binding energy of the envelope of the donor could
potentially have catastrophic effects on the formation rate of TŻOs due to the
substantial increase in systems where the envelope is unbound, but would
increase the formation of binaries with a NS and a stripped star. If we assume
that TŻOs with rejuvenated stellar donors would be unable to form through the
CE evolution channel due to the decrease in binding energy, the stripped
binary formation rate could increase by a factor of four (Fig. 2). For these
stripped binaries, the non-degenerate massive star in the binary would
continue to evolve and potentially form a compact object. If the objects are
close enough, they could merge within the age of the Universe, emitting
gravitational waves detectable by ground-based observatories. Alternatively,
the formation of wide DCO binaries in the Galaxy could also be detectable by
space-based gravitational wave detectors (see Lau et al. 2020, Wagg et al.
2022). In summary, rejuvenation could potentially increase DCO formation and
merger rates in the isolated binary formation scenario (but only by a factor
of a few), and could drastically decrease the TŻO formation rate.

###  4.3 Core-disruption and thin-envelope TŻOs

Recent studies by Everson et al. (2023) and Hutchinson-Smith et al. (2023)
have once again called into question the possibility of forming TŻOs (see
also, e.g., Tout et al., 2014; Maccarone & de Mink, 2016; Beasor et al.,
2018), specifically in field binaries undergoing CE evolution. In their 3D
hydrodynamical merger simulations, shocks generated by the inspiraling compact
object propagate through the companion star and spin up the envelope and core,
such that conditions are met to form a central accretion disk upon completion
of the merger. Accretion disk formation disallows the internal structure
required to form a stable TŻO (Podsiadlowski et al., 1995), so the parameter
space study of Everson et al. (2023) sought to determine cases in which
accretion disk formation might be avoided, first looking at populations in
which the stellar core of the companion remains intact through the merger.
They provide a criterion for core disruption based on the tidal radius of the
system ($r_{\rm tidal}\approx R_{\rm core}\,q_{\rm c}^{1/3}$;
$q_{\mathrm{c}}\equiv M_{\mathrm{NS}}/M_{\mathrm{core}}$), such that the core
will be destroyed and form an accretion disk about the NS upon merger if
$r_{\rm tidal}>R_{\rm core}$. This translates to a simple metric based on the
mass ratio, $q_{\mathrm{c}}\leq 1$, for the core to remain intact. Applying
this metric, we find that ≈20 of TŻO progenitors in our synthetic population
satisfy this criterion for core disruption (Fig. 6). These systems are
preferentially low-luminosity ($\lesssim 10^{4}\,\unit{L/L_{\odot}}$) post-MS
stars. More luminous progenitors have cores massive enough that they generally
avoid core disruption; this is primarily due to the binding energy of the core
being directly proportional to the square of mass.

Figure 6:  HRD showing an alternative outcome of the CE phase (cf. Fig. 4). We
analyse our synthetic population in the context of the evolutionary models
from Everson et al. (2023). We assess the outcome of the CE evolution based on
the criteria described by Everson et al. (2023) and show systems predicted to
have an intact core (light blue), disrupted core (purple), and those that
eject the CE (dark blue). The data from Everson et al. (2023) are overplotted
as the same symbols with black outlines. We show a 30 isoradial line, the
cutoff for TETŻO formation, as a solid line (gray). The $y$-axes of the
histograms are scaled by the number of observable SNe in the synthetic
population (see Sect. 2.1.2 for more details).

Using our synthetic population, we find that a non-negligible fraction of TŻO
progenitors are more massive and luminous, covering a wider temperature range
than those modeled by Everson et al. (2023). Their models are exclusively
post-MS with an upper mass limit of $\lesssim\qty{40}{}$, lacking the
resolution to coarsely fill the crossing of the HG.

A caveat within our method with respect to the simulations of Everson et al.
(2023) is that our stellar tracks are based on fitting formulae (Hurley et
al., 2000) which do not contain the details of the structure of the core.
These simplifying assumptions mean that we can only estimate $q_{\mathrm{c}}$
and thus the outcome of NS plunge-in. One effect of these differences can be
seen in Fig. 6, where the less luminous MS progenitors in our model retain an
intact core when detailed stellar models in hydrodynamic simulations clearly
suggest they will experience core disruption. For MS donors, we assume the
core mass of the star to be approximately equal to its total mass, which
naturally results in an overestimation of the core mass. Taking a fraction of
the star’s mass as the core would increase the $q_{\mathrm{c}}$ favoring an
outcome of core disruption.

Finally, Everson et al. (2023) find that all of their models, with both
disrupted and non-disrupted cores, will form an accretion disk due to the
angular momentum deposited by the NS during inspiral and merger, therefore
preventing the assembly of stable TŻOs. From these failed TŻOs, they suggest
instead that a thin-envelope TŻO (TETŻO) may be formed after the NS plunges in
and subsequently collapses into a BH due to mass accretion. For these systems,
if the binding energy of the core released during collapse or as accretion
feedback is not sufficient to fully eject the post-merger remnant envelope, a
small fraction ($\lesssim 1\%$) of the initial envelope mass may remain as an
optically-thick thin-envelope around the BH. The thin-envelope reprocesses
accretion luminosity from the central engine and emits in the X-ray band,
possibly contributing to the population of ultraluminous X-ray sources. We
explore their metrics within our population and provide an estimate of the
number of TŻO progenitors which could become possible TETŻOs. We estimate the
formation rate of these TETŻOs by considering the progenitors in which the
companion radius is $\leq\qty{30}{}$ (below which TETŻO formation may be most
likely; see Fig. 7 of Everson et al. 2023) and estimate their lifetime
following Equation 13 from Everson et al. (2023). From our population, and
following the metrics provided by Everson et al. (2023), we obtain $\approx
1.1\text{\times}{10}^{-4}$ TETŻOs per , which corresponds to
$\mathcal{R}_{\mathrm{TET\dot{Z}O}}\approx 0.015$ and results in $1.5(5)$
TETŻOs in the Milky Way, which would form instead of the $5(1)$ TŻOs
previously discussed.

###  4.4 Helium-rich TŻOs

We find a small fraction (≈0.1 of our TŻO subpopulation) of systems in which
the donor is a HeMS star (see Sect. 3.3). These HeMS donors originate via the
same formation pathway discussed throughout this work. Following Fig. 1: after
the successful ejection (step 6) of the CE, the newly stripped HeMS star in
the binary (step 7) might later fill its Roche lobe (step 8), beginning a new
MT episode. Typically, this episode of MT proceeds in a dynamically stable
manner (Hjellming & Webbink, 1987; Soberman et al., 1997; Vigna-Gómez et al.,
2018). The CE evolution leading to HeTŻOs occurs ∼e4∼e6 after the binary’s
first CE event, well after the stripped star has had time to thermally
readjust (Vigna-Gómez et al., 2022). In our population, ≈8 of the HeMS stars
with a NS companion experience dynamically unstable MT, resulting in a HeTŻO.
We find a formation rate of
$\mathcal{R}_{\mathrm{HeT\dot{Z}O}}\approx$4.0\text{\times}{10}^{-5}$$ for
these systems, which may constitute a lower limit for the following reasons.
First, the radii of these stars may be underestimated, although the effect is
small at high metallicity, and stability criteria in this regime are highly
uncertain and untested (Laplace et al., 2020). Second, COMPAS does not account
for tidal evolution nor circumbinary disks in post-CE binaries (see, e.g.,
Valli et al., 2024, and references therein) that could lower the separation of
post-CE systems, increasing the number of mergers. Therefore, both the
possible underestimation of the radius and the effect of tides could mean we
are underestimating the number of HeTŻOs.

HeMS TŻO progenitor donor stars in our synthetic population have luminosities
between $\approx 5\times 10^{3}$ and $\approx 5\times 10^{4}$ with masses in
the range between ≈2.5≈6. With respect to their parameters as star + NS
binaries, the mass ratio ($M_{\rm{HeMS}}/M_{\rm{NS}}$) of these systems is
from $2.35.3$ and the $\zeta_{\mathrm{RL}}$ is between approximately $28$. In
COMPAS, the radius-mass exponent of HeMS stars is simplistically defined as a
fixed number ($\zeta_{\rm{HeMS}}=2$) and the accretion rate onto a NS is
determined by the Eddington limit (Team-COMPAS et al., 2022), which is why the
condition for dynamical instability ($\zeta_{\rm{RL}}>\zeta_{\rm{HeMS}}$) is
met for a fraction of systems. It is not clear if these configurations should
actually result in a CE phase. Hjellming & Webbink (1987) used polytropic
models to assess rapid MT in binaries; for radiative MS stars, they do find
that a dynamical instability might arise if the mass ratio is above the
critical value of 2.14. However, the only way to better determine the outcome
for dynamical instability would be to depart from polytropic models and
consider the structure and response of HeMS stars, as well as the orbital
evolution this uncertain MT scenario, which is beyond the scope of this work.

In Fig. 5 we show that these systems have markedly different position in the
mass and $v_{\mathrm{CoM}}$ parameter space with respect to unbound stars.
TŻOs formed through the HeMS evolutionary channel might be identifiable
through their slow spatial velocities, unusually high masses, and lack of
hydrogen spectral lines. These donor stars have a mass range of 2.55.9, which
is distinguishable from their unbound counterparts. Even if rare, these
objects form an interesting avenue for future study, for example in the
context of Luminous Fast Blue Optical Transients (Metzger 2022; Soker 2022).

###  4.5 Nucleosynthetic signatures

While the potential nucleosynthetic signature of high mass TŻOs has been
discussed (see Levesque et al. 2014), the impact of rejuvenation and HeMS
donors on the nucleosynthetic signature has not. In the case of rejuvenation,
the accreted material is rich in He and N and poor in O and C, but these
differences would not impact the outcome of the irp-process. Therefore, we do
not expect nucleosynthesis from TŻOs to be significantly different due to
previous rejuvenation.

Helium content, however, could have a noticeable impact on the abundance
profile. The irp-process in high mass TŻOs is expected to lead to production
of rubidium, strontium, yttrium, and molybdenum (Cannon, 1993) as well as an
overabundance of lithium (e.g., Podsiadlowski et al. 1995). These spectral
features have been key tools in the search for TŻO candidates (Levesque et
al., 2014; O’Grady et al., 2020, 2023). However, Farmer et al. (2023) find
that evolutionary models with very high initial helium fractions can lack Rb.
Our synthetic population predicts that a small fraction of TŻOs with high
initial helium fractions will form, in particular those with post-terminal age
MS or HeMS donors (≈0.5 of our TŻO subpopulation). These TŻOs would likely
have a distinct abundance profile, differing from the TŻOs previously
considered in the literature (Farmer et al., 2023).

###  4.6 Collisions following a Supernova

While not the main focus of this paper, TŻOs have also been predicted to form
via collisions following a SN (SN-TŻOs). In this channel, the SN kick propels
the newly-formed NS directly into the companion star, which could potentially
form a TŻO. We calculated the rate of systems in our synthetic population
where one of the stars has experienced a SN resulting in a NS, and the
companion is a non-degenerate star, as $\mathcal{R}\approx 0.032$. We find
that ≈1 of these systems have a periastron $\leq$1 directly after the SN
explosion and thus could possibly result in a TŻO candidate from the SN
collision scenario; this order-of-magnitude estimate results in a formation
rate of
$\mathcal{R}_{\mathrm{SN-T\dot{Z}O}}\approx$3.4\text{\times}{10}^{-4}$$. While
this rate is low, these could still become relevant to the TŻO population and
are worth consideration.

Leonard et al. (1994) estimated a Galactic formation rate for this channel of
0.25kpc^-2 Myr ^-1 for TŻOs with an envelope mass over 5. Using a Galactic
radius of 15kpc and extrapolating their rate to the Galactic disk gives 2e-4/.
Our previously computed formation rate of ≈4e-4/ is higher than the above
extrapolated rate. However, Leonard et al. (1994) only considers dynamically-
formed TŻOs, and notes that including the CE evolution channel could increase
the formation rate.

Podsiadlowski et al. (1995) proposed that the SN colliding channel could be of
comparable importance to the CE channel, the latter with an estimated rate of
∼e-4/. Renzo et al. (2019) optimistically estimated that ${10}^{-4}$ TŻOs
would form per binary SN events at the time of the SN forming the NS, based on
the solid angle seen by the companion. Starting with the rate of binaries
where one is a NS and the other is a non-degenerate star, using the formation
rate of $\sim{10}^{-4}$ per core-collapse events in a binary from Renzo et al.
(2019) results in $\mathcal{R}\approx$3\text{\times}{10}^{-6}$$ TŻOs formed
via collision post-SN in our synthetic population. Translating this into a
rate per year (using the previously mentioned SFR), we obtain e-4/, the same
as Podsiadlowski et al. (1995) estimated for the collision channel and lower
than the rate estimated with the CE channel. Additionally, Hirai &
Podsiadlowski (2022) studied the scenario in which the NS becomes unbound and
then collides with its companion, but concluded that such events were rare
enough to not impact previous estimates of the TŻO rate.

##  5 Summary

In this study we analyzed the Galactic metallicity, CE efficiency parameter at
unity synthetic population from Grichener (2023) with a focus on potential
progenitors of TŻOs assembled in the field via CE evolution and SN kick-
induced mergers. We used this data to explore the potential progenitors of
TŻOs assembled in the field, particularly in the context of the CE channel.
The TŻO donor progenitors found in this population have an upper mass limit of
≈110, which is significantly more massive than existing TŻO models and suggest
that TŻOs should be explored at higher mass ranges.

##### Rates.

We found a formation rate of $\approx 2\text{\times}{10}^{-4}$ TŻOs per or
$\mathcal{R}_{\mathrm{T\dot{Z}O}}=N_{\mathrm{T\dot{Z}O}}/N_{\mathrm{SNe}}\approx
0.03$, in agreement with other studies in the literature (e.g., Cannon 1993;
Podsiadlowski et al. 1995; Hutilukejiang et al. 2018). We used TŻO lifetimes
as calculated in the one-dimensional detailed stellar simulations from Farmer
et al. (2023) which are of order e4e5. We predicted $\approx 5(1)$ TŻOs
presently in the Milky Way. We found that our formation rates are overall
consistent with those in the literature and that the predicted number of
Galactic TŻOs is largely sensitive to the TŻO predicted lifetime.

##### Rejuvenation.

We found that 92 (78) of TŻO (star + NS) progenitor stellar donors have
previously accreted mass during a MT episode and have become rejuvenated. If
the lower binding energy of rejuvenated stars results in easier CE ejection,
we speculate that many of the TŻO progenitors in our population could actually
result in star + NS systems, and potentially decrease (increase) the formation
rate of TŻOs (DCO formations and mergers).

##### Core disruption and TETŻOs.

We used the criterion of Everson et al. (2023) and found that the majority of
our TŻO progenitors avoid core disruption. We compared our results to the TŻO
models from Everson et al. (2023), which suggest that excess angular momentum
deposited during CE evolution likely prevents the formation of canonical TŻOs
in field binaries, identifying unexplored regions of the parameter space
described by our synthetic population. Everson et al. (2023) present TETŻOs as
a hypothetical merger product and X-ray source that may form from failed TŻOs;
we estimated the number of TETŻOs in the Milky Way to be $1.5(5)$.

##### He-rich TŻOs.

We found a rare outcome ($\mathcal{R}_{\mathrm{HeT\dot{Z}O}}\approx 4\times
10^{-5}$) in which a Helium-rich, stripped progenitor merges with a NS
companion, a configuration that is generally not accounted for in the
literature. These configurations would likely emit a distinct spectral
signature, and therefore pose an interesting topic and merit further study.

##### Kinematics.

We compared TŻO progenitors and unbound stars in mass and spatial velocity and
found nearly complete degeneracy for progenitors stars on the MS. However,
post-MS and HeMS TŻO progenitor stars are more massive and have lower spatial
velocities than their unbound counterparts. We find that TŻO progenitors can
be difficult to distinguish from stars which become unbound after the first
SN, while also being significantly less common.

##### SN colliding TŻOs.

We briefly explored TŻOs that form when the SN kicks the newly formed NS into
the companion star. We found these systems and filtered for those with a
periastron $\leq\qty{1}{}$ directly after the SN explosion in order to
estimate a formation rate of
$\mathcal{R}_{\mathrm{SN-T\dot{Z}O}}\sim${10}^{-4}$$. We consider that it is
possible that these systems collided and merged into TŻOs.

## Acknowledgements

We thank Martyna Chruslinska, Selma de Mink, Ryosuke Hirai, Cole Johnston,
Stephen Justham, Jing-Ze Ma, Ruggero Valli, and Reinhold Willcox for useful
discussions. We thank Richard O’Shaughnessy for his support. We thank the
Kavli Foundation and the Max Planck Institute for Astrophysics for supporting
the 2023 Kavli Summer Program during which much of this work was completed. KN
thanks the LSST-DA Data Science Fellowship Program, which is funded by LSST-
DA, the Brinson Foundation, and the Moore Foundation; her participation in the
program has benefited this work. AVG acknowledges funding from the Netherlands
Organisation for Scientific Research (NWO), as part of the Vidi research
program BinWaves (project number 639.042.728, PI: de Mink). AG acknowledges
support from the Miriam and Aaron Gutwirth Fellowship and from the Gruber
Fellowship funded by the IAU and The Gruber Foundation. RWE acknowledges
support from the Heising-Simons Foundation.

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Physical Review Letters, 119, 161101
  * Ablimit et al. (2022) Ablimit, I., Podsiadlowski, P., Hirai, R., & Wicker, J. 2022, Monthly Notices of the Royal Astronomical Society, 513, 4802
  * Abt (1983) Abt, H. A. 1983, Annual Review of Astronomy and Astrophysics, 21, 343
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, Annual Review of Astronomy and Astrophysics, 47, 481
  * Beasor et al. (2018) Beasor, E. R., Davies, B., Cabrera-Ziri, I., & Hurst, G. 2018, Monthly Notices of the Royal Astronomical Society, 479, 3101
  * Belczynski et al. (2016) Belczynski, K., Holz, D. E., Bulik, T., & O’shaughnessy, R. 2016, Nature, 534, 512
  * Biehle (1991) Biehle, G. T. 1991, The Astrophysical Journal, 380, 167
  * Biehle (1994) Biehle, G. T. 1994, The Astrophysical Journal, 420, 364
  * Blaauw (1961) Blaauw, A. 1961, Bulletin of the Astronomical Institutes of the Netherlands, 15, 265
  * Boersma (1961) Boersma, J. 1961, Bulletin of the Astronomical Institutes of the Netherlands, 15, 291
  * Broekgaarden et al. (2021) Broekgaarden, F. S., Berger, E., Neijssel, C. J., et al. 2021, Monthly Notices of the Royal Astronomical Society, 508, 5028
  * Cannon (1993) Cannon, R. C. 1993, Monthly Notices of the Royal Astronomical Society, 263, 817
  * Cannon et al. (1992) Cannon, R. C., Eggleton, P. P., Zytkow, A. N., & Podsiadlowski, P. 1992, The Astrophysical Journal, 386, 206
  * Chevalier (2012) Chevalier, R. A. 2012, The Astrophysical Journal, 752, L2
  * Chruslinska et al. (2018) Chruslinska, M., Belczynski, K., Klencki, J., & Benacquista, M. 2018, Monthly Notices of the Royal Astronomical Society, 474, 2937
  * Clopper & Pearson (1934) Clopper, C. J. & Pearson, E. S. 1934, Biometrika, 26, 404
  * Coe & Pightling (1998) Coe, M. J. & Pightling, S. L. 1998, Monthly Notices of the Royal Astronomical Society, 299, 223
  * DeMarchi et al. (2021) DeMarchi, L., Sanders, J. R., & Levesque, E. M. 2021, The Astrophysical Journal, 911, 101
  * Dong et al. (2021) Dong, D. Z., Hallinan, G., Nakar, E., et al. 2021, Science, 373, 1125
  * Eggleton (1983) Eggleton, P. P. 1983, The Astrophysical Journal, 268, 368
  * Eich et al. (1989) Eich, C., Zimmermann, M. E., Thorne, K. S., & Zytkow, A. N. 1989, The Astrophysical Journal, 346, 277
  * Eisner et al. (2022) Eisner, N. L., Johnston, C., Toonen, S., et al. 2022, Monthly Notices of the Royal Astronomical Society, 511, 4710
  * Eldridge et al. (2011) Eldridge, J. J., Langer, N., & Tout, C. A. 2011, Monthly Notices of the Royal Astronomical Society, 414, 3501
  * Elia et al. (2022) Elia, D., Molinari, S., Schisano, E., et al. 2022, The Astrophysical Journal, 941, 162
  * Everson et al. (2023) Everson, R. W., Hutchinson-Smith, T., Vigna-Gómez, A., & Ramirez-Ruiz, E. 2023, Rethinking Thorne-\\.{}Z{}ytkow Object Formation: Assembly via Common Envelope in Field Binaries
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, The Astrophysical Journal, 887, 53
  * Farmer et al. (2023) Farmer, R., Renzo, M., Götberg, Y., et al. 2023, Monthly Notices of the Royal Astronomical Society, 524, 1692
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, The Astrophysical Journal, 749, 91
  * Gallegos-Garcia et al. (2023) Gallegos-Garcia, M., Berry, C. P. L., & Kalogera, V. 2023, The Astrophysical Journal, 955, 133
  * Götberg et al. (2020) Götberg, Y., Korol, V., Lamberts, A., et al. 2020, The Astrophysical Journal, 904, 56
  * Grichener (2023) Grichener, A. 2023, Monthly Notices of the Royal Astronomical Society, 523, 221
  * Grichener (2024) Grichener, A. 2024, Population Synthesis Data of Massive Binaries: NS/BH-core Mergers and Double Compact Objects Formation and Mergers
  * Grichener & Soker (2023) Grichener, A. & Soker, N. 2023, Monthly Notices of the Royal Astronomical Society, 523, 6041
  * Hall & Tout (2014) Hall, P. D. & Tout, C. A. 2014, Monthly Notices of the Royal Astronomical Society, 444, 3209
  * Hellings (1983) Hellings, P. 1983, Astrophysics and Space Science, 96, 37
  * Hellings (1984) Hellings, P. 1984, Astrophysics and Space Science, 104, 83
  * Hirai & Podsiadlowski (2022) Hirai, R. & Podsiadlowski, P. 2022, Monthly Notices of the Royal Astronomical Society, 517, 4544
  * Hjellming & Webbink (1987) Hjellming, M. S. & Webbink, R. F. 1987, The Astrophysical Journal, 318, 794
  * Humphreys & Davidson (1994) Humphreys, R. M. & Davidson, K. 1994, PASP, 106, 1025
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, Monthly Notices of the Royal Astronomical Society, 315, 543
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, Monthly Notices of the Royal Astronomical Society, 329, 897
  * Hutchinson-Smith et al. (2023) Hutchinson-Smith, T., Everson, R. W., Batta, A., et al. 2023, Rethinking Thorne-Żytkow Object Formation: The Fate of X-ray Binary LMC X-4 and Implications for Ultra-long Gamma-ray Bursts
  * Hutilukejiang et al. (2018) Hutilukejiang, B., Zhu, C., Wang, Z., & Lü, G. 2018, Journal of Astrophysics and Astronomy, 39, 21
  * Iben & Livio (1993) Iben, Jr., I. & Livio, M. 1993, Publications of the Astronomical Society of the Pacific, 105, 1373
  * Ivanova et al. (2013) Ivanova, N., Justham, S., Chen, X., et al. 2013, The Astronomy and Astrophysics Review, 21, 59
  * Ivanova et al. (2020) Ivanova, N., Justham, S., & Ricker, P. 2020, Common Envelope Evolution (IOP Publishing)
  * Ivanova et al. (2024) Ivanova, N., Kundu, S., & Pourmand, A. 2024, Unified Rapid Mass Transfer
  * Klencki et al. (2021) Klencki, J., Nelemans, G., Istrate, A. G., & Chruslinska, M. 2021, Astronomy and Astrophysics, 645, A54
  * Kobulnicky et al. (2012) Kobulnicky, H. A., Smullen, R. A., Kiminki, D. C., et al. 2012, The Astrophysical Journal, 756, 50
  * Kroupa (2001) Kroupa, P. 2001, Monthly Notices of the Royal Astronomical Society, 322, 231
  * Kruckow et al. (2018) Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, Monthly Notices of the Royal Astronomical Society, 481, 1908
  * Landau (1938) Landau, L. 1938, Nature, 141, 333
  * Langer (2012) Langer, N. 2012, Annual Review of Astronomy and Astrophysics, 50, 107
  * Langer et al. (2007) Langer, N., Norman, C. A., de Koter, A., et al. 2007, Astronomy and Astrophysics, 475, L19
  * Laplace et al. (2020) Laplace, E., Götberg, Y., de Mink, S. E., Justham, S., & Farmer, R. 2020, A&A, 637, A6
  * Lau et al. (2020) Lau, M. Y. M., Mandel, I., Vigna-Gómez, A., et al. 2020, Monthly Notices of the Royal Astronomical Society, 492, 3061
  * Leonard et al. (1994) Leonard, P. J. T., Hills, J. G., & Dewey, R. J. 1994, The Astrophysical Journal, 423, L19
  * Levesque et al. (2014) Levesque, E. M., Massey, P., Zytkow, A. N., & Morrell, N. 2014, Monthly Notices of the Royal Astronomical Society, 443, L94
  * Livio & Soker (1988) Livio, M. & Soker, N. 1988, The Astrophysical Journal, 329, 764
  * Maccarone & de Mink (2016) Maccarone, T. J. & de Mink, S. E. 2016, Monthly Notices of the Royal Astronomical Society, 458, L1
  * Metzger (2022) Metzger, B. D. 2022, The Astrophysical Journal, 932, 84
  * Meurs & van den Heuvel (1989) Meurs, E. J. A. & van den Heuvel, E. P. J. 1989, Astronomy and Astrophysics, 226, 88
  * Neo et al. (1977) Neo, S., Miyaji, S., Nomoto, K., & Sugimoto, D. 1977, Publications of the Astronomical Society of Japan, 29, 249
  * Offner et al. (2023) Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 275
  * O’Grady et al. (2023) O’Grady, A. J. G., Drout, M. R., Gaensler, B. M., et al. 2023, The Astrophysical Journal, 943, 18
  * O’Grady et al. (2020) O’Grady, A. J. G., Drout, M. R., Shappee, B. J., et al. 2020, The Astrophysical Journal, 901, 135
  * Oppenheimer & Volkoff (1939) Oppenheimer, J. R. & Volkoff, G. M. 1939, Physical Review, 55, 374
  * Özel et al. (2010) Özel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010, The Astrophysical Journal, 725, 1918
  * Paczynski (1976) Paczynski, B. 1976, 73, 75
  * Pavlovskii & Ivanova (2015) Pavlovskii, K. & Ivanova, N. 2015, Monthly Notices of the Royal Astronomical Society, 449, 4415
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224
  * Podsiadlowski et al. (1995) Podsiadlowski, P., Cannon, R. C., & Rees, M. J. 1995, Monthly Notices of the Royal Astronomical Society, 274, 485
  * Pols et al. (1998) Pols, O. R., Schröder, K.-P., Hurley, J. R., Tout, C. A., & Eggleton, P. P. 1998, Monthly Notices of the Royal Astronomical Society, 298, 525
  * Ray et al. (1987) Ray, A., Kembhavi, A. K., & Antia, H. M. 1987, Astronomy and Astrophysics, 184, 164
  * Renzo et al. (2019) Renzo, M., Zapartas, E., de Mink, S. E., et al. 2019, Astronomy and Astrophysics, 624, A66
  * Renzo et al. (2023) Renzo, M., Zapartas, E., Justham, S., et al. 2023, The Astrophysical Journal, 942, L32
  * Romagnolo et al. (in prep.) Romagnolo, A., Klencki, J., Vigna-Gómez, A., & Belczynski, K. in prep.
  * Röpke & De Marco (2023) Röpke, F. K. & De Marco, O. 2023, Living Reviews in Computational Astrophysics, 9, 2
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444
  * Siegel et al. (2023) Siegel, J. C., Kiato, I., Kalogera, V., et al. 2023, The Astrophysical Journal, 954, 212
  * Soberman et al. (1997) Soberman, G. E., Phinney, E. S., & van den Heuvel, E. P. J. 1997, Astronomy and Astrophysics, 327, 620
  * Soker (2022) Soker, N. 2022, Research in Astronomy and Astrophysics, 22, 055010
  * Soker & Gilkis (2018) Soker, N. & Gilkis, A. 2018, Monthly Notices of the Royal Astronomical Society, 475, 1198
  * Soker et al. (2019) Soker, N., Grichener, A., & Gilkis, A. 2019, Monthly Notices of the Royal Astronomical Society, 484, 4972
  * Stevenson et al. (2017) Stevenson, S., Vigna-Gómez, A., Mandel, I., et al. 2017, Nature Communications, 8, 14906
  * Taam et al. (1978) Taam, R. E., Bodenheimer, P., & Ostriker, J. P. 1978, The Astrophysical Journal, 222, 269
  * Tabernero et al. (2021) Tabernero, H. M., Dorda, R., Negueruela, I., & Marfil, E. 2021, Astronomy and Astrophysics, 646, A98
  * Tauris et al. (2017) Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, The Astrophysical Journal, 846, 170
  * Team-COMPAS et al. (2022) Team-COMPAS, Riley, J., Agrawal, P., et al. 2022, The Journal of Open Source Software, 7, 3838
  * Thorne & Zytkow (1975) Thorne, K. S. & Zytkow, A. N. 1975, The Astrophysical Journal, 199, L19
  * Thorne & Zytkow (1977) Thorne, K. S. & Zytkow, A. N. 1977, The Astrophysical Journal, 212, 832
  * Tout et al. (2014) Tout, C. A., Zytkow, A. N., Church, R. P., et al. 2014, Monthly Notices of the Royal Astronomical Society, 445, L36
  * Valli et al. (2024) Valli, R., Tiede, C., Vigna-Gómez, A., et al. 2024, arXiv e-prints, arXiv:2401.17355
  * van der Meij et al. (2021) van der Meij, V., Guo, D., Kaper, L., & Renzo, M. 2021, A&A, 655, A31
  * Vanture et al. (1999) Vanture, A. D., Zucker, D., & Wallerstein, G. 1999, The Astrophysical Journal, 514, 932
  * Vigna-Gómez et al. (2020) Vigna-Gómez, A., MacLeod, M., Neijssel, C. J., et al. 2020, Publications of the Astronomical Society of Australia, 37, e038
  * Vigna-Gómez et al. (2018) Vigna-Gómez, A., Neijssel, C. J., Stevenson, S., et al. 2018, Monthly Notices of the Royal Astronomical Society, 481, 4009
  * Vigna-Gómez & Rebagliato (2020) Vigna-Gómez, A. & Rebagliato, T. 2020, Graphics: Evolutionary Pathways Leading to Double Neutron Star Formation
  * Vigna-Gómez et al. (2022) Vigna-Gómez, A., Wassink, M., Klencki, J., et al. 2022, MNRAS, 511, 2326
  * Wagg et al. (2022) Wagg, T., Broekgaarden, F. S., de Mink, S. E., et al. 2022, The Astrophysical Journal, 937, 118
  * Wagg et al. (2024) Wagg, T., Johnston, C., Bellinger, E. P., et al. 2024, The Asteroseismic Imprints of Mass Transfer: A Case Study of a Binary Mass Gainer in the SPB Instability Strip
  * Webbink (1984) Webbink, R. F. 1984, The Astrophysical Journal, 277, 355
  * Xin et al. (2022) Xin, C., Renzo, M., & Metzger, B. D. 2022, MNRAS, 516, 5816
  * Xu & Li (2010a) Xu, X.-J. & Li, X.-D. 2010a, The Astrophysical Journal, 722, 1985
  * Xu & Li (2010b) Xu, X.-J. & Li, X.-D. 2010b, The Astrophysical Journal, 716, 114
  * Zapartas et al. (2017) Zapartas, E., de Mink, S. E., Izzard, R. G., et al. 2017, Astronomy and Astrophysics, 601, A29
  * Zhu et al. (2024) Zhu, J.-P., Qin, Y., Wang, Z.-H.-T., et al. 2024, Monthly Notices of the Royal Astronomical Society, 529, 4554