On the Wind-driven Relaxation Cycle in Accretion Disks

Various astrophysical systems including X-ray binaries, young stellar objects, cataclysmic variables (CVs) and active galactic nuclei (AGNs) are powered by accretion disk processes. The relatively high luminosity of these objects is due to the efficient conversion of accretion power into radiation. The luminosity generated from such accretion processes tend to vary with time which provides us with important clues as to the nature of the accretion disk, the accretor, and also the object or source supplying matter to the disk.

For a constant rate of mass supply, time variability could be attributed to spatial and temporal variations in the disk structure or the strength and configuration of the disk magnetic field. Such variations can arise due to the extended nature of the disks, with the inner and outer radii differing by orders of magnitude (ranging between two orders in CVs, to about seven orders in AGNs). The surface properties at small radii differ from those at large radii and there is a host of different physical processes that may cause the accreting material to undergo time-dependent evolution. They include a variety of local instabilities such as convective, thermal or magnetorotational instability (Balbus & Hawley 1998; Fromang & Lesur 2019). Yet there are several nonlocal processes that could also lead to time variability.

One of the consequences of the large radial extent is the huge variation in escape velocity throughout the disk. Thus, as long as a disk has a slightly concave surface, the high-energy radiation that is emitted by the inner disk can irradiate the outer disk leading to the formation of a high-temperature surface layer in which thermal speeds can exceed the escape velocity. This can drive a strong wind from outer radii and cause a disruption in the accretion flow. When information about this disruption reaches the inner disk, it changes the local emission which, in turn, affects the disk irradiation. Thus, the radiation from the inner disk acts as a coupling between the inner and outer disks. The self-irradiated disk is an example of a “self-regulated accretion” process (Shields et al. 1986, S86 hereafter) with feedback.

The Compton-heated corona and eventual disk wind (Begelman et al. 1983) have widespread applications in understanding the absorption lines observed in AGNs (Woods et al. 1996) and X-ray binaries (Proga & Kallman 2002; Luketic et al. 2010; Miller et al. 2015; Waters & Proga 2018). However, it is unclear as to what degree a disk wind can destabilize the accretion disk and be responsible for the observed variability in the luminosity and spectral energy distribution (SED). To assess this role of the disk winds, we may define a variable ${\eta }_{w}\equiv {\dot{M}}_{w}/{\dot{M}}_{a}$, where ${\dot{M}}_{w}$ and ${\dot{M}}_{a}$ are the wind mass loss rate (at the outer disk region) and mass accretion rate onto an accretor, respectively. The ratio measures the efficiency of wind driving due to the accretion power and indicates how strongly the wind is coupled to the latter. The model of instantaneous response of wind-to-accretion and vice versa showed that a wind with η_w as low as one, destabilizes the disk (Begelman et al. 1983). However, taking the effect of viscosity into account, S86 found that the accretion at the inner disk edge responded much slower to the change in the disk surface density, Σ, at large radii. Viscosity stabilizes the disk by producing a “delay” or “relaxation time” to the propagation of perturbations throughout the disk. Therefore, a much higher η_w was needed to generate variability in the disk.

For systems with relatively high luminosities, the escape velocity from a disk at a given radius could be reduced by the radiation pressure on free electrons and due to opacity from spectral lines and bound-free processes. In X-ray binaries, the latter two are negligible because the gas is highly ionized and few lines as well as few bound transitions are present. Yet, as expected, and shown both numerically and theoretically, the radiation pressure on free electrons introduces a linear scaling between the launching radius of the thermal wind and the luminosity, L (Proga & Kallman 2002). In AGNs and CVs, the radiation pressure on lines (line driving) can produce a wind with ${\dot{M}}_{w}\propto {L}^{1/\alpha }$, where α is the force multiplier parameter. The value of α depends on the SED, the chemical composition, and the physical conditions in the gas but it generally ranges from  0.2 to 0.8 (e.g., Castor et al. 1975). The above scaling is quite universal as it holds for 1D stellar wind (Castor et al. 1975) as well as for 2D (e.g., Proga et al. 1998; Proga 1999) and even 3D disk winds (Dyda & Proga 2018).

