Neutron Star-Main-Sequence Star Collision: Formation and Long-Term Evolution of Thorne-Żytkow Objects

Introduction

Thorne–Żytkow objects (TZOs) are hypothetical stars consisting of a neutron-star (NS) or black-hole (BH) core embedded in a stellar envelope ([54, 55]), whose structure is supported by either the accretion power of the NS or nuclear burning ([5]). Observationally, massive TZOs are expected to have unique chemical signatures at the surface because of the rapid proton capture process at the base of the convective envelope ([9, 46]). Accordingly, some evolved stars with peculiar chemical abundance have been proposed to be the candidates of TZOs (e.g., HV2112; [56]), whose nature remains under debate ([4, 23]). Identifying a TZO from its surface properties relies on predictions concerning whether such an object can form and how it subsequently evolves. Recent studies have questioned whether the merger of an NS and an MS star can produce the classical TZO ([22, 28, 15 a]), and evolutionary calculations find that the surface properties of TZOs may differ from the classical expectations ([23]). Thus, establishing robust observational signatures requires models that connect the formation of a TZO and its long-term evolution.

Many scenarios for the formation of TZOs have been proposed (see A. O’Grady et al. [40], for a review), one of which involves supernova (SN) kicks in close binaries ([32]). The asymmetries in SN ejecta generate a momentum kick on the newborn NS, widening the orbit and pumping eccentricity if the binary is not disrupted. Meanwhile, the SN ejecta may heat up the companion star, causing envelope inflation ([33, 27, 11]). With these two effects, it is possible that the pericenter separation of the post-SN orbit is small enough that a direct collision between the NS and the companion occurs ([32]), creating the possibility for the NS to settle down to the center Email: [email protected] [email protected] [email protected] of the companion star. Previous works have suggested that such a collision can indeed lead to the formation of dynamically stable TZOs (R. Hirai & P. Podsiadlowski [26]; L. E. Williams et al. [59], see also Y. Shi et al. [52]).

However, rapid accretion of the NS within the stellar envelope may trigger NS collapse, preventing the formation of a long-lived TZO (P. C.-K. Cheong et al. [14]). The accretion rate within the envelope is usually expected to be Eddington-limited (N. Ivanova et al. [29]), but efficient neutrino cooling can enable hypercritical accretion (R. A. Chevalier [17]; R. Chevalier [16]; C. L. Fryer et al. [24]; R. A. Chevalier [18, 19]; P. Anninos et al. [1]). On the other hand, angular momentum can limit the accretion rate (M. MacLeod & E. Ramirez-Ruiz [34]), potentially allowing a long-lived TZO to form.

Independently, the basic structure of TZOs was established by assuming hydrostatic and thermal equilibrium (K. S. Thorne & A. N. Zytkow [54, 55]; G. T. Biehle [5]). We refer to this equilibrium configuration as the standard TZO structure. Later studies relied on these equilibrium configurations as starting models to explore the subsequent long-term evolution, nucleosynthesis signatures, and provide further predictions for identifying TZOs in observations (R. C. Cannon et al. [10]; R. C. Cannon [9]; P. Podsiadlowski et al. [46]; R. Farmer et al. [23]) and the final fate of TZOs (P. Podsiadlowski et al. [46]; A. O’Grady et al. [40]).

Nevertheless, studies that follow the evolution of TZOs from the post-SN orbit, hydrodynamical interactions, to the long-term evolution are still missing. As shown by L. E. Williams et al. [59], after the hydrodynamical interactions the stellar structure can deviate considerably from single stars. In particular, the core temperature is much higher than the typical temperature of the CNO cycle. So far, it is unclear whether and how this structure evolves to thermal equilibrium; this thermal relaxation process might affect the long-term evolution, and even the final fate of a TZO. We aim to bridge the dynamical formation of a TZO triggered by SN kick and NS-MS collision, and the long-term evolution to investigate how the stellar structure after the formation adjusts towards the standard TZO structure. In this work, we adopt an MS stellar model (unaffected by SN ejecta) as a reference framework to study the impact of the formation process on the long-term evolution, while a broader exploration of the parameter space and more realistic simulations are left for future studies. This paper is structured as follows. In § 2 we present our methodology of 3D hydrodynamic simulation, 1D stellar evolution and numerical setup of our

reference model. In § 3, we present our results. Finally, we discuss several issues related to our results in § 4 and conclude in § 5.

Method

We present a schematic to illustrate the scope of our work in Figure 1. We consider a TZO formed from the collision between an NS and an MS star. The companion is expected to be an MS star because the NS progenitor should be initially more massive and have a shorter lifetime than that of the secondary star. Consequently, when the primary star reaches core collapse, the secondary star is still on the MS. For simplicity, we use a zero-age main-sequence (ZAMS) model for the MS star, and our simulation starts from the post-SN eccentric orbits, assuming a suitable kick velocity and direction. We neglect the evolution of the MS star and the effects of the supernova ejecta on it, which we do not expect to significantly alter the main conclusions of this work (cf., § 4).

We choose our reference model to be a $1.4M_\odot$ newborn NS and a $10M_\odot$ ZAMS. The choice of $10M_\odot$ is based on past population synthesis predictions on Be X-ray binaries (e.g., Y. Shao & X.-D. Li [51]; X.-T. Xu et al. 2025; C. Schürmann et al. 2025). For the post-SN orbital configuration, we adopt a semi-major axis of $a = 9R_{\odot}$ and an eccentricity of $e = 0.8$ for the post-SN orbit3, corresponding to an initial orbital period $P_{\mathrm{orb,i}} \approx1$ d and a pericenter distance of $r_{\mathrm{p}} = 1.8R_{\odot}$, which is smaller than our initial companion star size $R_{\ast}^{\mathrm{ini}} \approx4R_{\odot}$, allowing for the NS–MS collision.

