Evolution of a Mode of Oscillation within Turbulent Accretion Disks

Consider ideal Newtonian hydrodynamics (Thorne & Blandford 2017). This is a useful first approximation in the following exploratory analysis of the interaction of a normal mode of oscillation with turbulence. Although the initial application will be to (geometrically thin) black hole accretion disks, the effects of general relativity should not change the general nature of our results (due to the Principle of Equivalence, within the small volume of the mode). The effects of the magnetic fields within such disks on the mode are unclear and will be discussed in the final section.

The changes that we find in the energy of the mode are produced mainly by the turbulent driving force. This amplification may be relevant to the high-frequency quasi-periodic oscillations (HFQPOs) observed in some of the stellar-mass to supermassive black hole (Figure 1) sources (Smith et al. 2018, 2021). They typically have stable frequencies but low duty cycles. The opposite is the case for the black hole low-frequency QPOs and most of the QPOs of accreting neutron stars (Remillard & McClintock 2006).

Zoom In Zoom Out Reset image size

Our major focus is the evolution of the amplitude of the mode. Previous investigations of oscillators subject to stochastic damping (Grue & Øksendal 1997; Ortega-Rodríguez et al. 2020) have mainly considered time averages and took the damping to have a positive mean value.

We take the background model of the accretion disk to be composed of two components: (a) the stationary and axisymmetric thin disk, and (b) the contribution of subsonic turbulence (continuously generated by the strong magnetorotational instability, which transfers the large free energy of the rate of shear (rdΩ/dr) in the disk into the turbulent cascades (Beckwith et al. 2011)). Thus, the background velocity, pressure, and mass density are decomposed as

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{u}} & = & r{\rm{\Omega }}(r){{\boldsymbol{e}}}_{\phi }+{\boldsymbol{v}}({turb}.),p={p}_{0}(r,z)+{p}_{1}({turb}.),\\ \rho & = & {\rho }_{0}(r,z)+{\rho }_{1}({turb}.),\end{array}\end{eqnarray} \tag{ 1 }$$

with ∣ ∇ · v ∣ ≪ ∣ ∇ × v ∣ and the turbulent Mach number ${ \mathcal M }\equiv v/{c}_{s}\lt 1$, where c_s is the speed of sound.

We consider an adiabatic perturbation of this background from a single normal mode and ignore mode−mode interactions. We also assume that the turbulence is unaffected by the mode. Then, the displacement vector of the mode is taken to be of the form

$$\begin{eqnarray}&&{\boldsymbol{\xi }}=A(t)\hat{{\boldsymbol{\xi }}}({\boldsymbol{x}}),\end{eqnarray} \tag{ 2 }$$

where the eigenfunction $\hat{{\boldsymbol{\xi }}}$ (taken to be axisymmetric, a property of many of the most observable modes) is approximated as that calculated from the zero-order steady background model, and the amplitude A is dimensionless. The density fluctuation produced by the mode is Δρ/ρ ≅ −∇ · ξ ∼ A. The perturbation theory that we shall now employ therefore requires that A² ≲ 1.

To analyze the dynamics of the mode, we employ the approach of Schenk et al. (2001) (especially Appendix I). With d/dt = ∂/∂t + u · ∇, and allowing for the time dependence of the unperturbed background model, Equations (I24) and (I25) of Schenk et al. (2001) give the components

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{\xi }_{i}}{{{dt}}^{2}}=\displaystyle \frac{{\partial }^{2}{\xi }_{i}}{\partial {t}^{2}}+2{\boldsymbol{u}}\cdot {\rm{\nabla }}\displaystyle \frac{\partial {\xi }_{i}}{\partial t}+\left(\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}\cdot {\rm{\nabla }}\right){\xi }_{i}+{\boldsymbol{u}}\cdot {\rm{\nabla }}({\boldsymbol{u}}\cdot {\rm{\nabla }}{\xi }_{i})\end{eqnarray} \tag{ 3 }$$

of the acceleration of the perturbation. We invoke the Cowling approximation (neglecting perturbations of the gravitational potential Φ, because the mass of the disk is much less than that of the black hole) and keep the first- and second-order perturbations of the pressure gradient restoring force. In addition, we add to the equation of motion (I37 of Schenk et al. 2001) the main external driving force per unit volume, the divergence of the turbulent Reynolds stress tensor.

This then gives the full equation of motion