In this paper, we study a more generalized model of disk wind coupling to accretion power. We verify our results against the classic case of linear dependence of wind on accretion rate (S86) and then extend our analysis to nonlinear dependencies of the wind. This allows us to test our model against different possibilities. Our main focus is to study the dynamic variability of the launching radius of the wind and explore our model for different free parameters. The outline of the paper is as follows: in Section 2, we describe in brief the mathematical and computational techniques used in our analysis. In Section 3, we introduce our calculation and verification of the S86 result (with increased numerical resolutions). Our generalized approach toward the problem allows us to look at the results of two different models holistically and identify several previously unexplored cases. Finally, in Section 4, we summarize the applications of disk oscillations that have been studied in the past and what future prospects they might hold.

2.1. Equations and Analytical Results

We assume azimuthal symmetry and perform 1D simulations on a geometrically thin and optically thick disk along the radial direction. We have adopted similar formulations and notations as used in S86. Their approach involves solving the diffusion equation that describes the disk evolution (Lynden-Bell & Pringle 1974; Pringle 1981). As in S86, we assume a constant rate of mass injection at the outermost disk radius, R_d, from an external source. This is true for all the cases studied henceforth. The loss of mass in the form of wind takes place at the very same radius as for the first model, which corresponds to the model described by Equation (3.5) in S86 (see Section III.(b) in S86). Figure 1 shows a schematic of this model. We also assume a simple radius-dependent viscosity ν ∝ R, similar to that used in Lynden-Bell & Pringle (1974), for all our models. The mass continuity and angular momentum conservation gives us the following diffusion equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}=\displaystyle \frac{1}{2\pi R}\displaystyle \frac{\partial \dot{M}}{\partial R}+S(R,t),\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&\dot{M}=6\pi {R}^{1/2}\displaystyle \frac{\partial }{\partial R}(\nu {\rm{\Sigma }}{R}^{1/2})\qquad \mathrm{and},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&S(R,t)={S}_{\mathrm{in}}(R)-{S}_{w}(R,t).\end{eqnarray} \tag{ 3 }$$

Here $\dot{M}$ is the mass accretion rate at a given radius, ν is the kinematic viscosity, and Σ is the local surface mass density. The net source term, S, accounts for both a steady mass input to the disk (S_in), as well as a mass loss due to wind ejection from the disk (S_w). The general mass input and output rates, respectively, are defined by

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{in}}=2\pi {\int }_{{R}_{i}}^{{R}_{d}}{S}_{\mathrm{in}}(R)R\quad {dR}\qquad \mathrm{and},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\dot{M}}_{w}(\lt R)=2\pi {\int }_{{R}_{i}}^{R}{S}_{w}(R^{\prime} )R^{\prime} \quad {dR}^{\prime} ,\end{eqnarray} \tag{ 5 }$$

where R_i and R_d are the innermost and outermost disk radii, respectively, and the mass loss is calculated up to an arbitrary radius R within which mass is being lost. The above equations are recast using dimensionless variables as in Equations (2.5)–(2.10) in S86. We present these equations below for clarity of our method description,

$$\begin{eqnarray}&&{R}_{* }\equiv \quad R/{R}_{d},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\nu }_{* }\equiv \quad \nu /{\nu }_{0},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\dot{M}}_{* }\equiv \quad \dot{M}/{\dot{M}}_{\mathrm{in}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{* }\equiv \quad {\rm{\Sigma }}/{{\rm{\Sigma }}}_{0},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{t}_{* }\equiv \quad t/{t}_{0},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{S}_{* }\equiv \quad {{SR}}_{d}^{2}/{\dot{M}}_{\mathrm{in}},\end{eqnarray} \tag{ 11 }$$

where ν₀ is the characteristic viscosity, ${{\rm{\Sigma }}}_{0}\equiv {\dot{M}}_{\mathrm{in}}/{\nu }_{0}$ is the characteristic surface density and ${t}_{0}\equiv {R}_{d}^{2}/{\nu }_{0}$ is the characteristic viscous timescale. Using the coordinate transformation, $x={R}_{* }^{1/2}$, Equations (1) and (2) can be rewritten as,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{* }}{\partial {t}_{* }}=\displaystyle \frac{1}{4\pi {x}^{3}}\displaystyle \frac{\partial {\dot{M}}_{* }}{\partial x}+{S}_{* }(x,t),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\dot{M}}_{* }=3\pi \displaystyle \frac{\partial }{\partial x}(x{\nu }_{* }{{\rm{\Sigma }}}_{* }).\end{eqnarray} \tag{ 13 }$$

