© The Authors 2023

1. Introduction

Close binary systems with a black hole (BH) component have been studied since the first detection of accreting BHs in binary systems by Roche lobe overflow (RLOF) in the 1960 decade, alongside with the first missions containing X-ray detectors (see, e.g., Giacconi et al. 1962; Lewin et al. 1967). The material lost by a normal companion of low mass, known as low mass X-Ray binary (LMXB) or of high mass, namely, a high mass X-Ray binary (HMXB) forms an accretion disk around the BH. Mass and angular momentum are thereby transferred to the BH, releasing an intense X-ray flux. This particular group of binary systems has been studied from both, an observational and a theoretical point of view (among the most recent works: e.g., Ivanova et al. 2017; Langer et al. 2020; Fukumura et al. 2021; Mata Sánchez et al. 2021; Mikołajewska et al. 2022; You et al. 2023).

V404 Cyg is a member of the LMXB family. It was discovered by the space satellite Ginga in May of 1989 as the transient X-ray source GS 2023+338 (Makino 1989). Its optic counterpart was identified as the variable star V404 Cyg (Wagner et al. 1989). Later, Charles et al. (1989) identified the source as an LMXB. The binary has an orbital period of P = 6.473 ± 0.001 d (Casares et al. 1992) and a mass function of f(M) = 6.08 ± 0.06 M_⊙ (Casares & Charles 1994). This high value of the mass function suggests the nature of the accretor is that of a BH. The mass ratio was determined by Casares et al. (1992) as $q\equiv M_{\mathrm{d}}/M_{\mathrm{BH}}=0.{06}_{-0.005}^{+0.004}$, with M_d and M_BH as the masses of the donor star and the BH, respectively. The companion was confirmed as a giant star when Khargharia et al. (2010) determined its spectral type as K3 III. They also found the binary’s inclination, $i=67{°}_{-1°}^{{}^{+}3°}$. Knowing all these parameters, the determination of the masses of the components is immediate, namely, M_d = 0.54 ± 0.05 M_⊙ and $M_{\mathrm{BH}}=9.0_{-0.6}^{+0.2}{M}_{\odot }$, for the donor and the BH, respectively. In 2009 a precise estimation of the distance was taken, giving a value of d = 2.39 ± 0.14 kpc by the measure of the parallax of the system on radio waves (Miller-Jones et al. 2009). Ziółkowski & Zdziarski (2018, henceforth ZZ18) presented, based on these observational data, the radius of the donor star, $R_{\mathrm{d}}=5.{50}_{-0.18}^{+0.17}$ R_⊙, the effective temperature (from the spectral type) $T_{\mathrm{eff}}={4274}_{-113}^{+116}$ K and the luminosity of the donor star, $L_{\mathrm{d}}=8.7_{-1.4}^{+1.7}{L}_{\odot }$.

Many studies have sought to calculate the accretion rate onto the compact object during the various outburst that had taken place in 1938, 1989, and 2015 (Chen et al. 1997; Życki et al. 1999; Motta et al. 2017). The system presented two episodes of outburst in 2015, namely: in July (Barthelmy et al. 2015) and in December (Martí et al. 2016; Motta et al. 2016). Due to the large absorption reported (Kimura et al. 2016), it was challenging to construct an X-ray luminosity curve this year and, thus, to estimate a mass accreted during this event. With the information provided by the above-mentioned authors, ZZ18 stated that the value ⟨Ṁ_BH⟩ = 4.0 × 10⁻¹⁰ M_⊙ yr⁻¹ is most likely an upper limit for the accretion rate onto the BH in V404 Cyg. This value is lower than the estimated mass loss rate for the donor star ⟨−Ṁ_d⟩ = 1.1 × 10⁻⁹ M_⊙ yr⁻¹ that was predicted using the Eq. (25a) from Webbink et al. (1983) for the system V404 Cyg. ZZ18 also re-obtained these values using an updated evolutionary model. The difference between the amount of mass lost by the donor and the mass accreted by the BH has been attributed to the mass that gets lost from the system, advecting angular momentum along with it. These mass and angular momentum losses make the system evolve in a non-conservative way. In addition, there are observational indications that V404 Cyg is currently losing mass (Muñoz-Darias et al. 2016).

In this work, we consider non-conservative close binary evolutionary models with the objective of reproducing observational data available for the main parameters of V404 Cyg, with the aim of obtaining a possible progenitor for the system. We also predict possible results for the evolution of the system’s donor and show some theoretical parameters that we expect to be observationally measured in the future, such as the time derivative of the orbital period. On the other hand, we also study the evolution of the spin parameter in the context of our models. We specify the numerical code used in Sect. 2, show the results of our models in Sect. 3, and present our conclusions in Sect. 4.

2. The binary evolution code