$$\begin{eqnarray}\begin{array}{rcl}\rho {d}^{2}{\xi }_{i}/{{dt}}^{2} & = & -({\rm{\nabla }}\cdot {\boldsymbol{\xi }}){{\rm{\nabla }}}_{i}p+{{\rm{\nabla }}}_{i}({\boldsymbol{\xi }}\cdot {\rm{\nabla }})p+{{\rm{\nabla }}}_{i}(p{{\rm{\Gamma }}}_{1}{\rm{\nabla }}\cdot {\boldsymbol{\xi }})\\ & & -\rho ({\boldsymbol{\xi }}\cdot {\rm{\nabla }}){{\rm{\nabla }}}_{i}{\rm{\Phi }}(r,z)-{{\rm{\nabla }}}^{j}({{pQ}}_{{ij}})\\ & & -(1/2)\rho {\xi }^{j}{\xi }^{k}{{\rm{\nabla }}}_{j}{{\rm{\nabla }}}_{k}{{\rm{\nabla }}}_{i}{\rm{\Phi }}(r,z)-{{\rm{\nabla }}}^{j}(\rho {v}_{i}{v}_{j}),\end{array}\end{eqnarray} \tag{ 4 }$$

where Γ₁ is the adiabatic index of the perturbed pressure. The quantity ${Q}_{j}^{i}={k}_{1}{{\rm{\Theta }}}_{j}^{i}+{k}_{2}{\rm{\Theta }}{\delta }_{j}^{i}+{k}_{3}{{\rm{\Xi }}}_{j}^{i}+{k}_{4}{\rm{\Xi }}{\delta }_{j}^{i}$, where ∣k_n ∣ ∼1, with ${{\rm{\Theta }}}_{j}^{i}={{\rm{\nabla }}}_{k}{\xi }^{k}{{\rm{\nabla }}}_{j}{\xi }^{i}$ and ${{\rm{\Xi }}}_{j}^{i}={{\rm{\nabla }}}_{k}{\xi }^{i}{{\rm{\nabla }}}_{j}{\xi }^{k}$ .

Operating with $\int {d}^{3}x\hat{{\xi }^{i}}$ on the above equation then gives

$$\begin{eqnarray}&&\begin{array}{l}{c}_{0}{d}^{2}A/{{dt}}^{2}+{f}_{1}(t){dA}/{dt}+[{g}_{0}+{g}_{1}(t)]A\\ \quad +[{n}_{0}+{n}_{1}(t)]{A}^{2}={f}_{2}(t),\end{array}\end{eqnarray} \tag{ 5 }$$

with the time-dependent functions generated by the turbulence. (Interactions with other modes would also contribute to the term g₁(t), proportional to their amplitudes.) The constant ${c}_{0}=\int {d}^{3}x\rho {\hat{\xi }}^{i}{\hat{\xi }}_{i}$ and the constant ${g}_{0}={c}_{0}{\omega }_{0}^{2}$, where ω₀(∼ Ω) is the eigenfrequency of the mode in the absence of turbulence and the nonlinearity in the restoring force. Employing the phase τ ≡ ω₀ t and dividing the above equation by c₀ gives our master equation

$$\begin{eqnarray}&&\begin{array}{l}{d}^{2}A/d{\tau }^{2}+{F}_{1}(\tau ){dA}/d\tau +[1+{G}_{1}(\tau )]A\\ \quad +[{N}_{0}+{N}_{1}(\tau )]{A}^{2}={F}_{2}(\tau ),\end{array}\end{eqnarray} \tag{ 6 }$$

in which all quantities are dimensionless. From Equation (4), it is seen that $| {N}_{0}| \sim | {{\rm{\nabla }}}_{i}{\hat{\xi }}^{j}| \sim 1$.

Multiplying Equation (6) by dA/d τ gives

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dE}}{d\tau } & = & {F}_{2}(\tau )\displaystyle \frac{{dA}}{d\tau }-{F}_{1}(\tau ){\left(\displaystyle \frac{{dA}}{d\tau }\right)}^{2}\\ & & -\displaystyle \frac{1}{2}{G}_{1}(\tau )\displaystyle \frac{{{dA}}^{2}}{d\tau }-\displaystyle \frac{1}{3}{N}_{1}(\tau )\displaystyle \frac{{{dA}}^{3}}{d\tau }.\end{array}\end{eqnarray} \tag{ 7 }$$

We have introduced the (dimensionless) energy