As in the wind model studied by S86, we treat mass loss using a delta function (launching of wind at a given radius), and S86 defined the main model parameter as

$$\begin{eqnarray}&&C\equiv {\dot{M}}_{w}/{\dot{M}}_{a},\end{eqnarray} \tag{ 14 }$$

where ${\dot{M}}_{a}\equiv \dot{M}({R}_{\mathrm{in}})$. In this paper, we alternatively refer to this ratio as the wind efficiency, η_w. While discussing or referring to the classic case studied in S86, we use the “C” notation for comparison. As demonstrated by S86, once the disk has attained a steady state, perturbations in disk surface density would either persist, grow, or decay with time. This forms the basis of the “self-regulated accretion” and η_w (or C) determines how strongly coupled the wind is to the central accretion rate and hence, to the luminosity of the central source. At steady state, the mass conservation relation applied to the disk, gives

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{in}}={\dot{M}}_{a}^{(s)}+{\dot{M}}_{w}^{(s)},\end{eqnarray} \tag{ 15 }$$

where the superscript (s) denotes the steady state value. The expression for steady state mass accretion rate normalized to the mass input rate ${\dot{M}}_{\mathrm{in}}$ reads in the following way,

$$\begin{eqnarray}&&{\dot{M}}_{a* }^{(s)}=\displaystyle \frac{1}{1+C}.\end{eqnarray} \tag{ 16 }$$

Using a straightforward radius-dependent viscosity law, S86 performed an analytical calculation to find the critical value of C that would lead to a perpetual oscillation in disk density about its steady state value. The amount of matter depleted through wind and accretion is continually replenished by the constant supply of matter. This results in a variability of mass accretion rate and luminosity at the inner edge of the disk. For the parameter C, S86 obtained the value required for sustained stable (critical) oscillations analytically as,

$$\begin{eqnarray}&&{C}_{\mathrm{crit}}=\cosh \pi \approx 11.6,\end{eqnarray} \tag{ 17 }$$

such that when C < C_crit, the oscillations decay, whereas for C > C_crit, they grow.

Zoom In Zoom Out Reset image size

The wind mass loss rate might not be a linear function of mass accretion rate. S86 considers such a possibility and they show based on their analytical treatment that for a power-law dependence, $C\propto {\dot{M}}_{a}^{{k}_{l}}$, where k_l is an arbitrary constant, the C_crit would be reduced by a factor of ${(1+{k}_{l})}^{-1}$ (see analogous Equation (4.5) in S86). The case with k_l = 0 corresponds to the case discussed above, where Equation (17) gives the value of C_crit.

We formally approach this possibility by generalizing our Equation (14) and writing it in the following way:

$$\begin{eqnarray}&&{\dot{M}}_{w* }\equiv C^{\prime} {\dot{M}}_{a* }^{p}.\end{eqnarray} \tag{ 18 }$$

where the subscript ∗ stands for mass loss rates normalized to the input mass rate ${\dot{M}}_{\mathrm{in}}$, where any remaining constant of proportionality has been absorbed into C′, and p is an arbitrary constant exponent. Note that Equation (18) reduces to Equation (14) when p = 1. The requirement for the steady state condition becomes

$$\begin{eqnarray}&&{\dot{M}}_{a* }^{(s)}+C^{\prime} {\left({\dot{M}}_{a* }^{(s)}\right)}^{p}-1=0,\end{eqnarray} \tag{ 19 }$$

but this equation needs to be solved numerically.

In the above analyses, a fixed wind launching radius has been assumed. However, in general, this may not be an ideal condition. Our main focus here is to examine the effects of relaxing this assumption by allowing the radiation from the inner disk to irradiate the entire disk and causing reduced local escape velocity. We can express luminosity L in units of the Eddington factor Γ, such that L = ΓL_Edd, where L_Edd is the Eddington luminosity. Here we assume that the irradiation luminosity L of the disk equals the total accretion luminosity. We express the coupling between the launching radius and ${\dot{M}}_{a}$ using the following expression:

$$\begin{eqnarray}&&{R}_{L* }=1-{\rm{\Gamma }}\displaystyle \frac{{\dot{M}}_{a* }}{{\dot{M}}_{a* }^{(s)}},\end{eqnarray} \tag{ 20 }$$