Following the method illustrated in Figure 1, we perform a 3D hydrodynamical simulation for the NS–MS collision by using ENZO ([8]), where we generate a non-rotating ZAMS model for the companion star by using the 1D stellar evolution code MESA r24.08.1 (Modules for Experiments in Stellar Astrophysics; [41, 42, 43, 44, 45, 30]). After the collision assembles a dynamically stable TZO, we map the configuration into the MESA code for long-term evolution. In Sect. 2.1, we outline our numerical setup for 3D hydrodynamical simulation. In Sect. 2.2, we describe the method on 3D-1D mappings. In Sect. 2.3.2, we give our setup for 1D stellar evolution.

NS–MS collision: 3D Hydrodynamical Simulation

The ENZO code solves the Eulerian hydrodynamic equations, and the self-gravity of gas is computed by the Poisson solver in ENZO. We adopt the ideal gas equation of state (e.o.s) with an adiabatic index of $\gamma= 5/3$. We neglect nuclear energy release during the hydrodynamical simulation, assuming that its effect on the dynamics is small. We also neglect chemical abundance evolution and assume a fixed composition (for each parcel of gas) during this phase. We evolve the NS and the MS star (i.e., the self gravitating gas) simultaneously within a grid described in Sect. 2.1.1. The NS is modeled as a Newtonian point source with a softening radius (i.e. the effective size of the NS $R_{\mathrm{core}}$) set to $0.05R_{\odot}$; it produces energy at the Eddington luminosity, which heats up the gas that the NS passes through (Sect. 2.1.2). We assume the accretion of the NS is limited by the Eddington rate (see below). Since the timescale of the NS-MS collision is too short for the NS to accrete a considerable amount of material, we can neglect the mass change of the NS.

Simulation grid

The computational domain is a Cartesian box centered on the binary center of mass (c.o.m.) with a half-size of $51.2R_{\odot}$. The base grid has a resolution of $0.4R_{\odot}$, with a finest resolution of $0.0125R_{\odot}$. To resolve the main interaction region during pericenter passage, we employ three nested fixed refinement boxes centered on the domain origin, with half-sizes of $9$, $9.5$, and $10R_{\odot}$. Together, they ensure that the global merger structure is resolved with a resolution of $0.05R_{\odot}$.

Additionally, since the NS is initially outside the fixed refinement region and the final merger product can acquire a recoil velocity due to asymmetric mass ejection during the NS–MS collision, we adopt adaptive mesh refinement following the NS, ensuring the majority of the mass of the resulting TZO stays at the center of the refined computational domain. To resolve both the global gas structure of the TZO and the local gravitational field of the NS, we introduce five nested refinement regions centered on the instantaneous NS position, with radii of $4.0R_{\odot}$, $3.5R_{\odot}$, $3.0R_{\odot}$, $0.5R_{\odot}$ and $0.3R_{\odot}$. The three outer levels cover the region of the TZO and maintain a resolution comparable to that of the central fixed-refinement grid. The two innermost levels resolve the effective numerical radius of the NS, $R_{\mathrm{core}}$, with four grid cells, so that the softened-gravity region is sampled by approximately $8^3 = 512$ cells.

Neutron star in hydrodynamical simulation

The gravitational potential of the NS is softened within an effective radius of NS as $R_{\mathrm{core}}$, using the cubic spline kernel introduced in D. J. Price & J. J. Monaghan [47]. We note that the adopted $R_{\mathrm{core}}$ is much larger than the radius of a realistic NS due to the resolution limit of our 3D simulation. We assume the NS isotropically radiates at the Eddington luminosity,

where $G$ is the gravitational constant, $c$ is the speed of light, and $\kappa= 0.2\left(1 + X\right)\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$ is the Thomson opacity for fully ionized hydrogen and helium with $X$ the hydrogen mass fraction. We take $\kappa= 0.34\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$ (i.e., $X = 0.7$), corresponding to $L_{\mathrm{Edd}} \approx5 \times10^{4}L_{\odot}$. We further assume that a fraction, $f_{\mathrm{th}}$, of the radiation energy is converted into local thermal energy at a rate $L_{\mathrm{th}}$. Such heating should be more efficient in optically thicker regions (i.e., higher density), and we therefore introduce the following prescription,

where $\rho_{0}$ is a characteristic density scale separating optically thin and thick regimes, and $\rho_{\mathrm{loc}} \equiv\rho(r_{\mathrm{NS}})$ is the local density at the NS position. This function ensures a smooth transition from inefficient coupling ($f_{\mathrm{th}} \to0$ for $\rho_{\mathrm{loc}} \ll\rho_{0}$) to full thermalization ($f_{\mathrm{th}} \to1$ for $\rho_{\mathrm{loc}} \gg\rho_{0}$). In the following, we describe the modified energy equation with the effects of heating. The thermal energy converted from radiation is then deposited into the surrounding gas within a finite spherical region of a volume $V_{\mathrm{th}} = 4/3\pi r_{\mathrm{th}}^3$, which is introduced to avoid numerical singularities associated with point-like energy injection. The deposited energy is spatially distributed using a normalized smoothing kernel $S(r,r_{\mathrm{th}})$ satisfying $\int_{|\boldsymbol{r}-\boldsymbol{r}_{\mathrm{NS}}|<r_{\mathrm{th}}} S(r,r_{\mathrm{th}})\,\mathrm{d}^3r = V_{\mathrm{th}}$. With this prescription, the gas energy equation is modified by an additional source term,

where $\rho$, $\boldsymbol{v}$, and $P$ denote the gas density, velocity, and pressure, respectively, and $\varepsilon\equiv u+|\boldsymbol{v}|^2/2$ is the specific total energy with $u$ being specific internal energy. The gravitational potential $\Phi$ includes contributions from both the self-gravity of the gas and the softened NS potential.

The injection scale is chosen as $r_{\mathrm{th}}=0.1R_\odot$, which is larger than the gravitational softening radius, ensuring that the energy deposition is numerically resolved and spatially smooth. The typical density scale is chosen as $\rho_0=10^{-6}\,\mathrm{g}/\mathrm{cm}^3$. We choose the typical density scale such that the heating happens within the stellar envelope. The density scale should then satisfy the condition that it is much larger than the minimum density limit of the stellar profile $10^{-9}\,\mathrm{g}/\mathrm{cm}^3$ and much lower than the mean density $1\,\mathrm{g}/\mathrm{cm}^3$. However, we test different $r_{\mathrm{th}}$ and $\rho_0$ and we find different values do not significantly change the profile of the final object.

