© The Authors 2025

1 Introduction

In the latest Decadal survey, one of the key themes is entitled “Worlds and suns in context” (National Academies of Sciences 2021). In particular, this title indicates that the next step for studies is to advance a joint description of the properties and evolution of stars and planets. One of the elements of this broad context is the problem of star-planet interactions.

There are many different possibilities when a star and a planet intensively interact with each other via such forces as tides, the magnetic field, or stellar wind (see Lanza 2015). However, all interactions require a planet to be located close to the star. Such a situation is not expected in naive scenarios of planet formation, as conditions at distances close to the star are not suitable for planet growth. Thus, an additional ingredient of a formation mechanism might be involved, namely migration.

Planetary migration is an actively studied topic (see, e.g., a brief review and references in Papaloizou 2021 and the recent review by Paardekooper et al. 2023). Migration is usually studied in protoplanetary disks (Armitage 2010). This involves many interlinked physical processes and results in numerous possibilities of a planetary orbit evolution. Planets can migrate both toward or away from a star in different parts of a disk and/or at different times. Additionally, trapping regions can appear. The gravitational interaction of a disk and a planet might be even more complicated due to planet-planet interactions and interactions with planetesimals (see, e.g., Raymond & Morbidelli 2022). Most importantly, migration is the main mechanism for the formation of “hot” and “warm” planets, especially gas and ice giants (hot jupiters, warm neptunes, etc.). Thus, migration can ultimately allow for intensive star-planet interactions, especially tidal interactions (see, e.g., Lazovik 2021, and references therein).

It is interesting and important to find and explore additional possibilities for planetary migration due to an interaction with a disk. In this paper, we advance the analysis initiated by Kulikova et al. (2019). The main idea is rather simple: In a binary system, the formation of an accretion disk around the secondary main sequence star is possible during the late evolutionary stages of its primary companion. The disk is fed by the stellar wind from the primary. If a planetary system preexists around the secondary component, then a new episode of migration due to a planet-disk interaction might appear. In this work, we study the circumstellar disks. However, the circumbinary disks at late stages of star evolution are also considered in the literature (Kluska et al. 2022), so a similar study is possible for planets in circumbinary disks.

About one-half of intermediate-mass stars are members of binary systems (Duchêne & Kraus 2013; Offner et al. 2023). It is well established that planets in binaries, including close systems with a semi-major axis a ≲ 100 AU, are quite frequent (Bonavita & Desidera 2020; also see the catalogs of planets in binary and multiple stellar systems in Schwarz et al. 2016 and Thebault & Haghighipour 2015)¹. The majority of systems are too wide to allow for the formation of a significant accretion disk, which can change the orbital distribution of planets. Nevertheless, this may not be the case in more than 10% of binaries in the appropriate mass range. In addition, in Hillen et al. (2016), the authors find that in systems with detected circumbinary disks, the secondary component shows the presence of a compact circumstellar accretion disk.

In Kulikova et al. (2019), the authors consider the stationary disk only. However, due to the evolution of the primary, the rate of accretion is variable, especially during periods of enhanced mass loss. In this paper, we present a model of a non-stationary disk around a solar-like main sequence star in a pair with an evolved and more massive companion (the donor, hereafter). The wind from the donor, which is enhanced during the evolution at the red giant and asymptotic giant branches (RGB and AGB, respectively), leads to the formation of a disk. The main part of this study is performing detailed calculations of the structure of such disks while accurately tracking the donor evolution. Our final results represent the migration of a single planet of a different mass in such a disk. Overall, we obtained that planets can experience significant migration in systems with a ≲ 100 AU, which leads to the appearance of “hot” and “warm” massive planets, which may experience intensive tidal interaction with the host star.

In the next section, we describe the scenario used in this study. It includes the binary evolution model, the two accretion disk models, and finally the model of planetary migration. In Sect. 3, the details about our numerical approach are given. Next, in Sect. 4, we present our results. These results and caveats regarding our modeling are discussed in Sect. 5. We draw our final conclusions in Sect. 6. We show additional figures in Appendix A. Details of our calculations are presented in Appendices B, C, D, E, and F.

2 Model

We considered binary systems with a major semi-axis a ≤ 100 AU. The primaries have initial masses, M₁, below 8 M_⊙, which guarantees the evolution of the star with smooth envelope loss at late stages (i.e., without a supernova explosion) and the formation of a white dwarf (WD). The secondary component in each system under consideration is formed as a Solar-like star: M₂ = 1 M_⊙ and evolves slower than the primary (M₂ < M₁). We neglect any possible interactions in a multi-planetary system, that is, only planet-disk interactions are analyzed. When the primary star evolves into a giant, the increasing stellar wind is captured by the secondary, and an accretion disk can be formed around it. The planet interacts with this disk which leads to a significant modification of the planetary orbit.

The disk is assumed to be coplanar with the orbital plane of the binary and with the equatorial plane of the secondary star. The eccentricity of the binary system is assumed to be zero, e = 0.

2.1 Binary evolution

The evolution of a binary system is mainly determined by the evolution of the primary which is calculated using the code MESA (Paxton et al. 2011)². Calculations are done for zero age main sequence (ZAMS) masses listed in Table 1. The formation and structure of an accretion disk around the secondary are mainly determined by the stellar wind of the primary. The profile of the mass loss rate by the primary, Ṁ_w, at late stages of its evolution is shown in Figs. 1, 2. It can be seen that the mass loss can be quite variable, such as the timescale of wind variations, $t_{\text{w}}\equiv {\dot{M}}_{\text{w}}/{\ddot{M}}_{\text{w}}$, can decrease down to 10³ years, while the characteristic viscous timescale for a disk with the size of several AU or more is typically higher than 10⁵ years. That is why a non- stationary model of the disk is an essential element of this study. It is also important that the primary component can lose up to more than half of its mass due to the stellar wind. Correspondingly, it is necessary to take into account the secular decrease of M₁, as this changes the orbital parameters: $$M_1(t)=M_1(0)-\underset{0}{\overset{t}{{\displaystyle \int }}}{\dot{M}}_{\text{w}}\left(t^{\prime }\right)\text{d}{t}^{\prime }.$$(1)

The mass of a secondary component is considered to be constant. For initial binary separations studied here a ∈ [10,100] AU, the evolving donor does not fill its Roche lobe. This allowed us to apply a simplified treatment using the Bondi– Hoyle accretion regime (Bondi & Hoyle 1944). The Bondi-Hoyle accretion rate is defined as $${\dot{M}}_{\text{acc}}={\dot{M}}_{\text{w}}\frac{r_{\text{a}}^2}{4a^2}.$$(2)

Here, r_a is the Bondi radius (aka the gravitational capture radius), which is $$r_{\text{a}}=\frac{2G{M}_2}{v_{rel}^2+c_{\text{w}}^2},$$(3)

where c_w is the sound speed in the stellar wind, $v_{\text{rel}}=\sqrt{v_{\text{w}}^2+v_{\text{s}}^2}$ is the relative velocity of the wind and the companion star, v_w is the wind velocity at the secondary star location, $v_{\text{s}}=\sqrt{G{M}_1^2/\left(M_1+M_2\right)/a}$ is the orbital velocity of the secondary, and G is the gravitational constant. In this work, we choose c_w = 10 km s⁻¹ and v_w = 20 km s⁻¹ (Habing 1996).

We make a rather simplified assumption about the wind velocity. However, this is justified by the fact that the formation of a sufficiently large disk as well as significant migration of planets are possible only at the RGB and AGB stages. These stages are characterized by the low velocity of the stellar wind (see, e.g., Soker & Rappaport 2000; Bennett 2010, and references therein). We discuss the impact of the wind velocity on the disk evolution in Sect. 5.1.

Binary separation changes with time due to the primary mass decrease and angular momentum loss by the system. Following the standard approach (e.g., see Postnov & Yungelson 2014), we obtained an equation for the binary separation a = a(t) evolution: $$\frac{\dot{a}}{a}=\frac{{\dot{M}}_{\text{w}}}{M_1}\left(1+\frac{r_{\text{a}}^2}{4a^2}-2\frac{r_{\text{a}}^2}{4a^2}\frac{1}{q}-\left(1-\frac{r_{\text{a}}^2}{4a^2}\right)\frac{q}{1+q}\right),$$(4)

where q ≡ M₂/M₁ is the stellar mass ratio, and r_a is defined from Eq. (3). Ṁ_w is defined positive. We note that M₁ (t), Ṁ_w (t), r_a (t), and q(t) are all functions of time and that Eq. (4) assumes that there is no stage with a common envelope (see, e.g., Iben & Livio 1993) in the system.See the initial and final binary separation values in Table 1. We note that we obtain Eq. (4) under the assumption of zero eccentricity of the system, e = 0. However, this may not be the case in the observed binaries (see, e.g., Wu et al. 2024). We relegate the study of non-stationary wind-fed disk in eccentric binaries for future work.

Fig. 1 Stellar wind rates for low mass donors with initial masses M1 = 1.7 M⊙ (left) and 3.0 M⊙ (right), calculated using MESA. The first peak of the stellar wind corresponds to the red giant branch (RGB) stage. The second peak associated with the envelope loss by the star corresponds to the stage of the AGB. The zero moment of time t = tZAMS = 0 corresponds to the ZAMS.

Fig. 2 Stellar wind rates for high mass donors with initial masses M1 = 5.0 M⊙ (left) and 7.0 M⊙ (right), calculated using MESA. The peaks of the stellar wind have the same nature as in Fig. 1. The zero moment of time t = tZAMS = 0 corresponds to ZAMS.

2.2 Accretion disk

An axisymmetric circumstellar disk is formed inside the Roche lobe of the planet-hosting main sequence secondary (see de Val-Borro et al. 2009). The outer boundary of the disk is limited by the tidal truncation radius. A stable disk cannot expand beyond this radius due to the tidal torque influence of the primary. The radius is estimated using Table 1 from Paczynski (1977) as $$r_{\text{out}}=r_{\text{tidal\hspace{0.17em}}}\approx 0.896{\left(\frac{1}{1+q}\right)}^{0.0467}{r}_1.$$(5)

Here, r₁ is the Roche lobe radius (Eggleton 1983), $$\frac{r_1}{a}=\frac{0.49q^{2/3}}{0.6q^{2/3}+\ln \left(1+q^{1/3}\right)}.$$(6)

We choose the value of the inner disk boundary as r_in = R₂. Actually, it can be defined, for example, by truncation by the magnetic field and/or wind of the secondary (see, e.g., Armitage 2017). However, the precise value of r_in has little effect on the overall evolution of the disk.

The disk is assumed to be geometrically thin. Since the disk mass is expected to be comparatively small, we neglect its selfgravity. The rotation is taken to be Keplerian: Ω = Ω_K, where ${\mathrm{\Omega }}_{\text{K}}\equiv \sqrt{G{M}_2/r^3}$ with r being the distance from the host star.

The viscous torque, F = F(r, t), is selected as an independent function describing the radial evolution of the disk (see, e.g., Lipunova & Shakura 2000 or Lipunova 2015). It is defined as $$F=-2\pi r\nu \mathrm{\Sigma }r\frac{\text{\hspace{0.17em}d}\mathrm{\Omega }}{\text{\hspace{0.17em}d}r}r=3\pi r^2\mathrm{\Omega }\nu \mathrm{\Sigma },$$(7)

where ν is a kinematic viscosity in the disk, and Σ is the disk surface density. The second equality is valid for the Keplerian rotation.

We introduce the specific angular momentum as $$l=r^2\mathrm{\Omega }=\sqrt{G{M}_{2}r}.$$(8)

Here on the right-hand side, it is assumed that the disk is Keplerian. This quantity simplifies the treatment of equations. We define the values l = l_in and l = l_out which correspond to r = r_in, r = r_out, respectively. Both of these values are used below.

With the notations at hand, the disk evolution is determined by the nonlinear diffusion equation $$\frac{\partial \mathrm{\Sigma }}{\partial t}=\frac{{\left(G{M}_2\right)}^2}{4\pi l^3}\frac{{\partial }^{2}F}{\partial l^2}+{\left(\frac{\partial \mathrm{\Sigma }}{\partial t}\right)}_{\text{ext}}$$(9)

obtained by Lyubarskij & Shakura (1987) with an additional source function from Perets & Kenyon (2013).

The initial condition for Eq. (9) is limited to the sufficiently small absolute value so that the initial spurious amount of matter does not affect further evolution (see Sect. 3). Boundary conditions necessary to solve Eq. (9) are specified as follows (e.g., Armitage 2017): $${F|}_{l_{\text{in}}}=0,{\frac{\partial F}{\partial l}|}_{l_{\text{out\hspace{0.17em}}}}=0,$$(10)

where the inner boundary condition corresponds to the zero viscous stress at r_in, while the outer one corresponds to the zero mass outflow rate at the outer boundary, r_out. This corresponds to the situation when the angular momentum is removed from the outer edge of the disk by the tidal action of the primary.

A different situation takes place when the disk does not expand up to r_out due to low optical thickness. In this case, the disk is optically thick in its own plane within r = r_τ < r_out. The edge r = r_τ can be estimated using the optical thickness in the radial direction introduced by Eq. (21) below (also see the introduction to Appendix D.2). Since the standard disk theory employed in this work is based on the local energy balance with respect to the whole disk plane, it becomes hardly reliable at r > r_τ. Therefore, instead of solving Eq. (9) beyond r_τ we assume that the ambient optically thin flow decrets due to the angular momentum coming from the disk with the velocity $${{\nu }_r|}_{l_{\tau }<l<l_{\text{out\hspace{0.17em}}}}=\frac{3}{2}\frac{v}{r},$$(11)

where v_r is the radial velocity of matter in a disk and l_τ is evaluated at r_τ. Equation (11) restated in terms of F reads $${\left(\frac{F}{l}+\frac{\partial F}{\partial l}\right)|}_{l_{\tau }<l<l_{\text{out\hspace{0.17em}}}}=0.$$(12)

Numerical tests show that such a boundary condition imposed at the stage of disk depletion is physically plausible since it is equivalent to a weakly evaporating tenuous flow outside r_τ. Hereafter, we refer to Eq. (12) as the floating boundary condition.

The source of matter from the stellar wind is calculated similarly to Perets & Kenyon (2013) as $${\left(\frac{\partial \mathrm{\Sigma }}{\partial t}\right)}_{\text{ext}}=\frac{{\dot{M}}_{\text{acc}}}{2\pi r{r}_{\text{a}}}=\frac{{\dot{M}}_{\text{w}}{r}_{\text{a}}}{8\pi r{a}^2}\theta \left(r_{\text{a}}-r\right),$$(13)

where θ (r_a − r) is the Heaviside function.

2.2.1 Quasi-stationary limit

One can obtain the stationary solution of Eq. (9) assuming ∂Σ/∂t = 0 as $$\nu \mathrm{\Sigma }=\frac{{\dot{M}}_{\text{acc}}(t)}{3\pi }f\left(l,l_{\text{a}}(t)\right).$$(14)