where R_L* is the launching radius normalized to R_d. Our Equation (20) is similar to Equation (22) in Proga & Kallman (2002) that was derived for the launching radius of a Compton-heated wind, corrected for radiation driving. When Γ = 0, R_L = R_d and we have the case studied by S86. The only numerical constraint is that R_L ≥ 0.

Using their analytical method for the simplest case (Γ = 0, p = 1), S86 derived an expression for the period of oscillations as

$$\begin{eqnarray}&&{P}_{* }\propto \displaystyle \frac{{R}_{L}}{{R}_{d}}.\end{eqnarray} \tag{ 21 }$$

We expect that in our variable R_L model, the disk stability condition and the variability period will be sensitive to Γ.

2.2. Numerical Methods

We have used numerical methods similar to that used in Bath & Pringle (1981). We have developed a Python code to study the effects of self-regulated accretion as discussed in the previous section. The wind launching zone is a delta function, with mass loss taking place from a single radial grid zone.

$$\begin{eqnarray}&&{S}_{w* }=\displaystyle \frac{C{\dot{M}}_{a* }}{\pi ({R}_{j* }^{{\prime} 2}-{R}_{j-1* }^{{\prime} 2})},\end{eqnarray} \tag{ 22 }$$

where ${R}_{j* }^{\prime} ={x}_{j}^{{\prime} 2}$ and ${x}_{j}^{\prime} =({x}_{j}+{x}_{j+1})/2$. Here R_j refers to the wind launching radius and ${R}_{j* }^{\prime}$ is the averaged launching radius value used to calculate S_w. The disk surface density Σ_* is assigned an initial value of 0. We take ${\dot{M}}_{\mathrm{in}}=1$, since all mass rates are normalized to ${\dot{M}}_{\mathrm{in}}$ and together with Equation (22), calculate the source term S_*. The diffusion Equation (12) is then solved using a forward difference method while updating Σ_*. The boundary condition is obtained by imposing the condition that mass flux is finite and conserved at the disk edges, i.e., $\partial (\nu {\rm{\Sigma }})/\partial R=0$.

Our resolution study shows that the calculated value of the accretion rate depends on the width of each radial zone. In particular, the variable Γ model is sensitive to smaller resolution in x, N_x (N_x < 200). Hence, for most of our models presented below, we use N_x = 200 (see Section 3 for more discussion).

To ensure that the value of Σ remains realistic at all times, it is quite common to impose the condition Σ_* = 0, whenever the numerical solution leads to Σ_* < 0. This condition is especially important when we deal with growing oscillations. Instead of allowing the disk to deplete completely, we introduce a floor value of 0.001 for Σ.

Three timescales that enter this problem are (1) the mass outflow timescale, t_out, (2) mass inflow timescale, t_in and, (3) diffusion timescales, t_diff. These scales can be defined as

$$\begin{eqnarray}&&{t}_{\mathrm{out}}=\dot{{\rm{\Sigma }}}({R}_{L})/{\dot{{\rm{\Sigma }}}}_{w}({R}_{L}),\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{in}}=\dot{{\rm{\Sigma }}}({R}_{\mathrm{in}})/{\dot{{\rm{\Sigma }}}}_{\mathrm{in}}({R}_{L})\,\mathrm{and},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{diff}}=0.25\displaystyle \frac{4{\rm{\Delta }}{x}^{2}}{3},\end{eqnarray} \tag{ 25 }$$

where Equation (25) expresses the stability criterion of Equation (12) which is a diffusion equation. In our numerical calculations, we choose a time step to be 20% of the shortest of the above three timescales.

The initial condition is to set the disk surface density to be 0 and let matter diffuse from the surrounding source until the accretion disk reaches a steady state. We then perturb it by switching on the wind. In practice, we allow the system to reach $(1-\epsilon ){\dot{M}}_{a* }^{(s)}$, where is a very small number of the order of 10⁻³. It takes roughly three viscous timescales to reach this near steady state. We use the fsolve method in SciPy’s¹ optimization library to find the root of Equation (19).