Initial model for the main-sequence companion

We use the MESA code to generate a $10M_\odot$ non-rotating ZAMS star with a metallicity of $Z=0.02$ for the companion star. The Eulerian hydrodynamical variables, initial density, velocity, and internal energy fields $\rho_{\mathrm{ini}}(r)$, $\boldsymbol{v}_{\mathrm{ini}}(r)$, $u_{\mathrm{ini}}(r)$ are assigned in a three dimensional grid based on a one dimensional solved stellar model. As our simulation is conducted in the c.o.m. frame, the MS initially is located at the apocenter of the orbit, and the initial velocity $\boldsymbol{v}_{\mathrm{ini}}(r)$ is set to be the Keplerian orbital velocity.

Additionally, since the Eulerian hydrodynamics solver in ENZO cannot handle vacuum, we add a low-density, high-temperature non-moving ambient medium outside the star. The medium is dilute enough not to affect our results.

1D stellar profile inferred from 3D structure

We evolve the system using 3D hydrodynamics until the NS is well settled at the center of the star and total bound gas mass converges to a nearly constant value, which we adopt as the criterion that the remnant has reached a quasi-steady configuration (see Eq. (6) below). We then map the 3D simulation output to 1D radial profiles for the subsequent stellar evolution by using the approach described in J. Schwab et al. [50].

For thermodynamic quantities, we perform angle-averaging on spherical shells in a spherical coordinate system $(r,\theta,\varphi)$ centered at the NS² with the orbital plane as its equatorial plane. For a given radius $r$, we calculate the angle-averaged quantities as

where $\mathrm{d}\Omega\equiv\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\varphi$ is the solid angle, and $\sigma(r)$ denotes the spherical surface at radius $r$. The “$*$” subscript denotes angle-averaged quantities. For the angular momentum profile, we compute a mass-weighted average specific angular momentum over a thin spherical shell with width $\Delta r$:

The lower and upper bounds of the radius of these 1D profiles are $R_{\mathrm{core}}$ and stellar radius $R_*$ respectively, where $R_*$ is defined as the boundary between the stellar material and the unbound gas. We estimate the mass of the bound gas $M_*$ through the Bernoulli criterion, which is

We then define the corresponding stellar surface $R_*$ with

This Bernoulli criterion provides an approximate estimate of the bound material, although it may not exactly represent the final stellar mass because the Bernoulli constant is not conserved given the presence of shock. We vary the adopted value of $M_*$ by up to 5% with corresponding $R_*$ and find that the resulting one-dimensional MESA models are insensitive to this choice.

Long-term evolution: 1D Stellar evolution

Our MESA calculation is initialized by an NS core with an effective core radius. For the stellar material above the NS core, we use the profiles mapped from the 3D hydrodynamical simulation. We describe our approach of generating the initial model for the MESA calculations in Sect. 2.3.1 and the input physics for MESA in Sect. 2.3.2.

Initial model of MESA

Consistent with the 3D simulation, the mass of the NS core $M_{\mathrm{core}}$ is initially taken to be $1.4\,M_{\odot}$, the effective core radius is equal to $R_{\mathrm{core}}$ (i.e., $0.05R_{\odot}$), and the accretion luminosity of the NS core is taken to be the Eddington luminosity (Eq. (1)), while the opacity uses the default opacity table in MESA. The mass growth of the core is governed by the Eddington accretion rate,

where $\eta$ is the accretion efficiency. For the envelope above the NS core, we use the mapped profiles to describe its structure (cf., Sect. 2.2). Since we do not expect significant nuclear reaction during the hydrodynamical phase, the metallicity $Z$ is taken to be the same as the initial stellar model used in the hydrodynamical simulation (i.e., $Z=0.02$). The total mass of the system is given by $M_{\mathrm{tot}}=M_{\mathrm{core}}+M_{\ast}$, where $M_{\ast}$ is determined in the 3D simulation (cf., Sect. 2.2).

Numerically, we use the relax_core method provided in MESA to set up the NS core and the relax_initial_entropy method to ensure the structure of the initial model matches the mapped 1D profiles. In addition, the angular momentum profile is included through a similar method, relax_initial_angular_momentum.

Input physics of MESA

We adopt the same input physics as that in C. Wang et al. [58] and X.-T. Xu et al. (2025). In the following, we briefly describe the most relevant inputs and the differences from X.-T. Xu et al. (2025). The boundary of convection zones is determined by the Ledoux criterion, which takes into account the effect of chemical gradient. We model convection by using the standard mixing-length theory with a parameter of $\alpha_{\mathrm{MLT}}=1.5$ (E. Böhm-Vitense [6]), and we consider the convective overshoot with a parameter of $\alpha_{\mathrm{OV}}=0.335$, as determined by I. Brott et al. [7]. While other types of mixing processes are also considered in our calculations, they are negligible in our reference model (cf., Figure 6). For the reference model, we include differential rotation. Following R. Farmer et al. [23], we adopt the approx21 nuclear network, which is suitable for the purpose of this work. We do not take into account the interrupted rapid-proton-capture process, as it occurs within $R_{\mathrm{core}}$, which is not resolved in our calculation, and it does not considerably contribute to the energy generation. For the same reason, the effects of general relativity are ignored. For simplicity, we turn off stellar wind for the reference model.

In addition, in previous studies, models were terminated when the NS core becomes a black-hole (BH) core (e.g., [10, 9]), which occurs when the core mass reaches $2M_{\odot}$. However, our model evolves across this BH formation point, with the same treatment for the core radius ($R_{\mathrm{core}}$), accretion luminosity, and mass accretion rate, and we neglect the additional energy release at the BH formation point (but see T. J. Moriya [37]). Then, we run the model until the MESA calculation stops, which we interpret as the breakdown of hydrostatic equilibrium. We discuss its implications in Sect. 4.

Result

Dynamical process: 3D Hydrodynamical simulation

NS-MS collision