$$\begin{eqnarray}&&E=\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{dA}}{d\tau }\right)}^{2}+\displaystyle \frac{1}{2}{A}^{2}+\displaystyle \frac{1}{3}{N}_{0}{A}^{3}\end{eqnarray} \tag{ 8 }$$

of the mode.

We find that

$$\begin{eqnarray}\begin{array}{rcl}{F}_{1} & = & \displaystyle \frac{2}{{c}_{0}{\omega }_{0}}\int {d}^{3}x\rho {\hat{\xi }}^{i}({\boldsymbol{u}}\cdot {\rm{\nabla }}){\hat{\xi }}_{i}\\ & = & \displaystyle \frac{1}{{c}_{0}}\displaystyle \frac{d}{d\tau }\int {d}^{3}x{\hat{\xi }}^{2}{\rho }_{1}\equiv {dH}/d\tau .\end{array}\end{eqnarray} \tag{ 9 }$$

To obtain the final expression, we have employed an integration by parts and the conservation of mass (∂ρ/∂t + ∇ · (ρ u ) = 0). We see that F₁ is proportional to the rate of change of the mass within the mode, so a decreasing mass produces growth of the mode. However, its long time average 〈F₁〉 = 0. Therefore, we now include a phenomenological source of damping, so that

$$\begin{eqnarray}&&{F}_{1}={dH}/d\tau +{D}_{* },\end{eqnarray} \tag{ 10 }$$

with the constant D_* > 0.

The function g₁(t) is generated by the last two terms in Equation (3) and the first three terms on the right-hand side of Equation (4) (which are larger than the fourth). The function n₁(t) is generated by the term in Equation (4) involving Qⁱ _j . Therefore,

$$\begin{eqnarray}&&| {G}_{1}| \sim | {N}_{1}| \sim | {p}_{1}| /{p}_{0}.\end{eqnarray} \tag{ 11 }$$

Finally, we obtain

$$\begin{eqnarray}\begin{array}{rcl}{F}_{2} & = & \displaystyle \frac{1}{{c}_{0}{\omega }_{0}^{2}}\int {d}^{3}x{\hat{\xi }}^{i}{{ \mathcal F }}_{i}\\ & = & \displaystyle \frac{1}{{c}_{0}{\omega }_{0}^{2}}\int {d}^{3}x{\rho }_{0}{v}^{i}({\boldsymbol{v}}\cdot {\rm{\nabla }}){\hat{\xi }}_{i},\end{array}\end{eqnarray} \tag{ 12 }$$

where ${{ \mathcal F }}_{i}=-{{\rm{\nabla }}}^{j}({\rho }_{0}{v}_{i}{v}_{j})$ is the turbulent driving force per unit volume.

Consider briefly the case when A² ≪ 1, so that we can neglect the nonlinear term in Equation (6). Employing the change of variable $A(\tau )=a(\tau )\exp [-0.5{\int }_{{\tau }_{0}}^{\tau }{F}_{1}(s){ds}]$ and the upper limits on the magnitudes of the stochastic functions (discussed in the next section) when ${ \mathcal M }\ll 1$ (giving $| {G}_{1}+(1/2){{dF}}_{1}/d\tau +(1/4){F}_{1}^{2}| \ll 1$), we obtain ${d}^{2}a/d{\tau }^{2}+a={F}_{2}(\tau )\exp [0.5{\int }_{{\tau }_{0}}^{\tau }{F}_{1}(s){ds}]$ from Equation (6). Employing the relevant Green's function, we then obtain the solution

$$\begin{eqnarray}\begin{array}{rcl}A(\tau ) & = & \cos \tau \left\{{A}_{0}\exp [-\tfrac{1}{2}{\int }_{0}^{\tau }{F}_{1}(s){ds}]\right.\\ & & \left.-{\int }_{0}^{\tau }\sin \tau ^{\prime} {F}_{2}(\tau ^{\prime} )\exp [-\tfrac{1}{2}{\int }_{\tau ^{\prime} }^{\tau }{F}_{1}(s){ds}]d\tau ^{\prime} \right\}\\ & & +\sin \tau \{{\dot{A}}_{0}\exp [-\tfrac{1}{2}{\int }_{0}^{\tau }{F}_{1}(s){ds}]\\ & & +{\int }_{0}^{\tau }\cos \tau ^{\prime} {F}_{2}(\tau ^{\prime} )\exp [-\tfrac{1}{2}{\int }_{\tau ^{\prime} }^{\tau }{F}_{1}(s){ds}]d\tau ^{\prime} \},\end{array}\end{eqnarray} \tag{ 13 }$$