First we checked the results from our simulations against the results presented in S86. Our resolution is higher than theirs by a factor of 10. This increased resolution does not considerably affect the key result although it does substantially reduce the amplitude of oscillation. This can be attributed to the fact that we have limited our analysis to a single zone wind regardless of the radial resolution. Namely, we still inject and eject the same amount of matter but now from a radial ring with smaller area. We find that C_crit ≈ 11.3 as opposed to 11.4 obtained by S86. Figure 2 illustrates our results in a similar manner to Figure 1 by S86. As expected, we observe damped oscillation for C < 11.3 and growing oscillations for C > 11.3. We note that we consider only those oscillations to be stable whose amplitudes vary by less than 10⁻⁴. This condition is consistently followed throughout our analysis, and is used to classify the oscillations into categories.

Zoom In Zoom Out Reset image size

In the growing oscillation cases, the mass accretion rate grows until it saturates after some time. This saturation of oscillation is caused by the local surface density reaching negligible values or, in other words, complete depletion of matter from that region of the disk. We expect that if S86 continued their calculations to longer times, they would likely find the same behavior. For C = 14.4 shown above, our simulations showed that during this phase, ${{\rm{\Sigma }}}_{* ,\min }$ approaches 0, whereas ${{\rm{\Sigma }}}_{* ,\max }$ approaches 2.5 times its steady state value. The outcome of this instability is a small amplitude ${\dot{M}}_{a* }$ oscillation but a large modification of the disk solution.

Table 1 summarizes the combination of p and C′ (for Γ = 0) values used in our simulations. The wind to accretion ratio sets the stability of the disk. The efficiency factor η_w is a function of time and hence, we consider the value of η_w at the beginning of the perturbation. The value of η_w decreases with increasing p with a slope of about −1 on a log–log plot (Figure 3). This confirms the inverse relationship between ${\eta }_{w,\mathrm{crit}}$ and p, that we discussed in Section 2.1 and demonstrates that the disk is easily destabilized if the wind is more strongly coupled to accretion (i.e., higher p).

Zoom In Zoom Out Reset image size

Table 1.  Summary of all Parameters (C′ and p) Used in Simulations for the Model Described by Equation (18)

p	${C}_{d}^{{\prime} }$	${C}_{\mathrm{crit}}^{{\prime} }$	${C}_{g}^{{\prime} }$
0.5	4.5	4.95	5
2/3	6	7.05	7.5
3/4	8	8.6	9
4/5	9	9.5	9.9
1	11	11.3	14.4
5/4	18.5	20	24
4/3	21	22.8	24
3/2	20.5	26	30.5
7/4	37	39	40.5
2	30.5	48.5	50.5

Note. The subscripts d, crit and g denote values of C′ for which the oscillations decay, persist and grow, respectively. The value of ${C}_{\mathrm{crit}}^{{\prime} }$ increases with increasing p. In Figure 3, we use the corresponding ${C}_{\mathrm{crit}}^{{\prime} }$ to plot the wind efficiency η_w,crit, which decreases with increasing p.

Download table as:  ASCIITypeset image

Table 2 contains our parameter survey for the variable R_L model. The Γ factor strongly controls and alters the outcome of disk evolution. In particular, we identified new cases for high Γ (Γ ≥ 0.7). In low Γ cases (Γ ∼ 0.2), we observe some deviation from the classical cases of stable oscillations.

Table 2.  Summary of Parameter Survey