We present the 3D hydrodynamical simulation of the NS-MS collision in Figure 2. The NS plunges into the MS star at the pericenter ($t=0.5P_{\mathrm{orb,i}}$), causing a strong shock near the stellar surface, which efficiently heats up and ejects the outer layers of the star. Meanwhile, orbital energy and angular momentum are dissipated by the gas friction, resulting in a strong orbital shrinkage and the deceleration of the NS. After the first interaction, the maximum separation between the NS and peak density location is about $5R_{\odot}$ (the lower panel of Figure 2), which is already comparable to the radius of the companion star before interaction (about $4R_{\odot}$), meaning that the NS has already been embedded into the star. The strong interaction then disrupts the outer part of the star at the second pericenter passage ($t=0.75P_{\mathrm{orb,i}}$). At the same time, the gas is focusing around the NS, making the NS position have the peak gas density (i.e., $\lvert r_{\mathrm{NS}}-r_{\mathrm{peak}}\rvert\approx0$). Finally, the gas falls back due to gravity. From $t=1.87P_{\mathrm{orb,i}}$, $M_{\ast}$ reaches a constant value, and the stellar structure does not evolve considerably (i.e. not contracting or expanding). The entire structure relaxed to a spherical configuration. Using the criterion introduced in Sect. 2.2, we find $M_{\ast}/M_{\ast}^{\mathrm{ini}}\approx90\%$, meaning that about $9M_{\odot}$ of stellar material is bound, giving a stellar radius of about $15R_{\odot}$. Due to the asymmetric envelope ejection, the resulting TZO has a small spatial velocity, and the remaining envelope material eventually co-moves with the NS. The corresponding stellar profile is taken to be the final profile given by the hydrodynamical simulation, which we present in the next subsection.

Final profile of hydrodynamical simulation

Figure 3 presents the top view (panel a) and side view (panel b) of the system at the end of the hydrodynamical simulation, and the profiles of the angle-averaged radial velocity (panel c) and angular velocity (panel d). The inner region of the system becomes nearly spherical, whereas the outer envelope remains highly asymmetric (top view of Figure 3); in particular, a tail formed by the passage of the NS on the left side. In order to measure how close the system is to the hydrostatic equilibrium, we present the ratio of the angle-averaged radial velocity ($v_r$) over Keplerian velocity at the equator ($v_{\mathrm{Kep}}$) in the panel (c) of Figure 3, where $v_r$ is given by

The ratio of $v_r/v_{\mathrm{Kep}}$ remains below 3% of the local Keplerian velocity throughout the structure, suggesting no significant radial motions, even in the outer region. Meanwhile, the envelope retains a certain amount of angular momentum transferred from the orbital motion of

the NS during the hydrodynamical inspiral phase. In the panel (d) of Figure 3, we present the ratio of the angular velocity ($\Omega_\ast$) over the Keplerian angular velocity at the equator ($\Omega_{\mathrm{Kep}}$), where $\Omega_\ast$ is obtained from the profile of specific angular momentum (Eq. (5)), which is

where the factor $3/2$ comes from the moment of inertia of a rotating shell. The inner region of the system, which contains most of the mass ($95\%M_\ast$ within $6R_\odot$), is slowly rotating. While $\Omega_\ast/\Omega_{\mathrm{Kep}}$ reaches about 0.4 at $r = 6R_\odot$, it is insufficient to form any disk-like structures (panel b of Figure 3). Therefore, the structure of the system is very close to the hydrostatic equilibrium and not considerably affected by the effects of rotation.

The 1D profiles obtained from this structure are shown in Figure 4. Compared with the initial profile, the post-dynamics profile has a much denser and hotter core and an extended envelope. The difference is mainly due to the gravity of the NS, which compresses the surrounding gas, forming a density spike at the inner region. At the same time, the compression significantly heats up the material near the NS. The temperature near the compact core ($2\times10^8\,\mathrm{K}$; Figure 4) can be understood by comparing the internal energy to the gravitational energy (i.e., virialized gas), which gives

where $T_{\mathrm{core}}$ is the temperature at the compact core surface, $\xi$ is a coefficient related to the e.o.s of the envelope gas, $k_{\mathrm{B}}$ is the Boltzmann constant, and $m_{\mathrm{p}}$ is the proton mass. This estimate roughly matches the innermost temperature shown in Figure 4. However, we note that the core temperature depends on $R_{\mathrm{core}}$ and can be higher in the smaller region (cf., Sect. 4). The temperature is high

near the compact core and remains at $\sim6 \times10^{6}\,\mathrm{K}$ at the radius enclosing 95% of the stellar mass (Figure 4). Consequently, intense nuclear burning extends from the compact core surface through most of the envelope by mass, significantly affecting the subsequent evolution in MESA (see Sect. 3.2).

In addition, we also compare the angle-averaged 1D profiles at the end of the hydrodynamical simulation and the profiles mapped into the MESA code (i.e., the starting point of the long-term evolution) in Figure 4. While the density profiles well match each other, there are small differences in the temperature profiles. The reason is that in the hydrodynamical simulation, the temperature is calculated by the ideal gas e.o.s in post-processing, while the mapped one by the realistic e.o.s in the MESA code. We do not expect this to have considerable effects on our results.

Stellar evolution

Figure 5 presents the long-term evolution of our TZO on the Hertzsprung–Russell (HR) diagram; we also show the track of a non-rotating single star model without stellar wind for comparison. The corresponding stellar structure is presented in Figure 6, and the evolution of the radius and surface luminosity in Figure 7. The evolution exhibits three major stages having different evolutionary timescales:

• Thermal relaxation phase: after the dynamical assembly, the resulting TZO is far from the thermal equilibrium. In particular, the energy generation rate of nuclear reaction is significantly higher than the luminosity at the surface $L$. The thermal relaxation phase lasts until the TZO reaches the global thermal equilibrium, i.e. when $L = L_{\mathrm{acc}} + L_{\mathrm{nuc}}$.

• Standard TZO phase: Once the system reaches the thermal equilibrium, the system remains stationary for about $10^{7}$ yr, and the compact core experiences negligible growth. Since the structure resembles the static solution explored in previous works (e.g., K. S. Thorne & A. N. Zytkow [54, 55]; G. T. Biehle [5]), we refer to this stage as the “Standard TZO” phase.