corresponding to the initial conditions A = A₀ and ${dA}/d\tau \,={\dot{A}}_{0}$, at τ = τ₀ = 0.

Let us consider the evolution of A(τ) at late times (〈F₁〉τ ≫ 1), with $\langle {F}_{1}\rangle \equiv {(\tau -\tau ^{\prime} )}^{-1}{\int }_{\tau ^{\prime} }^{\tau }{F}_{1}(s){ds}\cong {D}_{* }\ll 1$. From Equation (13) we then obtain (until effects of nonlinearity become important)

$$\begin{eqnarray}\begin{array}{rcl}A(\tau ) & \approx & {\int }_{\tau -{\tau }_{* }}^{\tau }{F}_{2}(\tau ^{\prime} )\sin (\tau -\tau ^{\prime} )d\tau ^{\prime} \\ & = & {\int }_{0}^{{\tau }_{* }}{F}_{2}(\tau -x)\sin (x){dx},\end{array}\end{eqnarray} \tag{ 14 }$$

where τ_* ≡ 1/D_* and τ ≫ τ_* ≫ 1.

Including the mild nonlinearity, the evolution of the energy E(τ) of the mode is governed by Equation (7). From Equation (8), the potential energy is PE = A²/2 + N₀ A³/3. ∣A∣ becomes unbounded if the energy reaches $1/(6{N}_{0}^{2})$, which occurs when A = − 1/N₀, as shown in Figure 2. Will equilibrium be achieved before this blowout can occur? (Of course, "blowout" only indicates that nonlinearities have become important, so the evolution cannot be accurately continued. The energy would remain bounded if the coefficient of an A⁴ term was positive.) Our observational predictions (in Section 4) are all based upon choices of parameters that do not produce a blowout.

Zoom In Zoom Out Reset image size

Averaging over a time interval δ τ (with 2π ≪ δ τ ≪ τ), Equation (7) gives

$$\begin{eqnarray}&&\langle {dE}/d\tau \rangle \cong \langle {F}_{2}({dA}/d\tau )\rangle -{D}_{* }\langle {\left({dA}/d\tau \right)}^{2}\rangle .\end{eqnarray} \tag{ 15 }$$

(We have found that the contributions of G₁ and N₁ are negligible). At blowout, $\langle {({dA}/d\tau )}^{2}\rangle \cong \langle {A}^{2}\rangle \cong 1/(2{N}_{0}^{2})$. Initially, the driving term will be larger than the damping term, because it is only proportional to A. Therefore, the condition for equilibrium to be achieved before a blowout could occur is

$$\begin{eqnarray}&&\langle {F}_{2}{dA}/d\tau \rangle \lesssim {D}_{* }/(2{N}_{0}^{2}).\end{eqnarray} \tag{ 16 }$$

We next consider how to characterize the turbulence. The ensemble and time averages of our four stochastic functions (H, F₂, G₁, N₁) should vanish. The major properties of a turbulent eddy of radius ℓ are its velocity v_ℓ and turnover time t_ℓ ∼ 2π ℓ/v_ℓ . The smallest dimension (radius) of the largest turbulent eddy will be of order the (half) thickness h(r) of the accretion disk. Vertical force balance gives hΩ ∼ c_s .

We shall consider a normal mode with a similar thickness (such as the fundamental g-mode (also called r-mode)), with period t_m ∼ 2π/Ω ∼ t_L (the period of the largest eddy). Goldreich & Kumar (1988) found that an acoustic mode in the Sun couples most strongly to the eddies of the convective turbulence with t_ℓ ≲ t_m , giving ${\ell }\lesssim { \mathcal M }h$ , where ${ \mathcal M }$ is now the Mach number of the largest eddies (radius L ∼ h). We shall assume that the same holds true for the modes in our more strongly turbulent accretion disks (Nowak & Wagoner 1995).

The pressure and density fluctuations produced by such eddies are

$$\begin{eqnarray}&&{p}_{1}({\boldsymbol{x}},t)\sim {\rho }_{0}{v}^{2},{\rho }_{1}({\boldsymbol{x}},t)\sim {\rho }_{0}{\left(v/{c}_{s}\right)}^{2}.\end{eqnarray} \tag{ 17 }$$