p = 0.5
Γ	${C}_{{sd}}^{{\prime} }$	${C}_{d}^{{\prime} }$	${C}_{\mathrm{crit}}^{{\prime} }$	${C}_{g}^{{\prime} }$	${C}_{{sg}}^{{\prime} }$
0.001	⋯	3.5 $\longrightarrow$ 4.7a	4.71	4.715 $\longrightarrow$ 4.75	4.8
0.01	⋯	$4\longrightarrow 4.71$	4.72	$4.8\longrightarrow 5$	⋯
0.05	⋯	$4.7\longrightarrow 4.75$	4.772	$4.8\longrightarrow 5$	⋯
0.1	⋯	$5\longrightarrow 5.05$	5.067	$5.1\longrightarrow 5.2$	⋯
0.2	$5.2\longrightarrow 5.5$	$5.8\longrightarrow 6$	$6.2\longrightarrow 12$	⋯
0.3	⋯	$4\longrightarrow 5.45$	$5.48\longrightarrow 6.2$	$6.5\longrightarrow 9$	$10\longrightarrow 15$
0.4	$4\longrightarrow 4.417$	⋯	4.418	⋯	$4.419\longrightarrow 5$
0.5	$3.4\longrightarrow 3.6$	⋯b	3.63	⋯	$3.7\longrightarrow 4$
p = 1
Γ	${C}_{{sd}}^{{\prime} }$	${C}_{d}^{{\prime} }$	${C}_{\mathrm{crit}}^{{\prime} }$	${C}_{g}^{{\prime} }$	${C}_{{sg}}^{{\prime} }$
0.001	⋯	11.4	11.6	⋯	11.65
0.01	⋯	11.6	11.65	11.7	⋯
0.05	⋯	11.8	$11.87\longrightarrow 11.9$	11.95	⋯
0.1	⋯	12.6	12.65	12.7	⋯
0.2	⋯	16.5	16.55	16.6	⋯
0.3	18	⋯	$20\longrightarrow 32$	⋯	45
0.4	15	⋯	15.93	⋯	$20\longrightarrow 25$
0.5	⋯	$11.8\longrightarrow 11.88$	11.89	⋯	$11.9\longrightarrow 11.95$
p = 2
Γ	${C}_{{sd}}^{{\prime} }$	${C}_{d}^{{\prime} }$	${C}_{\mathrm{crit}}^{{\prime} }$	${C}_{g}^{{\prime} }$	${C}_{{sg}}^{{\prime} }$
0.001	⋯	39.4	$35\longrightarrow 39.3$	$39.5\longrightarrow 40$	50
0.01	⋯	$11\longrightarrow 39.5$	39.7	40	⋯
0.05	⋯	$39.7\longrightarrow 41$	41.3	$41.5\longrightarrow 42$	⋯
0.1	⋯	$42\longrightarrow 43.5$	43.7	44	⋯
0.2	⋯	$43\longrightarrow 55$	60.35	$62\longrightarrow 65$	⋯
0.3	⋯	$60\longrightarrow 100$	107.8	$110\longrightarrow 120$	⋯
0.4	⋯	110	$130\longrightarrow 300$	500	⋯
0.5	78	$80\longrightarrow 83$	83.7	85	100

Notes. ${C}_{d}^{\prime} ,{C}_{\mathrm{crit}}^{{\prime} },{C}_{g}^{\prime}$ hold the same meaning as described in the Table 1. The two new types of outcome, ${C}_{{sd}}^{\prime}$ and ${C}_{{sg}}^{\prime}$, denote oscillations that decay/grow for a few timescales before stabilizing and oscillating with a constant amplitude. In Figures 6 and 7, we use the corresponding ${C}_{\mathrm{crit}}^{{\prime} }$ to plot the wind efficiency ${\eta }_{w,\mathrm{crit}}.$

^aThe $\longrightarrow$ denotes the range of tested values that fall under a certain category. The limits of the range are not absolute but give a more or less general idea of the behavior within those C′ values. ^bEmpty cells denote parameter combination for which we do not find the outcome in question. This does not rule out the possibility of such a case or of any other new cases. ${C}_{\mathrm{crit}}^{{\prime} }$ is the only value that was important for our analysis and the parameter survey was focused on finding the condition for stable (or critical) oscillations.

Download table as:  ASCIITypeset image

We limit our presentation to three cases, p = 0.5, 1, and 2. The case p = 1 corresponds to the special case studied in S86. The other two cases are representative of the lowest and highest p values considered (see Table 1). For 0.001 ≤ Γ ≤ 0.1, the disk behavior does not deviate much from the classical behavior depicted in Figure 2. However, for Γ ≥ 0.2, we find some new results. We describe these cases in detail below.

• 1.  
  Γ  = 0.2 and 0.3: We observe more than one ${C}_{\mathrm{crit}}^{{\prime} }$ value resulting in stable oscillations. For example, constant amplitude oscillations occur when ${C}_{\mathrm{crit}}^{{\prime} }$ is between 6.2 and 12. The lowermost panel of Figure 4 shows how the sinusoidal nature of the stable oscillations starts to change for C′ ≳ 13. We also find that for increasing C′, the frequency of oscillation starts to increase. For Γ = 0.3, the results are similar to those of Γ = 0.2.
