1. Introduction
Massive stars (
$M \gtrsim 8$
M
$_{\odot}$
) are almost always in binary and multiple systems with
$\sim$
70 % of them involved in various forms of interaction, including tidal interactions and mass transfer, leading eventually to close binaries and mergers (Moe & Di Stefano Reference Moe and Di Stefano2017). These interacting binary systems can give rise to high-energy phenomena, such as cataclysm variables (e.g. Warner Reference Warner1995), type Ia supernovae (e.g. Iben & Tutukov Reference Iben1984; Chevalier Reference Chevalier2012), short and long gamma-ray bursts (e.g. Fryer & Woosley Reference Fryer and Woosley1998; Brown, Lee, & Moreno Méndez Reference Brown, Lee and Moreno Méndez2007; Ramirez-Ruiz & Lee Reference Ramirez-Ruiz and Lee2009), and gravitational wave emission (e.g. Abbott et al. Reference Abbott2016).
For a certain range of binary and stellar parameters, the massive binary becomes a high-mass X-ray binary, where a red giant or red supergiant (RGS) feeds mass through the inner Lagrangian point,
$L_1$
, to a neutron star (or sometimes a black hole) companion in a phase of wind accretion and, possibly, Roche lobe overflow (RLOF). These systems form an accretion disk around the compact companion, which is X-ray bright and may develop jets. If the mass ratio is large, as is the case if the compact object is a neutron star, mass transfer reduces the orbital separation. In addition, the RGSs expand upon loss of mass, further accelerating the mass transfer. For neutron star–RGS systems, it is likely that unstable mass transfer and a common envelope (CE) phase may result (Ivanova et al. Reference Ivanova2013; Tauris et al. Reference Tauris2017). The outcome of this phase depends on whether or not the CE can be ejected before the neutron star merges with the helium core of the RGS.
Most observed X-ray binaries are undergoing a long-lived phase of stable wind accretion with timescales that depend on the parameters of the system. Some systems may be undergoing RLOF, in which case the evolution timescales are likely much shorter with a possibility of unstable mass transfer and CE (e.g. high mass X-ray binaries may remain in the RLOF state for only 10 000 yr; Savonije Reference Savonije1977, see also Tauris et al. Reference Tauris2017). Although, presumably, X-ray binaries in the stable wind accretion phase are more frequently observed (e.g. Cygnus X-1), it is possible that some may be caught in the faster phase of unstable mass transfer. Dickson (Reference Dickson2024) presented a model of X-ray binary M 33 X-7, that they believe to be in an unstable mass transfer phase based on a measurement by Ramachandran et al. (Reference Ramachandran2022) of the donor substantially overfilling its Roche lobe.
One open question in the study of massive CE interactions is the impact of the energetic feedback due to accretion of envelope gas onto the compact companion, before and during the in-spiral in the CE. For neutron stars and black hole companions, this could be critically important for the outcome of the interaction. Shiber et al. (Reference Shiber, Iaconi, De Marco and Soker2019) simulated an ad hoc jet, emanating from a companion in a low mass CE interaction to conclude that the jet would aid in unbinding envelope mass with consequences for the binary separation of the post-CE binary. On the other hand López-Cámara et al. (Reference López-Cámara, De Colle, Moreno Méndez, Shiber and Iaconi2022), using a self-regulated jet, powered by a fraction of the mass accretion rate that reached the inner boundary, showed that the jet would likely be quenched, even if it existed before the companion entered the CE.
A disk that forms during wind accretion and RLOF could survive inside the CE, be destroyed, or be destroyed and reform, as a result of accretion of envelope material on the secondary star (MacLeod & Ramirez-Ruiz Reference MacLeod and Ramirez-Ruiz2015; MacLeod et al. Reference MacLeod, Antoni, Murguia-Berthier, Macias and Ramirez-Ruiz2017; Moreno Méndez, López-Cámara, & De Colle Reference Moreno Méndez, López-Cámara and De Colle2017; Chamandy et al. Reference Chamandy2018; López-Cámara et al. Reference López-Cámara, De Colle and Moreno Méndez2019; Shiber et al. Reference Shiber, Iaconi, De Marco and Soker2019; López-Cámara, Moreno Méndez, & De Colle Reference López-Cámara, Moreno Méndez and De Colle2020; Moreno Méndez Reference Moreno Méndez2022). So far previous studies have found that the formation of the accretion disk inside a CE depends on the thermal properties of the envelope (adiabatic index of
$\lesssim\!1.2$
). Murguia-Berthier et al. (Reference Murguia-Berthier, MacLeod, Ramirez-Ruiz, Antoni and Macias2017) pointed out that this phase is only a transitory phase, due to the lack of stellar regions (zones of partial ionisation where
$\unicode{x03B3}$
is small enough) where the envelope is compressible enough to form a disk.
In this work, we study the formation of the disk around a
$1.4\ M_{\odot}$
neutron star, caused by RLOF mass transfer from a
$7\ M_{\odot}$
RGS, undergoing unstable mass transfer. We attempt to gain a quantitative idea of the parameters of an accretion disk to, eventually, determine its fate inside the CE. This work is also intended to contribute to the literature by studying accretion disks in 3D hydrodynamics. The goal is to determine when and how accretion disks form in response to mass accretion through
$L_1$
and as a function of a number of physical and numerical parameters to set this study in the broader context of disk formation (e.g. Makita, Miyawaki, & Matsuda Reference Makita, Miyawaki and Matsuda2000).
This paper is structured as follows. In Section 2, we outline the overall methodology with Sections 2.1 and 2.2 presenting the governing equations and simulation parameters. In Section 3 we give details of the formation and evolution of the disk and, in Section 3.1 the disk parameters. In Section 4, we discuss the sensitivity of our results to some physical and numerical parameters, while in Section 5 we present our conclusions.
2. Methods
To study the mass transfer phase through the
$L_1$
point, we consider a binary system consisting of a
$7\ M_{\odot}$
red supergiant as the donor star and a
$1.4\ M_{\odot}$
neutron star as the companion star (see Fig. 1 for a cartoon of the setup). We model the evolution of this binary system, between the RLOF phase and the CE phase, using the 1D implicit code mesa (Modules for Experiments in Stellar Astrophysics; version r21.12.1; Paxton et al. Reference Paxton2011, Reference Paxton2013, Reference Paxton2015, Reference Paxton2018, Reference Paxton2019).
Setup cartoon. The donor star is a red super giant of 7 M
$_{\odot}$
with a radius of 139 R
$_{\odot}$
, orbiting a compact object of 1.41 M
$_{\odot}$
, with a separation of 270 R
$_{\odot}$
. The dotted line (centred on the compact object) represents the computational domain in our 3D simulations.
We then use use the 3D hydrodynamic numerical code Mezcal (De Colle et al. Reference De Colle, Granot, López-Cámara and Ramirez-Ruiz2012) to simulate the formation of the accretion disk around the companion star using the mass transfer rate given by the mesa model as boundary condition to inject material into the computational domain. The 3D computational domain is represented by a dotted line in Fig. 1.
To only simulate the region around the accreting star, the 3D simulation is performed in a co-rotating system of reference, centred on the companion which is represented by a point mass particle encircled by an inflow boundary of radius
$\textrm{R}_\mathrm{in}$
. Mass is injected through a nozzle that represents the
$L_1$
point, located at the centre of the left boundary face. Below we motivate the method and the setup and explain the specific assumptions.
2.1 The hydrodynamic code and its governing equations
To study the formation and stability of accretion disks during the RLOF phase, we run a series of 3D numerical simulations employing the adaptive mesh refinement code Mezcal (De Colle et al. Reference De Colle, Granot, López-Cámara and Ramirez-Ruiz2012). The code integrates the hydrodynamic equations in a rotating frame. Self-gravity is not included in these simulations.
We solve the three-dimensional Euler equations for an inviscid gas (see, e.g. Makita et al. Reference Makita, Miyawaki and Matsuda2000), that is,
\begin{equation} \frac{\partial \unicode{x03C1}}{\partial t} + \nabla \cdot \left( \unicode{x03C1} \textbf{v} \right) = 0\;, \end{equation}
\begin{equation} \frac{\partial \!\left( \unicode{x03C1} \textbf{v} \right)}{\partial t} + \nabla \cdot \left(\unicode{x03C1} \textbf{v}\textbf{v} + P I \right) = -\unicode{x03C1} \textbf{f}\;, \end{equation}
\begin{equation} \frac{\partial e}{\partial t} + \nabla \cdot \left[\left( e+p \right)\textbf{v} \right] = -\unicode{x03C1} \textbf{v}\cdot \textbf{f} \;, \end{equation}
representing the evolution of mass density,
$\unicode{x03C1}$
, gas momentum,
$\unicode{x03C1} \textbf{v}$
, and gas energy density, e (for additional details see Pringle & King Reference Pringle and King2014). The variable P represents the pressure, I is the identity tensor,
$\textbf{v}$
is the velocity vector, and
$\textbf{f}$
is the specific force vector. The specific force vector
$\textbf{f}=\left(f_{x},f_{y},f_{z}\right)$
represents the gravity forces and fictitious forces associated with the rotating frame.
The binary system consists of a donor star of mass
$M_\mathrm{1}$
, and an accreting star of mass
$M_\mathrm{2}$
, with a mass ratio defined as
$q=M_\mathrm{1}/M_\mathrm{2}$
during their mass transfer phase. Our simulations are conducted in Cartesian coordinates in three dimensions with different levels of resolution. The binary system has an orbital separation a, with an orbital frequency
$\Omega=(G(M_\mathrm{1}+M_\mathrm{2})/a^3)^{1/2}$
. Then, the orbital period is given by
$P=2\unicode{x03C0}/\Omega$
. The origin of the rotating frame is the accreting star; the donor star is positioned on the left of the accretor (donor position is (
$-a$
,0,0)). We perform the simulation in dimensionless units using a as the length scale and
$a\Omega/2\unicode{x03C0}$
as the velocity scale. The time unit is the initial period of the binary. This means that quantities can be rescaled based on these units to physical units. In Table 1, we list the scaling factors between code and physical units.
Description of units, rescaling factors to physical systems, in terms of the binary separation (a), orbital frequency (
$\Omega$
), and total mass (
$M_\mathrm{tot}$
=
$M_\mathrm{{don}}$
+
$M_\mathrm{{acc}}$
, where
$M_\mathrm{{don}}$
is the mass of the donor star and
$M_\mathrm{{acc}}$
is the mass of the companion star).
With these assumptions, we can now write out the force vector in dimensionless form at each point in the computational domain, remembering that there is no self-gravity, but that we are operating in the rotating frame:
\begin{equation} f_\mathrm{x} = - 2v_\mathrm{y}\!-\!\left(x + \frac{q}{1+q}\right)\!+\!\frac{1/(1+q)}{|\textbf{r}_1|^3}(x - a)\!+\!\frac{q/(1+q)}{|\textbf{r}_2|^3}x, \end{equation}
\begin{equation} f_\mathrm{y} = 2v_x - y + \frac{1/(1+q)}{|\textbf{r}_1|^3}y + \frac{q/(1+q)}{|\textbf{r}_2|^3}y, \end{equation}
\begin{equation} f_\mathrm{z} =\frac{1/(1+q)}{|\textbf{r}_1|^3}z + \frac{q/(1+q)}{|\textbf{r}_2|^3}z, \end{equation}
where
$\textbf{r}_1$
and
$\textbf{r}_2$
are the distances from the point considered
$\textbf{r}=(x,y,z)$
to the centre of each star. In Equations (4) and (5), on the right-hand side, the first term represents the Coriolis force and the second term represents the centrifugal force. The last two terms in Equations (4) and (5), as well as the terms in Equation (6), represent the gravity force of each star.
2.2 Initial conditions
2.2.1 Calculation of the mass transfer rate through
$L_{{1}}$
using a 1D implicit code
To determine the mass transfer rate in a binary system between a massive red supergiant donor star and a compact accretor, we use mesa to simulate the evolution of a binary system comprising a 7 M
$_{\odot}$
, solar metallicity main sequence star with a radius of 57 R
$_{\odot}$
and a 1.4 M
$_{\odot}$
point mass companion, initially located at an orbital separation of a = 270 R
$_{\odot}$
. The Roche lobe radius of the primary star (
$\textrm{R}_{\textrm{L}_{1}}$
) is calculated according to the prescription by Eggleton (Reference Eggleton1983):
\begin{equation} \frac{\textrm{R}_{\textrm{L}_{1}}}{a} = \frac{0.49 q^{2/3}}{0.6 q ^{2/3}+ln{\left(1+q^{1/3}\right)}},\end{equation}
where a is the orbital separation, and
$q = M_1/M_2$
is the mass ratio between the donor star (
$M_1$
) and the accreting star (
$M_2$
).
We let the mesa simulation run for
$\sim 4.0 \times 10^{7}$
yr until the primary expands to fill its Roche lobe and starts to transfer mass to the companion. At this time, called time zero in Fig. 2, the red supergiant has a helium-burning core surrounded by a shell that experiences hydrogen burning (primarily through the CNO cycle), along with a massive convective hydrogen envelope. The orbital separation is 270.4 R
$_{\odot}$
, the red supergiant has a mass of 6.98 M
$_{\odot}$
, a radius of 134 R
$_{\odot}$
, an effective temperature of
$3\,981$
K, and a luminosity of
$3\,645$
L
$_\odot$
. The companion star has a mass of 1.4 M
$_\odot$
, and the mass transfer rate at time zero is
$1.4\times 10^{-8}$
M
$_\odot$
yr
$^{-1}$
.
Initial (start of RLOF) and final (the time of CE) binary parameters for 30 000 yr of Roche lobe mass transfer, modelled with the mesa code.
Temporal evolution of the mass transfer rate from the primary to the secondary star in our long-term, mesa simulation of the binary system. The initial masses of the stars are 7 M
$_{\odot}$
and 1.4 M
$_{\odot}$
; they have an initial orbital separation of 270 R
$_{\odot}$
and an orbital period of 177 days. We indicate with the orange box the mass transfer evolution that we are simulating in 3D.
We then continue the mesa simulation for an additional
$\sim$
30 000 yr, during which the mass transfer rate increases to a maximum value of 9.7
$\times10^{-2}$
M
$_\odot$
yr
$^{-1}$
(as we can see in Fig. 2). After 30 000 yr the mass of the red supergiant is 6.88 M
$_{\odot}$
its radius is 143 R
$_{\odot}$
and its effective temperature is
$\sim$
5 012 K, while the companion star mass has increased to 1.49 M
$_{\odot}$
, having accreted some mass. The orbital separation is 243 R
$_{\odot}$
. After this point we assume that the mass transfer leads to a CE in a short timescale. See Table 2 for a summary of the initial and final parameters of the mesa simulation.
The 3D simulation described next (Section 2.2.2) instead spans only 21 yr. This period of time is taken between 29 500 yr and 29 521 yr of the 30 000 yr stretch of the 1D simulation. Hence, the mass transfer rate in the 3D simulation, is prescribed from the 1D simulation to be between 2.3
$\times10^{-4}$
M
$_\odot$
yr
$^{-1}$
and 9.7
$\times10^{-2}$
M
$_\odot$
yr
$^{-1}$
. In Section 4.1, we will explore the impact that this mass-transfer rate choice has on the the accretion disk parameters.
Simulation summary: inputs are varied to understand the resilience of the accretion disk parameters to numerical and physical input parameter changes. We vary the initial mass transfer rate (sim-dot-
$\#$
), injection velocity through
$L_\mathrm{{1}}$
(sim-vel-
$\#$
), and background temperature and density (sim-bgT-
$\#$
). We compare each simulation to the reference simulation (sim-0). The majority of the simulations were carried out for 21 yr and start with a donor star of 7 M
$_{\odot}$
, a companion star of 1.4 M
$_{\odot}$
, an orbital separation of 266.34 R
$_{\odot}$
, and
$\Omega$
of 4.18
$\times10^{-7}$
s
$^{-1}$
.
2.2.2 Calculation of disk formation using a 3D explicit code
In the 3D simulations, a 1.4 M
$_{\odot}$
point mass particle represents the neutron star companion, at the centre of the domain. We assume that the donor star, with a mass of 7.0 M
$_{\odot}$
, is located outside the computational domain at a position
$({-}a,0,0)$
from the companion star, where
$a= 266$
R
$_{\odot}$
is the binary separation at the start of the hydrodynamic simulation. In this way, the position of
$L_\mathrm{1}$
is at 83.1 R
$_{\odot}$
(
$5.79\times 10^{12}$
cm) from the companion star. The gas is injected through the
$L_\mathrm{1}$
point, represented in the simulations by a small rectangular boundary with a variable size given by the cross section of the mass transfer stream, as explained in the following. The setup scheme is shown in Fig. 1.
The point-mass companion sits inside a spherical inflow boundary of radius
$\textrm{R}_\mathrm{in}$
=1.3 R
$_{\odot}$
(
$9.3\times10^{10}$
cm). Since cells cannot be empty, we set the density of gas within this inner boundary to
$\unicode{x03C1}_\mathrm{in}$
= 1.0
$\times10^{-20}$
g cm
$^{-3}$
and the pressure to 1.0
$\times10^{-8}$
dyne cm
$^{-2}$
. The cells inside the inner boundary are rewritten at every time step to have these values. Outside of the inner boundary we set a low background density
$\unicode{x03C1}_\mathrm{bg}$
= 2.6
$\times 10^{-22}$
g cm
$^{-3}$
, a temperature
$\textrm{T}_\mathrm{bg}$
= 1.0
$\times 10^{5}$
K, and a pressure
$P_\mathrm{bg}=8.9\times 10^{-11}$
dyne cm
$^{-2}$
. These initial conditions were chosen to have the lowest density possible before the simulation stops being able to calculate gas advection. The sensitivity of our simulations to different values of the background conditions is tested in Section 4.3. The gas adiabatic index is fixed to
$\unicode{x03B3}=1.1$
during the whole simulation, based on Makita et al. (Reference Makita, Miyawaki and Matsuda2000), MacLeod & Ramirez-Ruiz (Reference MacLeod and Ramirez-Ruiz2015) and Murguia-Berthier et al. (Reference Murguia-Berthier, MacLeod, Ramirez-Ruiz, Antoni and Macias2017). See Section 4.5 for further discussion about the adiabatic index in simulations of accretion disks.
Following the Jackson et al. (Reference Jackson, Arras, Penev, Peacock and Marchant2017) prescription for optically thin mass transfer, we define an elliptical area centred on
$L_\mathrm{{1}}$
, where most of the material escapes
$L_\mathrm{{1}}$
, referred from now on as the ‘nozzle’. The nozzle area is defined by
$S_{\textrm{L}_{1}}= \unicode{x03C0} \Delta z \Delta y$
, where the two dimensions are related to the pressure scale height in the y- and z-directions, and vary as
\begin{equation} \Delta y = \frac{\sqrt{2}c_T}{\Omega\sqrt{A-1}}, \;\;\; \Delta z = \frac{\sqrt{2}c_T}{\Omega\sqrt{A}}, \end{equation}
where
$c_\textrm{T} = \sqrt{\unicode{x03B3}k_\textrm{B} T / \unicode{x03BC}}$
is the isothermal sound speed,
$k_\textrm{B}$
is the Boltzmann constant,
$\unicode{x03BC}$
is the mean mass of a gas particle,
$\Omega$
is the orbital frequency and A is a dimensionless coefficient that depends on the mass ratio, q, (equivalently
$M_2/M_1$
or
$M_1/M_2$
in this equation), defined as
\begin{equation} A(q) = 4 + \frac{4.16}{-0.96 + q^{1/3} + q^{-1/3}}.\end{equation}
We set the length of the nozzle on the x-axis to be two cells thick at the coarsest level of refinement, or
$\Delta x$
= 1.30 R
$_{\odot}$
(
$9.05\times 10^{10}$
cm). The initial mass injection rate is
$\dot{M}_\mathrm{{ini}}=2.31\times10^{-4}$
M
$_\odot$
yr
$^{-1}$
, and it has a subsonic velocity in the x direction with value
$v_{\textrm{L}_{1}}=77.7\times 10 ^{3}$
cm s
$^{-1}$
(Lubow & Shu Reference Lubow and Shu1975; Jackson et al. Reference Jackson, Arras, Penev, Peacock and Marchant2017; Cehula & Pejcha Reference Cehula and Pejcha2023). The mass injection rate (
$\dot{M}$
) is interpolated at each time step (
$\Delta t$
), using the values given by the mesa simulation (see Section 2.2.1). Using the interpolated values, we calculate the nozzle volume (V
$_{{L_1}}$
), nozzle density as
$\unicode{x03C1}_{\textrm{L}_{1}}=\dot{M}{\Delta t}/V_{\textrm{L}_{1}}$
, and the pressure of the nozzle as
$P_{ L_1}=\unicode{x03C1}_{\textrm{L}_{1}}c_\textrm{T}^2$
, using the effective temperature of the donor star to calculate
$c_\textrm{T}$
.
We start the hydrodynamic simulation at
$t_\mathrm{{ini}}$
= 29 500 yr after the start of the mass transfer in the mesa simulation and let the hydrodynamic simulation run for 21 yr, with a final mass transfer rate of
$\dot{M}_\mathrm{ fin}=9.7\times10^{-2}$
M
$_\odot$
yr
$^{-1}$
. We also performed eight different simulations, altering different physical and numerical parameters (see Table 3). We will justify the need for these additional simulations in Section 4.
The size of the computational domain depends on the distance between the first Lagrange point and the accreting star, initially at
$-d_{\textrm{L}_{1}}$
= 83.1 R
$_{\odot}$
(5.79
$\times 10^{12}$
cm). Hence, the computational box has dimensions
$d_{\textrm{L}_{1}}\leq x\leq1.25 d_{\textrm{L}_{1}}$
, −1.5
$d_{\textrm{L}_{1}}\leq y\leq1.5 d_{\textrm{L}_{1}}$
and −0.5
$d_{\textrm{L}_{1}}\leq z\leq 0.5 d_{\textrm{L}_{1}}$
with outflow boundary conditions at each of the six faces (except for the location of the nozzle that is technically an inflow boundary). The grid was modified to have better resolution surrounding the companion star and around the nozzle, avoiding numerical problems. We employ (72, 96, 32) cells at the coarsest refinement level, with three levels of refinement (see Section 4.4 for convergence tests), corresponding to a maximum resolution of
$\Delta x = \Delta y = \Delta z=$
0.65 R
$_{\odot}$
(
$4.53\times 10^{10}$
cm).
3. Results: Formation and evolution of the disk
In what follows, we give details of the simulation whose input parameters are explained in Section 2.2. We refer to this simulation as the ‘Reference simulation’ (sim-0 in Table 3). Later, we test the results against changes in the input parameters.
In Fig. 3 we show a volumetric density rendering that shows the accretion disk forming over 21 yr. At the start of the simulation, the injected material with velocity in the x-direction (into the computational domain) effectively free falls towards the companion star under its own gravity, but due to its angular momentum, it sweeps around forming a disk. The injected material moves towards the companion, taking
$\approx$
0.12 yr to reach the internal boundary, and is deflected around the central potential toward the lower density medium. The interaction between the deflected material and the injection stream creates the accretion disk’s base structure, seen after 0.3 yr (upper left panel of Fig. 3). One year later, the accumulated deflected material forms two high density spiral arms on the orbital plane (see upper right panel).
Volumetric density rendering showing the accretion disk at four times (t =
$0.3$
yr,
$1.5$
yr,
$10.5$
yr, and 21 yrs). We can appreciate the 3D structure of the accretion disk, surrounding the companion star (represented by the white area in the middle of the accretion disk).
After 10.5 yr the shock between the spiral arms has created a high density structure with a disk-like shape, surrounding the companion star (see bottom left panel). The material in the disk is orbiting around the companion star, forming bow shocks with the injected material, showing a higher density on the left side of the computational domain than on the right side, consistent with the results of Makita et al. (Reference Makita, Miyawaki and Matsuda2000).Footnote
a
By the end of the simulation at 21 yr (see bottom right panel of Fig. 3) the accretion disk has maintained its approximate structure for
$10.5$
yr, albeit while growing somewhat in radius and scale height, due to the accumulation of mass.
In the next section, we quantify what we have just described qualitatively.
3.1 Accretion disk parameters
In Fig. 4, we present orbital and perpendicular density slices of the accretion disk at the same four times shown in Fig. 3. Along the orbital plane we see the formation of spiral arms, the injected mass leaves the nozzle and is deflected towards the companion star, due to its gravitational pull, Coriolis force, and centrifugal force. With an increase in the rate of mass transfer, a high density structure can be seen in the last panel of the top row of Fig. 4.
Density slices of the accretion disk on the orbital plane (top row) and perpendicular plane (bottom row) of sim-0 at times t = 0.3 yr, 1.5 yr, 10.5 yr, and 20.9 yr (from left to right). The letters in the third and fourth upper panels represent the positions used to calculate the density profiles (see text).
The edge-on panels in Fig. 4, bottom row, show the formation of the aforementioned bow shocks, due to the interaction between the injection stream and the accretion disk. Also, the edge-on view shows the presence of small, low density outflows along the polar axes, and may suggest the formation of hydrodynamically collimated jets in the future; however, our simulation ends too early to follow their development. These outflows are not visible in the rendered plots as a result of the choice of the colour bar limits.
In Fig. 5 we plot density profiles along the z-axis of the accretion disk, measured at 26.5 R
$_{\odot}$
(1.85
$\times10^{12}$
cm) from the companion star along the positive and negative side of the x- and y-axes, at
$t=10.5$
yr (blue line), and
$t=21$
yr (red line). The upper panels of Fig. 5 show the density profiles measured on the positive y-axis and x-axis, panels (a) and (b), respectively. Panels (c) and (d) show the density profiles measured on the negative sides of the y-axis and x-axis, respectively (these labels are also marked in Fig. 4).
Panel (d) shows the material to be more extended in the z-direction because the incoming material from the nozzle is constantly interacting with the disk material. Meanwhile, panel (b) presents a clear description of the disk thickness, with a sharply decreasing density above and below the midplane, just as we can see in the density maps (Fig. 4). At 10.5 yr the disk has an average thickness of 15.6 R
$_{\odot}$
(
$1.09\times10^{12}$
cm). By the end of the simulation (21 yr) the disk has reached an average thickness of 20.4 R
$_{\odot}$
(
$1.42\times10^{12}$
cm). We also measure the scale height of the disk at 10.5 yr to be H = 6.7 R
$_{\odot}$
(0.47
$\times10^{12}$
cm), while at the end of the simulation (21 yr) the disk scale height has slightly decreased to 4.9 R
$_{\odot}$
(
$0.34\times10^{12}$
cm).
Panels (b) and (d) in Fig. 5 show that the thickness of the disk perpendicular to the line that joins the two stars is thicker and less defined. Panel (d), includes the interaction of the deflected material that has left the nozzle and encounters the accretion disk, just as we can appreciate in the density maps of the orbital plane. The shape of the vertical density profiles at each of the four locations does not change significantly over 10 yr; only the density increases over time due to the constant injection of material.
We measured the mass inside one hundred radial points (Fig. 6), measured from the inner boundary (encircling the companion star) to the edge of the computational domain. We repeat this calculation for the same times as in Fig. 4. At
$t=0.3$
yr (orange line panel) the injected material has just reached the internal boundary, but the disk has not formed yet. Once the disk has formed, at
$1.5$
yr (maroon line panel in Fig. 6), the enclosed mass increases with radius until it reaches the edge of the disk. At this point, we see a flattening of the slope of the enclosed mass. The nozzle is located at a radius of 74.4 R
$_{\odot}$
(
$5.18\times10^{12}$
cm) from the inner boundary, and it manifests itself in the steepening of the gradient at the largest radius, particularly evident at the last time (bottom right panel in Fig. 6). We determine the radius of the disk by locating the first inflection point (black filled circle in all panels). The cumulative mass profile and hence the inferred mass and radius of the disk, depend somewhat on assumed parameters such as injection velocity. We will explore these dependencies in Sections 4.1 and 4.2. In Table 4, we summarise the disk characteristics discussed so far.
Density profile of the accretion disk versus height at a radius of 26.7 R
$_{\odot}$
(1.86
$\times10^{12}$
cm). Two times are shown:
$t=10.5$
yr (blue line) and 21 yr (red line). The readings were taken at four symmetric points around the companion indicated in Fig. 4 star along the x- and y-axis (panel (a):
$+y$
-axis, (b):
$+x$
-axis, (c):
$-y$
-axis, (d):
$-x$
-axis).
Disk properties at different moments in time in sim-0.
Cumulative mass as a function calculated from the companion star, at different times:
$t=0.3$
yr (top left panel, orange line),
$t=1.5$
yr (top right panel, maroon line),
$t=10.5$
yr (bottom left panel, purple line),
$t=21$
yr (bottom right panel, blue line). The black symbol in each panel indicates the adopted radius of the disk.
Using the inner inflow boundary, we measure the accretion rate onto the companion as a function of time (plotted in Fig. 7). We identify a cubic volume, 2
$\textrm{R}_\mathrm{in}$
on a side, which contains the inflow boundary sphere. The total accretion rate is measured by calculating the mass flux through the cubic boundary, by taking the projection of the velocity of each cell surrounding the boundary onto the normal direction to the boundary (solid line in Fig. 7). If the projected velocity of a cell is pointing inside the fixed box, we assume that the material in that cell is accreted and deleted from the simulation in the next timestep, or
\begin{equation} \dot{M} = \sum_{{i=1}}^{{N}} \unicode{x03C1}_{{i}} A_{{i}}\vec{v}_{{i}} \cdot \hat{n}, \end{equation}
where
$\unicode{x03C1}_{{i}}$
is the density of the cell,
$\vec{v}_{{i}}$
is the velocity vector of the cell,
$A_{{i}}$
represents the area of the cell face through which the material will cross in the next time step, and
$\hat{n}$
is the vector perpendicular to that face directed into the cubic boundary. We also carry out the same measurement with a slightly different approach: we project the velocity vector of each cell as well as the area of the face that the gas crosses, along the radial direction to the companion – this method results in somewhat smaller values for the accretion rate by a factor of 0.6 (dashed line in Fig. 7).
Mass accretion rate onto the companion (traversing the inner, inflow boundary) as a function of time. The solid line represents the mass accretion rate assuming the inner boundary is a cube with edge length 2R
$_\mathrm{{in}}$
and the dashed line shows the mass accretion rate calculated projecting the area and velocity of the cell over the spherical inner boundary.
The values of the accretion rate determined with the first method above are listed in Table 4. The accreted mass does not contribute significantly to the mass of the companion, which increases only by
$8.4\times10^{-3}$
M
$_{\odot}$
during 21 yr.
The accretion rate into the inner boundary in our simulations may be affected by numerical viscosity. However, the fact that the accretion rate values are approximately 100 times smaller than the injection rate at
$L_1$
, argues that numerical viscosity does not play a substantial role. We later also show that the accretion rate is well converged (Section 4.4), in agreement with this conclusion. It would therefore be reasonable to state that viscosity in our simulations is provided by disk turbulence and possibly the effect of spiral shocks. We can therefore compare the numerically-derived accretion rates to those theoretically predicted by Shakura & Sunyaev (Reference Shakura and Sunyaev1973):
\begin{equation} \dot{M}= 3 \unicode{x03C0} \Sigma \alpha c_{s} H_\mathrm{{disk}},\end{equation}
where
$\Sigma$
is the surface density of the disk,
$\alpha$
is an efficiency factor and where we assume
$\Sigma=M_\mathrm{{disk}}/\textrm{R}_\mathrm{{disk}}^2$
. For thin disks, we arbitrarily assume
$\alpha=0.1$
. The variable
$c_s$
is the speed of sound at
$\textrm{R}_\mathrm{{disk}}$
and
$H_\mathrm{{disk}}$
is the scale height of the disk. Using the aforementioned results at the 10.5 and 21 yr, we predict accretion rates of
$5.0\times10^{-4}$
M
$_\odot$
yr
$^{-1}$
and
$3.7\times10^{-3}$
M
$_\odot$
yr
$^{-1}$
, respectively, consistent with the accretion rate measured at the same two points in the simulation (Table 4).
The material in the nozzle has density (
$\unicode{x03C1}_{{L_1}}$
) and pressure (
$P_{{L_1}}$
) defined in terms of the mass injection rate, the volume of the nozzle, and the velocity, as explained in Section 2.2. Assuming an ideal gas, we show, in Fig. 8, temperature slices in the orbital and perpendicular planes at the end of the simulation (21 yr). The material inside the nozzle has an average temperature of
$\sim 5\,000$
K, once the material leaves the nozzle its temperature drops to
$\sim 3\,000$
K, due to the pressure difference between the nozzle and the material just outside the nozzle. Inside the inflow boundary around the accretor, the high temperature (
$\sim 10^{7}$
K) is due to the imposed low density. We observe a gradient of temperature decreasing radially away from the inner boundary, we can also observe on the perpendicular plane (x-z plane) the interaction between the injected mass and the mass that circles around the companion.
Temperature slices for sim-0 in the orbital (top panel) and perpendicular (bottom panel) planes of the disk at
$t=21$
yr.
In Fig. 9, we plot a temperature profile along the positive x-axis at 21 yr, with the internal boundary shown as a black vertical line and the radius of the disk as a dashed grey line. In the mid-plane along the positive x-direction, the cells closer to the inner boundary have a temperature of
$1.1\times10^{6}$
K, their location corresponds with the region where the pressure gradients is the highest; the cells located at the edge of the disk have lower temperatures of
$\sim$
60 000 K.
Temperature profile for sim-0 in the mid plane along the positive x-axis at
$t=21$
yr. The disk’s radius is indicated with a vertical grey dashed line, while the inner boundary’s radius is marked with a vertical black line. The solution for steady disks is indicated by the green curve (
$T\propto r^{-3/4}$
) and the best fit (
$\alpha = 1.12\pm0.02$
) is indicated by the red curve and red shaded area. The average percentage error between the temperature profile and the fit line is 15%.
Following the prescription of Shakura & Sunyaev (Reference Shakura and Sunyaev1973) for steady disks, the effective temperature profile should follow
$T\propto r^{-0.75}$
(green line in Fig. 9). This function fits our data with an average percentage error of 38%. The best fit is obtained using a steeper
$T\propto r^{-1.12\pm 0.02}$
(red line and shaded area in Fig. 9) with an average percentage error of 15%. Mineshige et al. (Reference Mineshige, Duschl, Frank, Meyer, Meyer-Hofmeister and Tscharnuter1994) stated that a value of
$\alpha \lt -0.75$
characterises X-ray binary systems in the flaring branch on the X-ray hardness-intensity diagram.
Slices of Mach number in the orbital plane (top panels) and perpendicular plane (bottom panels) of the disk at
$t=10.5$
yr (left panels) and
$t=21$
yr (right panels) for sim-0.
However, our simulation does not have explicit cooling, so we may expect it to diverge from the prediction of Shakura & Sunyaev (Reference Shakura and Sunyaev1973). Instead, our adopted equation of state implies that
$T \propto (\unicode{x03B3}-1)e$
, and as the energy dissipated per unit mass in an accretion disk scales as
$1/r$
, it is natural for the temperature to scale similarly.
The material is injected through the nozzle with a subsonic velocity of
$7.7\times10^4$
cm s
$^{-1}$
in the x direction (note that the sound speed value near the nozzle ranges between 6 and 8 kms
$^{-1}$
). In the upper panels of Fig. 10, we show the Mach number at 10.5 yr and at the end of the simulation (21 yr), with slices in the orbital and perpendicular planes. The nozzle is on the left of the domain, seen as a short vertical bar with subsonic velocity. Just outside the nozzle, the velocity becomes highly supersonic (
$\mathcal{M}$
= 27 at 21 yr). When the injected material approaches the centre of the domain, it circles it and collides with the material that was already in orbit around the companion star, slowing down.
In Fig. 11, we plot the gas velocity profile in the mid-plane along the positive and negative x-axes at
$t=10.5$
yr, normalised to the Keplerian velocity (
$v_{k}$
=
$\sqrt{GM/r}$
). We note that normalising velocities to the Keplerian values is appropriate because gas velocity vectors are effectively entirely in the azimuthal direction, except for a few cells near the centre that have a 20% radial component. The inner boundary radius,
$\textrm{R}_\mathrm{{in}}$
, is marked with a black line, while a grey line indicates the radius of the disk,
$\textrm{R}_\mathrm{{disk}}$
. The cell velocities around the inner boundary are Keplerian on average, though there is quite a bit of scatter in individual cells due to the gas there not moving entirely azimutally. Cells at the edge of the disk move with velocities approximately 80%–90% of the Keplerian value. This is likely due to the pressure support of the gas in the disk.
Velocity profile in the mid-plane along the positive and negative axes at
$10.5$
yr. The velocity is normalised to the Keplerian velocity (
$v_{k}$
). The inner boundary is marked by the black vertical solid line and the radius of the disk by the vertical grey dashed line.
We further study the disk by examining the locus of the gas cells in the specific energy vs. specific angular momentum plane (Hayashi et al. Reference Hayashi, Kawaguchi, Kiuchi, Kyutoku and Shibata2021, Fig. 12.Footnote b ) We calculate the specific orbital energy as
\begin{equation} e_\mathrm{{orb}} = {\frac{1}{2}} v^2 - {\frac{GM}{r}},\end{equation}
where v is the magnitude of the velocity in the centre of the cell, r is the position of each cell with respect to the companion star, and M is the mass of the companion star. The magnitude of the specific angular momentum for each cell is:
\begin{equation} j_\mathrm{{orb}} = r \cdot v.\end{equation}
Each cell is colour-coded by density (with the same colour bar as in Fig. 4), where the cells that compose the high density accretion disk are marked in red, and dark blue points represent the low density medium cells. The black line in Fig. 12 represents the analytical solution of a particle moving on a circular Keplerian orbit, while the black dashed line and the back dash-dotted line indicate orbits with eccentricity
$0.9$
and
$\sim 1$
, respectively.
Specific angular momentum versus specific orbital energy for every cell in the computational domain once the disk has formed (sim-0). The integration time corresponds to
$t=10.5$
yr. The density colour table is the same as in Fig. 4. The black line is the solution for a circular orbit. The black dashed line and the black dash-dotted line represent orbits with higher eccentricity (
$e=0.9$
and 0.9999, respectively). The time evolution of this plot is available at https://drive.google.com/file/d/10Kh7uXEJaRr8zY5GNmyQr24Vy9wC–G-/view?usp=sharing.
The majority of the high-density cells have negative orbital energy, which means that the material is bound to the companion star. On the other hand, the insert plot, showing a zoom-in near the origin, shows that some of the material has positive orbital energy, hence it is unbound. Some of the gas dropped by the nozzle onto the companion star is not initially bound to the system but slows down upon colliding with gas ahead of it. The cells in the high-density region are dispersed between the black line and the grey dashed line, which means that they move in closed orbits with eccentricities less than 0.9. The gas behaviour in this diagram is consistent with disk formation.
We finally determine whether the viscosity that enacts such large accretion rates is consistent with reasonable magnetic fields (noting that our simulations do not have magnetic fields and that in the simulations the viscosity is enacted by turbulent gas motion and possibly shocks). Using
$10^{15}$
G at the surface of the neutron star (Rea & De Grandis Reference Rea and De Grandis2025), one would find that the field at 1.3 R
$_{\odot}$
is of the order of a Gauss, following magnetic flux conservation for a dipole (
$B\sim r^{-3}$
). Such low magnetic flux is unlikely to be responsible for the accretion rates that we measure (or that we infer using the formalism of Shakura & Sunyaev Reference Shakura and Sunyaev1973). In fact using the expression of Wardle (Reference Wardle2007):
\begin{equation} \dot{M} \simeq 6\times10^{-12}\left(\frac{B}{1G} \right)^{2} \left(\frac{r_\mathrm{in}}{ 1.3 {\rm R}_\odot}\right)^{5/2} \left(\frac{M_\mathrm{NS}}{ 1.4\,{\rm M}_\odot}\right)^{-1/2} \mathrm{M}_\odot \mathrm{yr}^{-1},\end{equation}
where
$\dot{M}$
is the mass accretion rate, B is the magnetic field,
$r_\mathrm{in}$
is the accretion radius and
$M_\mathrm{NS}$
is the mass of the accretor, the deduced accretion rate would be very low indeed. To obtain instead a mass accretion rate of the order of
$4\times 10^{-3}$
M
$_{\odot}$
yr
$^{-1}$
, as is measured at the end of the simulation, would necessitate a field at the inner disk rim of 28 000 G. Such field could possibly be obtained by magneto-rotational instability amplification of a Gauss level field (Balbus & Hawley Reference Balbus and Hawley1991), over as little as
$\sim$
2 orbital periods of the gas at the inner rim (the amplification grows as
$\exp(3/4 \Omega t)$
, where
$\Omega$
is the orbital frequency and t is time). This would argue for a need to carry out these simulations in a magneto-hydrodynamics regime (e.g. Pjanka & Stone Reference Pjanka and Stone2020).
4. Sensitivity of results to the choice of some physical and numerical parameters
4.1 Mass injection rate
The full evolution of the mass transfer rate from the moment of Roche lobe overflow to the moment of CE as seen in the 1D mesa simulation lasts 30 000 yr (Fig. 2) and goes from
$1.6\times 10^{-9}$
M
$_{\odot}$
yr
$^{-1}$
to
$9.7\times 10^{-2}$
M
$_{\odot}$
yr
$^{-1}$
(Table 2). This entire period cannot be modelled in 3D which by necessity can only simulate a much shorter time (21 yr for us). The choice was therefore made to model a period of time towards the end of the 30 000 yr modelled in 1D between 29 500 and 29 521 yr, when the mass injection rate goes between
$2.3\times 10^{-4}$
and
$9.7\times10^{-2}$
M
$_\odot$
yr
$^{-1}$
(see sim-0 in Table 5).
Using sim-0 as reference, we computed three additional simulations with different initial mass transfer rates (called ‘sim-mdot-
$\#$
’ in Tables 3 and 5) and hence different start times in the context of the mesa simulation. Therefore, each 3D simulation samples a different part of the
$\dot{M}$
vs. time curve in Fig. 2. In Table 5, we present a summary of the initial and final parameters of these simulations, noting that, besides the mass transfer rate, all other parameters are the same as the Reference simulation sim-0 (see Section 2.2).
Parameters of a sequence of 3D simulations aiming to assess the resilience of disk parameters (
$M_\mathrm{ disk}$
and
$\textrm{R}_\mathrm{disk}$
) to the choice of mass transfer rate and simulation length. The mass-transfer rate, donor mass, accretor mass, and orbital separation are selected at the mesa start time. The disk mass and radius are measured at the end of the simulation.
In Fig. 13, we present density slices in the orbital plane of these 3D simulations at their final time step (see the fifth column in Table 5). In the upper left panel, we show the simulation sim-mdot-1, with the smallest initial mass transfer rate of
$1.4\times10^{-8}$
M
$_\odot$
yr
$^{-1}$
, running for the first 128 yr of the 30 000 yr mesa simulation. As expected, only a small, low mass disk forms at such low mass transfer rate, even if the run time is relatively long.
Densities slices on the orbital plane at the end of the 3D simulations. The models that are shown are: sim-0 (lower left panel), sim-mdot-1 (upper left panel), sim-mdot-2 (upper right panel) and sim-mdot-3 (lower right panel). Note how the density colour bar may have different maximum limits.
In sim-mdot-2 with initial mass transfer rate of
$1.1\times10^{-5}$
M
$_\odot$
yr
$^{-1}$
(top right panel), starting towards the end of the 30 000-yr long mesa simulation (only 120 yr before sim-0) and running for 29 yr, a disk forms with a radius of just under 50 R
$_{\odot}$
and a mass that is not that dissimilar to that of sim-mdot-1, despite the mass injection rate being larger by three orders of magnitude (the run time was four times smaller).
In the lower left panel of Fig. 13, we present sim-0, our reference simulation. The initial mass accretion rate is only 20 times higher at the start than the previous simulation, sim-0, but it increases far more steeply than for sim-mdot-2. The disk size is of the same order as the one in sim-mdot-2, but the disk mass is
$\sim$
100 times larger.
Finally, sim-mdot-3 is very similar to sim-0 but starts 8 yr later with an initial mass transfer rate that is again,
$\sim$
20 times larger. It runs only 13 yr to the same end point as sim-0 (see the lower right panel, Fig. 13). The disk reaches a similar radius (
$\sim 43$
R
$_{\odot}$
) and a mass that is 1.6 times larger than the mass of the accretion disk of sim-0, despite the fact that the simulation started 8 yr later and therefore ran for only 13 yr compared to 21 of the reference simulation and injected slightly less mass. This is due to the fact that disk growth is not only dependent on the mass transfer rate and length of simulation, but also on the specific geometry of the flow, which dictates how much mass accretes through the inner boundary as well as the shape of the disk at the time it is measured, whereby the ‘edge’ of the disk as defined by our criteria (Section 3.1) can vary slightly.
With hindsight, these tests could be performed more systematically so as to gain a better idea of whether the disk parameters as stated in Table 5 are close to what we might expect to be the disk just before CE for a systems such as ours. With these tests as they are we can only state that the disk’s mass and radius are likely reasonable within a factor of
$\lesssim 2$
for the mass, and 15% for the radius. Given other sources of uncertainty this is a reasonable and sufficient statement for now.
4.2 Velocity through the nozzle
The distribution and kinematics of the gas at
$L_1$
(the ‘nozzle’) is calculated according to the prescription of Lubow & Shu (Reference Lubow and Shu1975), Ritter (Reference Ritter1988), and Jackson et al. (Reference Jackson, Arras, Penev, Peacock and Marchant2017). They used Bernoulli’s principle to describe the evolution of the gas moving from the donor’s surface toward
$L_\mathrm{{1}}$
. The gas above the donor photosphere moves with a velocity
$v_{{L_1}}$
$\ll$
c
$_{{T}}$
, where
$c_{{T}}$
is the isothermal sound speed, while near
$L_\mathrm{{1}}$
, the gas is assumed to reach a velocity comparable to the isothermal sound speed,
$v_{{L_1}}$
$\lesssim$
$c_{{T}}$
. After the gas passes
$L_\mathrm{{1}}$
, due to the pressure gradient, the gas free falls supersonically into the companion’s Roche Lobe.
In our simulation, we need to set an injection velocity because otherwise gas placed in the nozzle does not enter the computational domain at the prescribed rate. This initial nozzle velocity is therefore arbitrary, and we thus need to ensure that changing its value does not affect the disk parameters.
We test the dependency of the simulation on the velocity of injection (
$v_{\textrm{L}_{1}}$
) at the nozzle, by executing three different simulations using the same parameters as the reference simulation, but setting different nozzle velocities in the x-direction: a subsonic velocity of 7.75
$\times 10 ^{4}$
cm s
$^{-1}$
(sim-0; 0.01 in code units), the isothermal sound speed at the donor’s photosphere, or 5.74
$\times 10 ^{5}$
cm s
$^{-1}$
(sim-vel-1), and a supersonic velocity of 6.16
$\times 10 ^{5}$
cm s
$^{-1}$
(sim-vel-2). See Table 3 for a summary of all the simulations parameters.
The top panel of Fig. 14 shows the cumulative mass as a function of radius at 21 yr, similar to Fig. 6, the nozzle is located at a radius of 74 R
$_{\odot}$
(5.2
$\times10^{12}$
cm) from the inner boundary. For the reference simulation, sim-0, with the lowest injection velocity, the formation of the disk takes longer since the injected mass needs more time to leave the nozzle and to fall into the companion’s Roche lobe. The accretion disk in this simulation has 11 percent more mass than simulations sim-vel-1 and sim-vel-2. We conclude that relatively small differences in the injection velocity around the isothermal velocity value
$c_{{T}}$
, do not greatly affect the disk parameters.
The cumulative mass as a function of radius (top panel) and the vertical density profile (bottom panel) for models with different injection velocities at
$t=21$
yr. Solid dark green line:
$v_{\textrm{L}_{1}}$
= 7.75
$\times 10 ^{4}$
cm s
$^{-1}$
(model sim-0), dotted green line:
$v_{ L_1}$
= 5.74
$\times 10 ^{5}$
cm s
$^{-1}$
(sim-vel-1), dashed emerald line:
$v_{\textrm{L}_{1}}$
= 6.16
$\times 10 ^{5}$
cm s
$^{-1}$
(sim-vel-2, which has the necessary velocity to leave the nozzle every time step).
The density profile along the z-axis at 21 yr (Fig. 14, bottom panel) is similar in the three simulations at z-values close to the inner boundary, the disk scale height in the reference simulation at 21 yr (
$H_{sim-0}$
= 4.9 R
$_{\odot}$
) is similar to the disk scale height for sim-vel-1 (
$H_\mathrm{{sim-vel-1}}$
= 4.3 R
$_{\odot}$
) and the disk scale height for sim-vel-2 (
$H_\mathrm{{sim-vel-2}}$
= 5.0 R
$_{\odot}$
). The density at the edges of the box is an unbound low-density gas interacting with the even lower-density background; this gas does not affect the dynamics/structure of the accretion disk. The arbitrarily assumed velocity injection does not affect the final results of the simulation.
4.3 Sensitivity of the results to background density and temperature
Filling the background with low-density gas is an expedient to ensure that no cell is empty, which would cause the inability to calculate pressure gradients. The value of the background density is
$\unicode{x03C1}_\mathrm{bg}=2.60 \times\ 10^{-22}$
g cm
$^{-3}$
(1
$\times10^{-16}$
code units), the lowest viable value, below which the code does not run. As consequence of this background density choice, at the beginning of the simulation a low-density rarefaction wave propagates from the discontinuity between the inner boundary and the background moving outward and out of the computational domain. The rarefaction wave leaves the computational box before the injected material reaches the inner boundary. This rarefaction wave does not affect the movement of the injected material, and the evolution of the disk formation is shown in Fig. 4.
We repeat the simulation using the same background density but decreasing and increasing the background temperature to
$1\times10^{4}$
K (sim-bgT-1) and
$1\times10^{6}$
K (sim-bgT-2), respectively (note that sim-0 has a temperature of
$1\times10^5$
K). The adiabatic index is
$\unicode{x03B3} = 1.1$
in all simulations (see Section 4.5 for a discussion about the adiabatic index). In Fig. 15, we show the total cumulative mass as a function of radius (top panel) and the density profile of the disk along the z-axis (bottom panel) at the end of the simulation. The vertical density profile of the disk is similar for all three simulations, but the highest temperature (and pressure) simulation results in approximately half of the mass in the disk at 21 yr.
The cumulative mass as a function of radius (top panel) and the vertical density profile (bottom panel) for models with different background temperatures at
$t=21$
yr. Solid purple line:
$\textrm{T}_\mathrm{bg}$
=
$10^{5}$
cm s
$^{-1}$
(model sim-0), dotted red line:
$\textrm{T}_\mathrm{bg}$
=
$10^{4}$
cm s
$^{-1}$
(sim-bgT-1), dashed orange line:
$\textrm{T}_\mathrm{bg}$
=
$10^{6}$
cm s
$^{-1}$
(sim-bgT-2).
Increasing the background temperature by an order of magnitude increases the pressure of the background such that it is higher than the pressure inside the nozzle (see fifth column in Table 3), reducing the amount of mass that can flow out from the nozzle. In the case of sim-bgT-2 only half of the injected material manage to leave the nozzle compared to the reference simulation.
4.4 Convergence tests
We finally test the sensitivity of the simulation to spatial and temporal resolution. Our comparison models have the same parameters as sim-0 (see Table 3) with (
$72 \times\ 96 \times\ 32$
) cells at the coarsest level; we also preserved the size of the computational box (
$-d_{\textrm{L}_{1}}\leq x \leq 1.25d_{\textrm{L}_{1}}$
,
$-1.5 d_\mathrm{L_1} \leq y \leq 1.5 d_\mathrm{L_1}$
and
$-0.5 d_\mathrm{L_1} \leq z \leq 0.5 d_\mathrm{L_1}$
, where
$d_{ L_1}=83$
R
$_{\odot}$
). We repeated the simulation with increasing levels of refinements: 2 levels (as the low-resolution simulation), 3, 4, and 5 levels (high-resolution simulation). The inner boundary around the companion star has a constant radius of
$\textrm{R}_\mathrm{{in}}=1.3$
R
$_{\odot}$
in all four simulations, and the size of the nozzle is also not resolution dependent. Refinement occurs primarily in the higher density regions close to the accretor, as concentric circles on the orbital plane, and around the nozzle area; this high-resolution refinement expands along the line that connects the nozzle with the inner boundary on the XZ plane.
The cumulative mass as a function of radius (top panel), the vertical density profile (middle panel) of the accretion disk at
$t=21$
yr for different resolutions, and the mass accretion rate onto the companion (bottom panel) as function of time for different resolutions. sim-0 has 3 levels of refinement.
In Fig. 16, we show the convergent behaviour of the disk mass and scale height (top and middle panels). We measure the properties of the disk at the end of the simulation (21 yr; see Table 6). The results show that the disk properties converge across different resolution levels. For the highest resolution (5 levels) simulation, the disk mass shows a 12% discrepancy compared to the Reference simulation (sim-0; 3 levels of refinement), the disk radius and scale height as measured using the criteria established in Section 3.1 are consistent across the three highest resolution simulations.
Disk properties for the reference simulation (sim-0), at 21 yr, with different resolutions.
We test the convergence of the mass accretion rate onto the neutron star as we have done for other quantities. In the bottom panel of Fig. 16, it is clear that the accretion rate converges. At the highest resolution (5 levels of refinement), the total accreted mass is 7% smaller than for the reference simulation (3 levels of refinement).
4.5 Simulation dimensionality and adiabatic index
The equation of state, given by
$P=(\unicode{x03B3}-1)\unicode{x03C1} e$
, relates the gas pressure P to the internal specific energy e and adiabatic index
$\unicode{x03B3}$
. In the absence of dissipation, a fluid parcel’s internal energy changes only due to the
$P\,\mathrm{d}V$
work that it does on the surrounding fluid or vice-versa, yielding the adiabatic relationship
$P \unicode{x03C1}^{-\unicode{x03B3}} = \mathrm{constant}$
. A higher adiabatic index implies that a higher fraction of the internal energy is in the kinetic energy of the particles and manifests as pressure, resulting in lower gas compressibility. On the other hand, with a lower adiabatic index the majority of the work done during compression is ‘hidden’ in the internal excitation of the particles, so the pressure does not increase as significantly. The adopted adiabatic index also plays an analogous role in shock-heating, in which bulk kinetic energy is converted to internal energy.
The adoption of an isothermal gas was proposed by Lubow & Shu (Reference Lubow and Shu1975). Armitage & Livio (Reference Armitage and Livio2000) then showed, using an adiabatic index of 4/3 in their 2D simulations of a mass transfer into the Roche lobe of a neutron star, that an accretion disk could form. Subsequently, Makita et al. (Reference Makita, Miyawaki and Matsuda2000), MacLeod et al. (Reference MacLeod, Antoni, Murguia-Berthier, Macias and Ramirez-Ruiz2017), and Murguia-Berthier et al. (Reference Murguia-Berthier, MacLeod, Ramirez-Ruiz, Antoni and Macias2017) showed, using 3D simulations of mass transfer through a Roche lobe with various system parameters, that the maximum adiabatic index that would allow a disk to form was
$\unicode{x03B3} \leq 1.2$
, for simulations without radiation or cooling function implemented.
We tested these claims by first performing a 2D comparison simulation of the reference simulation (sim-0, Table 3), and three additional 2D simulations with larger adiabatic indices and then carrying out two 3D simulations with higher adiabatic indices than sim-0. For the 2D simulations, we assume the gas is restricted to the orbital plane. We drop the z-coordinate from the grid,
$\Delta z = 1$
for all cells, assuming symmetry along the z-axis, and the gas cannot be deflected in the z-direction. The code solved the two-dimensional Euler equations for inviscid gas with external forces
$\textbf{f}=\left(f_{x},f_{y}\right)$
. The simulations have the same physical parameters as the reference simulation (sim-0).
In Fig. 17, we present density slices of the four, 2D simulations with
$\unicode{x03B3}=1.1$
(top left panel, a 2D version of sim-0),
$1.2$
(top right panel),
$4/3$
(bottom left panel) and
$5/3$
(bottom right panel), after 21 yr. In all simulations, the injected material revolves around the companion and forms a high-density disk structure.
Density in the orbital plane for 2D simulations with different adiabatic indexes. Adiabatic index
$\unicode{x03B3}= 1.1$
is on the upper-left panel (a 2D version of the 3D sim-0),
$\unicode{x03B3}= 1.2$
in the upper-right panel,
$\unicode{x03B3}= 4/3$
in the lower-left panel, and
$\unicode{x03B3}= 5/3$
in the lower-right panel). All simulations are plotted at
$t=21$
yr.
Density slices in the orbital plane (upper panels) and perpendicular plane (lower panels) for 3D simulations with different adiabatic indexes (
$\unicode{x03B3}= 4/3$
left, and
$\unicode{x03B3}= 5/3$
right). The simulations are plotted at
$t=21$
yr, and can be compared with last panel of Fig. 4.
In Fig. 18, we show density slices of two 3D simulations performed with adiabatic indexes
$\unicode{x03B3}=4/3$
and
$\unicode{x03B3} = 5/3$
to investigate the reluctance to form a disk for higher values of
$\unicode{x03B3}$
observed by other authors. This figure should be compared with sim-0 (
$\unicode{x03B3} = 1.1$
) shown in the last column of Fig. 4. For
$\unicode{x03B3} = 4/3$
, only slightly higher than the value used in sim-0, some of the material gets deflected around the companion, which may prevent the formation of a stable disk. For
$\unicode{x03B3}=5/3$
, the pressure gradient dominates over the gravitational force of the companion so that a fraction of the material is deflected away from the centre and may inhibit the formation of a disk.
We quantify the difference between 3D simulations with different adiabatic indices (see last row of Table 4). The top panel of Fig. 19 shows the cumulative mass as a function of radius, from which we measure the radius and mass of the disk-like structure in each simulation. Even though the simulations have the same mass transfer rates as the reference simulation, the high-density structure formed on the simulation with an adiabatic index
$\unicode{x03B3}=4/3$
is five times smaller (15.2 R
$_{\odot}$
or 1.06
$\times 10^{12}$
cm) and two orders of magnitude less massive (5.37
$\times 10^{-5}$
M
$_\odot$
). In the case of the simulation with
$\unicode{x03B3}=5/3$
a radius and a mass value cannot be determined from the cumulative mass plot.
In the bottom panels of Fig. 19, we compare the density profile along the z-axis of the simulations with different thermal properties. It is clear that when the adiabatic index is closer to the isothermal value the vertical structure is more defined, as expected from the bottom panels of Fig. 18, while for the higher adiabatic indices there is only a very marginal equatorial compression.
These results are consistent with those of Makita et al. (Reference Makita, Miyawaki and Matsuda2000), MacLeod & Ramirez-Ruiz (Reference MacLeod and Ramirez-Ruiz2015), and Murguia-Berthier et al. (Reference Murguia-Berthier, MacLeod, Ramirez-Ruiz, Antoni and Macias2017), who showed that disk formation depends primarily on its thermal properties: a disk forms only with an adiabatic index lower than
$\unicode{x03B3}\leq 1.2$
when the gas is cooler and has more compressibility. Although in the case of
$\unicode{x03B3}=4/3$
, a thick disk seems to be formed in our simulations.
5. Discussion and conclusion
The goal of this paper was to study the formation of an accretion disk around a neutron star (1.4 M
$_{\odot}$
) due to unstable mass transfer from an intermediate mass red supergiant (7 M
$_{\odot}$
) through Roche lobe overflow (RLOF). This phase likely immediately precedes a CE in-spiral with concomitant accretion onto the companion and possibly the formation of a jet that may affect the pre-in-spiral, as well as the in-spiral phases. Such feedback may lead to a different outcome than what has been modelled thus far by 3D hydrodynamic simulations without feedback (e.g. Lau et al. Reference Lau2022) or by 3D hydrodynamic simulations with a simplified version of feedback (e.g. Hillel, Schreier, & Soker Reference Hillel, Schreier and Soker2022).
By necessity, with an explicit 3D simulation, we only model a short time (21 yr) of the evolution of the RLOF phase predicted to last 30 000 yr by a 1D, implicit model. By carefully choosing to model in 3D the last phases of the mass-transfer before the putative CE in-spiral, we show that the 3D disk mass is likely only a factor of a few smaller than it might be if the entire phase were modelled, and very similar in radius.
We show that the accretion disk in our system grows to a mass of
$\sim 5 \times 10^{-3}$
M
$_\odot$
, a radius of
$\sim 40\textit{R}_\odot$
, and a scale height of
$\sim 5\textit{R}_\odot$
, just before it presumably goes into the CE in-spiral phase. This disk has approximately Keplerian rotation near the inner boundary, while the outer regions have a rotation velocity slightly slower than Keplerian due to pressure support. The temperature profile between the inner boundary and the disk’s outer radius can be fit with an exponential law with index
$-1.1$
steeper than
$-0.75$
predicted for accretion disks. An immediate improvement to understanding the disk’s temperature (and structure) would be to include explicit cooling. At the one significant figure level, these results are resilient with respect to various physical and numerical choices. The results are well converged with respect to spatial and temporal resolution.
The accretion rate through the inner boundary that surrounds the companion reaches
$4\times 10^{-3}$
M
$_\odot$
yr
$^{-1}$
, at the end of the simulation, which is consistent with the accretion rate predicted using a Shakura & Sunyaev (Reference Shakura and Sunyaev1973) formalism with
$\alpha=0.1$
. We show that the high accretion rate is driven by turbulent gas motion and possibly shocks, rather than by numerical viscosity. The influence of magnetic fields on the accretion rates cannot be measured in our non-magnetic code, but it remains a possibility in nature. While disk magnetic fields at the inner boundary of the disk should be Gauss-level, such magnetic field strength could be amplified rapidly by the magnetic-rotational instability (Balbus & Hawley Reference Balbus and Hawley1991) – within a couple of orbital periods and their effect on disk accretion could therefore be substantial.
The Eddington limit for mass accretion onto a neutron star of 1.4 M
$_{\odot}$
is on the order of
$10^{-8}$
M
$_\odot$
yr
$^{-1}$
. X-ray binaries are known to accrete at rates that can be 100 times the Eddington rate. That said the rate modelled for systems accreting in an unstable mass transfer regime can be larger than that (e.g. 1 000 times, see Dickson Reference Dickson2024). Extreme mass accretion rates predicted in the last phases of mass transfer before the CE in-spiral are likely far more complex physical phenomena, with likely extreme X-ray feedback. Moreno Méndez (Reference Moreno Méndez2022) explained the formation of high mass X-ray binaries, with a black hole as a companion, with hyper accretion rates (
$\sim 1 \times 10^{8}\dot{M}_\mathrm{Edd}$
) before the CE phase.
Also, the measured maximum accretion rate of
$4\times 10^{-4}$
M
$_\odot$
yr
$^{-1}$
could result in a jet mass loss rate of
$4\times 10^{-5}$
M
$_\odot$
yr
$^{-1}$
. An escape velocity at the inner boundary of 644 km s
$^{-1}$
, would result in a mechanical luminosity of
$\sim3\times 10^3$
L
$_\odot$
. On the other hand, if the jet is launched close to the surface of the neutron star (which is not resolved in our numerical simulations), the escape velocity would be
$\sim 0.3$
c, and the mechanical luminosity would be
$\sim 10^8$
L
$_\odot$
. Understanding the potential for jet formation from disks such as this one would be the next critical step because it may impact the early CE, even if the accretion disk and jet do not survive entering the envelope (Murguia-Berthier et al. Reference Murguia-Berthier, MacLeod, Ramirez-Ruiz, Antoni and Macias2017). Another possibility in this case, could be that the companion jets can prevent the companion from entering the CE. Instead, the companion undergoes ‘grazing envelop evolution’, where the companion slowly enters the donor’s envelope while its jets help unbind the envelope (Shiber & Iaconi Reference Shiber and Iaconi2024). Further investigation on the survival of the accretion disk during the CE phase for this system is intended on a follow up study.
Acknowledgements
We acknowledge the computing time granted by DGTIC UNAM on the supercomputer Miztli (projects LANCAD-UNAM-DGTIC-281 and LANCAD-UNAM-DGTIC-321). AJG and FDC acknowledge support from the UNAM-PAPIIT grant IN113424. AJG also acknowledges funding support from Macquarie University through the International Macquarie University Research Excellence Scholarship (‘iMQRES’) and the Postgraduate Research Fund (’PGRF’).
Data availability statement
The data underlying this article will be shared on reasonable request to the corresponding author.
Appendix A. Numerical considerations
Appendix A.1 Conservation of mass, energy and angular momentum
As we explain in Section 2.2.2, the Reference simulation has constant mass injection through the nozzle, we have external outflow boundaries and the inner inflow boundary. To test for conservation, we run the reference simulation, for 9.5 yr after which we switch off the nozzle, change the outflow boundary condition into reflective boundaries, and remove the outflow boundary around the companion star in the middle of the domain. We then run the simulation for an additional 11.5 yr, during which we evaluate the conservation of mass, angular momentum, and energy.
Total mass (top panel), total angular momentum (middle panel) within the computational domain as a function of time. The bottom panel shows total kinetic energy (Kin energy), potential energy (Pot energy), thermal energy (Therm energy), and total energy (Tot energy) as a function of time. During the grey shadow area, the mass injection is on, and internal inflow boundary is present on the computational box. The dashed line in each panel represents the comparative value after the injection is over.
For each time step, we calculate the total mass, angular momentum, and energy in the inertial frame of reference. For this we integrate the mass, angular momentum in each cell at every time step. In the case of the total energy, we take into account the kinetic energy of each cell, the potential energy between the fluid and the companion star, and the thermal energy of each cell. The integrated values of mass, angular momentum and energy are shown in Fig. A1, where the grey area represents the time when the nozzle and mass injection are on, causing the total mass, energy and angular momentum to increase over time. The dashed light line in each panel represents the first value of mass, energy angular momentum after the nozzle is switched off.
As shown in the top panel of Fig. A1, the mass within the reflective boundaries remains constant throughout the simulation, with a maximum variation of
$0.02$
% at 21 yr. As expected, due to numerical effects such as numerical viscosity, and the way the conservation of the angular momentum equation is discretised, in the middle panel of Fig. A1, the total angular momentum grows up to 6% over the 10 yr the simulation has reflective boundaries. For the measure of the total energy within the computational box (bottom panel in Fig. A1), after turning off the nozzle, there is a decrease in kinetic energy over 0.15 yr. This change is reflected in a step increase in thermal energy of 2.5
$\times10^{18}$
erg s
$^{-1}$
. The energy has a maximum change of 13% after the mass injection finish and remains constant until the end of the simulation. The extent of the non-conservation of the preceding values justifies the approximations made in the simulations.
Open access