The main tool for this work is the binary evolutionary code described in Benvenuto & De Vito (2003), De Vito & Benvenuto (2012), Benvenuto et al. (2012). When components remain detached, it works as a standard evolutionary code for isolated stars. In the case of semi-detached configurations, our code includes the mass transfer rate, Ṁ₁, as a new variable in the difference equations. Then, the mass of the donor is M₁ = M^1,prev + Ṁ₁Δt, where M_1, prev is the mass of the donor in the previous stage and Δt is the time step. As M₁ appears in the equations of the structure of the entire model, Ṁ₁ is treated as a global variable to be solved. This is in contrast with all the other variables that are local and meant to be relaxed. When handling the corresponding generalized Henyey matrix, this treatment involves a non-zero column. The resulting matrix equation can be solved with a slight modification of the standard algebra. This solves the structure of the donor star, the orbital evolution, and the value of the mass transfer rate simultaneously in a fully implicit way, which makes the algorithm numerically stable. A detailed explanation of the procedure is given in Benvenuto & De Vito (2003). We assume that the mass is only transferred via RLOF. As for opacities, we used OPAL libraries (Iglesias & Rogers 1996) for temperatures of T ≥ 10⁴ K and molecular opacities computed by Ferguson et al. (2005) for lower values of T. A detailed description of how the code works may be found in Benvenuto et al. (2012).

Regarding the abundances assumed for our models, in the first step of our calculations, we followed ZZ18 to employ the solar metallicities. We have set them to X = 0.710, Y = 0.276, and Z = 0.014, whereas the mixing length parameter has been set to α_MLT = 1.50. With these values, our code is able to compute a solar structure compatible with observations at its present age. We remark here that these abundances are slightly different from those given in Asplund et al. (2021) who measured values of X = 0.7438 ± 0.0054, Y = 0.2423 ± 0.0054, and Z = 0.0139 ± 0.0006 at the surface of the Sun. If we set the abundances to these values, the Sun would be slightly under-luminous by 0.05 dex, which is a small discrepancy since the physical ingredients employed by Asplund et al. (2021) are different from ours. So, we decided to slightly adjust the initial abundances to produce a Sun compatible with observations. In the second step of the model calculations, we took into account the determination of abundances for the donor star in V404 Cyg, presented in González Hernández et al. (2011), and we employed X = 0.71, Y = 0.262 and Z = 0.028.

2.1. Non-conservative mass transfer and orbital evolution

For cases of conservative binary evolution calculations, total mass, and orbital angular momentum remain as constants. However, in analyzing the difference between the estimated mass loss rate from the donor component and the estimated accretion rate on the BH from V404 Cyg, it is commonly assumed that mass gets lost in a non-conservative mass transfer episode advecting angular momentum away from the system (Webbink et al. 1983; Chen et al. 1997; Życki et al. 1999; Motta et al. 2017; ZZ18). This is expected to occur in some astrophysical scenarios of interest, so this phenomenon was included in the calculations.

We employed the usual equation to compute the evolution of the orbital semi-axis. This can be obtained using the definition of the angular moment combined with Kepler’s Third Law. The episode of non-conservative mass transfer is specified by two free parameters, as in Rappaport et al. (1982, 1983): (1) the fraction β of mass lost by the primary star¹ that is accreted by the secondary star, and (2) the specific angular momentum of matter lost away from the system α in units of the same quantity for the compact object. We assume that the orbit is always well approximated by a circle of radius r_orb (where r_orb is a function of time) and we have neglected the rotational angular momentum of the components. In the case where the angular momentum is lost only by mass ejection from the system, Eq. (1) can be rewritten using Kepler’s third law and the expression of total angular as, as the following differential equation we find:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{J}_{\mathrm{ME}}}{\mathrm{d}t}=\alpha (1-\beta )\sqrt{GM{r}_{\mathrm{orb}}}(\frac{M_1}{M}{)}^2\dot{M_1},\end{array}$$(1)

where G is the gravitational constant and M = M₁ + M₂ is the total mass of the system, with M₁ and M₂ the masses for the donor star and the BH, respectively.

Angular momentum can also be lost from the system by gravitational radiation and it is calculated according to the standard formula (Landau & Lifshitz 1971):

$$\begin{array}{c}\hfill \frac{\mathrm{d}\ln (J_{\mathrm{RG}})}{\mathrm{d}t}=-\frac{32G^3\mu }{5c^5}\frac{M^2}{r_{\mathrm{orb}}^4},\end{array}$$(2)

where c is the vacuum speed of light and $\mu =\frac{M_1{M}_2}{M}$.

The code also considers angular momentum loss due to magnetic braking, using the prescription of Rappaport et al. (1983), based on the magnetic-braking law of Verbunt & Zwaan (1981),

$$\begin{array}{c}\hfill \frac{\mathrm{d}{J}_{\mathrm{MB}}}{\mathrm{d}t}=-3.8\times {10}^{-30}{M}_1{R}_1^4{\omega }^3\mathrm{dyn}\mathrm{cm},\end{array}$$(3)

where ω is the angular rotation frequency of the donor star, assumed to be synchronized with the orbit, and R₁ is the donor’s radius. The code includes full magnetic braking when the star has a sizable convective envelope embracing a mass fraction ≥0.02.