We assume that these fluctuations arise mainly from the conversion of the eddies into acoustic modes (Goldreich & Kumar 1988) and take ∣ ∇ ∣ ≲ 1/h. A limit on the turbulent velocity comes from the effects of a "turbulent viscosity" on the structure and evolution of accretion disks (Kato et al. 2008). Beckwith et al. (2011) find that ${ \mathcal M }\approx 0.06\mbox{--}0.10$.

We generate the (now dimensionless) eddy lifetimes Δτ (assumed equal to the turnover time τ_ℓ ) by random sampling a distribution function with mean 〈Δτ〉 ≲ π, and we shall assume no delay between eddies. (Note that the unperturbed period of the mode is τ_m = 2π.) We choose a Rayleigh distribution $f{({\rm{\Delta }}\tau )=({\rm{\Delta }}\tau /{\sigma }_{* }^{2})\exp [-({\rm{\Delta }}\tau )}^{2}/(2{\sigma }_{* }^{2})]$. Its mean $\langle {\rm{\Delta }}\tau \rangle \,={\sigma }_{* }\sqrt{\pi /2}$, and its variance is $(2-\pi /2){\sigma }_{* }^{2}$.

From the above equations, it is seen that

$$\begin{eqnarray}&&| H| \sim | {G}_{1}| \sim | {N}_{1}| \sim | {F}_{2}| \sim \eta {{ \mathcal M }}^{2},\end{eqnarray} \tag{ 18 }$$

where η is the ratio of the volume of the dominant correlated eddies within the mode to the volume of the mode. Thus, because $\eta {{ \mathcal M }}^{2}\ll 1$, we neglect G₁ and N₁ (compared to unity) in our master Equation (6). Also, if F₁ is given by Equation (10), we again see that H(τ) does not affect the long-term solution (Equation (13)) to our linear equation, because ${\int }_{{\tau }_{1}}^{{\tau }_{2}}{F}_{1}(s){ds}\,=[H({\tau }_{2})-H({\tau }_{1})]+{D}_{* }({\tau }_{2}-{\tau }_{1})\cong {D}_{* }({\tau }_{2}-{\tau }_{1})$.

We approximate the functions H(τ) and F₂(τ) as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Psi }}}_{k}(\tau ) & = & 16{Y}_{k,n}{x}^{2}{\left(1-x\right)}^{2},\\ x & \equiv & (\tau -{\tau }_{n})/{\rm{\Delta }}{\tau }_{n}\ \ (k=1-2,0\leqslant x\leqslant 1)\end{array}\end{eqnarray} \tag{ 19 }$$

during each eddy lifetime Δτ_n . Note that ∣Y_k,n ∣ is the maximum value of ∣Ψ_k ∣ in eddy n, and the value and first derivative of Ψ_k vanish at the beginning and end of the eddy.

The values of Y_k,n are generated by random sampling a Gaussian–Markov conditional probability function $P({Y}_{k,n},{\bar{\tau }}_{n}| {Y}_{k,n-1})$, where ${\bar{\tau }}_{n}=({\rm{\Delta }}{\tau }_{n}+{\rm{\Delta }}{\tau }_{n-1})/2$. Its mean $\langle {Y}_{k,n}\rangle =\bar{Y}+{e}^{-{\bar{\tau }}_{n}/{\tau }_{r}}({Y}_{k,n-1}-\bar{Y})$, and its variance ${\sigma }_{k,n}^{2}\,=(1-{e}^{-2{\bar{\tau }}_{n}/{\tau }_{r}}){\sigma }_{k}^{2}$. We choose the equilibrium mean $\bar{Y}=0$, consistent with our assumption that 〈Ψ_k 〉 = 0 for averaging times δ τ ≫ Δτ. We choose the relaxation time τ_r = K〈Δτ〉, with the value of K allowing for correlations between subsequent eddies. We usually choose values of ${\sigma }_{k}=\eta {{ \mathcal M }}^{2}\lesssim {10}^{-3}$, consistent with the expectation that η ≪ 1 and ${{ \mathcal M }}^{2}\lesssim {10}^{-2}$.

4.1. Evolutions

For a calculation characterized by particular values of the physical parameters ($\eta {{ \mathcal M }}^{2}$, D_*, 〈Δτ〉, and K), we employ a second-order integrator, Huen's method. We generate M samples of Δτ and Ψ_k (τ). Each evolution has a duration ${\tau }_{\max }\approx M\langle {\rm{\Delta }}\tau \rangle$ , with different random realizations of the distributions of Δτ and Y_k during each eddy n. The initial conditions are chosen to be ${A}_{0}={\dot{A}}_{0}=0$ at τ = 0, but the long-term evolution does not depend on them.