Equation (14) has been obtained by Kulikova et al. (2019) using the zero torque condition at the inner boundary along with the explicit form of the source of matter (see Eq. (13)). They showed that the additional factor f can be written as³ $$f\left(l,l_{\text{a}}\right)=\{\begin{array}{ll}1-\frac{l_{\text{in}}}{l}-\frac{l^2}{3l_{\text{a}}^2}+\frac{l_{\text{in}}^3}{3l_{\text{a}}^{2}l},\hfill & l\le l_{\text{a}},\hfill \\ \frac{2l_{\text{a}}}{3l}-\frac{l_{\text{in}}}{l}+\frac{l_{\text{in}}^3}{3l_{\text{a}}^{2}l},\hfill & l>l_{\text{a}},\hfill \end{array}$$(15)

where l_in and l_a are values of the specific angular momentum at r_in and r_a, respectively, Eq. (8). Here, we additionally assume that the local accretion rate through the disk vanishes beyond r_a.

2.2.2 Radial velocity of gas in a disk

Viscous torque, Eq. (7), is related to the local accretion rate in a disk according to the azimuthal projection of the Euler equation for the accreting matter, $$\dot{M}(l)=\frac{\text{d}F}{\text{\hspace{0.17em}d}l},$$(16)

where by definition $$\dot{M}(l)\equiv -2\pi r(l)\mathrm{\Sigma }(F,l)v_{\text{r}}>0.$$(17)

Using Eq. (7), we obtained $$v_{\text{r}}=-\frac{1}{2\pi r(l)\mathrm{\Sigma }(F,l)}\dot{M}(l)=-v\frac{3}{F}\frac{\text{\hspace{0.17em}d}F}{\text{\hspace{0.17em}d}r}.$$(18)

Equation (18) is used below to get the rate of the type II planetary migration in Sect. 2.4 and to consider the disk structure in Sect. 4.1.

2.3 Energy balance

Equation (9) must be solved taking into account details of internal heating and cooling of the disk. Geometrically thin disks are known to have a local energy balance: thermal energy extracted from the mean shear motion of matter, for example, via turbulent cascade, at some r is emitted away at the same distance from the host star. A number of further simplifying assumptions on the turbulent dissipation of energy as well as on the transfer of thermal energy across the disk leads us to the energy balance equation, which determines kinematic viscosity standing in Eq. (7). Thus, these assumptions allowed us to establish the missing relationship between F and Σ, which is necessary to solve Eq. (9).

Following Shakura (1972) and Shakura & Sunyaev (1973) we assume that the turbulent viscosity is expressed through the dimensionless parameter ɑ, $$\nu =\alpha c_{\text{s}}H=\alpha c_{\text{s}}^2{\mathrm{\Omega }}^{-1}=\alpha \frac{\gamma k_{\text{B}}}{\mu m_{\text{H}}}{T}_{\text{c}}\frac{r^{3/2}}{\sqrt{G{M}_2}},$$(19)

where H is the disk vertical scale height and H ≈ c_s/Ω due to the vertical hydrostatic equilibrium. The sound speed, c_s, is related to the temperature according to the ideal gas equation of state, $c_{\text{s}}^2=\frac{\gamma k_{\text{B}}}{\mu m_{\text{H}}}{T}_{\text{c}}$, where γ is a specific heat ratio, k_B is the Boltzmann constant, µ is the mean molecular weight in the disk matter, m_H is a hydrogen atom mass, and T_c is the disk temperature at its midplane. In our simulations we set γ = 1.4, µ = 2.3, which corresponds to a standard hydrogen-helium composition of the gas.

We set α = 0.01 for all our calculations unless otherwise stated. This is a conservative estimate of dimensionless viscosity based on spectroscopic and interferometric observations of protoplanetary disks combined with each other within the alphadisk paradigm, see Hartmann et al. (1998). The more recent observations indicate that α varies from disk to disk or/and may be different in the inner and the outer parts of the disk approximately in the range 10⁻⁴−10⁻² (see, e.g., Rosotti et al. 2020, Trapman et al. 2020, Flaherty et al. 2024, and the review by Rosotti 2023). For this reason, we also perform additional calculations taking the lower α = 0.001 discussed in Sect. 5.4.

While considering the radiative energy transfer in a disk we do not restrict the model with the optically thick case only. We estimated the Rosseland optical thickness in the ɀ-direction as $${\tau }_{\text{z},\text{R}}(r,t)=\frac{{\kappa }_{\text{R}}(r,t)\mathrm{\Sigma }(r,t)}{2},$$(20)

where κ_R(r, t) is the Rosseland mean opacity taken from Semenov et al. (2003)⁴. Everywhere below it is assumed that κ_R(r, t) is taken at the disk midplane. Following Nakamoto & Nakagawa (1994) we also consider disks which are optically thin in the ɀ-direction (i.e., τ_z,R < 1).

Additionally, we introduced an estimate of optical thickness in the radial direction, $${\tau }_{\text{r},\text{R}}(r,t)=\underset{r}{\overset{r_{\text{out}}}{{\displaystyle \int }}}\frac{{\tau }_{\text{z},\text{R}}(r,t)}{H(r,t)}\text{d}r.$$(21)

We note that the integration limits chosen in the definition (21) yield the lower estimate of optical thickness in the radial direction at the given r as far as we assume the normal behavior of surface density decreasing outward.

The local energy balance is valid in regions with τ_r,R > 1 only. We assume that the condition τ_r,R > 1 provides a reasonable restriction for the validity of our model to low-mass disks that are optically thin in the ɀ-direction.

The midplane temperature that determines the viscosity of a disk is obtained from the energy balance equation (for the details of its derivation, see Appendix B), $$\begin{array}{ll}{T}_{\text{c}}^4=\hfill & \frac{9\nu \left(T_{\text{c}}\right){\text{ΣΩ}}_{\text{K}}^2}{8{\sigma }_{\text{B}}}\left(1+\frac{3{\kappa }_{\text{R}}\left(T_{\text{c}}\right)\text{Σ}}{8}+\frac{1}{2{\kappa }_{\text{p}}\left(T_{\text{c}}\right)\text{Σ}}\right)\hfill \\ \hfill & +\left(T_{\text{I},\text{A}}^4+T_{\text{I},\text{D}}^4+T_{\text{w}}^4\right)\left(1+\frac{1}{2{\kappa }_{\text{p}}\left(T_{\text{c}}\right)\text{Σ}}\right),\hfill \end{array}$$(22)

where σ_B is the Stefan-Boltzmann constant, κ_P is the Planck mean opacity essential for optically thin regions where τ_z,R < 1 (but τ_r,R > 1, see Appendix C) obtained by Semenov et al. (2003). We note that Eq. (22) contains terms due to irradiation flux from the host star, $F_{\text{I},\text{\hspace{0.17em}A}}~{\sigma }_{\text{B}}{T}_{\text{I},\text{\hspace{0.17em}A}}^4$, the irradiation flux from the donor, $F_{\text{I},\text{D}}~{\sigma }_{\text{B}}{T}_{\text{I},\text{D}}^4$ (see Appendix D for the details), and additional heating of the disk due to dissipation of the kinetic energy of settling matter, $F_{\text{w}}~{\sigma }_{\text{B}}{T}_{\text{w}}^4$ and the rest “viscous” flux is labeled as F_v.

2.4 Planetary migration

In this work, we consider a planet on a circular orbit with size r_p. The eccentricity of a planetary orbit is assumed to be zero, e = 0. The orbit inclination with respect to the disk is equal to zero as well. The largest dynamically stable circular planetary orbit can be estimated as $$\frac{r_{\text{p}}^{\max }}{a}=0.464-0.380\frac{M_1}{M_1+M_2}$$(23)

(see Holman & Wiegert 1999). Normally, a planet transfers its angular momentum to the disk, so that it migrates inward. The planet continues to migrate until reaching the orbit where tidal interactions with the central star become significant. We estimate this critical distance, $r_{\text{p}}^{\text{fin}}$, as the radius from which a planet can migrate due to tidal forces down to the stellar surface by t_life. For calculations, we use the rate of tidal migration (TM hereafter) from Jackson et al. (2008) defined as $$\frac{r_{\text{p}}^{\text{fin}}\left(t_{\text{life}}\right)}{R_2}={\left(\frac{117}{4Q}\left(\frac{m_{\text{p}}}{M_2}\right)\sqrt{\frac{G{M}_2}{R_2^3}}{t}_{\text{life}}\right)}^{2/13}.$$(24)

Here, Q is the stellar tidal dissipation parameter. Its expected values for solar-like stars are Q ~ 10⁵–10⁶ (we select Q = 10⁶ for the estimation), but higher values of Q up to 10⁸ are also possible. In Eq. (24), m_p is the planetary mass which is assumed to be constant (i.e., we neglect mass accretion by the planets), and t_life is the disk characteristic lifetime, or, equally, total migration time of the planet. A numerical estimate of the distance (24) for the given values of parameters is: $r_{\text{p}}^{\text{fin}}\left(t_{\text{life\hspace{0.17em}}}\right)/R_2\approx 2.6{\left(m_{\text{p}}/300{\text{m}}_{\oplus }\right)}^{2/13}{\left(t_{\text{life\hspace{0.17em}}}/\text{Myr}\right)}^{2/13}$. Obviously, this value is usually close to unity. For example, the largest value in our simulations is achieved for a planet with m_p = 300 m_⊕ in a system with $t_{\text{life\hspace{0.17em}}}~120\text{Myr}:r_{\text{p}}^{\text{fin\hspace{0.17em}}}\left(t_{\text{life\hspace{0.17em}}}\right)\approx 5.4R_2$.

We considered two types of migration: type I migration for less massive planets and type II migration for more massive planets which open a gap in the disk. Type I migration rate is estimated according to Tanaka et al. (2002) as $$\frac{\text{d}{r}_{\text{p}}}{\text{\hspace{0.17em}d}t}=-(2.72+1.08\beta )\frac{m_{\text{p}}}{M_2^2}{h}^{-2}\mathrm{\Sigma }{r}_{\text{p}}^3\mathrm{\Omega }\left(r_{\text{p}}\right),$$(25)

where h = H/r is the aspect ratio of the disk at the planetary orbit, β is the power law index of the surface density profile Σ ∝ r^−β.

The transition from the type I migration to the type II migration is determined by the planet-to-star mass ratio, q ≡ m_p/M₂ (see, e.g., Baruteau et al. 2014). Its critical value is $$q_{\text{crit\hspace{0.17em}}}=\frac{100}{ℛ}{\left({(X+1)}^{1/3}-{(X-1)}^{1/3}\right)}^{-3},$$(26)

where $X=\sqrt{1+\frac{3}{800}{h}^3ℛ}$ is the numerical coefficient and $ℛ=r_{\text{p}}^2{\mathrm{\Omega }}_{\text{p}}/\nu$ is the Reynolds number defined by the viscosity of the disk. If the planet-to-star mass ratio q is less than the critical value, then the type I migration takes place. The type II migration occurs in the other case.

We compared q with the critical mass ratio all the time the planet migrates embedded in the evolving disk. So, the planet can change the type of migration in the course of the disk evolution. The transition from type II to type I happens not instantaneously. We assume that the gap can be closed on a characteristic timescale which we estimate according to Armitage & Rice (2005) as $$\mathrm{\Delta }{t}_{\text{close\hspace{0.17em}}}~\frac{H^2}{\nu }=\frac{1}{\alpha {\mathrm{\Omega }}_{\text{p}}}.$$(27)

For example, if a planet at first evolves according to type II migration and later occurs in circumstances that correspond to q < q_crit, then the transition to the type I happens only after Δt_close.

The type II migration rate is estimated according to $$\frac{\text{d}{r}_{\text{p}}}{\text{\hspace{0.17em}d}t}=v_{\text{r}}\frac{M_{\text{d}}}{M_{\text{d}}+m_{\text{p}}}=\frac{3}{2}\alpha h^2{r}_{\text{p}}{\mathrm{\Omega }}_{\text{p}}\frac{2r_{\text{p}}}{F}\frac{\text{\hspace{0.17em}d}F}{\text{\hspace{0.17em}d}r}\frac{M_{\text{d}}}{M_{\text{d}}+m_{\text{p}}},$$(28)

where $M_{\text{d}}=2\pi r_{\text{p}}^2{\mathrm{\Sigma }}_{\text{p}}$ is the characteristic mass of the disk at the planetary orbit, v_r is the radial velocity of the accreting matter (18), and F is the viscous torque taken from Eq. (7). We note that at a given r the value of v_r can be both positive or negative. Therefore, the outward type II migration is possible. Equation (28) takes into account the decrease of the migration rate in the case when m_p > M_d (see Ivanov et al. 1999; Scardoni et al. 2020).

In order to improve a qualitative presentation of our results, we also define a characteristic migration timescale $${\widehat{t}}_{\text{p}}=\frac{r_{\text{p}}}{\left|\text{\hspace{0.17em}d}{r}_{\text{p}}/\text{d}t\right|}.$$(29)

For the type II migration, the timescale given by Eq. (29) corresponds to the local viscous time in a disk. The shorter the timescale, the faster the migration of the planet, and vice versa.

In this work, we would like to avoid a bulk of complexity associated with the description of migration developed recently, see Paardekooper et al. (2023). Instead, our objective is to carry out a comprehensive study of simple migration in a complex disk in order to obtain a reliable basis for future work. This simplified approach allowed us to study a parameterized problem of migration. We discuss the complexity of migration in Sect. 5.5.

3 Numerical setup of the model

Equations governing the evolution of the accretion disk have been described in the previous Section. These are Eq. (9) which determines Σ = Σ(F, l, t) and Eq. (22) which determines T_c = T_c(Σ, l, t). These equations are related to each other through the constraint given by Eq. (7), which defines F = F(Σ, T_c, l, t). We use Eq. (7) to express Σ via F and T_c. The set of Eqs. (7), (9), (22) allowed us to determine F, Σ, and T_c as functions of the radial coordinate (r and l are interchangeable variables) and time. The scheme employed to solve this nonlinear problem numerically is described in detail in Appendix E.

We additionally introduce the semi-analytical quasi- stationary disk model by replacing Eq. (9) with Eq. (14). The corresponding solution has been obtained in Kulikova et al. (2019). We compare our results with this solution. It is denoted as QSD (quasi-stationary disk), while the full non-stationary model described by Eq. (9) is denoted as NSD (non-stationary disk) hereafter. The non-stationary disk is referred to as the NS disk hereafter.

The numerical grid for disk evolution is defined as follows. We set a uniform spatial grid in terms of the specific angular momentum, $l:\left\{l_n=l_{\text{in\hspace{0.17em}}}+n\mathrm{\Delta },\mathrm{\Delta }=\left(l_{\text{out\hspace{0.17em}}}-l_{\text{in\hspace{0.17em}}}\right)/N,n=\overline{0,N}\right\}$. For all simulations, we set N = 100 spatial cells distributed uniformly in the range from l_in to l_out. We set the non-uniform numerical grid over time $\left\{t_k:t_0<t_1<\cdots <t_K,k=\overline{0,K}\right\}$, where the slices t_k are taken from the MESA track grid with optional additional grid thickening performed using a linear interpolation of MESA parameters. The MESA time grid that we use for post-MS stars concentrates cells to the RGB and AGB peaks and has from roughly 3000 time cells for the most massive donor to more than 11000 time cells for the least massive donor (see Figs. 1 and 2 for the time cells used). We also add K′ = 500 uniformly distributed time cells to have a more dense coverage of evolution stages far from RGB and AGB peaks. The implicit unconditionally stable numerical scheme (see Appendix E.1 for details) allowed us to use an arbitrary ratio of grid steps by time and coordinate.

Besides the mathematical stability of the scheme, we ensure the physical convergence of the numerical solution. Regarding the time resolution, our time grid has a time step of about or less than 1 year close to the RGB and AGB peaks. This is at least two orders of magnitude less than the characteristic wind variability timescale, which has its minima at RGB and AGB peaks. At the same time, the viscous timescale is much larger.

The sharp gradients of the opacity may occur due to evaporation of dust species, which leads to significant gradients of temperature and surface density (or, interchangeably, viscous torque). This implies the limitation of the size of the spatial grid. For N = 100 and typical (varying) outer radius, we get the grid cell size of ~0.05 AU at the inner edge of the disk. The grid cell size increases with radius $\propto \sqrt{r}$. We do not expect sharp gradients of the disk variables at its outer parts, where the size of the grid cell gradually attains values up to one AU or to a few AU depending on the location of the outer boundary which can approach 100 AU. We ensure that the numerical solution does not depend on the number of the spatial cells if N ≳ 100.

The spurious initial value of the disk viscous torque is set to $F_{\text{init\hspace{0.17em}}}={10}^{20}\sin \left(\frac{\pi }{2}\frac{l-l_{\text{in\hspace{0.17em}}}}{l_{\text{out\hspace{0.17em}}}-l_{\text{in\hspace{0.17em}}}}\right)$ dyn cm, which satisfies the boundary conditions. However, we note that the exact choice of the initial condition is of small significance for the following solution, provided that the absolute value of this condition is sufficiently small. The disk loses its initial state over a viscous timescale, which is always much smaller than the time between ZAMS and RGB or AGB. The stability of the numerical solution has been tested for varying amplitude of the initial condition. It has been varied from 10²⁰ dyn cm to zero. So, the zero initial condition is also acceptable in our scheme.

At the initial slice and for all subsequent slices the surface density Σ is expressed through F and T_c using Eq. (7). The temperature T_c is obtained as the solution of Eq. (22). This is performed rearranging this equation to the form f (T_c) = 0 and looking for the roots by the Brent algorithm (Brent 1973), for which we use the eponymous method BRENTQ from PYTHON library SCIPY.OPTIMIZE.

The boundary condition is imposed according to a numerical form of Eq. (12). We note that the outer boundary of a disk defined in Eq. (5) depends on time, l_out = l_out(a) = l_out(t_k). Thus, the spatial step Δ = Δ(t_k) = Δ_k and spatial nodes l_n = l_n (t_k) = l_k,n also depend on time. They must be updated at every slice accounting for the binary orbital evolution (see Eq. (4) which is also integrated numerically). For the numerical simulations, we select the range of initial separations a from 10AU to 100 AU with the step of 10AU for all donor masses (see Table 1).

After the calculation of the accretion disk evolution, we calculate the planet migration. The details of the scheme are described in Appendix E.2. We interpolate the numerical variables describing the accretion disk evolution to a better spatially and time-resolved grid. The spatial resolution of the new grid becomes less than 0.02 AU. The time resolution of the new grid is about 1 year during the rapidly varying wind epochs (RGB and AGB) and not worse than 100 years during the whole evolution. Using this grid, we simulate planet migration according to equations described in Sect. 2.4. We apply the described scheme to a range of planet masses from 1 m_⊕ to 3000 m_⊕ and a range of initial planet orbit separations from 0.25 AU to the largest dynamically stable circular planetary orbit, defined by Eq. (23). For the visual representation of results over the wide range of free parameters, we use bilinear interpolation to plot the contours of migration as described in Appendix E.3.

4 Results

4.1 Accretion disk

For further analysis, we define several time moments that are the same for all disks with the specified donor:

• t_RGB is the moment of the maximum wind rate at the red giant stage (see the first peak in Figs. 1, 2);

• t_AGB is the moment of the maximum wind rate at the stage of the asymptotic giant branch (see the second (main) peak in Figs. 1, 2);

• t_env is the moment of the envelope loss by the donor when the disk is not supplied with matter anymore.

Additionally, we introduce the time moments that specify the disk in the system with particular initial separation:

• t_frm is the moment of disk formation when the accretion rate becomes high enough to produce a disk with considerable size. We formally take t_frm corresponding to the disk with size 2 AU when the accretion rate for the first time attains ${\dot{M}}_{\text{acc}}\simeq {\dot{M}}_{\text{acc}}^{\lim }$ introduced below in Sect. 4.1.3.

• t_dpl is the moment after the envelope loss by the donor, which is in the midst of disk depletion during its closed evolution;

• t_dec is the moment of the disk decay at the end of its closed evolution when r_τ → r_in.

The time moments t_frm, t_AGB, and t_dpl are used to show the disk profiles in Fig. 3.

The universal sequence of the time moments is t_ZAMS < $t_{\text{frm}}^{*}$ < t_RGB < $t_{\text{frm}}^{*}$ < t_AGB < t_env < t_dpl < t_dec, where the asterisk at t_frm means, that this time moment can be either before t_RGB or after t_RGB depending on the initial system separation and initial donor mass (see also Table 2, where these time moments are presented for different donors in a system with initial binary separation a = 30 AU).

Finally, we select several spatial regions in the disk. The regions are related to the optical thickness determined by the wind rate or, alternatively, the accretion rate. The optical thickness substantially varies across the disk.

• Region 1 (R1) is the standard optically thick region where τ_z,R > 1 (i.e., τ_r,R >> 1). The disk contains an extended R1 at high accretion rates. Typically, this corresponds to Ṁ_acc ≳ 10⁻⁸ M_⊙ yr⁻¹. This region is the most important for planetary migration because it is characterized by the largest disk mass up to ~0.01 M_⊙, the largest surface density up to ~ 10⁴ g cm⁻² and the largest viscous torques. Below, in Sections 4.1.4 and 4.1.5, the disk with R1 extending up to a considerable value which is not less than 1 AU is represented by the curves taken at t = t_AGB for all donors and by the curves at t_dpl only for donors M₁ = 7.0,5.0 M_⊙ (see Fig. 3).

• Region 2 (R2) is transitional region where τ_z,R ≲ 1 but τ_r,R > 1. Usually, a disk with an outer R2 having a size ∼2 AU forms in epochs when Ṁ_acc ≳ 10⁻¹¹ M_⊙ yr⁻¹ (see the derivation of the corresponding estimate, Eq. (33), in Section 4.1.3)⁵. In this case, a limited R1 appears in the inner part of the disk. As far as Ṁ_acc ≲ 10⁻¹¹ M_⊙ yr⁻¹, disk becomes fully optically thin in the vertical direction. Planet migration in R2 becomes slower in general. Nevertheless, it can still be non-negligible for high-mass planets experiencing type I migration in comparison with the migration of low-mass planets in R1. In Sections 4.1.4 and 4.1.5, the disk with R2 extending enough is represented by the curves at t_frm for all donors and the curves at t_dpl only for donors M₁ = 3.0,1.7 M_⊙ (see Fig. 3).

• Region 3 (R3) is the optically thin region. In this region, the disk becomes optically thin in the radial direction: τ_r,R < 1. An assumption of the local energy balance is far from being valid and thus, our model predictions must be inaccurate in this region. Therefore, we use a floating bound in our calculations (see Eq. (12)). It separates R2 from R3. Also, we apply the condition of decretion of R3 (see Eq. (12)). The low value of optical thickness in the radial direction specific to R3 prevails at low accretion rates: disk usually becomes fully optically thin for Ṁ_acc ≪ 10^–11 M_⊙ yr^–1. The disk has extremely low mass if R3 is dominating (see Sections 4.1.4 and 4.1.5 for details). The corresponding planet migration is also negligible (for details see Sect. 4.2.2). We note that these results were also checked within the model with the outer boundary imposed at the tidal truncation radius only.

The evolution of an accretion disk starts from the formation of R2 close to the host star. As the increasing stellar wind rate approaches the value corresponding to the accretion rate Ṁ_acc ≪ 10^–11 M_⊙ yr^–1, R2 expands into the surrounding R3 where the accreted matter stays to be optically thin. The further increase of the wind rate leads to the formation of R1 close to the host star. Subsequently, R1 and R2 both spread outward. Once the wind rate starts to decrease, the reverse process takes place. In the case when the wind ceases after the envelope loss by the donor, the reverse process continues until the disk decays. By the disk decay we mean the moment when both R1 and R2 vanish. In the systems with low mass donors, the wind rate may decrease so much that the disk decays temporarily between t_RGB and t_AGB.

In Sect. 4.1.4, we present a detailed description of disk evolution in the system with the most massive donor M₁ = 7.0 M_⊙. After that, in Sect. 4.1.5, we highlight the variations of disk evolution in the case of less massive donors, M₁ = 5.0, 3.0, and 1.7 M_⊙. We start by demonstrating the importance of the non- stationary approach and then describe the heating of the disk by the donor.

We also remind that the abbreviations NSD and QSD mean non-stationary and quasi-stationary disk models, respectively. When referring to a physical non-stationary disk, we use the notation “NS disk”.

Fig. 3 Radial profiles of the accretion disk variables. Each column corresponds to a different initial donor mass, M1, of our sample. The initial binary separation is a = 30 AU for all columns. The rows from 1 to 6 show the disk profiles of surface density Σ, midplane temperature Tc, radial velocity vr, aspect ratio h = H/r, optical thickness in vertical and radial directions τz,R, τr,R, respectively. Blue, red, and green curves represent a disk at the time t = tfrm, tAGB, tdpl, respectively (for the values of the time moments see Table 2). Solid and dashed curves show, respectively, NSD and QSD. The markers of the solid curves correspond to the spatial grid of the solution. QSD curves are presented only on plots with Σ, τz,R and τr,R, otherwise, for the moment tAGB only (see Table 4 for the values of the Bondi radius and the snow line radius).

4.1.1 Relevance of the non-stationary evolution model

The results of Kulikova et al. (2019) are based on QSD. That model is satisfactory, while the characteristic timescale of the wind rate variation, t_w, is much longer than the timescale of the disk non-stationary evolution. The latter scale is usually approximated by the viscous timescale of a disk, $$t_v=1/\left(\alpha h^2\mathrm{\Omega }\right)\ll t_{\text{w}},$$(30)

where Ω and h are taken at the outer boundary of the disk. Employing NSD, we compare t_ν with t_w in the systems with the lowest mass donor and the most massive one for several initial separations (see Fig. 4).

This comparison suggests that QSD should describe well the evolution of a disk supplied from the wind in epochs far from RGB and AGB peaks. Moreover, t_ν decreases in wider systems since the smaller accretion rate leads to the formation of smaller disks. This picture is complicated by the truncation of the disk at the tidal radius in the closest systems (see the curve for a = 10 AU in Fig. 4). As one can see, the residual disk truncated at the tidal radius lives longer than the disk with the floating outer boundary represented by the curves for a = 30, 100 AU, although the latter is initially larger and thus initially has a larger t_ν. Such a peculiarity reflects the underestimation of the timescale of the truncated disk non-stationary evolution in comparison with its viscous timescale. In general, we find that an external long-term source of matter allows the disk to exist for a much longer time than t_ν, which stays within ≲10⁵ years, except the time around AGB peak, where it attains ∼10⁶ years. In the latter case, the viscous timescale corresponds to the disk with the size ~10 AU and h ~ 0.05, which is in accordance with the residual disk profiles in Fig. 3 (see below for more details).

However, the condition of quasi-stationarity, Eq. (30), is not valid during the whole period of the system evolution. It can be seen that t_w becomes comparable or even less than t_ν around t_AGB for all donors and initial separations. This is also the case around t_RGB for initially close systems with high mass donors. The predictions of QSD are inaccurate at these times. This is demonstrated below in Fig. 3 by the large difference between the disk profiles of QSD and NSD at t = t_AGB. Additionally, the disk profiles of QSD and NSD significantly differ from each other at t = t_frm for the donor with M₁ = 3.0 M_⊙ (see column 3 of Fig. 3), because in this case the sufficiently large disk with the size ~2 AU is accumulated right by the beginning of prominent oscillations of the wind which is typical for AGB phase of the lower mass donors. These wind oscillations can be seen on the bottom panel in Fig. 4.

We find that in all these cases NSD provides an order of magnitude smaller surface density and somewhat smaller aspect ratio of the disk than is predicted by QSD. Looking at the disk profiles at the moment of the increasing accretion rate, it is clear that NSD always lags behind QSD. Thus, disks in NSD are found to be more (less) massive than in QSD at the same moments of time during the decrease (increase) of the wind rate. In other words, the NSD solution tends to the QSD solution on a rather long timescale corresponding to a few viscous timescales of the disk, which has been checked in our test calculations. The lag between NSD and QSD is negligible as long as the condition t_ν ≪ t_w is valid, but the lag grows along with the ratio t_ν/t_w.

We also note that, unlike QSD, NSD enabled us to study the residual disk left after the envelope loss by the donor. NSD provides accurate initial disk profiles that determine the isolated evolution of the residual disk without the matter inflow. We note that these initial data are the outcome of essentially non-stationary dynamics of wind-fed disks during the preceding AGB phase. Also, note that the radial velocity of the residual disk acquires a smooth profile and becomes negative throughout the disk up to its floating boundary. This is in accordance with its quasi-stationary state (see the corresponding row 3 in Fig. 3).

It is important to emphasize the peculiar profile of the radial velocity in the essentially NS disk. NSD demonstrates an extended zone of decretion, v_r > 0, in the outer parts of the disk at t = t_AGB. The profile of the radial velocity significantly deviates from the known self-similar solution first obtained by (Lynden-Bell & Pringle 1974). The zone of decretion starts from ≳2AU (≲2 AU) for the least (most) massive donor. Thus, it goes far inside the Bondi radius, which is within r_a = 3–4 AU for all systems around t = t_AGB (see also Table 4 for the exact values). It is also instructive to compare the zone of decretion at t = t_AGB with that at t = t_dpl when the disk is close to a quasi- stationary state. Despite the latter case, the disk is even smaller than its progenitor at t = t_AGB, and its own zone of decretion starts considerably farther from the host star.

Decretion takes place due to the rapid accumulation of matter in the middle part of the disk. This part is limited by the Bondi radius from the outside and the radius where the local viscous timescale becomes longer than the timescale of wind variation from the inside. The corresponding abnormal distribution of surface density causes the gradient of the viscous stress to change sign. This leads to the inverse flux of an angular momentum through the middle part of the disk. We note that the absolute value of vr in the middle of the zone of decretion can exceed the standard quasi-stationary value inferred from the instant disk temperature using Eq. (11) up to several times when taken at the same location in a disk and up to an order of magnitude when taken at the outer boundary. This is visible, for instance, by the profile of vr shown at t = t_frm for the donor with M₁ = 3.0 M_⊙ (see the range 1.5 ≲ r ≲ 2.0 AU in Fig. 3, column 3). We check the validity of an enhanced decretion in the essentially NS disk performing an additional numerical test with a simple disk and wind models (see Appendix F for details).

For the lower mass donors, M₁ = 1.7 M_⊙ and M₁ = 3.0 M_⊙, the zone of decretion exists during the whole AGB phase. That manifests itself in the form of prominent growing oscillations of the wind rate. For example, it lasts up to ~0.5 Myr in systems with the donor with M₁ = 3.0 M_⊙. The zone of decretion shrinks outward for a while, which occurs right after the sharp wind rate minimum. Just before that, for a short period ~10³–10⁴ yr, the zone may split into two parts. Thus, an additional layer of accretion may emerge for a while in the range 2–4 AU. We check that this feature occurs due to the heating of the outer parts of the disk by the donor. The higher mass donors, M₁ = 5.0 M_⊙ and M₁ = 7.0 M_⊙ have a single peak of the wind rate which occurs right at t = t_AGB, by definition. The corresponding zone of decretion arises already ~1 Myr before t = t_AGB being attached to R2. However, it becomes significant only ~0.1 Myr before t = t_AGB and vanishes shortly after the peak of the wind rate. We note that an additional layer of accretion appears for ~5 kyr shortly after the wind rate attains its maximum. However, this layer is attached to R2, which is adjacent to the floating boundary of the disk. At the same time, R2 widens due to a prominent expansion of the system as the donor intensively loses its mass. Similar to the systems with the lower mass donors, an additional layer of accretion emerges due to the heating of R2 by the donor. As the binary becomes wider, the zone of decretion becomes less extended for all donors. For the binary initial separation as large as 100 AU, it stays attached to R2 adjacent to the floating boundary of the disk.

Fig. 4 Disk viscous and wind variability timescales for two donors and three initial separations over the disk evolution. Black curve: characteristic timescale of wind variability, $t_{\text{w}}={\dot{M}}_{\text{w}}/{\ddot{M}}_{\text{w}}$. Blue, green, and red curves: characteristic viscous timescale of the disk, tν = r2/ν, taken at the outer boundary of NS disk in the system with initial separations a = 10, 30, 100 AU, respectively. The initial donor mass M1 = 7.0 M⊙ and M1 = 1.7 M⊙ on the top and the bottom panels, respectively.

4.1.2 Donor impact on the disk structure

The luminosity of the donor changes significantly in the course of donor evolution after the Main Sequence. The dramatic increase in the donor size leads to essential heating of the outer parts of the disk. This kind of heating requires the following condition, obtained for the conical disk geometry assumption (see Fig. D.1 for a visual representation of the calculation): $$R_1>H_{\tau }\approx 1\text{AU}\left(\frac{r_{\tau }}{10\text{AU}}\right)\left(\frac{h_{\tau }}{0.1}\right),$$(31)

where H_τ is the vertical scale height of the disk at r = r_τ. H_τ approximates the largest vertical scale height in a disk. Let us remind that r_τ is the radial coordinate corresponding to the radial optical thickness τ_r,R = 1, while h_τ = H_τ/r_τ.

The condition (31) suggests that the disk heating by the donor takes place for the donor size above one to several AU depending on the wind rate and the binary separation. Such a value is typical for the considered donors at RGB and AGB stages, t_RGB ≲ t ≤ t_AGB. Our calculations confirm this expectation, (see the typical case in Fig. 5 where a “hump” of T_c near the outer edge of the disk is caused by the donor impact). The temperature hump is seen on the solid curves, which are compared to the dashed curves obtained in the absence of the donor heating. Numerical tests show that the temperature hump exists in a disk for the whole range of the initial orbital separation studied here, a ∈ (10, 100) AU, as well as for all considered initial donor masses. We note that this feature can be found also in Fig. 3.

The rise of disk temperature produced by donor heating is usually several times larger than that produced solely by the host star. It can exceed the temperature of the disk obtained in the absence of the donor heating up to an order of magnitude for the most extreme combination we studied here: the high mass donor with M₁ = 7.0 M_⊙ and the closest binary with initial separation a = 10 AU. At the same time, the donor heating becomes weaker as the initial separation increases and/or we consider the less massive donor. The radial extent of the temperature hump caused by donor heating usually attains several AU in the outer parts of the disk, while the inner parts of the disk within ~1 AU remain unaffected due to the disk self-shadowing. The temperature hump also leads to the corresponding decrease of surface density, since T_c and Σ are approximately related through the condition F ~ νΣ ~ T_cΣ ~ const valid at least in the quasi-stationary regime. As a consequence, an ordinary decrease of surface density toward the disk outskirts is replaced by an approximately flat radial profile Σ(r) ≈ const, or in some cases, even by the inverse radial profile Σ(r) ∝ r^β, β > 0. Such a radial profile of surface density is a distinctive feature of a disk heated by a luminous donor. Moreover, it can be accompanied by the inverse radial profile of temperature (see the curves for a cooler disk at t_frm in Fig. 5). The detailed analytical derivation of donor heating is relegated to Appendix D.

Fig. 5 Donor heating of a disk sometime after the disk formation at t = tfrm = 42.966 Myr (blue curve; see Table 2) and at t = tAGB = 48.111 Myr (red curve). The NSD profiles of surface density (top panel) and midplane temperature (bottom panel) are presented. The solid lines show the results with donor heating taken into account. The dashed lines show the results without accounting for the donor heating, i.e., for $T_{\text{I},\text{D}}^4$ set to zero in Eq. (22).

4.1.3 The condition of disk formation

Our simulations apply to disks that are optically thick in the radial direction. This corresponds to some lower limit on Ṁ_acc. To obtain the condition of the optically thick disk formation, we put rough limits on the optical thickness in both ɀ- and r- directions. We derive an estimate for τ_z,R from Eqs. (20) and (22) assuming that T_c ≈ T_I,A which is true for a low-mass disk, and also assuming that the disk can be described within the QSD approach. Thus, employing Eq. (14), we obtained $${\tau }_{\text{z},\text{R}}\approx 0.03\left(\frac{{\kappa }_0}{3.0{\text{cm}}^2{\text{g}}^{-1}}\right){\left(\frac{\alpha }{0.01}\right)}^{-1}{\left(\frac{r}{2\text{AU}}\right)}^{-1}\left(\frac{{\dot{M}}_{\text{acc}}}{{10}^{-11}{\text{M}}_{\odot }{\text{yr}}^{-1}}\right),$$(32)

where κ₀ is the value of Rosseland opacity for a given temperature: at low temperatures, it varies within the range ~1–5 cm² g^–1. We note that the normalization of distance at r = 2 AU is a matter of convention. We chose this value as it is close to the location of the snow line in the preceding protoplanetary disk. The corresponding part of the protoplanetary disk has been a nursery of planets which later appeared to be embedded in the secondary disk studied in this work. Thus, such a choice of the secondary disk size can be justified by the general idea that the region of the snow line should be most populated by planets.

A similar rough estimate is derived for τr,R from Eqs. (21), (22), and (32) using the assumptions already made in order to obtain Eq. (32). In addition, it is necessary to employ the radial power-law scalings of the optical thickness in the vertical direction as well as the vertical scale height which are derived within QSD. We obtain the following dependencies: τ_z,R ~ r^–1 and H ~ r^5/4. Finally, we found $$\begin{array}{l}{\tau }_{\text{r},\text{R}}=\underset{r}{\overset{r_{\text{out}}}{{\displaystyle \int }}}\frac{{\tau }_{\text{z},\text{R}}}{H}\text{d}{r}^{\prime }\approx \frac{{\tau }_{\text{z},\text{R}}(r)}{H(r)}\underset{r}{\overset{r_{\text{out}}}{{\displaystyle \int }}}{\left(\frac{r^{\prime }}{r}\right)}^{-9/4}\text{d}{r}^{\prime }={\tau }_{\text{z},\text{R}}(r)\frac{r}{H(r)}\\ \begin{array}{l}\times \frac{4}{5}\left[1-{\left(\frac{r}{r_{\text{out}}}\right)}^{5/4}\right]\approx \frac{{\tau }_{\text{z},\text{R}}(r)}{h}\frac{r_{\text{out}}-r}{r_{\text{out}}}\approx 1.0\left(\frac{r_{\text{out}}-r}{r_{\text{out}}}\right)\hfill \\ \times {\left(\frac{h}{0.05}\right)}^{-1}\left(\frac{{\kappa }_0}{3.0{\text{cm}}^2{\text{g}}^{-1}}\right){\left(\frac{\alpha }{0.01}\right)}^{-1}{\left(\frac{r}{2AU}\right)}^{-1}\left(\frac{{\dot{M}}_{\text{acc}}}{{\dot{M}}_{\text{acc}}^{\lim }}\right),\hfill \end{array}\end{array}$$(33)

where ${\dot{M}}_{\text{acc}}^{\lim }$ is the previously mentioned threshold value for the formation of a disk with a radial size of 2 AU, which is optically thick in the radial direction. For a given normalization, this value is ${\dot{M}}_{\text{acc}}^{\lim }\approx 1.7\times {10}^{-11}{\text{M}}_{\odot }{\text{yr}}^{-1}$.

Equation (33) implies that disk of size r ~ 2 AU sufficient to cover the substantial part of preexisting planetary system is formed provided that an accretion rate ${\dot{M}}_{\text{acc}}\gtrsim {\dot{M}}_{\text{acc}}^{\lim }$. Clearly, if we considered the formation of a larger (smaller) disk, the value of ${\dot{M}}_{\text{acc}}^{\lim }$ would become larger (smaller). By our calculations of NSD, we check that the estimates (32), (33) are in reasonable accordance with an accurate value of optical thickness, so that, indeed, τ_r,R(r) ≳ 1 in the significant part of the area where the disk can exist, r ∈ (r_in, r_out), as soon as the accretion rate exceeds its threshold value. Equation (33) also shows that τ_r,R exceeds τ_z,R by a factor of ~1/h in a geometrically thin disk.

The value of ${\dot{M}}_{\text{acc}}^{\lim }$ estimated above is used to define the moment of disk formation, t_frm. The corresponding disk profiles are shown in Fig. 3.

4.1.4 Disk evolution for the most massive donor

Here we describe general features of NSD obtained in our simulations. Fig. 3, column 1 demonstrates the disk radial structure in the binary with the donor initial mass M₁ = 7.0 M_⊙ and intermediate initial separation a = 30 AU. The selected time moments are given in Table 2, while some informative values describing the disk at t = t_RGB and t = t_AGB are given in Table 3. We note that these Tables provide results for all donors, whereas we relegate the description of disks produced by other donors to the next section.

In Fig. 3, the curves taken at t_frm show the appearance of the disk when the accretion rate for the first time attains ${\dot{M}}_{\text{acc}}\simeq {\dot{M}}_{\text{acc}}^{\lim }$ due to the growing wind at the beginning of the donor RGB stage. We note that t_frm is sufficiently close to t_RGB. Thus, we do not consider the structure of the disk at earlier stages when our crude assumption about the velocity of the wind is not valid. The curves generated at t = t_AGB represent the most massive disk ever accumulated during the evolution of the binary. The curves at t_dpl show the depleting disk ~0.5 Myr after the envelope loss by the donor.

We find that a sufficiently large disk is formed already at the RGB stage. Despite a decrease of the wind rate by order of magnitude between t_RGB and t_AGB, the disk persists for ~5 Myr until the envelope loss by the donor because ${\dot{M}}_{\text{acc}}\gtrsim {\dot{M}}_{\text{acc}}^{\lim }$ stays all the time. Generally, the AGB stage of the donor is much shorter than its RGB stage. At the same time, the disk at the RGB stage is much smaller than its follower at the AGB stage, check the disk size values in Table 3. This is obviously explained by a much stronger wind at the AGB stage. The ratio of the peak wind rates at t_AGB and t_RGB is of order ~10⁵. The peak intensity wind at t = t_AGB gives birth to the largest disk with a size exceeding ~11 AU. Such a value substantially exceeds the corresponding Bondi radius, that is, the distance from the host star limiting the wind accretion, which we find to be always within the range 2– 4 AU, with an exact value weakly depending on the initial binary separation. At the same time, the disk never reaches the tidal truncation radius, which takes the instant value r_out = 28.23 AU at t = t_AGB in the system with a = 30 AU. Shortly after the wind attains the peak intensity, the donor loses the envelope and the wind ceases: the disk starts its closed evolution characterized by an accretion timescale. An exact way of this closed evolution is determined by the structure of disk at t = t_env, which we finally obtain simulating NS disk during the whole binary evolution. It takes approximately 1.1 Myr for the disk matter to accrete onto the host star (see the difference between t_dec and t_env in Table 2). In principle, a new generation of planets can start to grow in long-lived massive disks, but analysis of this possibility is beyond the scope of our study. We check that the lifetime of the residual disk remaining after the envelope loss by the donor is in reasonable accordance with the viscous timescale estimated at the outer edge of the disk at t = t_env.

The low-mass disk accumulated by t = t_frm has the surface density varying from Σ ~ 1–10 g cm^–2 in the inner parts down to Σ ~ 10^–1 g cm^–2 at the outskirts. The midplane temperature is found in the range T_c ~10²–10³ K corresponding to a cool disk heated mostly by the host star and having a small aspect ratio approaching h ~ 0.01. We note that temperature does not decrease even lower due to the donor heating: the profiles of T_c and h have a distinctive hump at the outskirts of the disk (see curve at t_frm in column 1 of Fig. 3 and see Sect. 4.1.2 for the general description of the donor influence). The donor heating is particularly efficient in R2 since the disk radiative cooling is suppressed by a factor of low vertical optical thickness. At the same time, dust starts to evaporate in the innermost part of the disk ~0.1 AU, as the temperature rises. This results in the transition of the disk to the optically thin state. Thus, the temperature rises further and exceeds ~10³ K. In this case, the innermost part of the disk which is mostly heated by the host star, becomes puffed for a reason similar to that described just above for the disk outskirts: its radiative cooling is suppressed as compared to the optically thick middle part of the disk. We find that the low-mass disk described above has a mass as small as ~10^–6 M_⊙.

The largest disk accumulated by t = t_AGB has a surface density up to Σ ~ 10³–10⁴ g cm^–2 in the inner parts down to Σ ~ 10 g cm^–2 at the outskirts. The midplane temperature of the largest disk attains ~10³ K much farther away from the host star as compared with the low-mass disk, compare the corresponding curves at t_frm and t_AGB in Fig. 3. The snow line shifts beyond r = 5 AU. At the same time, the inner parts of the disk where dust is evaporated approach the size of ~0.7 AU. The corresponding optical thickness substantially decreases; however, the disk remains optically thick. Hence, the enhanced radiative cooling makes the inner part of the disk dense and geometrically thin, note the decrease of h by a few times as changing from the outskirts to the inner parts of the disk shown in column 1 Fig. 3. This feature becomes more distinct due to the donor heating, which makes the disk outskirts geometrically thicker similar to that of the low-mass disk. The disk acquires a distinctly flaring shape. We find that the innermost part of the largest disk is heated up to T_c ≳ 10⁴ K. The pure gas opacity gets sharply higher at such temperature due to ionization processes that inflate the disk close to its inner boundary (see the corresponding hump on the profile of h at r ≲ 0.05 AU in Fig. 3, column 1). In this case, the surface density is found to have decreased considerably. We also note that the disk accumulates mass up to 0.0013 M_⊙ at t_AGB.

The residual disk left after the envelope loss by the donor is shown at half of its lifetime (see the t_dpl curves in column 1 of Fig. 3). It is seen that despite the considerable time ~0.5 Myr has passed since the envelope loss, the disk still holds a substantial amount of matter. Its mass and size are ~10^–4 M_⊙ and ~6 AU, respectively. We note that the disk radial structure becomes quite smooth in comparison with the epoch of wind accretion. The surface density profile is close to a power law, which indicates that the disk evolution proceeds well within a quasi-stationary regime (Lynden-Bell & Pringle 1974; Lyubarskij & Shakura 1987). The disk is distinctly flaring with an aspect ratio substantially decreased as compared to the moment t = t_AGB. An aspect ratio varies from ~0.03 at the outskirts of the disk down to ~0.015 in its innermost part. One can see that, similar to the low-mass disk produced at the wind rate corresponding to ${\dot{M}}_{\text{acc}}\simeq {\dot{M}}_{\text{acc}}^{\lim }$, such a disk is already light enough to become optically thin due to evaporation of dust in its innermost part. The residual disk has a temperature and aspect ratio rather close to that of the low-mass disk; however, it is much larger and heavier than the low-mass disk generated by the weak wind from the donor. The results concerning the residual disk presented here allowed us to predict the appearance of extended cool accretion disks around MS stars in wide binaries with young WDs.

As the initial separation of the binary becomes larger, the amount of wind matter captured by the disk decreases. This shortens its lifetime and decreases its mass and size. As above in this Section, we describe here general results with an emphasis on the sufficiently large disk with ${\dot{M}}_{\text{acc}}\gtrsim {\dot{M}}_{\text{acc}}^{\lim }$, or equivalently, the disk size not less than ~2 AU. For a increasing up to 100 AU, such a disk does not exist far away from t = t_RGB anymore. Nevertheless, for the widest system with a = 100 AU, disk approaches the size (the mass) r ~ 2 AU (10^–6 M_⊙) right at t = t_RGB. At the same time, the existence of a disk with a size not less than ~2 AU before t = t_AGB shortens up to 0.5 Myr. However, disk holds the size (mass) around ~1 AU (10^–7 M_⊙) for a considerably longer time ~5 Myr. For the initial binary separation increasing from 30 AU to 100 AU, the mass of the largest disk accumulated by t = t_AGB decreases by half an order of magnitude taking the values down to 0.0003 M_⊙. The corresponding disk size takes the values down to 8 AU (compare this value with the corresponding value in Table 3). After the envelope loss by the donor, the disk decays approximately twice as fast as the binary initial separation increases from 30 AU up to 100 AU. Thus, the residual disk lifetime in the system with initial a = 100 AU is around ~0.5 Myr.

In the binaries closer than initial a = 30 AU, the disk becomes larger as compared to the case demonstrated in Fig. 3 and Table 3. For example, the system with the initial separation a = 20AU gives birth to the largest disk we meet in calculations. It attains the size ~18 AU about 50 kyr after t = t_AGB. However, the disk produced in close binaries with initial separation approaching 10 AU extends up to the tidal truncation radius around t = t_RGB as well as t = t_AGB. In this case, we obtain the most massive but more compact disk. At t = t_RGB, its size stays within ~2 AU. Despite that, the disk contains much more matter as compared to a wider binary: the disk mass attains 3 × 10^–5 M_⊙, which should be compared to an order of magnitude smaller value given in Table 3. Similarly, at t = t_AGB, we obtain the disk a few times more massive as compared to the binary with initial a = 30 AU. Its mass and size are, respectively, 0.0071 M_⊙ and 12.7 AU.

4.1.5 Disk evolution for less massive donors

The results outlined in the previous Section basically apply to less massive donors. The disk radial structure is shown for donors with initial masses M₁ = 5.0, 3.0, 1.7 M_⊙ in Fig. 3, Cols. 2–4, respectively. As in the previous Section, we first analyze the disk properties in binaries with the initial separation a = 30 AU and then explain the changes corresponding to the changes of the initial binary separation. Being generated at their own t = t_frm, t_AGB, and t_dpl, the corresponding radial profiles of disk look qualitatively similar to those for the most massive donor, compare Cols. 2–4 to Col. 1 in Fig. 3. At the same time, the number of significant differences can be identified by the values contained in Tables 2 and 3. Both qualitative and quantitative differences between the disks produced by different donors are determined by variations of the wind rate profile depending on the donor mass (see Figs. 1–2). Below we highlight the main deviations from the disk formation and evolution obtained for donor M₁ = 7.0 M_⊙ while proceeding with less massive donors. Additionally, we point out some important features of the disk that are common to all donors.

Contrary to the case with the most massive donor, sufficiently large disk with ${\dot{M}}_{\text{acc}}\gtrsim {\dot{M}}_{\text{acc}}^{\lim }$ is not formed at the RGB stage of intermediate-mass donors with M₁ = 3.0 M_⊙ and 5.0 M_⊙. However, such disk is formed in the case of low mass donor with M₁ = 1.7 M_⊙. In this case, the disk exists for 2.5 Myr before t = t_RGB and attains the size of 3.5 AU which is even larger than the corresponding disk produced by the most massive donor (see Table 3). We note that its mass is twice that of the disk produced in the binary with M₁ = 7.0 M_⊙ and is greater by almost two orders of magnitude than the disk produced in the binary with M₁ = 3.0 M_⊙. Such a situation takes place due to the peculiarities of stellar evolution. Indeed, the lowest mass red giant demonstrates a more powerful and long-standing wind than heavier donors, comparing the plots in Figs. 1, 2. As soon as the wind from the low-mass red giant abruptly weakens, the disk depletes by ~0.1 Myr after t = t_RGB and shrinks back to a small size. The intermediate-mass donors with M₁ = 3.0M_⊙ and 5.0 M_⊙ produce a sufficiently large disk for the first time, respectively, ~0.4 Myr and ~2.4 Myr before the peak intensity wind at t = t_AGB. Whereas the lowest mass donor does the same once again ~1 Myr before t = t_AGB (see Table 2).

Since the ratio of the peak wind rates at t_AGB and t_RGB stays in the range 10⁴–10⁶ for less massive donors, their AGB stage also produces the largest and the most massive disk throughout the whole evolution of the system (see Table 3). The largest disk size weakly depends on the donor mass and attains ~13 AU for the intermediate mass donor with M₁ = 5.0 M_⊙. As in the case of the most massive donor, the largest disk size substantially exceeds the corresponding Bondi radius for less massive donors. We checked that for a = 30 AU no donor produces the disk at the AGB stage expanding up to the tidal truncation radius. The latter takes values from r_out = 26.07 down to 14.08 AU at t = t_AGB while changing from M₁ = 5.0M_⊙ down to M₁ = 1.7 M_⊙, respectively. We note that despite the same initial separation of the binary, r_out increases by t = t_AGB more for the heavier donors, which is explained by the growing expansion of the binary due to the more intensive mass loss by the heavier donors. We also check that the closed evolution of the residual disk produced by less massive donors lasts not less than in the case of the most massive donor. The residual disk exists even longer for the intermediate-mass donors decaying 1.1–1.5 Myr after they lose their envelopes, check Table 2.

We note that the low-mass disks shown in Fig. 3 for t = t_frm are nearly identical to each other having the same mass ~10^–6 M_⊙ for all donors. This is not the case at t = t_AGB when the surface density gradually increases by a factor of a few in the middle of the disk while changing from the lowest mass donor to the most massive one. Accordingly, the largest disk accumulated by t = t_AGB is more massive in binaries with heavier donors (see Table 3). Contrary to the cases of the most massive donor and the donor with 5.0 M_⊙, the innermost part of this disk produced by the less massive donors, 1.7 M_⊙ and 3.0 M_⊙, is not puffed (see the corresponding values at t_AGB in Fig. 3). It takes the smallest value of an aspect ratio slightly above 0.01, while its temperature stays well below 10⁴ K. Additionally, the snow line typically shifts closer to the host star in a disk produced by the less massive donors. For example, it is located at r = 3.3 AU and 5.8 AU for donors 3.0 M_⊙ and 5.0 M_⊙, respectively (see Table 4).

It is notable that the peak values of the donor wind rate at t = t_AGB for the lowest mass donor, 1.7 M_⊙, and the most massive one, 7.0 M_⊙, show a considerably larger difference than masses accumulated by disk in these two cases (see Figs. 1, 2 and Table 3). This is explained by the variation of an expansion rate of the binary due to the decrease of the donor mass in the systems with different donors. Indeed, the different expansion rates of the binaries with different donors change the ratio of the cross sections for the gravitational capture of wind in systems with the same initial separation, and finally, the ratio of accretion rates in the course of binary evolution. Using Eq. (4), we make sure that by t = t_AGB, the binaries with initial separation a = 30 AU expand by a factor of ~2 and ~4 for the donors with M₁ = 1.7 M_⊙ and 7.0 M_⊙, respectively. Since the corresponding values of Bondi radius in these two systems are almost the same (see Table 4), we find that the growth of the disk in the systems with heavier donors is inhibited mostly by the relatively faster expansion of the binary. For the same reason, the largest disk size weakly depends on the donor mass (this has already been pointed out above in this section). Moreover, the size of such a disk in the system with donor M₁ = 5.0 M_⊙ slightly exceeds that of the disk in the system with the most massive donor, M₁ = 7.0 M_⊙ (see Table 3).

The residual disk left after the envelope loss by the donor evolves similarly for all donors. As we change to the lowest mass donor, its mass and size by t = t_dpl decreases down to ~10^–5 M_⊙ and ~3 AU, respectively. We note that the humps of temperature and aspect ratio arising in the optically thin innermost part of the disk become more significant for the less massive donor, compare the corresponding curves at t_dpl shown in Fig. 3.

The less massive donors in wider systems with initial separation increasing up to 100 AU do not produce a sufficiently large disk at t = t_RGB. In this case, disk size (mass) rises up to ~1 AU (10^–7M_⊙) for the lowest mass donor. As a increases up to 100 AU, the existence of a disk with a size not less than ~2 AU before t = t_AGB shortens up to ~0.3, and 0.4 Myr for donors with M₁ = 1.7, and 5.0 M_⊙, respectively, while it remains almost unchanged for donor with M₁ = 3.0 M_⊙. Nevertheless, the disk holds size (mass) around ~1 AU (10^–7 M_⊙) for a considerably longer time ~0.7 and 2.4 Myr for donors with M₁ = 1.7 and 5.0 M_⊙, respectively. A different trend in the formation of the disk in wider systems with donor M₁ = 3.0 M_⊙ is explained by the too-low wind rate right up to the beginning of the AGB stage. We check that for the initial binary separation increasing from 30 AU up to 100 AU, the mass of the largest disk accumulated by t = t_AGB decreases by half an order of magnitude for all less massive donors. For example, the size (mass) of such disk becomes as small as 6 AU (7 × 10^–5 M_⊙) for a = 100 AU in the binary with M₁ = 1.7 M_⊙.

In the closer binaries with a down to 10 AU, the disk attains the tidal truncation radius around t = t_RGB and t = t_AGB for all less massive donors. In this case, each given donor produces the most massive but compact disk. At t = t_RGB, the disk size stays within ~3 AU for the lowest mass donor with M₁ = 1.7 M_⊙ and even smaller for the intermediate mass donors. The corresponding disk mass takes values 3 × 10^–5, 5 × 10^–7, and 8 × 10^–7 M_⊙ for the donors M₁ = 1.7, 3.0, and 5.0 M_⊙, respectively. These values should be compared with the corresponding data in Table 3. Similarly to the case of the most massive donor, the less massive donors produce at t = t_AGB a few times more massive but more compact disk as compared to the corresponding binary with initial a = 30 AU. The disk mass (size) is 0.0011, 0.0036, 0.0058 M_⊙ (6.4, 8.9, 10.5 AU) for donors 1.7, 3.0, 5.0 M_⊙, respectively.

4.2 Planetary migration

4.2.1 Analytical estimates concerning planetary migration

We consider migration of a single constant mass planet embedded in the NS disk. The type I migration timescale can be estimated analytically in the low-mass disk within QSD under the assumption T_c ≈ T_I,A. The same assumptions were used to obtain Eq. (32). With these assumptions and the given values of parameters from Eq. (32), that is, ${\dot{M}}_{\text{acc}}={\dot{M}}_{\text{acc}}^{\text{lim}}$, r = 2 AU, we obtained the following upper limit on the migration timescale, $$\begin{array}{l}{{\widehat{t}}_{\text{p}}|}_{\text{I}}\approx 2.8\times {10}^8\text{yr}{\left(\frac{2.72+1.08\beta }{4.34}\right)}^{-1}\left(\frac{\alpha }{0.01}\right){\left(\frac{m_{\text{p}}}{{\text{m}}_{\oplus }}\right)}^{-1}\\ \times {\left(\frac{r}{2\text{AU}}\right)}^2{\left(\frac{{\dot{M}}_{\text{acc}}}{{\dot{M}}_{\text{acc}}^{\lim }}\right)}^{-1}.\end{array}$$(34)

Equation (34) is an estimate of the migration rate taking place in a low-mass disk, which we conventionally defined in Sect. 4.1.3 as the one with ${\dot{M}}_{\text{acc}}\simeq {\dot{M}}_{\text{acc}}^{\lim }$. Therefore, the migration of the Earth-like planets is expected to be negligible in such a disk. At the same time, the Earth- to the Neptune-like planets on the close orbits within 1 AU are expected to migrate weakly given the long-term existence of a low-mass disk with ${\dot{M}}_{\text{acc}}\simeq {\dot{M}}_{\text{acc}}^{\lim }$ before t = t_RGB. According to Eq. (34), the type I migration should be significant for giant planets provided that the low-mass disk exists at least for ~1 Myr. We further expect that giant planets embedded in the largest disk around t = t_AGB should reach the distance of TM, $r_{\text{p}}^{\text{fin}}$. On the other hand, the Earth- to the Neptune-like planets can migrate far from their original orbits at the AGB stage since the accretion rate is enhanced by a few orders of magnitude for hundreds of kyr. These conclusions on planetary migration at the AGB stage should be the case for planets that preexisted on orbits within ~10 AU which is constrained by the disk size (see Table 3).

Giant planets can open the gap and undergo the type II migration. We do not consider the other types of planet influence on the disk. As soon as the viscous timescale of a disk, t_ν, is significantly less than its lifetime (see Fig. 4), we expect that type II migration alone can shift the planets up to $r_{\text{p}}^{\text{fin}}$. It is instructive to derive the approximate equation for q_crit valid for values of disk viscosity and aspect ratio met in this work. As far as the parameter X in Eq. (26) is close to unity even for h highly exceeding α, we have $$q_{\text{crit}}\approx \frac{100}{ℛ}{\left(2^{1/3}-\frac{1}{8}{\left(\frac{h}{\alpha }\right)}^{1/3}\right)}^{-3},$$

where the second term in brackets can also be safely omitted for our case for α = 0.01 and h ≲ 0.1. We finally obtain $$q_{\text{crit}}\approx 1.25\times {10}^{-3}\left(\frac{\alpha }{0.01}\right){\left(\frac{h}{0.05}\right)}^2.$$(35)

According to this estimate, we expect that the opening gap mass of the planet varies in the range $m_{\text{p}}^{\text{crit\hspace{0.17em}}}=15-1500{\text{m}}_{\oplus }$ as we find that the thickness of the NS disk corresponds to h ~ 0.01 close to the time of the disk formation or decay and attains h ~ 0.1 at t = t_AGB (see Sects. 4.1.4 and 4.1.5). As soon as the opening gap mass strongly depends on the aspect ratio, $m_{\text{p}}^{\text{crit\hspace{0.17em}}}\propto h^2$, we expect that the outcome of migration of the giants embedded in NS disk is determined by a sequence of type I and type II migration periods changing each other. The change between the migration types can be caused by either the radial gradient or the time evolution of the disk aspect ratio.

We estimated the ratio between the type I and II migration rates for the gap opening mass of the planet, q = q_crit. To do so, we compared the velocities (25) and (28) in order to obtain an estimate: $$\begin{array}{l}\frac{\text{d}{r}_{\text{p}}/{\text{d}t|}_2}{\text{\hspace{0.17em}d}{r}_{\text{p}}/{\text{d}t|}_1}=\frac{v_2}{v_1}=\frac{1.5\alpha h^2{f}^{-1}{M}_{\text{d}}/\left(M_{\text{d}}+m_{\text{p}}\right)}{(2.72+1.08\beta )h^{-2}{m}_{\text{p}}{M}_{\text{d}}/\left(2\pi M_2^2\right)}\\ =\frac{3\pi \alpha h^4}{(2.72+1.08\beta )f}\frac{M_2^2}{m_{\text{p}}^2}\frac{1}{1+M_{\text{d}}/m_{\text{p}}}=0.98\left(\frac{\alpha }{0.01}\right){\left(\frac{h}{0.08}\right)}^4\\ \times {\left(\frac{f}{1}\right)}^{-1}{\left(\frac{2.72+1.08\beta }{4.34}\right)}^{-1}{\left(\frac{m_{\text{p}}}{300{\text{m}}_{\oplus }}\right)}^{-2}{\left(\frac{M_2}{{\text{M}}_{\odot }}\right)}^2\frac{1}{1+M_{\text{d}}/m_{\text{p}}},\end{array}$$(36)

where the notations from Eqs. (25), (28) are used. It is also assumed that migration occurs in a quasi-stationary disk, which allowed us to rewrite the term (2r/F)(dF/dr) = 1/f, where f is the factor from Eq. (14). Equation (36) reproduces the known result that the type I migration proceeds faster than the type II migration for sufficiently massive planets in cold and low-mass disks (h ≲ 0.08 and M_d ≲ m_p) (see Ward 1997). We note the strong dependence of this feature on the disk aspect ratio.

Therefore, the fast type I migration takes place as soon as the ratio given by Eq. (36) becomes less than unity for m_p given by Eq. (35). Explicitly, $$\frac{v_2}{v_1}=0.078{\left(\frac{\alpha }{0.01}\right)}^{-1}{\left(\frac{f}{1}\right)}^{-1}{\left(\frac{2.72+1.08\beta }{4.34}\right)}^{-1}\frac{1}{1+M_{\text{d}}/m_{\text{p}}}.$$(37)

This estimate shows that the fast type I migration is expected in any disk with sufficient viscosity. We note that this is not the case in Kulikova et al. (2019), where the type II migration is considered in the limit of a heavy disk only, while the disk characteristic mass is actually smaller than the opening gap mass of the planet. The latter takes place because Kulikova et al. (2019) do not consider the relatively dense part of the disk beyond the snow line.

4.2.2 General description of migration

The examples of migration tracks are shown in Fig. 6 and in Fig. A.1 for planets m_p = 10 m_⊕ and 100 m_⊕ and for all donors in the system with intermediate initial separation. The dependence of the planetary orbit on time is accompanied by the migration timescale as well as the disk aspect ratio and the critical mass of the gap formation in the vicinity of the migrating planet. These additional quantities help one follow the condition of gap formation around the giants. The initial planetary orbits are set to 1 and 4 AU, with the latter being close to the largest stable orbit, $r_{\text{p}}^{\max }$, in the particular case of high-mass donors in the systems with initial a = 30 AU. We show migration around the AGB stage only. However, it can be seen that the planet with m_p = 100 m_⊕ in the system with a donor with M₁ = 7.0 M_⊙ and planets m_p = 10, 100 m_⊕ in the system with a donor with M₁ = 1.7 M_⊙ have already considerably changed their initial 1 AU orbit by this time. This is explained by migration in a disk accumulated at the RGB stage. In turn, this makes a lighter planet reach the distance of TM before the donor loses its envelope. The latter is not the case in systems with donors with M₁ = 3.0 M_⊙ and 5.0 M_⊙. At the same time, a heavier planet migrates up to $r_{\text{p}}^{\text{fin}}$ irrespective of the donor mass and initial orbit size.

It is seen that distant planets start migrating not long before t = t_AGB when the disk spreads beyond their initial orbit. Also, a distant giant with m_p = 100 mθ always starts migrating according to the type I. Then, in systems with low mass donors, migration can switch to the type II in the innermost part of a disk due to oscillations of the accretion rate and the disk parameters.

Once a low-orbit giant appears in a disk, it usually opens the gap for the following reasons. First, $m_{\text{p}}^{\text{crit}}$ depends mainly on the disk aspect ratio according to Eq. (26) or (35). Second, the disk aspect ratio is smaller closer to the host star and, additionally, at an earlier time compared to the AGB peak. If it were a stationary system, the planet would undergo the type II migration further on. However, the growth of the accretion rate on the way to the AGB peak leads to the growth of the disk aspect ratio, which closes the gap and migration changes to the type I. We find that this always leads to an increase in the migration rate in our model. As far as the temporal growth rate of h exceeds the rate of its decrease around the planet migrating inward in a flaring disk, the type I migration holds further on. On the contrary, as soon as the temporal growth of h slows down or ceases, the decrease of h around the planet migrating inward in a flaring disk opens the gap back again. Numerical analysis reveals an additional type of migration which may last for a substantial time. In our simplified model of migration, this new type of migration exhibits jumps between the migration of the types I and II (see Figs. 6 and Fig. A.1, left plot). This situation takes place provided that h of a disk grows up faster (slower) than it decreases around the planet changing its orbit with the type II (I) migration rate (see Sect. 4.2.3 for details).

The profile of h measured in the vicinity of a planet migrating in NS disk may be quite complex already in the case of a simple shape of the accretion rate at the AGB peak and without the change between different types of migration (see the migration track of the low-mass planet starting from 1 AU in the systems with high-mass donors). The prominent peak of such a h closely follows the peak of the accretion rate in the neighboring case of a distant low-mass planet initially situated at 4 AU. However, in the former case of the smaller initial orbit, such a peak of h is diminished by the substantial spatial gradient of h in the NS disk down from 0.7 AU. This behavior is observed in the model close to t = t_AGB (see Fig. 3). The same effect is seen for the giant planet, where it is additionally complicated by the change between the types of migration in the low-orbit case (see Fig. 6). The shaping of h measured in the vicinity of the migrating planet in the systems with low mass donors is even more complicated by the oscillations of the accretion rate. Nevertheless, as one proceeds from the distant low-mass planet to the low-orbit giant planet, a similar trend can be found for these systems as well (see Fig. A.1, center and right plots).

We note that a relatively slow migration of low-mass planets takes place in the residual disk. Nevertheless, given the long lifetime of the residual disk, it shifts the planets a few times closer to the host star as compared to that at the AGB stage. Additionally, the type I migration timescale decreases outward due to the profiles of Σ and h specific for residual disk (see Fig. A.1, left and center plots). This facilitates the relocation of distant planets toward the host star.

With regards to giant planets, the migration tracks show that as the planet opens (closes) the gap, its migration slows down (speeds up). This feature is confirmed analytically (see Eq. (37) above). The estimated ratio of migration rates is in accordance with the numerical calculations of migration of the gas giants (see Figs. 6 and A.1, left plot).

In order to obtain a picture of an overall transformation of a planetary system caused by a planetary migration in a wind- fed disk, we plot contours of constant final orbits. First, they are shown on the plane of initial binary separations vs. initial planetary orbits (see Fig. 7). We choose the case m_p = 30 m_⊕ as here the contribution of the type II migration is still negligible for all combinations of a and r_p. At the same time, such a planet already undergoes the substantial type I migration. We find a prominent transformation of planetary orbits even for the widest binaries. For example, the planet at initial orbit 4 AU migrates down to 1 AU by the time of the disk decay in the system with the initial separation 100 AU. Moreover, in systems with initial separations less than ~35 AU, the planets at all dynamically stable initial orbits finally reach the distance of TM. It is notable that as the initial orbit decreases, planets approach the host star progressively tighter. Thus, a large number of planets on the plane (a, r_p) end up within the distance ≲ 1 AU from the host star. For example, this is the case for the planets at all dynamically stable initial orbits in the systems with initial separations less than ~65 AU.

Further, the contours of constant final orbits are shown on the plane of planetary masses, and initial orbits (see Fig. 8). We choose the case of binary with intermediate initial separation, a = 40 AU, and a massive donor with M₁ = 5.0 M_⊙. The left part of the plot should be considered as an extension to the data in Fig. 7. It shows the transformation of orbits belonging to the low-mass planets 1–40 m_⊕ undergoing the type I migration. It can be seen that the Earth-like planets shift by ~10% toward the host star. At the same time, all planets with masses ~10 m_⊕ on the dynamically stable initial orbits relocate to orbits ≲2 AU with a small number of them reaching the distance of TM in such a system, $r_{\text{p}}^{\text{fin}}$. Hereafter, the situation when a planet reaches $r_{\text{p}}^{\text{fin}}$ during the existence of the wind-fed disk is referred to as the complete migration. Planets with masses ~40 m_⊕ are subject to complete migration from all dynamically stable initial orbits. We check numerically that such a great relocation of distant planets toward the star is also caused by the significant migration in a residual disk where the migration timescale decreases outward, see the migration tracks above.

The range of planetary masses includes also the super-Jovian values ≲3000 m_⊕. We demonstrate that planets undergoing complete migration are limited from above by the value ~300 m_⊕. This number just weakly depends on the initial planetary orbit (see the right red line in Fig. 8). The migration tracks of the low- orbit super-Jupiters, r_p ≲ 2 AU and m_p ≳ 10³ m_⊕, show that their migration (being the type II) is slowed down on the account of the small ratio of the disk characteristic mass enclosed within the low planetary orbit and the planetary mass.

For planets of even smaller mass, two effects contribute to the enhancement of the overall migration:

1. the type II migration becomes less reduced due to the finite disk mass;

2. the planet no longer opens the gap in the largest disk accumulated around t = t_AGB. This causes an intermediate period of the fast type I migration. Such behavior can be seen from the contour describing the final orbit 1.0 AU in Fig. 8.

In the second case, the final orbit is almost independent of the initial orbit. This can be explained by a longer intermediate period of the fast type I migration of the planet with the same mass located farther from the host star, thus having a greater h and so a greater $m_{\text{p}}^{\text{crit}}$ in its vicinity. It is also notable that the sufficiently distant super-Jupiters with m_p ≳ 10³ m_⊕ migrate outward (to see it, follow the contour describing the final orbits 3.0 and 4.0 AU in Fig. 8). We consider this novel effect below in Sect. 4.2.4.

Fig. 6 Orbital evolution (top panel) and the corresponding characteristic migration time (second panel from top) for planets with masses mp = 10 mθ and 100 mθ at initial orbits rp = 1 AU and 4 AU. Additional panels from top to bottom show the total accretion rate as well as the disk aspect ratio and the critical planet mass of the gap formation (Eq. (26)) in the vicinity of the planet for the system with initial M1 = 7.0 M⊙, a = 30 AU. We note that the critical planetary mass of the gap formation is shown only for giants.

Fig. 7 Size of the final planetary orbit as a function of its initial value, rp, and initial binary separation, a. A planet migrates all the time it is embedded in an NS disk with α = 10−2. The curves show contours of the same size of the final planetary orbit marked by its value in AU. The masses are: M1 = 5.0 M⊙ and mp = 30 m⊕. The gray area at the upper left part of the figure shows the region of dynamically unstable planetary orbits calculated using Eq. (23). Migration contours are obtained by interpolation using the nodes taken with the step ∆a = 10 AU and ∆rp = 1 AU with additional values rp = 0.25,0.5,0.75,1.5 AU closer to the host star. The curve marked as “0.0” (red color) shows the planets that reached the TM distance, $r_{\text{p}}^{\text{fin}}$, by the moment of the disk decay.

Fig. 8 Size of the final planetary orbit as a function of the initial planetary orbit, rp, and planetary mass, mp. A planet migrates all the time it is embedded in an NSD with α = 10−2. The initial parameters are M1 = 5.0 M⊙ and a = 40 AU. The notations are the same as in Fig. 7. The vertical dashed lines additionally highlight the planets reaching the region of tidal migration.

4.2.3 The rolling migration

In this section, we discuss a special type of migration that can only take place in the NS disk. Let us consider the general case of a planet embedded in the disk with an aspect ratio variable over the radial coordinate as well as changing over time, h = h(r, t). The rate of evolution of h at the planetary position can be written as $$\frac{\text{d}h}{\text{\hspace{0.17em}d}t}=\dot{h}+h^{\prime }{v}_{\text{mig}},$$(38)

where ḣ and h′ denote the partial derivatives of h(r, t) over t and r, respectively. v_mig is the migration rate of the planet having mass m_p. We assume that the value of v_mig is equal to v1 in case of type I migration, and to v₂ for the type II migration. The migration rate is negative for the rolling migration. Hereafter, we consider the situation of the transition between the types of migration which means that h measured by the observer on the planet takes the critical value, h = h_crit, related to m_p via Eq. (35). As discussed in Sect. 4.2.2, we usually have |v₂| < |v₁| in this situation (see Eq. (37)). Therefore, if dh/dt > 0 (< 0), the gap is closed (opened) and the planet accelerates (decelerates) toward the host star changing to the migration type I (II). As far as dh/dt saves its sign during such a transition, we find an ordinary single change of the migration type. Obviously, this is always the case for the stationary disk when ḣ = 0. However, the NS disk adds a new option. We find that the sign of dh/dt changes while planet crosses the point where h = h_crit, provided that $$\begin{array}{l}\dot{h}>0,\\ h^{\prime }>0,\\ v_1<v_{\text{roll}}<v_2<0,\end{array}$$(39)

where $$v_{\text{roll}}\equiv -{\frac{\dot{h}}{h^{\prime }}|}_{h=h_{\text{crit}}}$$(40)

is the velocity of the constant disk aspect ratio. If the conditions (39) are satisfied, the planet starts migrating at an average rate equal to v_roll. The detailed hydrodynamical consideration of gravitational disk-planet interaction in a more general case of NS disk is required in order to predict the exact behavior of the migrating planet in this situation. Since the timescales of the gap opening and closure along with the timescale of crossing the gap by the line h = h_crit are all of the same order, it is possible that the migration rate is subject to oscillations around |v_roll|. We refer to this new type of planetary migration as the rolling migration, since it is determined by the propagation of line h = h_crit along the disk growing at the AGB stage. Let us also note that, generally, the other variants of the rolling migration are possible for other combinations of conditions entering Eq. (39) (see Table 5).

The conditions (39) are usually met in the NS disk considered in this work. Employing the numerical tests, we find that the rolling migration takes place for the giants 100–300 m_⊕ orbiting the host star within ~ 1 AU in the binaries with high-mass donors and intermediate initial separations. In this case, the planet is embedded in the inner part of the sufficiently large disk growing at the AGB stage, which corresponds to the zone of evaporation of dust (see Fig. 3). We note that such planets undergo complete migration. Oscillations of the wind rate taking place in binaries with low mass donors hinder the rolling migration. An example of a rather long period of the rolling migration is shown in Fig. A.2, left plot. It can be seen that, at first, the planet undergoes the type I migration which leads to a gradual decrease of $m_{\text{p}}^{\text{crit}}$ due to the decrease of an aspect ratio around the planet approaching the inner part of a flaring disk. Further, the gap opens and the planet slows down changing to the type II migration. This period of migration lasts for ~400 kyr. As soon as the system approaches t = t_AGB, the disk aspect ratio grows faster. Once the condition |v_roll| > |v₂| becomes satisfied, the period of the rolling migration starts, which lasts for ~300 kyr and ends with a return back to the type II migration. Finally, the planet enters the innermost hot part of the disk puffed due to opacity enhanced by the partial ionization of the gas. This causes the new change to the fast type I migration and the complete migration before the donor loses its envelope.

4.2.4 Outward migration in the zone of decretion

Here we show examples of the type II migration in the zone of decretion. Details about this essentially NS variant of the disk evolution are presented in Sect. 4.1.1. Since the critical mass of the planet opening the gap is large in the outer parts of the disk, we are obliged to consider a super-Jupiter with m_p = 2000 m_⊕ (see Fig. A.2, center and right plots). Therefore, we obtain a moderate outward migration reduced by the ratio of the characteristic disk mass and the planet mass. We note that it takes place for both high and low mass donors. As soon as the accretion rate vanishes, outward migration is replaced by the usual inward migration in a residual disk, since it turns into a quasi- stationary state. If we took a smaller planet, it would undergo the fast type I migration close to the peak of the accretion rate associated with the AGB stage. Since the type I of migration is directed inward in our simple model, it overwhelms a relatively weak outward migration due to the decretion of the disk matter. Nevertheless, an outward migration can be found in contours of constant final planetary orbits for the largest planets (see Fig. 8).

There is a caveat that the disk structure in the zone of decre- tion is far from that of a simplified isothermal disk used to obtain the type I migration rate in Tanaka et al. (2002). An abnormal negative surface density gradient, which makes the net viscous torque positive despite the negative specific net viscous torque, may reverse the type I migration. If this is the case, we should expect a different shape of contours of the constant final planetary orbits for distant planets embedded in the zone of decretion: such planets may avoid the complete migration in the whole range of planet masses (see Fig. 8).

To get a feeling of the significance of possible outward type I migration, we show the characteristic timescale of decretion and the corresponding migration track of the planet embedded in the zone of decretion. This corresponds to the situation when the planet migrates at the same rate as the gas around it (see the upper panel of Fig. A.2, right plot). Clearly, the rate of migration significantly exceeds the one obtained for the outward migration in our simplified model and leads to fast escape of the planet to the zone of dynamically unstable orbits. The possible outward type I migration of distant planets may change substantially the overall transformation of a planetary system caused by migration in NS wind-fed disk, as well as the picture of the complete migration of giant planets (see Sect. 4.2.5 below). The gravitational interaction between a planet and an NS disk with the zone of decretion should be the subject of a future study. The zone of decretion with sufficiently large surface density may probably lead to other effects such as the formation of planets in density traps as discussed in Alessi et al. (2020).

4.2.5 Complete migration

In this Section, we consider the planets subject to complete migration (CM). We define the CM as follows: migration from the initial orbit to the innermost orbit close to the star where additional processes or the orbital evolution (e.g., the tides) are important.

We represent this type of migration in the plane (a, r_p) separately for each donor (see Figs. 9 and Fig. A.3). There are contours for planets of several given masses. The bottom left part of the plane (a, r_p) concerning each contour introduces the CM of planets with the corresponding mass. We find these plots to be qualitatively similar to each other. Most importantly, there is always a limiting contour of CM which shows the largest initial orbit of the planet reaching $r_{\text{p}}^{\text{fin}}$ for the binary with a specified initial separation regardless of the planetary mass. Conversely, all planets preexisting in the binary on orbits above this limiting contour are not subject to CM in our definition. Thus, either they survive in an evolved binary or they are ejected beyond the Roche lobe of the host star due to a dynamical instability of their new orbits.

The planetary mass representing the limiting contour of CM belongs approximately to the range 60–80 m_⊕ for all donors. We check that all lower mass planets representing their contours of CM reach $r_{\text{p}}^{\text{fin}}$ undergoing the type I migration. Therefore, the shape of the limiting contour is similar to that of CM for lower mass planets (see the contours for planets 10 m_⊕ and 30 m_⊕ in Figs. 9 and Fig. A.3). In the case of the high-mass donors, 5.0 M_⊙ and 7.0 M_⊙, the limiting contour of CM is close to that of CM of a planet with m_p = 100 m_⊕. For example, in the binaries with intermediate initial separation a = 40 AU, the largest initial orbit of CM exceeds the largest dynamically stable orbit for all donors except for the lowest one, whereas it becomes smaller than 4 AU or even 3 AU in the binaries with initial a = 60 AU (see Figs. 9 and A.3). The planets with m_p ≲ 10 m_⊕ undergo CM in close binaries starting with initial orbits r_p ≲ 2 AU. The intermediatemass planets, 30 m_⊕, undergo CM for the whole range of binary initial separations including the distant planets with r_p ≳ 2 AU in sufficiently close binaries with a ≲ 40 AU. CM of the lower mass planets becomes a little weaker as one switches from high-mass to low mass donor.

The existence of the limiting contour of CM is explained as follows. For larger planetary masses, the increasing rate of the type I migration makes the planet reach $r_{\text{p}}^{\text{fin}}$ from a larger initial orbit. This is no longer the case as soon as episodes of the type II migration appear. For relevant values of r_p, a, and m_p the planet opens the gap for the first time in a residual disk. We check that the critical mass of the gap formation decreases down to the planet’s mass as soon as the disk cools down in the absence of wind, and additionally, the planet reaches the inner part of the disk which has the least aspect ratio due to the evaporation of dust. Under these conditions, the further increase of the planetary mass leads to longer episodes of relatively slow type II migration. Therefore, the largest initial orbit of CM stops growing for an increasing planetary mass in the given binary. We note that the final period of migration yielding the limiting contour of CM always occurs with a closed gap because it takes place in the innermost hot region of the young residual disk. This hot region is puffed due to opacity enhanced by partial ionization of the gas. Therefore, the final period of migration in this case occurs at a comparatively high rate. As the planet becomes even more massive, it holds the gap opened all the time it stays in the residual disk. This slows down the migration and shifts the contour of CM back to more close binaries and lower initial planetary orbits on the plane (a, r_p). This proceeds unevenly along the contour. Thus, we find that CM ceases in the binaries with the largest initial separations a ≳ 80 AU as the planet mass increases from 100 m_⊕ up to 150 m_⊕ regardless of the planet’s initial location (see Figs. 9 and A.3, left plot for the case of high-mass donors). At the same time, CM ceases at initial separations a ≳ 80 AU already for the planetary mass ~100 m_⊕ and at a ≳ 60 AU for the planetary mass ~150 m_⊕ (see Fig. A.3, center plot for the case of donor 3.0 M_⊙).

The least massive donor 1.7 M_⊙ cannot provide CM of planets with m_p = ~100 m_⊕ already at initial separations a ≳ 70 AU regardless of the planetary initial location (see Fig. A.3, right plot). Hence, a part of the contour describing CM for the largest initial binary separations is almost independent on r_p and shifts to the closer binaries as the mass of the giant further increases.

The contour of CM for Jupiters, m_p = 300 m_⊕, becomes independent on r_p for all dynamically stable orbits regardless of the donor. It corresponds to binaries with initial separation slightly exceeding 40, 35, 30, and 20 AU from the largest to the smallest donor, respectively. Such a feature of CM contour is caused by the weak dependence of the final orbit of the giant on its initial orbit. This is also in accordance with contours of partial migration represented on the plane (m_p, r_p) (see Fig. 8). In other words, distant giants catch up with close giants on the way to the host star. As discussed at the end of Sect. 4.2.2, this is caused by the different ratio of the time periods of either opened or closed gap spent by the giant as migrating across the wind-fed disk from different initial orbits. Indeed, the further away the planet is from the host star initially, the earlier it switches to the type I migration due to the growth of NS flaring disk while approaching t = t_AGB. Since the type I migration under the conditions near the gap formation is always faster than the type II migration, the total time of migration up to the distance of TM hardly increases with increasing initial orbit. We check that the structure of the NS disk in the vicinity of the planet undergoing CM does not facilitate the transition to the rolling migration. Thus, we usually find a single transition from the type II to the I migration in such a situation. Finally, the CM of Jupiters and super-Jupiters is additionally suppressed by the decrease of the type II migration rate on account of the finite characteristic mass of the disk.

Fig. 9 The complete migration of planets. The contours show planets of a given mass that reached the tidal migration distance, $r_{\text{p}}^{\text{fin}}$, by the time of the disk decay. In other words, all planets of the same mass at the given contour and below will also reach the tidal migration distance. The yellow line for mp = 1.0 m⊕ is invisible because it is below the plotted region. The donor initial mass is M1 = 7.0 M⊙. The notations are the same as in Fig. 7.

5 Discussion

5.1 Stellar wind properties

In this study, we make simplifying assumptions on the stellar wind properties. This is motivated by the main goal, which is to consider the significant planetary migration in a relatively massive disk formed during the relatively short stages of an extensive mass loss by the primary. In the first place, we are interested in the AGB stage when the stellar wind is intensive and it has a low velocity. Typically, the corresponding value is taken as ~10– 20 km s⁻¹. Such a value has been used in previous studies of accretion disks in binaries with evolved donors (e.g., Perets & Kenyon 2013; Kulikova et al. 2019).

In our model, we do not study the dependence of the disk structure on the stellar wind velocity. Instead, we use the same value 20 km s⁻¹ for the whole evolutionary track, which is unrealistic. Obviously, stars have much larger wind velocities at the main sequence stage (e.g., Cassinelli 1979). This should result in a far smaller Bondi radius and lower accretion rate. Thus, we significantly overestimate the accretion rate for a long-term formation of the disk before the RGB stage, as well as between the RGB and AGB stages of the low mass donors. Still, we find that such an overestimation does not lead to any considerable planetary migration before the RGB stage of the primary. Additionally, the wind velocity might depend on the stellar parameters (mass, radius, temperature, metallicity; see, e.g., Vink et al. 2000, and references therein). Analysis of these dependencies is out of the scope of this paper, which is to evaluate the upper range of the disk growth and planetary migration. More detailed calculations are a subject for future study.

Our computations of the NS disk show that a small and very light disk with a size less than 1 AU consisting entirely of R2 appears at an accretion rate much below the formal threshold value, ${\dot{M}}_{\text{acc}}^{\lim }$. Within our simplified model of wind, this formally corresponds to the epochs far earlier than the RGB stage, or additionally, far later than the RGB and earlier than the AGB stage in the case of the low mass donors, for the considered close binaries with a ≲ 100 AU. Of course, it is unlikely that such a disk will exist in the model with a variable wind velocity. However, we expect that such a disk may appear in wide binaries with a > 100 AU, which provides a low accretion rate close to RGB and AGB peaks of the low-velocity wind. If confirmed, the formation of a light wind-fed disk may cause moderate migration of giants in wide binaries.

There is also a possibility that the appearance of such a disk can be prevented by the counteraction of the wind from the accretor. The corresponding condition can be roughly formulated by comparing the wind ram pressure from both the donor and the accretor at the Bondi radius of the accretor (see Eq. (3)). Given the spherical symmetry of winds, the corresponding estimate shows that counteraction of the solar-like accretor is not enough to prevent the disk formation at the wind rate of the donor above ~10⁻¹⁰–10⁻¹¹ M_⊙ yr⁻¹, which corresponds to the accretion rate well below $M_{\text{acc}}^{\lim }$.

5.2 Circularization radius versus Bondi radius and process of matter settling

By circularization radius we mean a specific distance from the rotation axis of the accretion disk where the bulk of the captured matter gets into circular streamlines with a velocity close to the Keplerian rotation of the disk, and afterward, the matter cools down and settles into the disk. In realistic accretion disks the circularization radius can be much smaller than the radius of gravitational capture, r_a, which is identified in our model with a radius under which all matter settles into the disk. Thus, we assume that r_a defines the efficiency of disk accretion η: Ṁ_acc ~ ηṀ_w (see Eq. (2)). However, the situation may differ.

An exact picture of the captured wind matter settling into the disk might be quite complicated (see, e.g., Huarte-Espinosa et al. 2013). Depending on the orbital, sound, and wind velocities, the structure of the flow inside the Bondi radius can vary. Additional deviations in comparison with a simplified Bondi–Hoyle–Lyttleton accretion (see, e.g., Edgar 2004, for a review and references to early studies) may be due to the properties of the outflow from the donor, especially in the case of AGB stars (see, e.g., El Mellah et al. 2020, and references therein). The outflow can be spherically non-symmetric and the wind velocity may have a complicated profile depending on the detailed parameters of the donor (e.g., rotation and chemical composition). Additionally, the pulsations of the AGB star may influence the outflow as well as the accretion flow (Chen et al. 2017). Some simplifications are necessary for the purposes of our study as we want to study the long-term evolution of an accretion disk. Hence, we follow the approach presented in Perets & Kenyon (2013)⁶ making no difference between the circularization radius and the Bondi radius.

5.3 Photoevaporation of a residual disk

In this study, we follow the disk evolution after the donor forms a WD. The newborn WD is characterized by high effective temperature in the range from ~10^4.5 K for low mass donors up to ~10⁵ K for high mass donors (see Blouin 2024, for a recent review on WD physics). That leads to irradiation of the disk by high-energy photons. Such irradiation affects the outer parts of the disk via photoevaporation.

We closely follow Armitage (2017) and estimate the photoevaporation rate of the disk considering the irradiation of the outer edge of the disk by photons with energy sufficient to ionize and evaporate the disk matter. For simplicity, we assume that the disk consists of neutral hydrogen with threshold ionization energy 13.6 eV. This value corresponds to the high-energy tail of the spectrum of WD. We estimate the photon illumination of the disk as $${\mathrm{\Phi }}_{\text{disk\hspace{0.17em}}}={\mathrm{\Phi }}_{\text{total\hspace{0.17em}}}\frac{S_{\text{disk\hspace{0.17em}}}}{S_{\text{total\hspace{0.17em}}}}\approx {\mathrm{\Phi }}_{\text{total\hspace{0.17em}}}\frac{H_{\max }{r}_{\max }}{\pi a^2},$$

where ф_total is the luminosity of WD in the band ≥13.6 eV, S_total is the surface of the sphere with the radius a centered at WD, which is the final separation of the system after the AGB stage, S_disk is the cross-section of the disk edge as seen from WD, roughly estimated as the rectangle 2H_max × 2r_max corresponding to the largest radial and vertical sizes of the disk.

We estimate the high-energy fraction of illumination by WD integrating the black-body energy distribution over the energies from 13.6 eV to infinity. The result of integration gives the photon illumination in the range 10⁴⁴–10⁴⁶ photons s⁻¹ for WD temperatures in the range 10^4.5–10⁵ K. Finally, we obtain an estimate of the photoevaporation rate of the disk by multiplying ф_disk by the hydrogen atom mass, m_H, as $$\begin{array}{l}{\dot{M}}_{\text{evap}}={\text{m}}_{\text{H}}{\text{Φ}}_{\text{disk}}\\ \approx {10}^{-7}\frac{{\text{M}}_{\odot }}{\text{yr}}\left(\frac{{\text{Φ}}_{\text{total}}}{{10}^{46}{\text{s}}^{-1}}\right)\left(\frac{H_{\max }}{1\text{AU}}\right)\left(\frac{r_{\text{max}}}{10\text{AU}}\right){\left(\frac{a}{60\text{AU}}\right)}^{-2}.\end{array}$$(41)

The estimate (41) gives us an upper limit on the evaporation rate of the accretion disk.

Photoevaporation at rate ≲10⁻⁷M_⊙ yr⁻¹ depletes the disk much faster than its viscous evolution. Indeed, the timescale of photoevaporation is estimated as t_evap ~ M_disk ${\dot{M}}_{\text{evap}}^{-1}\gtrsim {10}^4$ yr for typical disk mass ~10⁻³M_⊙ after AGB, which is about two orders of magnitude below the viscous timescale approaching t_ν ~ 10⁶ yr. However, it is known that thermal wind from the disk surface produced by external high-energy radiation stalls under the radius of gravitational capture of matter heated in the vicinity of the host star, see Hollenbach et al. (1994). The temperature of the corresponding ionized atmosphere above the disk surface attains 10⁴ K, which yields gravitational capture of the heated matter up to r ~ 9 AU around a Solar-like star. Here, we obtain that the largest residual disk only slightly exceeds this value (see Table 3). Thus, we conclude that photoevaporation of a residual disk has little effect on its evolution and the corresponding planetary migration.

5.4 Lower disk viscosity

In this work, we have assumed that α = 0.01. However, current observations do not exclude the much lower values of α. Such uncertainty obliges us to track how the results can change. Here we briefly report additional calculations of the accretion disk evolution and the corresponding planetary migration for α = 0.001.

The lower viscosity leads to a substantial decrease of the radial velocity in a disk and the respective increase of its viscous timescale. As far as the disk evolution stays in a quasi-stationary regime, we expect the respective increase of surface density provided that the accretion rate remains the same. We note that according to the analytical stationary solutions specified by the opacity typical in protoplanetary disks, the temperature increases only slightly in this case. We find that the NSD solution is in accordance with these conclusions at the stages of donor evolution with slowly varying wind. However, as the wind rate rapidly grows toward RGB and AGB peaks, the surface density increases more slowly than it might be expected from the comparison with the corresponding NSD solution obtained for α = 0.01. The same applies to the disk mass. We attribute this to the non- stationary lag behind the corresponding QSD solution, which becomes even greater than that in the case of α = 0.01. Thus, the lower viscosity disk attains only twice as the mass of the disk described in Sect. 4.1.4 by t = t_AGB for the most massive donor, while its size does not exceed ~15 AU. We find that at that moment the zone of decretion expands inside the disk considerably farther than that in the case of α = 0.01, while the abnormal amplitude of positive radial velocity decreases considerably less than it might be expected from the difference between the viscosities. We note that the residual disk with α = 0.001 lives for more than 10 Myr, which is by order of magnitude longer as compared with the case of α = 0.01.

Due to the larger disk volume density the type I migration becomes faster. In contrast, the type II migration becomes slower due to the slower accretion. The transition between the migration types occurs for less massive planets according to Eq. (35). Thus, the distant sub-Jupiters at several AU and low-orbit Neptunes at one AU become the gap-opening planets, which is not the case for α = 0·01. We expect planetary migration to become extensive also due to the long-lived residual disk.

Let us introduce the comparative results with regards to the complete migration (CM here below) for the most massive donor. We find that the limiting contour of CM shifts to the lower planetary mass, 40–50 m_⊕. At the same time, CM covers the much larger part of the plane (a, r_p) due to the prominent type I migration. For example, planets m_p = 10 m_⊕ are subject to CM from any initial dynamically stable orbit provided that the initial separation of the binary is a ≤ 60 AU, whereas in systems with a ≤ 100 AU the same is the case for all initial orbits r_p < 6 AU. Moreover, the Earth-mass planets m_p = 1 m_⊕ undergo CM from any initial dynamically stable orbit for initial binary separations a ≤ 30 AU. In contrast, the giants migrate much slower than that in the disk with a = 0·01. So, planets m_p = 300 m_⊕ undergo CM from any initial dynamically stable orbit for initial binary separation as small as a ≤ 20 AU.

5.5 Advanced migration

During the last decades, the description of migration has been developed far beyond the basic variant introduced above (see the review by Paardekooper et al. 2023). Foremost, it was shown that the unsaturated corotation torque arising in real disks with thermal and viscous diffusion slows down type I migration of super-Earths (see, e.g., Paardekooper & Mellema 2006; Kley & Crida 2008). The generalized migration rates derived by Paardekooper et al. (2010), Paardekooper et al. (2011) and by Jiménez & Masset (2017) later on have been used to upgrade, for example, simulations of planetary populations (Dittkrist et al. 2014) and planet formation models (Guilera et al. 2019). It has become clear that the generalized type I migration is quite sensitive to the disk structure. For example, it may trap the planet in the vicinity of the opacity transitions (see, e.g., Kretke & Lin 2012). The original type II migration built on the assumption of an empty gap has been also refined. Numerical studies of the giants interacting with the surrounding accreting matter showed that the gap is never swept out completely. The remaining gap-crossing flow allows the planet to migrate independently on the accreting gas so that the actual type II migration rate may exceed its original value by a few times (see, e.g., Dürmann & Kley 2015). An advanced analytical model of type II migration suggests that gravitational torque on the planet is determined by the surface density at the bottom of the gap (see Kanagawa et al. 2018). The dependence of this surface density on the viscosity and aspect ratio of the disk as well as the mass of the planet determines the gradual transition between the two types of migration. We relegate the application of the advanced migration theory to the migration of planets preexisting in evolved wide binaries for future work.

6 Conclusions

We have studied planetary migration in a wind-fed accretion disk around a secondary component of a binary system. For that, we constructed a detailed model of the NS disk within the paradigm of the standard disk accretion. For the first time, the energy balance across the disk includes both irradiation by the host star and by the donor. We find the donor heating to be significant at the outer parts of the disk for the whole range of binary separations as well as for all the considered donor masses. This is the case due to the combination of the sufficiently small disk size and aspect ratio along with the proximity of the donor and its large size at RGB and AGB stages. At the same time, the disk energy balance is generalized onto the optically thin medium. The latter makes it possible to take into account the formation of very light disks, on the order of ~10⁻⁷ M_⊙ and even smaller, which takes place in the sufficiently wide binaries and/or for the low mass donors (see Table 3, the details in Section 4.1.5, and the discussion in Section 5.1).

We have shown that the disk exists up to 5 Myr before the envelope loss of the most massive donor in binaries with an intermediate initial separation of ~30 AU. The disk accumulates mass up to 0.01 M_⊙, reaching a size of up to 18 AU by the AGB peak from the wind of the most massive donor in closer binaries with an initial separation of ~20 AU. The further decrease of initial separation leads to the tidal truncation of the disk. The disk is non-stationary because it is supplied by a variable source of matter. This variability is provided by two occurrences. First, the rapid evolution of a donor at RGB and AGB stages. Second, the accompanying change of binary separation and mass ratio.

We find that the disk behaves approximately quasi-stationary only when sufficiently far from RGB and AGB peaks of the donor wind. At the same time, it becomes essentially non- stationary close to the maxima of the wind rate, which is expressed in the appearance of an extended zone of decretion in its outer parts beyond ~2 AU. The absolute value of the radial velocity of decreting matter substantially exceeds the standard quasi-stationary value. The zone of decretion exists up to several hundred thousand years before the AGB peak of the wind rate. The full NS disk model allowed us to determine the residual disk left after the envelope loss by the donor. It depletes for ~ 1 Myr in the quasi-stationary regime.

We found the new physical features that planetary migration acquires in the NS disk. In a growing disk, a planet undergoing the type II migration can be accelerated up to the rate that is intermediate between the rates of the type II and the faster type I migrations. It appears that the planet is trapped by the line of the critical value of the disk aspect ratio determining the gap formation, which propagates toward the host star through a flaring disk. That is why we refer to this new type of migration as rolling migration. We suppose that rolling migration is activated by the partially closed gap, though the corresponding torque on the planet should differ from that in the advanced type II migration developed by Kanagawa et al. (2018) and others in the case of a stationary disk with a constant aspect ratio. This will be the subject of future research in order to give a comprehensive hydrodynamic description of rolling migration. In the zone of decretion, the planet that opens the gap should migrate outward. Moreover, the corresponding abnormal distribution of matter in the zone of decretion should presumably provide an outward type I migration. Additionally, there is a general issue with the transition between the two types of migration in the regions of a disk with substantial radial and/or time derivatives of an aspect ratio and the surface density when considering them jointly with the radial motion of a migrating planet. This may be important to address since the characteristic time of variation of these quantities in the vicinity of the migrating planet can be comparable to the characteristic time of the gap opening and closure. However, these conjectures should be verified through the solution of the corresponding hydrodynamical problems, which are far beyond the scope of this study.

We find that planets with masses up to ~60 m_⊕ are mostly subject to the type I migration regardless of the donor mass. Already, the Neptune-like planets from initial orbits at a few AU significantly approach the host star until the disk decays – even in systems with an initial a ~ l00 AU. Moreover, the Neptune-like planets from all dynamically stable orbits reach the tidal migration distance in the close systems with an initial a ≲ 20 AU for all donors. We also checked that migration in the residual disk substantially contributes to that.

As planetary mass increases up to ~100m_⊕, the most distant planets reaching the tidal migration distance begin to open the gap in the innermost parts of the decaying residual disk. The corresponding episode of the relatively slow type II migration inhibits the further increase of the largest initial orbit of the planet subject to complete migration. Therefore, for every system with a specified initial separation, there is a limiting initial planetary orbit beyond which no one planet approaches the tidal migration distance. Even heavier planets undergo the fast type I migration in the vicinity of the AGB peak only, while the slow type II migration takes place the rest of the time. The limiting initial orbit of the planet subject to complete migration equals the largest dynamically stable planetary orbit in the systems with an initial separation of ~30 AU and ~50 AU for the low and the high mass donors, respectively. It decreases down to ~2 AU in the wider systems with the initial separation of ~80 to 100 AU for all donors.

Data availability

The data underlying this article will be shared on reasonable request to the corresponding author.

Acknowledgements

We thank the anonymous referee for helpful and insightful suggestions which allowed us to improve the paper. The authors thank Galina V. Lipunova and Nikolai I. Shakura for the insightful discussions. The authors thank Alexander S. Andryushin for providing the stellar evolution tracks calculated with MESA. This research is supported by the DFG research unit FOR 5195 ‘Relativistic Jets in Active Galaxies’ (project number 443220636, grant number WI 1860/20-1). VVZ acknowledges the support of Roscosmos. This research makes use of the open-source PYTHON libraries NUMPY, SCIPY, ASTROPY and MATPLOTLIB.

Appendix A Additional figures

Here we present our additional figures used in the corresponding sections of the main text.

Fig. A.1 Same as Fig. 6 but for the systems with M1 = 5.0 M⊙ (left plot), M1 = 3.0 M⊙ (center plot), M1 = 1.7 M⊙ (right plot).

Fig. A.2 Same as Fig. 6 but the left plot is for the planet with mp = 100 m⊕ and initial orbits rp = 2AU and 3 AU in the system with initial M1 = 7.0 M⊙ and a = 30 AU. The center plot is for the planet with mp = 2000 m⊕ and initial orbits rp = 3 AU and 4 AU in the system with initial M1 = 5.0 M⊙ and a = 30 AU. Additionally, the dashed curve shows the migration track determined by the rate of decretion. The dotted and the dot-dashed curves show, respectively, the timescale of decretion and the location of the largest dynamically stable planetary orbit. The right plot is as the center plot but for M1 = 1.7 M⊙ and for planets with initial orbits rp = 4 AU and 5 AU.

Fig. A.3 Same as Fig. 9 but for M1 = 5.0 M⊙ (left plot), M1 = 3.0 M⊙ (center plot), M1 = 1.7 M⊙ (right plot). The yellow line for mp = 1.0 m⊕ is invisible (except for the right plot) because it is below the plotted region.

Appendix B Derivation of energy balance equation

Following Nakamoto & Nakagawa (1994), the local energy balance for a geometrically thin accretion disk reads $${\dot{E}}_{\text{V}}=2F_{\text{z}},$$(B.1)

where Ė_v is the viscous energy dissipation rate per unit surface area, while F_z is net energy flux from either side of the disk.

The viscous energy dissipation rate is, explicitly, $${\dot{E}}_{\text{V}}=\nu \mathrm{\Sigma }{\left(r\frac{\partial \text{Ω}}{\partial r}\right)}^2=\frac{9}{4}\nu \mathrm{\Sigma }{\text{Ω}}_{\text{K}}^2.$$(B.2)

The net energy flux F_z is the sum of components representing various physical processes, $$F_{\text{Z}}=F_{\text{d}}-F_{\text{I},\text{\hspace{0.17em}A}}-F_{\text{I},\text{D}}-F_{\text{w}},$$(B.3)

where F_d is the radiation flux from the disk surface, F_I,A is the disk surface irradiation flux from the host star, F_I,D is the disk surface irradiation flux from the primary, and F_w is the heating flux due to the kinetic energy dissipation of the settling wind material. Below we consider each component of Eq. (B.3) in detail.

B.1 Radiation flux from the disk surface

B.1.1 Optically thick regime, τ_z,R > 1

Introducing the effective temperature of the disk surface, T_d, we have: $$F_{\text{d}}={\sigma }_{\text{B}}{T}_{\text{d}}^4,$$(B.4)

where σ_B is the Stefan-Boltzmann constant. The surface temperature of the disk is generally not equal to the midplane disk temperature, T_c. The diffusive approximation of radiation transfer yields the corresponding relation between T_d and T_c (see, e.g., Ruden & Lin 1986), $${\sigma }_{\text{B}}{T}_{\text{d}}^4={\sigma }_{\text{B}}{T}_{\text{c}}^4-\frac{3}{8}{\tau }_{\text{z},\text{R}}{\dot{E}}_{\text{ν}},$$(B.5)

where τ_ɀ,R = κ_RΣ/2 is the Planck mean optical thickness in ɀ-direction.

B.1.2 Optically thin regime, τ_z,R < 1

Following Nakamoto & Nakagawa (1994), we obtain (see Appendix C for the details) $$F_{\text{d}}=4{\tau }_{\text{z},\text{P}}{\sigma }_{\text{B}}{T}_{\text{d}}^4,$$(B.6)

where τ_ɀ,P = κ_PΣ/2 is the Planck mean optical thickness in ɀ- direction estimated by the Plank mean opacity ĸ_P taken at the disk midplane. In this case, we assume the disk to be isothermal in vertical direction, $$T_{\text{d}}=T_{\text{c}}.$$(B.7)

B.2 Irradiation flux from the host star

In this case, the disk is heated by the radiation absorbed at the surface layers of the disk. The corresponding energy flux can be obtained similarly to Chiang & Goldreich (1997). Assuming that the disk is conical we obtain $$F_{\text{I},\text{\hspace{0.17em}A}}=\frac{2}{3\pi }{\left(\frac{R_2}{r}\right)}^3{\sigma }_{\text{B}}{T}_2^4,$$(B.8)

where R₂ is the host star radius and T₂ is its effective temperature.

B.3 Irradiation flux from the primary

Under certain conditions, the disk surface can be additionally irradiated by a donor. It is possible to obtain analytically the corresponding irradiation flux for a conical disk (see Appendix D for detailed calculations). Here we give the final result, $$\begin{array}{l}{F}_{\text{I},\text{D}}=I_1\cos h\frac{\sqrt{R_1^2-H^2}}{a}\times \\ \times [\frac{{\left(R_1-H\right)}^2}{a^2}-\frac{{\left(H_{\tau }-H\right)}^2}{r_{\tau }^2}\left(1+\frac{r^2}{2r_{\tau }^2}+\frac{r}{a}\frac{r}{r_{\tau }}\right)+\\ +\tan h\left(\frac{r}{r_{\tau }}\frac{H_{\tau }-H}{r_{\tau }}-\frac{r}{a}\frac{R_1-H}{a}\right)],\end{array}$$(B.9)

where R₁ is the radius of the donor, T₁ is its effective temperature, and r_τ is the distance from the host star where the disk scale height attains a maximum, H_τ. In all cases r_τ = r_max due to conical geometry.

B.4 Heating flux due to the kinetic energy dissipation of the settling wind material

This contribution to the energy balance is taken following Perets & Kenyon (2013) as $$F_{\text{w}}=\frac{G{M}_2{\dot{M}}_{\text{acc}}}{2\pi r^2{r}_{\text{a}}}=\frac{G{M}_2{\dot{M}}_{\text{w}}{r}_{\text{a}}}{8\pi a^2{r}^2}\theta \left(r_{\text{a}}-r\right),\forall {\tau }_{\text{z},\text{R}},$$(B.10)

where all the values are defined in Sect. 2.1.

B.5 Final equation

Putting Eq.(B.2) along with Eqs. (B.4)–(B.10) together into Eq. (B.1) we obtain energy balance equations for the two cases, $$T_{\text{c}}^4=\frac{9\nu \mathrm{\Sigma }{\text{Ω}}_{\text{K}}^2}{8{\sigma }_{\text{B}}}\left(1+\frac{3{\tau }_{\text{z},\text{R}}}{4}\right)+\frac{F_{\text{I},\text{\hspace{0.17em}A}}}{{\sigma }_{\text{B}}}+\frac{F_{\text{I},\text{D}}}{{\sigma }_{\text{B}}}+\frac{G{M}_2{\dot{M}}_{\text{w}}{r}_{\text{a}}}{8{\sigma }_{\text{B}}\pi a^2{r}^2}\theta \left(r_{\text{a}}-r\right),$$(B.11)

for the optically thick case, and $$T_{\text{c}}^4=\frac{9\nu \mathrm{\Sigma }{\text{Ω}}_{\text{K}}^2}{32{\sigma }_{\text{B}}{\tau }_{\text{z},\text{P}}}+\frac{F_{\text{I},\text{\hspace{0.17em}A}}}{4{\sigma }_{\text{B}}{\tau }_{\text{z},\text{P}}}+\frac{F_{\text{I},\text{D}}}{4{\sigma }_{\text{B}}{\tau }_{\text{z},\text{P}}}+\frac{G{M}_2{\dot{M}}_{\text{w}}{r}_{\text{a}}}{32{\sigma }_{\text{B}}{\tau }_{\text{z},\text{p}}\pi a^2{r}^2}\theta \left(r_{\text{a}}-r\right),$$(B.12)

for the optically thin case.

Combining Eqs. (B.11)–(B.12) we obtain the energy balance equation which is applicable for any value of the optical thickness, τ_ɀ,R, $$T_{\text{c}}^4=\frac{9\nu \mathrm{\Sigma }{\text{Ω}}_{\text{K}}^2}{8{\sigma }_{\text{B}}}\left(1+\frac{3{\tau }_{\text{z},\text{R}}}{4}+\frac{1}{4{\tau }_{\text{z},\text{P}}}\right)+\left(T_{\text{I},\text{\hspace{0.17em}A}}^4+T_{\text{I},\text{D}}^4+T_{\text{w}}^4\right)\left(1+\frac{1}{4{\tau }_{\text{z},\text{P}}}\right),$$(B.13)

where the terms $T_{\text{I},\text{\hspace{0.17em}A}}^4$, $T_{\text{I},\text{D}}^4$, $T_{\text{w}}^4$ are related to Eqs. (B.8)–(B.10) according to the Stefan-Boltzmann law.

Appendix C Optically thin disk irradiation by the host star

Following Nakamoto & Nakagawa (1994) we assume that an optically thin disk is locally represented by a plane-parallel layer of finite thickness exposed to the ambient radiation. It is useful to define the optical thickness of a disk per unit frequency as $${\tau }_{\nu }(ɀ)=\underset{-H}{\overset{ɀ}{{\displaystyle \int }}}{\kappa }_{\nu }\rho (ɀ)\text{d}ɀ,$$(C.1)

where κ_ν is opacity at frequency ν.

The radiation transfer equation for the intensity of radiation with frequency ν yields the solution, $$I_{\nu }\left({\tau }_{\nu }(ɀ),\alpha \right)=I_{\nu }(0,\theta )e^{-{\tau }_{\nu }(ɀ)/\cos \alpha }+\left(1-e^{-{\tau }_{\nu }(ɀ)/\cos \alpha }\right)B_{\nu }\left(T_{\text{c}}\right),$$(C.2)

where B_ν (T_c) is the Planck function taken for the midplane temperature of disk, T_c, and α is the angle between the direction of propagating radiation and the normal to the disk surface. We note that the disk is assumed to be isothermal here. Since radiation from the host star is assumed to be a blackbody, the total frequency-integrated intensity at ɀ = −H, where τ_ν(ɀ) ≡ 0 for the incoming radiation, along the line of sight, intersecting the star photosphere simplifies and becomes $$I(0,\alpha )=\underset{0}{\overset{\infty }{{\displaystyle \int }}}{B}_{\nu }\left(T_2\right)d\nu =\frac{{\sigma }_{\text{B}}}{\pi }{T}_2^4,$$(C.3)

where T₂ is the effective temperature of the host star.

The first term of Eq. (C.2) becomes negligible during the following propagation of radiation through the disk, since in this case τ_ν(ɀ)/cos α ≫ 1 because the source of radiation is located at α ≈ π/2. On the contrary, for the second term the extraction by τ_ν (ɀ)/cos α ≪ 1 can be used since the disk itself irradiates mostly at α ≈ 0. Hence, Eq. (C.2) becomes $$I_{\nu }\left({\tau }_{\nu }(ɀ),\alpha \right)\approx \frac{{\tau }_{\nu }(ɀ)}{\cos \alpha }{B}_{\nu }\left(T_{\text{c}}\right).$$(C.4)

The total frequency-integrated intensity at ɀ = + H reads $$I(H,\alpha )=\underset{0}{\overset{\infty }{{\displaystyle \int }}}{I}_{\nu }\left({\tau }_{\nu },\alpha \right)\text{d}\nu =\underset{0}{\overset{\infty }{{\displaystyle \int }}}\frac{{\tau }_{\nu }}{\cos \alpha }{B}_{\nu }\left(T_{\text{c}}\right)\text{d}\nu =\frac{{\tau }_{\text{z},\text{P}}}{\cos \alpha }\frac{{\sigma }_{\text{B}}}{\pi }{T}_{\text{c}}^4,$$(C.5)

where τ_ɀ,P is the Planck mean opacity given by $${\tau }_{\text{z,P}}={\displaystyle \underset{0}{\overset{\infty }{\int }}{\tau }_{\nu }{B}_{\nu }\text{d}\nu }/{\displaystyle \underset{0}{\overset{\infty }{\int }}{B}_{\nu }\text{d}\nu }.$$(C.6)

Generally, the radiation flux is related to the intensity as $$F={{\displaystyle \int }}^{\text{}}I\cos \alpha \text{dΩ},$$(C.7)

where dΩ = sin θdθdϕ is an element of the solid angle (θ and ϕ are polar and azimuthal angles of the spherical coordinates, which we select in this case with ɀ-direction to the center of the star of the linked Cartesian system), α is the angle between the direction of propagating radiation and the normal to the disk surface. In the considered case, the angle α is defined as follows: cos α = sin θ cos ϕ. Then, for the upward flux, we obtain $$\begin{array}{l}{F}^{+}(H)={{\displaystyle \int }}^{\text{}}\text{\hspace{0.17em}d}\varphi {{\displaystyle \int }}^{\text{}}I(H,\alpha )\cos \alpha \sin \theta \text{d}\theta =\\ =\left({\varphi }_2-{\varphi }_1\right)\underset{{\theta }_2}{\overset{{\theta }_1}{{\displaystyle \int }}}\text{\hspace{0.17em}d}(\cos \theta )\frac{{\tau }_{\text{z},\text{P}}}{\cos \alpha }\frac{{\sigma }_{\text{B}}{T}_{\text{c}}^4}{\pi }\cos \alpha =\\ =\frac{{\varphi }_2-{\varphi }_1}{\pi }{\tau }_{\text{z},\text{P}}{\sigma }_{\text{B}}{T}_{\text{c}}^4\left(\cos {\theta }_1-\cos {\theta }_2\right)=2{\tau }_{\text{z},\text{P}}{\sigma }_{\text{B}}{T}_{\text{c}}^4,\end{array}$$(C.8)

where ϕ₁ = −π/2, ϕ₂ = π/2, θ₁ = 0 and θ₂ = π sets the range of angles when the disk radiates from the surface. Due to the symmetry of the disk, the downward intensity is $$I^{-}(H,\alpha )=-I^{+}(0,\alpha )=-\frac{{\sigma }_{\text{B}}{T}_2^4}{\pi },$$(C.9)

while the corresponding flux is $$\begin{array}{l}{F}^{-}(H)=\underset{{\varphi }_1^{\prime }}{\overset{{\varphi }_2^{\prime }}{{\displaystyle \int }}}\text{\hspace{0.17em}d}\varphi \underset{{\theta }_1^{\prime }}{\overset{{\theta }_2^{\prime }}{{\displaystyle \int }}}\text{\hspace{0.17em}d}\theta I^{-}(H,\alpha ){\sin }^2\theta \sin \varphi =\\ =-{\frac{{\sigma }_{\text{B}}{T}_2^4}{\pi }\left({\varphi }_2^{\prime }-{\varphi }_1^{\prime }\right)\frac{\theta -\sin \theta \cos \theta }{2}|}_{{\theta }_1^{\prime }}^{{\theta }_2^{\prime }}\approx -\frac{2}{3\pi }{\left(\frac{R_2}{r}\right)}^3{\sigma }_{\text{B}}{T}_2^4,\end{array}$$(C.10)