Replacing Eqs. (1)–(3) in the expression for the evolution of the angular momentum and considering Ṁ₂ = −βṀ₁ from the definition of β, we obtain a differential equation for the orbital separation, which has no analytical solution.

2.2. Eddington limit and black hole spin parameter

This is the first work in which we employ our code to calculate the evolution of a binary system with a BH; in previous papers, the companion was a neutron star or a normal star. Thus, we had to change the accretion efficiency of the compact object. Furthermore, we are interested in calculating the evolution of the spin angular momentum of the BH as it receives mass and angular momentum from its companion. For that purpose, we followed the prescriptions given in Podsiadlowski et al. (2003, henceforth PRH03).

The luminosity released due to accretion onto the BH is:

$$\begin{array}{c}\hfill L_2=\eta {\dot{M}}_2{c}^2,\end{array}$$(4)

where Ṁ₂ is the BH accretion rate, c is the vacuum speed of light, and η is the efficiency with which the BH radiates, determined by the last stable particle orbit. This parameter can be expressed as:

$$\begin{array}{c}\hfill \eta =1-\sqrt{1-\left(\frac{M_2}{3M_{\mathrm{BH}}^0}\right)},\end{array}$$(5)

where the quantities $M_{\text{BH}}^0$ and M₂ are the initial and present mass of the BH, respectively.

Equaling the BH luminosity L₂ with Eddington’s luminosity and assuming spherical accretion, an expression for the maximum accretion rate onto de BH can be obtained (see PRH03) as:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{Edd}}\simeq 2.6\times {10}^{-7}{M}_{\odot }{\mathrm{yr}}^{-1}\left(\frac{M_2}{10M_{\odot }}\right){\left(\frac{\eta }{0.1}\right)}^{-1}{\left(\frac{1+X}{1.7}\right)}^{-1},\end{array}$$(6)

where X is the hydrogen mass fraction. The BH accretion rate is limited by this value through its evolution.

On the other hand, the accretion phenomena not only affects the mass of the BH. As it accretes matter, and since that matter carries angular momentum with it, the spin parameter of the BH defined as $a^{*}\equiv c{J}_2/G{M}_2^2$ also increases, according to:

$$\begin{array}{c}\hfill a^{\ast }={\left(\frac{2}{3}\right)}^{\frac{1}{2}}\frac{M_{\mathrm{BH}}^0}{M_2}\{4-{[18{\left(\frac{M_{\mathrm{BH}}^0}{M_2}\right)}^2-2]}^{\frac{1}{2}}\}.\end{array}$$(7)

It is important to remark that these expressions are adequate for an initially non-rotating BH (see PRZ03) and are valid when $M_2<\sqrt{6}{M}_{\mathrm{BH}}^0$, which is the case for all our calculations (for a detailed treatment see Bardeen 1970; King & Kolb 1999).

3. Models and results

Our primary objective in this work is to obtain possible progenitors for the binary system V404 Cyg. With this goal, we analyzed the results generated by different sets of initial parameters: masses of the donor star and the BH ($M_{\text{d}}^0$ and $M_{\text{BH}}^0$, respectively), the orbital period $P_{\text{orb}}^0$ and the fraction β of the mass lost by the donor that is accreted by the BH (as defined in Sect. 2.1). We fixed the free parameter that describes the specific angular momentum of matter lost as α = 1.

An important consideration that stands out is that this code calculates the evolution of the donor star with an already existing BH. That is to say, we are not discussing how the compact object formed and we are also avoiding the common envelope phase. The code assumes that the orbit of the system components is circularized, so the circular restricted three-body problem can be applied. This is a very reasonable assumption for V404 Cyg, where we can get the orbital eccentricity using the definition of the mass function and the observational estimations for the orbital period, the K semi-amplitude, and the mass function given by Casares et al. (1992) obtaining a value of e ∼ 0.024.

3.1. Models with solar abundances

We computed 72 evolutionary sequences with solar abundances exploring the combinations of initial donor masses of 1.5 and 2.0 M_⊙, initial BH masses of 7, 8, and 9 M_⊙, initial orbital period of 0.75, 1.00, and 1.25 days and values of 0.1, 0.3, 0.7 and 1.0 for the β fraction. The identification of each model was built using the first character associated with the initial masses of the model (see Table 1), followed by labels that are related to the initial orbital period and the β parameter. For example, the model C_075_07 was calculated with initial masses of 1.5 M_⊙ and 8 M_⊙ for the donor and the BH, respectively, as well as an initial orbital period of 0.75 d and β = 0.7.

We introduce the function ϵ² that helps us to determine how close the parameters given by our models are from the ones observationally acquired. This quantity is defined as:

$$\begin{array}{c}\hfill {ϵ}^2=\sum _i{ϵ}_i^2,{ϵ}_i=\frac{E_i-E_i^{\mathrm{obs}}}{E_i^{\mathrm{obs}}}\end{array}$$(8)