• Black-hole growth phase: the formation of a BH core is expected when the core mass reaches $2 M_{\odot}$. After that, the compact core grows at an accelerated rate, which we call the “BH growth” phase.

In the following subsections, we describe the evolution during these stages in detail.

Thermal relaxation phase

The entire thermal relaxation takes about 1000 yr (A to B in Figure 5 and 0–1000 yr in Figure 6-Figure 7). In the early stage, the outer stellar envelope is radiative, and the luminosity and surface temperature increase as the star evolves (Figure 5). This is because the nuclear power is initially about 6 orders of magnitude higher than the surface luminosity (Figure 7), driving the expansion of the entire star. As seen in Figure 6, the nuclear burning efficiency then drops accordingly, making the burning region shrink toward the surface of the compact object core. Meanwhile, the cooling makes the envelope more opaque, leading to a shrinking radiative zone. Once the star cools down to a certain extent, the nuclear power can no longer support the entire structure, which leads to a slight contraction of the star. Consequently, the accretion power of the NS core dominates the evolution of stellar structure making the surface luminosity adjust to the accretion power ($t = 100$ yr in Figure 7). The star restarts the expansion at the same time, with the envelope becoming fully convective. Thermal equilibrium (i.e., $L_{\mathrm{acc}} \approx L$) is achieved at about 1000 yr. This is consistent with the thermal timescale ($\tau_{\mathrm{th}}$), given by

where $U$ is the total internal energy of the star estimated using the virial theorem.

As shown in Figure 6, intense nuclear burning initially occurs in a large fraction of the mass, and convection transports the CNO-cycle processed material to the surface, producing carbon depletion and nitrogen enhancement with little change in ${}^{16}\mathrm{O}$ (Figure 8). Since the temperature in these nuclear burning regions is significantly higher than that in a normal MS star (Figure 4), the corresponding abundance pattern is different from the canonical CNO cycle³. Once the envelope becomes fully convective, the processed material is well mixed with the unburnt envelope, resulting in the final surface abundance. These spatially extended, formation-induced nuclear burning regions differ from standard TŻO models, which may provide a complementary chemical signature for identifying TZOs through observations (Sect. 4).

In addition, we also present the evolution of the surface angular velocity (in units of the critical angular velocity) in Figure 8. The expansion during the thermal relaxation stage substantially reduces its rotational rate. By the time thermal equilibrium is reached, the rotation rate has fallen below 5% of the critical value and has little effect on the subsequent evolution (see App. § 7.§ 1). We therefore ignore the effect of rotation in the following discussion.

Standard TŻO phase

This phase corresponds to the track from B to C on the HR diagram shown in Figure 5. Since our reference model does not take into account stellar wind, the evolution after reaching the thermal equilibrium is driven by the growth of the compact core. The growth rate is determined by the Eddington accretion rate (Eq. (8)), and it takes about $10^{7}$ yr for the NS core to reach $2M_{\odot}$ and become a BH. Consistent with previous studies (e.g., R. C. Cannon et al. [10]; C. Cannon [9]), the envelope supported by accretion power is fully convective, with a stellar radius of about $10^{3}R_{\odot}$. The surface luminosity is $6.3 \times10^{4}L_{\odot}$, matching the scaling law given in Cannon [9].

Black-hole growth phase

After the NS becomes a BH, the mass of the BH grows on a timescale of about $10^{8}$ yr. The evolutionary track on the HR diagram remains near the Hayashi line due to its strongly convective envelope (C to D in Figure 5).

Different from a normal post-MS star, the structure of our model is supported by the accretion power, allowing it to be brighter than the maximal luminosity expected by the single star model of the same mass. With the growth of the BH core, the envelope mass is reduced, and the Eddington factor in the envelope increases. At the end of the simulation ($t \approx4 \times10^{8}\,\mathrm{yr}$), the radiation pressure dominates over the gas pressure with the Eddington factor largely exceeding 1, which the MESA code cannot handle. We discuss the implications of this configuration in Sect. 4.3.

Model uncertainties

In this subsection, we investigate how the uncertainties in the input physics affect our result. We vary the size of the compact core ($R_{\mathrm{core}}$), the mass accretion rate of the compact core (controlled by $\eta$), wind mass loss, and nuclear network. In the following computations, we turn off stellar rotation for simplicity, since the evolution without rotation on the HR diagram (fiducial model without rotation; see the upper left panel in Figure 9) does not show considerable difference with the one with rotation (fiducial model; see Figure 5). All the models are computed until the break of hydrostatic equilibrium.

Size of the compact core

Our calculation cannot resolve the structure around the NS (radius about 10 km $\sim10^{-5}R_{\odot}$). Here we address this uncertainty by varying the core radius $R_{\mathrm{core}}$ from $0.025\,R_{\odot}$ to $0.1\,R_{\odot}$. As shown in Figure 9, the choice of $R_{\mathrm{core}}$ only affects the evolution during thermal relaxation, while the evolution after reaching thermal equilibrium remains unchanged. The value of $R_{\mathrm{core}}$ determines the temperature around the compact object core (Eq. (11)), affecting the nuclear power of the star. In the $R_{\mathrm{core}} = 0.025R_{\odot}$ case, the nuclear power is strong enough to drive a continuous expansion, while the other cases show a brief contraction due to the rapid decay of the nuclear power. Stronger nuclear burning may therefore process more envelope material and produce greater surface enrichment. We emphasize that enrichment is a robust qualitative outcome of our models, without attempting to predict its magnitude quantitatively.

Accretion rate

In the reference model, we compute the mass accretion rate by taking $\eta= 1$ in Eq. (8). However, the accretion rate can be much higher if the accretion flow is cooled by neutrino emission (R. A. Chevalier [17, 19]). To address this uncertainty, we vary $\eta$ from 1 down to 0.1 and 0.01 (only the $\eta= 0.01$ case is presented in Figure 9). The evolutionary tracks appear to be identical to the reference case. The difference lies in the evolutionary timescale, as the compact core grows at different rates