• 2.  
  Γ = 0.4: For p = 0.5 and p = 1, we find a single value of ${C}_{\mathrm{crit}}^{{\prime} }$. For p = 2, we see that ${C}_{\mathrm{crit}}^{{\prime} }$ lies between 130 and 300. The range of ${C}_{\mathrm{crit}}^{{\prime} }$ and the critical value of $C^{\prime}$ rise considerably.
• 3.  
  ${\rm{\Gamma }}=0.5$: For all p values, we obtain single-valued ${C}_{\mathrm{crit}}^{{\prime} }$. We note that ${C}_{\mathrm{crit}}^{{\prime} }$ for Γ = 0.5 is less than that for Γ = 0.4. This result holds true for all three p values investigated.
• 4.  
  Γ = 0.7: Figure 5 summarizes the distinct cases obtained for p = 1. We see that up to a certain value of C′, the accretion rate initially increases and at the same time oscillates with relatively high frequency and small amplitude. In the case of C′ = 6.5, we see that eventually the fundamental oscillation dominates and steadies around 1. As C′ increases, we see a growing oscillation that eventually saturates.
• 5.  
  Γ = 0.9: The behavior is similar to the Γ = 0.7 case. For higher C′, the frequency increases rapidly with time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The amplitude of oscillations in ${\dot{M}}_{a* }$ remain constrained to ∼50% of its steady state value in all our simulations. This can be mostly attributed to the amount of matter available to be launched as wind and also the constraint on the lowest possible launching radius. In all of the cases studied, the amplitude of oscillation increases for increasing η_w.

As stated in Section 2.2, we have chosen a spatial resolution of N_x = 200 for the above cases. Our resolution study showed that there are quantitative changes in most cases as well as qualitative changes in some of the extreme cases. For example, for $p=1,C^{\prime} =6,{\rm{\Gamma }}=0.7$, the oscillations cease earlier for N_x = 400, while for N_x = 800, they persist for a longer number of time steps. In addition, the amplitude of oscillations also decrease with increasing resolution. For higher N_x, there is a clear tendency toward convergence. This resolution study was conducted for several other cases in our parameter survey. Convergence was evident for most cases with Γ ≤ 0.7. For higher Γ, the behavior was more erratic and unpredictable for different N_x. We restricted our N_x to 200 despite this fact, since the nature of oscillations and the feedback on ${\dot{M}}_{w}$ remained unaffected from a qualitative perspective.

In Figure 6, we show the steady state critical η_w–p relation for various Γ. We find that the slope of this relation is nearly constant for Γ ≤ 0.2. For higher Γ, the slope changes and more than one ${\eta }_{w,\mathrm{crit}}$ exists. For such cases we have shaded the region between all the possible straight line fits. For Γ > 0.5, oscillations are distinctly different from that of the classical cases and we cannot group them under simple categories (see above and Figure 5). We do not plot these points in Figure 6 or list them in our table of classification of disk oscillations.

Zoom In Zoom Out Reset image size

To visualize our results for the wind efficiency in a different way, we also plot steady state critical η_w as a function of Γ for different p values (the upper panel of Figure 7). The curves for different p values generally resemble each other. For Γ ≤ 0.2, ${\eta }_{w,\mathrm{crit}}$ increases with increasing Γ. For higher Γ, η_w decreases with increasing Γ, in all three p cases. To more directly compare the results for various p, in the bottom panel of Figure 7, we plot $p{\eta }_{w,\mathrm{crit}}$ versus Γ. For Γ < 0.1, we see that $p{\eta }_{w,\mathrm{crit}}$ is nearly constant, which is what we concluded from Figure 3.

Zoom In Zoom Out Reset image size

The study of disk winds in the context of state transition has evoked a lot of interest over the past decades (e.g., Fender et al. 2005; Körding et al. 2006). The wind launching mechanism and location may lead to time variability, the effects of which could be coupled to the signatures of mass accretion by the accretor. We explore this aspect through a self-regulated accretion disk. We found that our model is unlikely to explain state transitions in accretion disk spectra. However, it could be responsible for persistent small amplitude regular single-mode oscillations in the mass accretion rate. Admittedly, our treatment is very simple. We use one viscosity law in all our analyses, ν_* = R_* = x². We did not consider thermodynamic effects. We also did not incorporate magnetic fields in our analysis. Moreover, we assumed a δ-function model for the wind, which exaggerates the role of mass removal.