where E_i and $E_i^{\text{obs}}$ are the model and observational data for the parameter i, and i = 1, …, 5 correspond to the BH mass (M_BH), the donor’s mass (M_d), the orbital period (P_orb), effective temperature (T_eff) and luminosity (L_d) or radius (R_d) of the donor star, respectively. For example, ${ϵ}_{\text{1}}=(E_{\text{1}}-E_1^{\text{obs}})/E_1^{\text{obs}}$ means ϵ_MBH = (M₂ − M_BH)/M_BH. The time dependence for ${ϵ}_i^2$ will be given by the change of the quantities along the evolutionary sequences. on the modeled parameter on the evolution code. All values used for the observational measures can be found in Table 2. For ϵ² = 0, we have a case where all the i parameters calculated from the models are equal to the ones obtained observationally. This makes it easier to see that this quantity helps us to determine when and how the models represent the system’s observational data simultaneously.

Computing this function along the calculated sequences allows us to obtain a minimum value for each of them. This minimum value corresponds to the time when the quantities modeled are closest to the ones observationally estimated. We consider that the model represents the characteristics of the observed system well enough when the minimum value of the ϵ² function is lower than 0.0485, which is the value of the sum obtained when $E_i=E_i^{\mathrm{obs}}+{\sigma }_i^{\mathrm{obs}}$ for each i, where ${\sigma }_i^{\mathrm{obs}}$ is the observational error for each parameter that can be found on Table 2. This restriction guarantees that the modeled parameters are not far from their observational uncertainties. Only two of our models calculated with solar abundances satisfy this condition, namely: E_100_01 with ${ϵ}_{\text{min}}^2=0.0449$ and E_100_03 with ${ϵ}_{\text{min}}^2=0.0482$, as seen on Fig. 1. In Sects. 3.1.1 and 3.1.2, we present our results for these models.

Fig. 1. Quantity ϵ2 as a function of time for our best models with solar composition: E_100_01 and E_100_03. With a dashed line and a grey area is indicated the value of ϵ2 = 0.0485 and the acceptance region we considered for our models. With a black solid line, the value ϵ2 = 0 represents the situation where all the parameters modeled are equal to the ones observed simultaneously.

3.1.1. Mass transfer

As we already said, the observational data for the mass transfer episode for V404 Cyg seems to imply that a part of the mass gets lost from the system, advecting angular momentum away.

Among other studies, ZZ18 have studied the non-conservative mass transfer episode for this system. In their work, they estimated the accretion rate over intervals between outbursts and got a likely upper limit for it, with a value of ⟨Ṁ_BH⟩ = 4 × 10⁻¹⁰ M_⊙ yr⁻¹. They also stated that their models with β ≲ 0.33 were the ones that are in agreement with this limit. For the donor mass loss rate, the existing estimation is of Ṁ_d = 1.1 × 10⁻⁹ M_⊙ yr⁻¹, value taken from Webbink et al. (1983) using the equation 25a with the V404 Cyg system’s parameters.

Our results for the mass transfer episode are resumed in Fig. 2, where the donor mass loss rate (Ṁ₁) for the best models is shown. In Fig. 3, we show the accretion rate on the compact object (Ṁ₂) for the same models.

Fig. 2. Donor star mass loss rate for each of the best models with solar composition: E_100_01 and E_100_03. As the only variation between these models is on the β value, and this parameter does not affect strongly the mass loss episode, their plots are mostly overlapping. The black dashed horizontal lines represent the estimation for the mass loss rate for V404 Cyg of Ṁd = 1.1 × 10−9 M⊙ yr−1 is denoted with a dashed horizontal line. As for the dashed vertical lines, they represent the times of the minimum value of the epsilon squared function.