with different $\eta$. For $\eta= 1$, the final age of the model is $\sim100$ Myr, but $\sim1$ Myr for $\eta= 0.01$.

Stellar wind

In the fiducial model, we do not consider the effects of stellar wind for simplicity, which may describe the long-term evolution of TZOs at low metallicity. At high metallicity, strong stellar winds can significantly alter the evolution of massive stars. We account for stellar winds by recomputing the reference model (no rotation) with the wind prescription determined by I. Brott et al. [7]. The thermal relaxation phase persists, since the timescale is too short for a significant wind mass loss. In the following long-term evolution, the stellar envelope is stripped by winds within about 3 Myr, leaving a naked neutron star, and the mass of the compact core remains roughly unchanged during the evolution. In Figure 10, we provide a zoom-in for the evolution after reaching the thermal equilibrium. The luminosity remains roughly constant, since the structure is supported by the accretion power of the compact core.

Nuclear Network

Our fiducial model adopts the simplified nuclear network approx21.net, which is designed to capture the main hydrogen-burning processes using approximate reaction pathways, while omitting many intermediate reactions and species. It does not provide a complete treatment of high-temperature burning processes, such as the Ne–Na and Mg–Al cycles. As shown in Eq. (11) and Figure 4, the inner region reaches temperatures well above those required for standard hydrogen burning, allowing these additional reactions to occur. Although this hot-burning region contains only a small fraction of the stellar mass and has limited impact on the overall stellar structure, it may still leave behind a chemical signature through the evolution process.

We recompute the reference (no rotation) model with a larger network, stern.net tested in H. Jin et al. [31], which contains 77 reactions, including not only the CNO cycle but also the Ne–Na and Mg–Al cycles, and the burnings of helium, carbon, neon, and oxygen. The initial metallicity remains unchanged but with the metal distribution by M. Asplund et al. [2] to initialize extra elements. We see that the evolutionary track computed using stern.net is similar to the fiducial model, but the resulting chemical abundances are slightly different.

Discussion

Caveats of our assumptions

We start our hydrodynamical simulation with a non-rotating $10M_{\odot}$ ZAMS star. However, after the SN explosion, the companion would first interact with the SN ejecta. The impact can strip a small amount of envelope material and deposit thermal energy into the star, causing the envelope to inflate ([33, 27, 11]), enhancing the probability of the NS-MS collision (e.g., A. Ercolino et al. [21]). However, this process is sensitive to many factors, like explosion energy and ejecta mass, which is beyond the scope of this work. Importantly, the impact of the SN ejecta mainly affects the outer envelope, leaving the inner stellar structure nearly unchanged (M. Ogata et al. [39]). Since the inflated envelope only contains a small fraction of the stellar mass, we do not expect this uncertainty to considerably affect our results.

In the classical TŻO models, energy transport in the innermost region around the NS core is dominated by radiative diffusion instead of convection. The transition region between the radiative and convective zones is known as the “knee”, located at a radius of $\sim100\,\mathrm{km}$ (K. S. Thorne & A. N. Zytkow [55]). This knee feature plays a crucial role in determining the structure and energy source of TŻOs, and in massive TŻOs, the interrupted rapid-proton reaction occurring around the knee produces unique chemical signatures ([5, 9]). While our MESA simulation cannot resolve this knee region due to the resolution limit of the hydrodynamical simulation, the properties of the standard TŻO stage are broadly consistent with previous studies (cf., Sect. 3.2), indicating that our treatment is adequate for investigating the structure of the envelope.

Comparison with classical TŻO Evolution

By construction, the energy generation of classical TŻO models is dominated by the innermost region around the NS core (K. S. Thorne & A. N. Zytkow [55]; G. T. Biehle [5]; R. C. Cannon [9]). However, we demonstrate that, for the TŻOs emerging from NS-MS collisions, a large fraction of the mass is initially hot enough to initiate hydrogen burning (Figure 6), and the corresponding nuclear power can be orders of magnitude higher than the Eddington luminosity of the NS core (Figure 7). This condition leads to a stage of thermal relaxation, which is not explored in previous studies, and the evolution during this stage is sensitive to the choice of $R_{\mathrm{core}}$ (cf., Sect. 3.3). In particular, the contribution of nuclear burning increases with decreasing $R_{\mathrm{core}}$ (Figure 12). While the structure of our TŻO model is supported by the accretion power after the thermal relaxation stage, it is unclear whether the structure becomes nuclear-power-supported if our calculation can resolve deeper regions.

In classical TŻO models, there is no stable solution for the TŻO mass in the range $8$–$14M_{\odot}$ (e.g., R. C. Cannon [9]). This is because, in this mass range, pair production around the knee region becomes important, altering the opacity and therefore reducing the Eddington luminosity (Eq. (1)), and the nuclear power near the NS core is not strong enough to compensate this decrease in the accretion luminosity. However, we show that the core region is not the only significant energy source and that nuclear burning in the overlying envelope also contributes to supporting the envelope. Consequently, the thermal relaxation process may adjust the nuclear power in the envelope to compensate the decrease in the accretion luminosity, filling the mass gap expected by the classical TZO models.

Final fate of TŻOs

As described in Sect. 2.3.2, our long-term evolution runs until the crash of the MESA calculation, indicating the breakdown of the hydrostatic equilibrium. At the end of the calculation, the Eddington ratio in the stellar envelope typically exceeds unity. For example, in the reference model, the Eddington factor is over 50 at a mass coordinate near $9\,M_{\odot}$ (Figure 11). Such a configuration can also occur in massive MS stars but with a lower maximal Eddington factor ([25, 48]). In our model, the radiation energy stored beneath the outer envelope can exceed the gravitational binding energy of the overlying material, providing a potential energy source for eruptive mass ejections (e.g., [13]). Radiative heat exchange can excite pulsations and convert part of this energy into the kinetic energy of the gas. Our stability analysis identifies a growing radial mode, suggesting that these pulsations may become strong enough to unbind and eject the outer layers (Appendix § 7.§ 2). Therefore, we expect that, before becoming an isolated NS or BH, a TZO may experience mass ejections similar to luminous blue variables near the end of its life.