Despite the simplifications, this analysis might be a stepping stone toward developing models that can account for the existence of state transition signatures in accretion disks. Time variability in accretion rates may be explained by physical processes considered here. They are rich in features and may hold the key to understanding the coupling forces operating in an accretion disk. For example, the study of GRS 1915+105 by Neilsen et al. (2011) is a very detailed analysis in that direction. They nicknamed the oscillations in X-ray spectra as the “heartbeat” state and conducted a study of the geometry of the accretion disk using the X-ray continuum and emission lines in the optical spectrum. They demonstrate a strong correlation between mass loss in the form of wind and oscillations in the accretion rate that would explain the long-term effects in the disk.

A study of the same source, GRS 1915+105, by Zoghbi et al. (2016) considered in detail the reflection spectrum during the oscillatory phase of the source. Their calculations indicate winds being launched from very small disk radius which remains unchanged during this phase. Our variable R_L model allows the launching of wind from as close as the innermost disk region for high Γ cases. The time evolution of the accretion rate from our simulation does not indicate any distinct spike or sharp flare, which leads us to conclude that the system does not produce outbursts. Higginbottom et al. (2017) used their photoionization modeling of the SEDs to argue that thermally driven winds may hold the key to explaining the state changes in such systems. On the other hand, Neilsen (2013) goes on to demonstrate how heavy outflows not only quench the disk, thus affecting the formation of jets and causing state transition, but also influence the further production of winds. Disk variability and the state changes could also be caused or affected by instabilities of the disk itself. For instance, Janiuk et al. (2002) study the effect of radiation instabilities leading to limit-cycle behavior and modulations in accretion rates. As opposed to our model, this radiation-driven instability may result in sharp spikes indicating a high outflow from the system.

There are a number of studies related to fast outflows with high mass loss rates leading to state transitions and subsequent detection of jets in the system (Neilsen & Lee 2009; King et al. 2013; Gatuzz et al. 2019). A similar situation has been discussed in Casares et al. (2019), where a particular case of V404 Cyg indicates the presence of massive outflows, almost 2 orders of magnitude higher than central mass accretion rate ${\dot{M}}_{a}$. They speculated that these outflows are produced by radiation-driven winds coupled with classical thermal winds. This work stands out as providing direct observational evidence of powerful outflows leading to a quenching of accretion. The optical H_α line profile clearly indicates disk contraction following the massive outflow phase, consistent with what we would expect happens when irradiation is reduced. Our analysis of the Σ radial profile showed that Σ approaching 0 is responsible for saturating the oscillations.

Recently, Tomaru et al. (2019) showed that the detection of a blueshifted line from a black hole binary source H1743-322 strongly suggests a thermal-radiative disk wind. Their work indicates a disappearing wind in the hard state which could be attributed to the shadowing of the outer disk region by the inner corona. They also went on to state that the absorption features in other black hole sources such as GRS 1915+105 and GRO J1655-40 are most likely due to thermal-radiative winds as opposed to previously speculated magnetic effects. Another recent paper (Dubus et al. 2019) studies the effect of thermal-viscous instability on the light curves and stability diagrams associated with black hole binary systems. Additionally, they consider a fraction of the X-ray irradiation to be scattered by the wind and partially impinge on the outer disk regions. This has a stabilizing effect and can explain the shortened outbursts but cannot explain the rapid decay of outbursts. They studied a particular BHXB, GRO J1655-40, for which their model was able to reproduce the observed features of the light curve. They speculate that magnetic fields would need to be considered for a more promising explanation for the outbursts.

These studies indicate several possible aspects of our simplistic approach toward the study of wind-accretion coupling. The criteria for instability derived by S86 is often invoked when discussing consequences of observed or model disk winds (e.g., Luketic et al. 2010). The main conclusion of our work is that upon satisfying this criteria, a disk wind might not be responsible for large scale variations in luminosity because instability saturates at a relatively low level in terms of ${\dot{M}}_{a}$. However, it could result in a harder-to-detect significant reduction of Σ at large radii.

This work was supported by NASA under ATP grant 80NSSC18K1011. We thank Drs. Tim Waters and Rebecca Martin for their valuable comments and discussions.