In Figure 3, we see how short-term changes in 〈F₂ dA/d τ〉_{δ τ} correlate with changes in the energy E. The damping function 〈F₁(dA/d τ)²〉_{δ τ} is usually subdominant on this averaging timescale, but notice that antidamping becomes strong just before the blowout. As predicted in Section 2, the energy reaches E = 1/6 at blowout for our choice N₀ = 1.

Zoom In Zoom Out Reset image size

In Figure 4, we see how the negative extent of the amplitude A grows as the blowout is approached and reaches A = −1/N₀ = −1 at blowout. It is also seen that the period of the mode has increased by a factor of about 1.7 since the beginning of the evolution, as also expected from Section 2. Also note the behavior of the driving and damping functions during the eddies.

Zoom In Zoom Out Reset image size

In Figure 5, the damping constant has been increased by a factor of 10, which is sufficient to prevent a blowout (at least over the extended time interval indicated). This is our fiducial evolution. The maximum value of the energy (about 1/2 that required for a blowout) is about 10 times larger than its average (〈E〉 = 8 × 10⁻³ ).

Zoom In Zoom Out Reset image size

In Figure 6, the parameter choices that produce a blowout or not are indicated, with the corresponding rms amplitudes. Consider times τ ≫ τ_* = 1/D_*, so we can employ Equation (14) while the effects of nonlinearity are small. Then we obtain from Equation (15), when averaging over times δ τ = 10⁴,

$$\begin{eqnarray}&&\begin{array}{l}\langle E\rangle (\tau )\cong {\int }_{0}^{\tau }\langle S\rangle (\tau ^{\prime} )\exp [{D}_{* }(\tau ^{\prime} -\tau )]d\tau ^{\prime} \\ \ \ \approx \ {\int }_{\tau -{\tau }_{* }}^{\tau }\langle S\rangle (\tau ^{\prime} )d\tau ^{\prime} ,\end{array}\end{eqnarray} \tag{ 20 }$$

where $S(\tau )={F}_{2}(\tau ){\int }_{0}^{{\tau }_{* }}[{{dF}}_{2}(\tau -x)/d\tau ]\sin (x){dx}$. As $S\,\propto {(\eta {{ \mathcal M }}^{2})}^{2}$ and the integral spans τ_* = 1/D_*, $\langle E\rangle (\tau )\,\propto {(\eta {{ \mathcal M }}^{2})}^{2}/{D}_{* }$ . Because the blowout occurs at a fixed value of E = 1/6, the boundary of the blowout region does appear to approximately agree with the dependence ${D}_{* }\sim 10{(\eta {{ \mathcal M }}^{2})}^{2}$.

Zoom In Zoom Out Reset image size

From Figure 7, we see how 〈A²〉 and $\langle {F}_{2}^{2}{\rangle }^{1/2}$ depend on our four parameters. The dependence $\langle {A}^{2}\rangle \approxeq \langle E\rangle \sim 0.1{D}_{* }^{-1}{(\eta {{ \mathcal M }}^{2})}^{2}$ is also obtained from this figure.

Zoom In Zoom Out Reset image size

How does the average energy in the mode compare to the average energy in the effective turbulent eddies? The energy in the mode is approximately ∫d³ x ρ(∂ξⁱ /∂t)(∂ξ_i /∂t) ≈ c₀ ω² A². The energy in a dominant eddy is approximately 0.5∫d³ x ρ v² ≈ c₀ ω² F₂. Thus, equipartition would imply that $\langle {A}^{2}\rangle \approx \langle {F}_{2}^{2}{\rangle }^{1/2}$. The results shown in Figure 7 indicate that this is approximately true for small values of $\eta {{ \mathcal M }}^{2}$.

4.2. Contribution of the Mode to an Observed Power Spectral Density

Consider an axisymmetric (the most observable) mode of oscillation at radius r_m with volume v = 2π r_m ζ h²(r_m ). Its observed normalized power spectral density (PSD) P_x (f) (f = ω/2π) is related to the the photon count rate x by

$$\begin{eqnarray}&&\int {P}_{x}(f){df}=\langle {\left({\rm{\Delta }}x\right)}^{2}\rangle /\langle x{\rangle }^{2}\cong {P}_{x}(\max ){\rm{\Delta }}f,\end{eqnarray} \tag{ 21 }$$