Our calculation does not include any possible energy release associated with the NS-to-BH transition. This is reasonable when most material and energy are trapped by the newly formed event horizon (e.g., [3]). However, if the progenitor NS has a large amount of rotational angular momentum, a disk can form during the collapse of the NS, producing jets ([36]). It is unclear whether such jets can break out from the envelope and produce transient events.

In addition, if nuclear reaction is the dominant energy source, TZOs may end up with a violent explosion ([46, 37, 38]). When such TZOs exhaust their nuclear fuel, gravitational collapse is expected, as the accretion power is not strong enough to support their envelopes.

Observation signatures

TŻOs can be observationally identified by their unusual surface abundance. Previous searches focus on lithium enrichment and the products of the interrupted proton-capture process ([56, 4, 53]). However, lithium enrichment through the Cameron-Fowler process is not unique for TŻOs (e.g., [4, 53]). In addition to these signatures, our results imply that the elements related to the CNO cycle (even Ne-Na, and Mg-Al cycles) may also help with identifying TŻOs. During the NS-MS collision, stellar material is significantly heated by shocks. Thus, intense nuclear reaction occurs in a large fraction of the mass (cf., Figure 6), causing considerable nitrogen and $^{13}\mathrm{C}/^{12}\mathrm{C}$ enrichment at the surface.

While the interrupted proton-capture process is expected to occur in massive TŻOs, the enrichment in nitrogen and $^{13}\mathrm{C}/^{12}\mathrm{C}$ can occur in less massive TŻOs as long as the stellar material is significantly heated during the formation process. However, these signatures are not unique to TŻOs either. For example, rotational mixing can also transport CNO-cycle processed material to the surface ([7]), and the hot-bottom burning in super-asymptotic giant branch (sAGB) stars can produce carbon depletion and nitrogen enhancement at the surface ([20]). In addition, the chemical signatures emerging from our models may also occur in a higher-mass regime, potentially appearing as chemically peculiar red supergiants.

Conclusion

In this work, we have investigated the formation of a TŻO through NS-MS collision and its long-term evolution. We focus on a $10\,M_{\odot}$ MS star. By combining three-dimensional hydrodynamical simulations with one-dimensional stellar-evolution calculations, we connect the dynamical formation of a TŻO to its subsequent thermal and nuclear evolution. The whole evolution can be summarized as follows:

• The 3D hydrodynamical simulations show that an NS colliding with an MS star can lead to the formation of a bound remnant that rapidly settles into a quasi-spherical configuration. In our simulation, roughly 90% of the original stellar mass remains bound within $R_{*} \approx15R_{\odot}$, with the NS located at the density maximum. This dynamically formed remnant provides the initial structure for the subsequent stellar-evolution calculations.

• The newly formed TŻO is initially far from thermal equilibrium and undergoes an early thermal-

relaxation phase lasting approximately $10^{3}$ yr. Compression and shocks during the collision phase heat up a substantial fraction of the star, producing an initial nuclear luminosity much larger than the surface luminosity. The envelope consequently expands and cools, becomes fully convective, and eventually approaches thermal equilibrium with the surface luminosity $L$ equal to the accretion luminosity $L_{\mathrm{acc}}$. The products of the nuclear reaction are transported to the surface by convection, producing carbon depletion and nitrogen enhancement.

• After thermal relaxation, the system enters a long-lived “standard TZO” phase. With weak winds and accretion at the Eddington rate, the object has a cool, fully convective envelope with $R_{\ast} \approx10^{3} R_{\odot}$ and $L \simeq6.3 \times10^{4} L_{\odot}$, and remains near the Hayashi line for $\gtrsim10^{7}$ yr. The global structure during this stage broadly resembles the equilibrium TZO solutions explored in previous studies.

• At later times, continued growth of the compact core and the transition from an NS to a BH make the envelope increasingly radiation-pressure dominated. Our fiducial model terminates when part of the envelope becomes strongly super-Eddington, suggesting the onset of substantial or eruptive mass loss. The present calculations do not directly determine the subsequent envelope disruption. When stellar winds are included, the envelope is instead stripped within approximately 3 Myr, leaving a nearly naked NS before substantial core growth, and the envelope also becomes super-Eddington in the end.

• Our exploration with various model parameters shows that the early thermal-relaxation phase and its associated nucleosynthesis are sensitive to the unresolved compact-core boundary. By contrast, the global evolution after thermal equilibrium is comparatively insensitive to the inner-boundary radii and nuclear-reaction networks. The accretion efficiency primarily changes the evolutionary timescale, whereas stellar-wind mass loss can qualitatively alter the lifetime and final fate of the system.

The key distinction between our model and classical TZO models concerns the origin of the surface chemical anomalies. In standard models, these anomalies are produced by nuclear burning confined to the innermost region surrounding the NS core. In our dynamically formed TZO, the NS–MS collision heats up a much larger fraction of the star through shocks and compression, leading to widespread nuclear burning during the thermal-relaxation phase. Convection subsequently transports the processed material to the surface, producing carbon depletion and nitrogen enhancement.

After thermal relaxation, the stellar structure approaches a classical accretion-supported TZO configuration, while the surface composition can retain a chemical imprint of the earlier nuclear burning. The resulting carbon and nitrogen anomalies may help identify TZO candidates when combined with other observational properties, although similar abundance patterns are likely to occur in other stars. More detailed nuclear networks and better-resolved calculations of the compact-core region will be required to obtain quantitative abundance predictions.

Acknowledgements

We thank Philipp Podsiadlowski and Liang Dai for insightful discussions. We also thank Charlie Gibson and Qinghui Sun for helpful technical advice and assistance. The simulations presented in this work were performed on the TDLI-Astro cluster at the Tsung-Dao Lee Institute. X.T.X. was supported by the National Natural Science Foundation of China (grant No. 12603055).

Further details of the Reference model

Rotation