Fig. 3. BH accretion rate for the models with solar abundances of the E group with $P_{\text{orb}}^0$ = 1 d. The qualitative form of the functions is related to the mass loss rate of the donor star by the relation Ṁ2 = −βṀ1. This is, for a higher value for β larger would be the accretion rate onto the compact object. The shaded area represents the zone below the likely upper limit for the mass accretion rate (⟨Ṁ2⟩ = 4 × 10−10 M⊙ yr−1 given by ZZ18. Models with β ≤ 0.3 agree with this upper limit and the results reached by the mentioned authors.

For the first case (shown in Fig. 2), the existing estimation has been featured with a horizontal dashed line. We have highlighted the age predicted by our models when the epsilon squared function reaches its minimum value (see Fig. 1) with a vertical dashed line. At this time, the mass loss rate of the donor star for models E_100_01 and E_100_03 nicely agrees with the above-quoted estimation while for the conservative model C_125_10 we find an appreciable difference (of ∼37%).

As for the mass accretion rate onto the BH (Fig. 3), we added to our two best models the ones from the same group and initial orbital period so the effect of the variation of β becomes evident. As is typical in the computation of close binary models, this parameter relates the mass loss by the donor with the one that is accreted by the compact object as Ṁ₂ = −βṀ₁. This relation makes the mass accretion rate function similar in form to the mass loss rate function of the previous figure, but a change in the β parameter does not provoke a strong variation to Ṁ₁ as to Ṁ₂. The shaded area on this figure represents the values below the likely upper limit suggested by ZZ18. We found that both of our models with β ≤ 0.3 remains under this limit, confirming previous results.

3.1.2. Orbital period and BH spin parameter

A non-conservative mass transfer episode in a binary system essentially consists of a fraction of the mass transferred by the donor being accreted by the companion, while the rest is lost from the system. The mass accreted by the BH accelerates its rotation. This phenomenon can be studied by analyzing the evolution of the BH spin parameter $a^{*}=c{J}_{\text{BH}}/G{M}_{\text{BH}}^2$, with J_BH as the rotational angular momentum of the BH.

The only available data for the spin parameter of V404 Cyg BHs is given in Walton et al. (2017). This study used a spectral analysis of NuSTAR X-ray observations of V404 Cyg for its 2015 outburst and proposed different models that fit the observations using the reflection method. Although these authors obtained multiple solutions for the BH spin parameter, they stated that the most robust one to be a^* > 0.92 with a 99% of statistical uncertainty.

For the computation of this parameter, assuming that the BH is initially not rotating, we employed Eq. (7). We obtained the evolution for the BH spin parameter over time, as shown in Fig. 4. In this graph (and Table 3), we show our results for all the models from group E with an initial orbital period of 1 d, so it becomes evident that the larger the β value, the faster the BH’s final rotation.

Fig. 4. BH spin parameter for the V404 Cyg compact object for each of the best models, according to Eq. (7). As the BH accretes matter it spins up, provoking the increase of $a^{*}=c{J}_2/G{M}_2^2$. Thus, the lower the β, the lower the BH spin. The vertical dashed lines are the predicted ages for V404 Cyg. The values of the BH spin for these times are given in Table 3.

The solutions for the BH spin parameter at the presumed ages for the system (see Table 3) correspond to a slow rotation regime when the existing estimations (even the ones for slow rotators obtained by Walton et al. 2017) are much higher.

As matter leaves the system it carries away angular momentum, as described by Eq. (1). Other effects that also modify the orbital period are gravitational radiation (Eq. (2)) and magnetic braking (Eq. (3)). In this work, we consider all three effects together, ultimately finding that the orbital period mostly increases with time, reaching values of P_orb = 46–48 d at the end of our calculations (age of the donor star of 14 Gyr). The evolution of this quantity for each of our best models is shown on the top panel of Fig. 5, where the initial orbital period is $P_{\text{orb}}^0$ = 1.00 d for both. Our results are in good agreement with the well-determined value of P_orb = 6.47 d, when ϵ² has its minimum value. As for the bottom panel, we show the time derivative of the orbital period of our models, calculated for each time as an approximation of an incremental quotient. The results for this quantity and the characteristic timescale (P_orb/Ṗ_orb) at the expected ages of the system can be found in Table 4. We found an increased timescale of ∼2 × 10⁸ yr, which is in good agreement with the one predicted by ZZ18.

Fig. 5. Orbital period as a function of time for our best models with solar metallicity (solid lines), observational estimation for the present orbital period Porb = 6.47 d (dashed horizontal line) and age of the system for the system predicted for each model (dashed vertical lines) shown at the top. The bottom panel shows the time derivative for the orbital period for the same models and ages.

Although the time derivative of the orbital period would be very useful to test the evolutionary scenario, this quantity is not yet known and there are no prospects for deriving it any time soon. For measuring such quantity with an adequate degree of certainty, we would need a time basis that is far longer than what is presently available. King & Lasota (2021) stated that this time basis should be long enough for the radii of the donor star and its Roche lobe to vary at least on a density scale height, H_ρ (H_ρ ≡ −dr/dlnρ), which corresponds to thousands of years. Any measurement in the near future will surely reflect the occurrence of short-timescale phenomena, neglected in our calculations. In this sense, the values of Ṗ_orb we have presented above are related to the ingredients considered for a modeling of the evolution of the whole system (magnetic braking, gravitational radiation, and mass loss from the system).

3.2. Models with higher metallicity

As V404 Cyg is an object that belongs to the field of our Galaxy, we begin our exploration using solar metallicity. As described above and as done in ZZ18, we initially considered solar abundances. Nevertheless, González Hernández et al. (2011) presented a chemical abundance analysis for the donor star, and obtained [Fe/H] = 0.23 ± 0.19. This value is well addressed with a metallicity of Z = 0.028, two times the value corresponding to the Sun. Therefore, we calculated 30 additional models with this new value of Z, fixing the hydrogen abundance on X = 0.71. For this instance, we fixed the value of the initial donor’s mass at 1.5 M_⊙ and explored only the values of β = 0.3 and 0.1, based on the results obtained from the solar metallicity analysis. For the values of the initial BH mass, we still considered $M_{\text{BH}}^0$ = 8, 9, and 10 M_⊙, and we explored the same interval of initial orbital periods, $P_{\text{orb}}^0$ = 0.75, 1.00, and 1.25 d, adding 1.50 and 1.75 d values to the analysis. These models are identified similarly to those corresponding to solar metallicity, but adding “_Z028” at the end of the name.

The best models we obtained are: E_150_01_Z028 and E_150_03_Z028, with minimum ϵ² values of 0.0119 and 0.0132, respectively (see Fig. 6). These models not only reach lower values than our acceptance one, but also each of their parameters gets within their respective observational uncertainty listed on Table 2 at ages of 5.170 and 5.176 Gyr².

Fig. 6. Quantity ϵ2 as a function of time for our best models with the change on metallicity: E_150_01_Z028 and E_150_03_Z028. The elements represented are the same as in Fig. 1. With a dashed line and a grey area is indicated the value of ϵ2 = 0.0485 and the acceptance region we considered for our models. With a black solid line, the value ϵ2 = 0 represents the situation where all the parameters modeled are equal to the ones observed simultaneously.

3.2.1. Mass transfer

A change in the metallicity of the donor star implies a change in the donor’s outer opacities and, thus, in its structure as well. As for the mass transfer episode, the results can be seen in Figs. 7 and 8. We estimate the present mass loss rate as 1.24 × 10⁻⁹ M_⊙ yr⁻¹, which is also in good accordance with the estimation given by ZZ18.

Fig. 7. Donor star mass loss rate for each of the best models with the metallicity change: E_150_01_Z028 and E_150_03_Z028. The elements represented are the same as in Fig. 2.

Fig. 8. BH accretion rate for the two best models with double solar metallicity. The elements represented are the same as in Fig. 3. The qualitative form of the functions is related to the mass loss rate of the donor star by the relation Ṁ2 = −βṀ1. This is, for a higher value for β larger would be the accretion rate onto the compact object. The shaded area represents the zone below the likely upper limit for the mass accretion rate (⟨Ṁ2⟩ = 4 × 10−10) given by ZZ18.

Considering the accretion rate, the model computed with β = 0.3 seems to exceed the upper limit of 4.0 × 10⁻¹⁰ M_⊙ yr⁻¹ at some parts of its evolution, but is below this limit near the present age. The other model, computed with β = 0.1, is still fully under the limit.

3.2.2. Orbital period and BH spin parameter

As the BH accretion episode has not changed quantitatively, a significant amount from the models with lower metallicity, it is expected that the results for the BH spin parameter do not differ too much from the ones presented above. This parameter evolution for the two models that have been taken into account is shown in Fig. 9. Once again, our models do not reach the observational estimation given by Walton et al. (2017).

Fig. 9. BH spin parameter evolution for the V404 Cyg compact object according to Eq. (7). The values of the BH spin for the predicted ages for V404 Cyg (vertical dashed lines) are very similar to the two best models computed with solar composition in Table 3.

The orbital period evolution considering these models can be resumed in Fig. 10 with some results given in Table 5. Models with this metallicity and initial orbital period reach the observed one for V404 Cyg at the same moment when the other quantities are still within their observational uncertainties while predicting a present value of Ṗ_orb ∼ 9.2 × 10⁻¹¹. This value deduces a characteristic increase time consistent with the one obtained by ZZ18. Our models with this metallicity deduced final orbital periods between 68 and 70 d.

Fig. 10. Orbital period and its time derivative as a function of time for models with a metallicity of Z = 0.028. The represented elements are the same as the ones shown in Fig. 5 but adequate to these models.

3.3. Donor star evolution and proposal of our bests progenitors

This work aims to model the characteristics of V404 Cyg to get a predecessor system and also to analyze its present and future evolution. With this in mind, we calculated our models up to an age of t = 14 Gyr. This allowed us to make some estimations for the complete evolution of the system.

Figure 11 shows the different evolutionary tracks for the models in the Hertzsprung–Russell diagram for some of our best progenitors. Here, we can offer some remarks. As for the tracks that are calculated with solar abundances, the mass-transfer episode occurs in three parts. The first one begins on the main sequence for both models (mass transfer episode in case A; Kippenhahn & Weigert 1967) and ends abruptly due to the contraction of the donor when the central hydrogen is exhausted. After a very short time, the second mass transfer episode begins. This behavior cannot be appreciated in the evolutionary tracks because of the small portion of the diagram covered while the donor star is detached from its Roche lobe. Once the donor evolves blueward and gets dimmer, a hydrogen thermonuclear flash occurs making the star raise its luminosity and expand. Then, it starts the third mass transfer episode. The whole thermonuclear flash event is very fast, making the third mass transfer episode look like a Dirac’s delta distribution.

Fig. 11. Evolutionary track in an HR diagram for two of our models. Top: model E_100_01 for the case of solar metallicity. Bottom: model E_150_01_Z028 for the case of higher metallicity. The cases of the best models with higher β value are practically overlapped with the models shown. On the left panels, the whole track is shown with a shaded area when the mass loss episode occurs. There can be seen the observationally estimated position of the donor star and the position for the age of the system predicted by our models (with circles on the evolutionary track). Once the star is getting near to becoming a white dwarf it suffers a thermonuclear flash, making it expand and lose mass again for a very short period of time. The end of these tracks are very similar, the donor becomes a low-mass white dwarf with a hydrogen-rich surface. A dashed line denoting R = 0.02 R⊙ was included; it approximately corresponds to the asymptotic radius for the white dwarf. In the right panels, the zoom is on the nearness of this estimation, where the shaded area considers the observational error for the effective temperature and luminosity.

Because the donor star is not massive enough to start the 3α reactions finishing the red giant branch, no helium burning occurs. Once the thermonuclear flash event has occurred, the star has lost almost 20% of the initial superficial hydrogen abundances. Our prediction for the final fate of the V404 Cyg donor star is to become a low-mass helium white dwarf with a mass of 0.28 M_⊙ and a radius of ∼0.02 R_⊙ with a hydrogen-rich surface.

As for the present state of the donor star, observational data can be seen in Fig. 11. These data were placed on the HR diagram using the estimations of L_d = 8.7 L_⊙ and T_eff ∼ 4200 K and its respective error bars (see Table 2). Also, our theoretical models predict that the donor star is currently on the red giant branch and getting near the end of the first mass transfer episode (remaining ∼0.2 Gyr).

The evolutionary tracks corresponding to a metallicity of Z = 0.028, shown in Fig. 11 (bottom panel), shares lots of characteristics in common with our previous analysis for the models with solar abundances. However, we do note that the mass transfer episode occurs in two parts, where the first one begins after core hydrogen exhaustion (case B of mass transfer episode; Kippenhahn & Weigert 1967). The remnant compact object is still a helium white dwarf of mass M ∼ 0.29 M_⊙. The estimated age of the system predicts that the donor star is losing mass on the main mass transfer episode and there are ∼0.3 Gyr remaining for the end of it. The entirety of the mass transfer episodes takes place within ∼1 Gyr.

In Fig. 12, we present the contribution for every parameter considered in the construction of the epsilon squared function for one of the best models with each metallicity. For the case of solar abundances (left panel), we can see that near the minimum value, the parameters that dominate the total epsilon squared function are the system orbital period and the donor’s luminosity with a great contribution of the donor’s mass, while with Z = 0.028 (right panel) the best models get a better estimation for this last parameter but it’s behavior is still dominated by the luminosity and orbital period. The evolution of binaries is very sensitive to the variation of these parameters since the orbital evolution of the system takes a primary role in the RLOF episodes and the donor’s initial mass determines the initial position on the ZAMS and the way it evolves. So, even with a thin grid on these parameters, we would have to get very specific for these initial parameters to find better models. As for the donor’s luminosity, this is the observed parameter that has the largest relative error. This is due to the fact that the system is located in the bulge of our galaxy, which is a very obscured area, so we do not rely much on this quantity.

Fig. 12. Epsilon squared function decomposed on the contributions of every parameter analyzed for two of our models. Left: model with Z = 0.014. Right: model with Z = 0.028. The models computed with the same initial parameters except for β have similar behavior. The shaded area represents the acceptance zone, where the function takes values lower than 0.0485 and the black horizontal line at ϵ2 = 0 represents where the modeled parameter equals the observed estimation for it.

4. Conclusions

Based on the calculations and analysis we performed in this work, we are in a position to propose a plausible progenitor for V404 Cyg. We considered, in the first step of calculations, solar abundances, and in the second step, a more metallic donor star. From the results given by our models and the analysis performed in the previous section, we found that two models with each metallicity considered could represent adequately the current state of V404 Cyg.

Considering the epsilon squared (ϵ²) function (see Eq. (8)), we selected the two best models for each metallicity: E_100_01 reaching ${ϵ}_{\text{min}}^2=0.0449$ and E_100_03 with ${ϵ}_{\text{min}}^2=0.0482$ (solar abundances, Z = 0.014), and E_150_01_Z028 with ${ϵ}_{\text{min}}^2=0.0119$ and E_150_03_Z028, with ${ϵ}_{\text{min}}^2=0.0132$ (Z = 0.028). The last two, not only reached minimum values lower than the accepted one but also each quantity is simultaneously within its observational uncertainty.

Then, we consider that the best progenitor for V404 Cyg is a system formed by a BH of 9 M_⊙ together with a normal star of 1.5 M_⊙, with a metal content of Z = 0.028. The orbital period of this progenitor was 1.5 d, and the BH accretes between 10 and 30% of the mass lost by its companion.

This model predicts that the donor star of this system may have a hydrogen thermonuclear flash event leading to a short mass loss episode. The remnant of the evolution for the donor star is predicted to be a low-mass helium white dwarf with a hydrogen-rich envelope of mass M = 0.29 M_⊙ and radii R = 0.02 R_⊙.

Although most of the main characteristics of the V404 Cyg system are accounted for by our models, there is one specifically that is not. This is the BH spin parameter, for which we obtained values that are far below the only observational available one, presented in Walton et al. (2017). It seems natural to consider that this discrepancy is due to the assumption that the BH is initially not rotating. Thus, our results may be interpreted as giving some evidence that in the context of close binary systems, stellar mass BHs may be born with appreciable angular momentum.

Acknowledgments

The authors want to thank our anonymous referee for his/her comments and suggestions that helped us to improve our work. They also thank Florencia Vieyro and Federico García, whose comments were of very good use.

References

All Tables

All Figures

	Fig. 1. Quantity ϵ2 as a function of time for our best models with solar composition: E_100_01 and E_100_03. With a dashed line and a grey area is indicated the value of ϵ2 = 0.0485 and the acceptance region we considered for our models. With a black solid line, the value ϵ2 = 0 represents the situation where all the parameters modeled are equal to the ones observed simultaneously.
In the text

Fig. 2. Donor star mass loss rate for each of the best models with solar composition: E_100_01 and E_100_03. As the only variation between these models is on the β value, and this parameter does not affect strongly the mass loss episode, their plots are mostly overlapping. The black dashed horizontal lines represent the estimation for the mass loss rate for V404 Cyg of Ṁd = 1.1 × 10−9 M⊙ yr−1 is denoted with a dashed horizontal line. As for the dashed vertical lines, they represent the times of the minimum value of the epsilon squared function.

In the text

Fig. 3. BH accretion rate for the models with solar abundances of the E group with $P_{\text{orb}}^0$ = 1 d. The qualitative form of the functions is related to the mass loss rate of the donor star by the relation Ṁ2 = −βṀ1. This is, for a higher value for β larger would be the accretion rate onto the compact object. The shaded area represents the zone below the likely upper limit for the mass accretion rate (⟨Ṁ2⟩ = 4 × 10−10 M⊙ yr−1 given by ZZ18. Models with β ≤ 0.3 agree with this upper limit and the results reached by the mentioned authors.

In the text

	Fig. 4. BH spin parameter for the V404 Cyg compact object for each of the best models, according to Eq. (7). As the BH accretes matter it spins up, provoking the increase of $a^{*}=c{J}_2/G{M}_2^2$. Thus, the lower the β, the lower the BH spin. The vertical dashed lines are the predicted ages for V404 Cyg. The values of the BH spin for these times are given in Table 3.
In the text

	Fig. 5. Orbital period as a function of time for our best models with solar metallicity (solid lines), observational estimation for the present orbital period Porb = 6.47 d (dashed horizontal line) and age of the system for the system predicted for each model (dashed vertical lines) shown at the top. The bottom panel shows the time derivative for the orbital period for the same models and ages.
In the text

Fig. 6. Quantity ϵ2 as a function of time for our best models with the change on metallicity: E_150_01_Z028 and E_150_03_Z028. The elements represented are the same as in Fig. 1. With a dashed line and a grey area is indicated the value of ϵ2 = 0.0485 and the acceptance region we considered for our models. With a black solid line, the value ϵ2 = 0 represents the situation where all the parameters modeled are equal to the ones observed simultaneously.

In the text

	Fig. 7. Donor star mass loss rate for each of the best models with the metallicity change: E_150_01_Z028 and E_150_03_Z028. The elements represented are the same as in Fig. 2.
In the text

Fig. 8. BH accretion rate for the two best models with double solar metallicity. The elements represented are the same as in Fig. 3. The qualitative form of the functions is related to the mass loss rate of the donor star by the relation Ṁ2 = −βṀ1. This is, for a higher value for β larger would be the accretion rate onto the compact object. The shaded area represents the zone below the likely upper limit for the mass accretion rate (⟨Ṁ2⟩ = 4 × 10−10) given by ZZ18.

In the text

	Fig. 9. BH spin parameter evolution for the V404 Cyg compact object according to Eq. (7). The values of the BH spin for the predicted ages for V404 Cyg (vertical dashed lines) are very similar to the two best models computed with solar composition in Table 3.
In the text

	Fig. 10. Orbital period and its time derivative as a function of time for models with a metallicity of Z = 0.028. The represented elements are the same as the ones shown in Fig. 5 but adequate to these models.
In the text

Fig. 11. Evolutionary track in an HR diagram for two of our models. Top: model E_100_01 for the case of solar metallicity. Bottom: model E_150_01_Z028 for the case of higher metallicity. The cases of the best models with higher β value are practically overlapped with the models shown. On the left panels, the whole track is shown with a shaded area when the mass loss episode occurs. There can be seen the observationally estimated position of the donor star and the position for the age of the system predicted by our models (with circles on the evolutionary track). Once the star is getting near to becoming a white dwarf it suffers a thermonuclear flash, making it expand and lose mass again for a very short period of time. The end of these tracks are very similar, the donor becomes a low-mass white dwarf with a hydrogen-rich surface. A dashed line denoting R = 0.02 R⊙ was included; it approximately corresponds to the asymptotic radius for the white dwarf. In the right panels, the zoom is on the nearness of this estimation, where the shaded area considers the observational error for the effective temperature and luminosity.

In the text

Fig. 12. Epsilon squared function decomposed on the contributions of every parameter analyzed for two of our models. Left: model with Z = 0.014. Right: model with Z = 0.028. The models computed with the same initial parameters except for β have similar behavior. The shaded area represents the acceptance zone, where the function takes values lower than 0.0485 and the black horizontal line at ϵ2 = 0 represents where the modeled parameter equals the observed estimation for it.

In the text