where Δf is its FWHM. Adopting the approach of Nowak & Wagoner (1995),

$$\begin{eqnarray}&&x-\langle x\rangle ={\rm{\Delta }}x=({ \mathcal F }v/h)({\rm{\Delta }}\rho /\rho )\end{eqnarray} \tag{ 22 }$$

at r = r_m . The average photon number flux from the disk is ${ \mathcal F }(r)={\int }_{{\nu }_{1}}^{{\nu }_{2}}{{ \mathcal F }}_{\nu }(\nu ,r)\xi (\nu )d\nu$ for a detector efficiency ξ(ν). We shall employ the approximation ${ \mathcal F }{(r)={ \mathcal F }({r}_{m})(r/{r}_{m})}^{-{C}_{1}}$.

In addition, recall that the density fluctuation Δρ/ρ ≈ A(τ), related to its (not normalized) PSD by

$$\begin{eqnarray}&&\int {P}_{A}(q){dq}=2\pi \langle {A}^{2}\rangle .\end{eqnarray} \tag{ 23 }$$

The average count rate is (assuming that C₁ > 2)

$$\begin{eqnarray}\begin{array}{rcl}\langle x\rangle & = & 2\pi {\int }_{{r}_{\min }}^{{r}_{\max }}{ \mathcal F }(r){rdr}\\ & \approx & 2\pi { \mathcal F }({r}_{m}){\left({C}_{1}-2\right)}^{-1}{r}_{m}^{2}{\left({r}_{\min }/{r}_{m}\right)}^{(2-{C}_{1})},\end{array}\end{eqnarray} \tag{ 24 }$$

giving

$$\begin{eqnarray}\begin{array}{rcl}\int {P}_{x}(f){df} & \approx & {\left({C}_{1}-2\right)}^{2}{\zeta }^{2}{\left[h({r}_{m})/{r}_{m}\right]}^{2}{\left({r}_{\min }/{r}_{m}\right)}^{2({C}_{1}-2)}\langle {A}^{2}\rangle \\ & \equiv & \chi \langle {A}^{2}\rangle .\end{array}\end{eqnarray} \tag{ 25 }$$

Now consider the observed continuum PSD P_x (cont.) of black hole sources near the frequency f₀ of an HFQPO, with (M/10M_⊙)f₀ ∼ 200 Hz for black hole binaries (BHBs) and ∼50 Hz for NLS1 AGNs (Smith et al. 2021). We find that f₀ P_x (cont.) ∼ (2–20) × 10⁻⁴ for BHBs and ∼ (3 − 30) × 10⁻³ for NLS1s. The corresponding quantity for our mode is given (from Equation (25)) by

$$\begin{eqnarray}&&{f}_{0}{P}_{x}({\rm{max.,}}\,{\rm{QPO}})\approx \chi Q\langle {A}^{2}\rangle .\end{eqnarray} \tag{ 26 }$$

The quality factor Q ≡ f₀/Δf and $\chi \sim {\left(\zeta h/r\right)}^{2}$ .

During a quasi-equilibrium, when E(τ) is changing slowly but is not far below its blowout value of 1/6 (for ${N}_{0}^{2}=1$), the effects of our lowest-order nonlinearity control the width of the PSD peak, with the damping and driving forces subdominant. Then the standard analysis (Landau & Lifshitz 1976) predicts an amplitude-dependent shift of the anharmonic oscillator frequency of

$$\begin{eqnarray}&&\delta \omega /{\omega }_{0}=\delta q\approx -(5/12){N}_{0}^{2}E.\end{eqnarray} \tag{ 27 }$$

If the energy varies over a range ΔE during an evolution, it produces a corresponding width Δq of the PSD peak. Comparing the range of energies (that occupy a significant fraction of of the time) seen in Figure 5 with the width of the PSD peak corresponding to the same (fiducial) parameters in Figures 8 and 9, we see that this relation is approximately valid. It is only indicative, as the effects of higher-order nonlinearities have been neglected. From the results shown in Figure 9, we obtain the approximate relation Δq ∼ 1.0〈A²〉^0.9.

However, it is also seen from Figures 8 and 9 (and their extension for other values of the parameters) that the maximum value of each of the PSDs of the fundamental mode is P₀ ∼ 3, where ∫P_A (q)dq ≈ 1.1P₀Δq. Therefore, employing Equations (23) and (26), we see that visibility of an HFQPO (P_x (max.,QPO) > P_x (cont.)) requires that