In this appendix, we examine the rotational evolution of the TZO and its impact on stellar evolution. Tidal interactions before the SN may reduce the MS star’s spin toward synchronization; for simplicity, we adopt an initially non-rotating MS star, so the newly formed TZO acquires its angular momentum primarily from the post-SN orbit. During thermal relaxation, expansion slows the rotation, interrupted by a brief contraction-induced spin-up, while internal transport redistributes angular momentum within the envelope. By the time thermal equilibrium is reached, the surface angular velocity falls below $0.05\Omega_{\mathrm{crit}}$ (Figure 8).

Following MESA, we define $\Omega_{\mathrm{crit}} = \Omega_{\mathrm{Kep}}(1-\Gamma_{\mathrm{Edd}})^{1/2}$, where $\Omega_{\mathrm{Kep}} = (GM_{\mathrm{tot}}/R_{\ast}^{3})^{1/2}$ and $M_{\mathrm{tot}} = M_{\mathrm{core}} + M_{\ast}$ (Sect. 2.3.1). Here, $\Gamma_{\mathrm{Edd}} = \kappa L/(4\pi cGM_{\mathrm{tot}})$ is the surface Eddington factor, evaluated using the surface-layer opacity and the total stellar mass, rather than the compact-core mass used in Eq. (1). We use this definition during the hydrostatic, sub-Eddington phase.

The rotation rate after thermal relaxation can be estimated from the initial orbital angular momentum, $J_{\mathrm{orb}} = \mu[G(M_{\ast}^{\mathrm{ini}} + M_{\mathrm{NS}})a(1-e^{2})]^{1/2}$, where $\mu= M_{\ast}^{\mathrm{ini}}M_{\mathrm{NS}}/(M_{\ast}^{\mathrm{ini}} + M_{\mathrm{NS}})$, and $a$ and $e$ are the post-SN orbital parameters.

Approximating the TZO as a rigidly rotating body gives $I\Omega= f_JJ_{\mathrm{orb}}$, where $f_J$ is the retained fraction of orbital angular momentum and $I = k_IM_{\mathrm{tot}}R_{\ast}^{2}$, with $k_I$ the dimensionless moment-of-inertia coefficient. Thus,

For the reference system, $M_{\ast}^{\mathrm{ini}} = 10\,M_{\odot}$ and $M_{\mathrm{NS}} = 1.4\,M_{\odot}$, while the bound envelope has $M_{\ast} \simeq9\,M_{\odot}$ after the collision, giving $M_{\mathrm{tot}} \simeq10.4\,M_{\odot}$. Equation (A1) becomes

Taking $f_J = 1$ and $k_I = 2/5$ (a uniform-density sphere), the reference values $a = 9R_{\odot}$, $e = 0.8$, and $R_{\ast} = 10^{3}R_{\odot}$ give $\Omega/\Omega_{\mathrm{Kep}} \simeq0.018$. The normalization used in Figure 8 is instead $\Omega/\Omega_{\mathrm{crit}} = (\Omega/\Omega_{\mathrm{Kep}})(1-\Gamma_{\mathrm{Edd}})^{-1/2}$. This estimate illustrates the reduction in rotation as the star expands, but neglects the actual density profile, angular-momentum loss, and differential rotation.

For the same illustrative parameters, the centrifugal acceleration relative to gravity at the stellar surface is

The small influence of rotation is also supported by the non-rotating calculation, whose HR-diagram track is nearly indistinguishable from that of the rotating reference model (Sect. 3.3; Figure 9). Rotation therefore has little effect on the subsequent evolution in our models or the principal conclusions of this work.

The final profile of the Eddington factor

Figure 11 shows the local Eddington-factor profile of the last converged reference model during the BH-growth phase. Here, $\Gamma_{\mathrm{Edd}} = L/L_{\mathrm{Edd}}$ uses the total local luminosity, with $L_{\mathrm{Edd}}$ evaluated from the local opacity and enclosed mass. The vertical dashed line marks the BH-core boundary at $M_{\mathrm{core}}$. Throughout most of the envelope, $\Gamma_{\mathrm{Edd}} > 1$, reaching a maximum of $\simeq55$ at an enclosed mass of $\simeq8.8\,M_{\odot}$, before decreasing sharply in the outermost layers.

To investigate why our MESA calculation terminates, we first compare the radiation energy beneath the outer envelope with the gravitational binding energy of the overlying material. We divide the final profile at $m_b \simeq8.80\,M_{\odot}$, the peak of the total-luminosity ratio $L/L_{\mathrm{Edd}}$, and define the radiation energy of the inner envelope as $E_{\mathrm{rad,in}} = \int_{M_{\mathrm{BH}}}^{m_b}(3P_{\mathrm{rad}}/\rho)\,dm$ and the gravitational binding energy of the outer envelope as $E_{\mathrm{grav,out}} = \int_{m_b}^{M_{\mathrm{tot}}}(Gm/r)\,dm$. The ratio $E_{\mathrm{rad,in}}/E_{\mathrm{grav,out}} \simeq1.34$ shows that the stored radiation provides an energy reservoir large enough in principle to unbind these outer layers, provided sufficient energy is transferred to them.

To further assess the stability of our final model, we perform a linear perturbation analysis. GYRE is a numerical tool for stellar pulsation and stability analysis, integrated with MESA ([57, 43]). We use version 8.1 to examine the stability of non-adiabatic radial perturbations throughout the envelope of our final model, assuming that perturbations to heat transport are dominated by radiation. GYRE retains radial acceleration in the momentum equation,

where $P$ includes gas and radiation pressure. We adopt a time dependence $e^{-i\omega t/t_{\mathrm{dyn}}}$, where $\omega$ is the dimensionless complex frequency and $t_{\mathrm{dyn}} = [R_{\ast}^{3}/(GM_{\mathrm{tot}})]^{1/2}$. GYRE solves the linearized stellar evolution equations for $\omega$ in our final MESA model.

We find a growing radial mode with $\omega\simeq0.741 + 0.115i$, indicating that the final profile is pulsationally unstable on a timescale of approximately 10 yr. Together with the energy budget, this supports our interpretation that the model is approaching mass ejection, which would prevent further hydrostatic evolution.

Models with different $R_{\mathrm{core}}$