$$\begin{eqnarray}&&\chi Q\langle {A}^{2}\rangle \gtrsim {f}_{0}{P}_{x}({\rm{cont.}})\sim (2\mbox{--}300)\times {10}^{-4}.\end{eqnarray} \tag{ 28 }$$

Recalling that χ ∼ (ζ h/r)², typical values of ζ ∼ 1 and h/r ∼ 10⁻² and the requirement A² < 1 could produce a visible HFQPO if Q ≫ 1 or χ ≫ χ(typical) ∼ 10⁻⁴. For instance, O'Neill et al. (2009) found that the trapped g-mode leaked into outgoing p-waves in their viscous hydrodynamic simulation. This would increase the value of ζ = Δr/h.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Nowak & Wagoner (1995) estimated the contribution of the turbulence to the continuum PSD. Referring to their Figure 1, it is seen that near the frequencies of the HFQPOs in BHBs,

$$\begin{eqnarray}&&{f}_{0}{P}_{x}({\rm{turb.}})\approx 0.1(L/{L}_{\mathrm{Edd}}){{ \mathcal M }}^{4}.\end{eqnarray} \tag{ 29 }$$

Then with L ≲ L_Edd and ${ \mathcal M }\lesssim 0.1$, we obtain f₀ P_x (turb.) ≲ 10⁻⁵. Comparing with the observed values of f₀ P_x (cont.) above, we see that it is unlikely that turbulence is the major contributor to the observed continuum PSD. This result is consistent with the fact that at lower frequencies, it is also seen from their Figure 1 that the turbulence in the accretion disk at the correspondingly larger radii provides very little of the observed power in the X-ray band.

A critical feature of the physical conditions that we have investigated is the fundamental difference between the nature of our turbulent damping function dH/d τ and the "turbulent viscosity" ν_t ≈ ℓ v_ℓ /3 (analogous to molecular viscosity) employed in the analysis of the effects of turbulence acting on the rate of shear of a quasi-steady flow (Thorne & Blandford 2017). See in particular Kato et al. (2008, Section 11.5). Because the contribution of turbulence to F₁(τ) is a total time derivative, there is an equal probability of short-term positive and negative damping. There should be no long-term temporal correlations between the mode and the turbulent eddies.

The physical origin of the damping of the p-modes in the Sun is uncertain (Basu 2016), but it should involve the coupling to higher frequency p-waves, which are damped when their wavelength becomes less than the scale height at the photosphere. This may also be true for our accretion disks. The generation of Alfvén waves could also contribute to the damping. Nowak et al. (1997) found that changes in entropy produced by radiative transfer effects led to the growth of modes within the types of accretion disks considered here.

We hope to consider magnetic forces in the future. In relevant numerical MHD simulations, although the ratio of magnetic to gas pressure is typically a few percent, the ratio of magnetic pressure to the Reynolds stress can be greater than unity (Dewberry et al. 2020). In addition, the magnetic forces can reduce the trapping of the g-modes (Fu & Lai 2009), but the amount depends on the relative magnitude of the poloidal and toroidal components (Ortega-Rodríguez et al. 2015). Reynolds & Miller (2009) found that the g-mode did not appear in their GRMHD simulations, although the number of orbits may have been insufficient to see growth of the mode. However, Dewberry et al. (2020) found that a small eccentricity in accretion disk orbits can excite r (g)-modes to large amplitudes in their MHD simulations. Warps can also excite such modes (Kato 2004, 2008).

In the future, we also hope to refine the modeling of the turbulence. One issue is the collective correlated effect of the eddies within the mode. Another is the effects of anisotropy, in particular the stretching of eddies in the e _ϕ direction. In addition, could (a) the infrequent significant increases of the mode energy or (b) a temporary reduction in the damping or (c) intermittency in the turbulent cascade be relevant to the low duty cycle exhibited by the HFQPOs?

The approach that we have taken to this problem may be applicable to other physical systems in which an oscillator is coupled to a turbulent environment.

We thank Nicole Lloyd-Ronning, Jeff Scargle, and James Stone for helpful comments. C.R.T. acknowledges financial support from the Stanford Physics Department Summer Undergraduate Research Program and NASA award 80NSSC20K0591 to Krista Lynne Smith. We would like to thank Stanford University and the Stanford Research Computing Center for providing computational resources on the Sherlock cluster and support. Support was also provided by the Center for Space Science and Astrophysics at Stanford. We thank the referee for helpful comments.