Viscous Evolution of Magnetized Clumps: A Source for X-Ray Flares in Gamma-Ray Bursts

X-ray flares are detected in about one third of Swift gamma-ray bursts (GRBs) (Burrows et al. 2005; Falcone et al. 2006; Nousek et al. 2006; Zhang et al. 2006), and marked as one of the most common phenomena that reveal the late-time erratic behavior of GRBs' central engines. Flares are detected in both long and short GRBs (Campana et al. 2006; Falcone et al. 2006; Romano et al. 2006; Margutti et al. 2011), and appear mostly in a time window of 10²–10⁵ s (Chincarini et al. 2010; Yi et al. 2016), which overlaps the afterglow timescale.

Several models have been proposed over the years to describe this late-time flaring activity. The leading external shock scenario (e.g., Meszaros & Rees 1997) failed to justify the temporal properties of X-ray flares. In particular, Lazzati & Perna (2007) argued that inhomogeneities in the external shock prevent this model from reproducing the observational properties of the flares. On the other hand, making more detailed analogies between the temporal and spectral properties of the flares and those of prompt emissions (e.g., Margutti et al. 2010) favors the idea that flares and prompt pulses have a common origin, that is to say they both trace the activity of the central engine (for more discussions see, e.g., Ioka et al. 2005; Beniamini & Kumar 2016). Therefore, studying the origin of the flares might open a new window on GRBs' central engines.

A variety of efforts have been made to attribute the flares to evolution of the central engine. King et al. (2005) suggested that fragmentation of a collapsing star and its subsequent accretion can lead to X-ray flare production. Moreover, ring-like fragmentation of the outer regions in a hyperaccreting disk might cause such a flaring activity, as proposed by Perna et al. (2006). From the magnetic point of view, some scenarios have been proposed that support the long-term activity of the central engine, e.g., late-time episodic accretion caused by a magnetic flux accumulation in the inner disk regimes, i.e., a magnetic barrier (Proga & Zhang 2006), and magnetic reconnection in a differentially rotating millisecond pulsar (Dai et al. 2006).

Motivated by the observed correlation between a flare's duration and its arrival time (Cusumano et al. 2006; O'Brien et al. 2006), and the similarity in the distributions of waiting times of X-ray flares and prompt gamma-ray emissions, which might lead to the consideration of a similar physical source for both events (Guidorzi et al. 2015), Dall'Osso et al. (2017) developed the idea of disk fragmentation in a more quantitative manner. They suggest that various flares with different arrival times (including both early and late-time ones) might be attributed to the viscous spreading of different clumps created during the early or late-time phase of central engine evolution. More specifically, they introduced an offset time as the time delay between the GRB trigger and the time at which X-ray flares set in, and argued that different offset times are subject to the starting point of various clumps' viscous evolution. So the early flares are considered to be caused by the prompt accretion of clumps generated during the early activity phase, while the delayed flares are due to clumps accreted later, or to fragments generated early that need to migrate toward the inner regions in order for viscous stress to dominate. They then provided a semi-analytical solution to this viscously evolved clump, and reproduced light curves of X-ray flares that, interestingly, matched the data well, after implementing some spectral corrections.

On the other hand, Shahamat & Abbassi (2017) studied the effects of magnetic field on both the vertical structure and neutrino luminosity of a self-gravitating neutrino-dominated accretion disk, as a plausible candidate for early-phase activity of GRB's central engine. The strongly magnetized nature of the disk, with a magnetic field of about 10^15–16 G, made both small- and large-scale effects of the magnetic field worth studying (the former was studied through magnetic viscosity, the latter via the magnetic braking process). They mainly found that such a consideration may shrink the gravitationally unstable regions toward the outer initial disk, and they highlight a magnetic barrier as a plausible mechanism for the lower accretion rates. All these findings, together with the fact that GRBs' central engines are supposed to be highly magnetized, lead us to develop the scenario of Dall'Osso et al. (2017) through taking the effects of magnetic field (small and large scale) into account, and to probe how the combination of magnetic barrier and fragmentation might affect this intriguing model to describe the light curves of X-ray flares.

Different approaches in the literature have addressed the magnetic barrier mechanism. For instance, D'Angelo & Spruit (2010, 2011, 2012) studied a magnetically truncated disk around a star (such as a neutron star). As they pointed out, the strong magnetic field near the star truncates the accretion disk. Due to angular momentum transfer between magnetic field and disk, one plausible state is that the disk is truncated outside the corotation radius ${R}_{c}={({{GM}}_{\ast }/{{\rm{\Omega }}}_{\ast }^{2})}^{1/3}$ (at which the frequency of the star's rotation equals the Keplerian frequency), where M_* and Ω_* are the mass and rotational frequency of the star. In this case a centrifugal barrier prevents the disk from accreting. Furthermore, the magnetic field lines form a "propeller regime." Such a topology of the field lines causes the gas inside the magnetospheric radius, R_m (the radius at which the accumulated magnetic field disrupts the accretion flow), to flow in freely along field lines toward the star's magnetic poles (Narayan et al. 2003). On the other hand, Narayan et al. (2003) discussed the physical conditions that govern the magnetically arrested disks (MADs) around black holes. They argue that MADs are considered to be accreting structures that are disrupted at R_m by the accumulated poloidal magnetic field in the vicinity of black holes. For R > R_m the flow is axisymmetric, while for R < R_m the flow breaks up into blobs or streams. Such streams or blobs can flow in slowly (with a velocity much less than the freefall velocity) toward the black hole via magnetic interchanges and reconnection. Regarding physical conditions imposed by a magnetic barrier around black holes, we adopt the strategy of Narayan et al. (2003), which leads us to consider two limiting boundary conditions. The first one is zero mass flux rate at the inner radius, and the second is taking the magnetospheric radius (R_m) as the inner boundary of the disk, namely R_in.

In the present paper, Section 2 manifests our model framework with a clarification of all the assumptions and approximations we made in order to smooth our way in extracting a semi-analytical solution, following the approach of Dall'Osso et al. (2017). The correlation between our model parameters and those related to the shape of the light curve will be considered in Section 3, together with some key analogies between our model predictions and observations. We then summarize and discuss our main conclusions in Section 4.

Hyperaccreting disks are gravitationally unstable in their outer regions, which might result in disk fragmentation and be a source of a disk's late-time activity (i.e., X-ray flares), as suggested by Perna et al. (2006), and discussed by Liu et al. (2014) and Shahamat & Abbassi (2017). A disk becomes gravitationally unstable once the Toomre parameter, Q, satisfies (Toomre 1964)

$$\begin{eqnarray}&&Q=\displaystyle \frac{{c}_{s}{\rm{\Omega }}}{\pi G{\rm{\Sigma }}}\lt 1\end{eqnarray} \tag{ 1 }$$

where c_s is the sound speed, Ω is the local Keplerian angular velocity, G is the gravitational constant, and Σ is the disk surface density. In the magnetized case, this criterion takes the form (Shu 1992)

$$\begin{eqnarray}&&{Q}_{\mathrm{mag}}=\displaystyle \frac{\sqrt{{c}_{s}^{2}+{v}_{{\rm{A}}}^{2}}{\rm{\Omega }}}{\pi {\rm{\Sigma }}G}\lt 1\end{eqnarray} \tag{ 2 }$$

where ${v}_{{\rm{A}}}=\tfrac{B}{\sqrt{4\pi \rho }}$ is the Alfvén velocity. However, in these unstable regimes, fragmentation into bound objects will occur if (Gammie 2001; Perna et al. 2006)

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}\lt {t}_{\mathrm{cirt}}\approx 3{{\rm{\Omega }}}^{-1},\end{eqnarray} \tag{ 3 }$$

where the cooling timescale is denoted by t_cool ≈ (H/R)²t_ν (Pringle 1991), with ${t}_{\nu }=\tfrac{2}{3}\tfrac{{R}^{2}}{\nu }$ regarded as the viscous timescale and t_cirt is the critical time scale. Fragments merge and/or accrete until their tidal influence gets strong enough to open a gap in the disk. This happens when the mass of the clump has increased to

$$\begin{eqnarray}&&{M}_{{\rm{frag}}}\simeq {\left(\displaystyle \frac{H}{R}\right)}^{2}{\alpha }^{1/2}{M}_{{\rm{BH}}}.\end{eqnarray} \tag{ 4 }$$

where M_BH is the mass of the central black hole (Takeuchi et al. 1996; Dall'Osso et al. 2017), H is the half thickness of the disk, and α is the viscous parameter.

In what follows, regarding the model proposed by Perna et al. (2006), we suppose that such a ring-like clump has been created in the outer regions as a sharp accumulated mass, more specifically a delta function, at radius R₀. Computing the viscous evolution of the clump, we also consider the standard equation for axisymmetric accretion disks (Kato et al. 2008):

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}t}{\rm{\Sigma }}(R,t)=\displaystyle \frac{1}{R}\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}R}\left({R}^{1/2}\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}R}(3\nu {\rm{\Sigma }}{R}^{1/2})\right).\end{eqnarray} \tag{ 5 }$$

Although this equation seems to ignore magnetic field effects, the small-scale effect of magnetic field can be taken into account through the magnetic viscosity term, i.e., ν. Hence, being interested in considering the large-scale effect of magnetic field, we add a term corresponding to the torque exerted by the Lorentz force. To do so, we make use of the approach of Lee et al. (2000) and Shahamat & Abbassi (2017) to modify Equation (5) to

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\rm{\Sigma }}(R,t)=-\displaystyle \frac{1}{R}\displaystyle \frac{\partial }{\partial R}\left[\displaystyle \frac{\tfrac{\partial }{\partial R}\left({R}^{3}\nu {\rm{\Sigma }}\tfrac{d{\rm{\Omega }}}{{dr}}\right)+\tfrac{{R}^{2}{B}_{\phi }{B}_{z}}{2\pi }}{\tfrac{d}{{dR}}({R}^{2}{\rm{\Omega }})}\right].\end{eqnarray} \tag{ 6 }$$

Considering the common assumption B_z ≈ B_R (e.g., Proga & Zhang 2006; Shahamat & Abbassi 2017), in order to get rid of the magnetic field components in the second term of the numerator, we can make use of the magnetic viscosity equation

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{R}{B}_{\phi }}{4\pi }=-\displaystyle \frac{3}{2}\alpha P,\end{eqnarray} \tag{ 7 }$$

where P is the total pressure of the accretion flow. On the other hand, we know that

$$\begin{eqnarray}&&P=-\displaystyle \frac{1}{\alpha }\rho \nu R\displaystyle \frac{d{\rm{\Omega }}}{{dR}},\end{eqnarray} \tag{ 8 }$$

and we consider Hρ ≈ Σ as a vertically averaged approximation (Kato et al. 2008). Finally, one can rewrite Equation (6) as

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\rm{\Sigma }}(R,t)=-\displaystyle \frac{1}{R}\displaystyle \frac{\partial }{\partial R}\left[\displaystyle \frac{\tfrac{\partial }{\partial R}\left({R}^{3}\nu {\rm{\Sigma }}\tfrac{d{\rm{\Omega }}}{{dr}}\right)+3\tfrac{d{\rm{\Omega }}}{{dR}}\tfrac{{\rm{\Sigma }}}{h}\nu {R}^{2}}{\tfrac{d}{{dR}}({R}^{2}{\rm{\Omega }})}\right],\end{eqnarray} \tag{ 9 }$$

where $h=\tfrac{H}{R}$.

In general, the viscosity ν depends on the surface density and Equation (9) is nonlinear. If, however, ν is only a function of radius, then the equation is linear and much more amenable to analytic methods. Therefore, following Dall'Osso et al. (2017) and Tanaka (2011), we assume the viscosity to follow a radial power law, ν ∝ Rⁿ, to achieve an exact solution for Σ(R, t) using a Green's function G,

$$\begin{eqnarray}&&{\rm{\Sigma }}(R,t)={\int }_{{R}_{{\rm{in}}}}^{\infty }G(R,{R}^{{\rm{{\prime} }}},t){\rm{\Sigma }}({R}^{{\rm{{\prime} }}},t=0)\,{{dR}}^{{\rm{{\prime} }}},\end{eqnarray} \tag{ 10 }$$

in which Σ(R, t = 0) is a given arbitrary profile at t = 0. Having Σ(R, t), the accretion power L_acc, due to the viscous spreading of the clump, is approximately

$$\begin{eqnarray}&&{L}_{{\rm{acc}}}\simeq {\int }_{{R}_{{\rm{in}}}}^{\infty }\displaystyle \frac{9}{4}{\rm{\Sigma }}(R,t)\nu (R){{\rm{\Omega }}}^{2}(R)2\pi R\,{dR}.\end{eqnarray} \tag{ 11 }$$

To compute the Green's function, we need to determine our desired boundary condition. In this regard, we encounter two strategies. The first is a boundary condition of zero central torque, which is of astrophysical interest especially in the case of accretion onto a black hole or a slowly rotating star, at radii larger than the radius of the innermost circular orbit or the stellar surface, respectively. The other is a boundary condition of zero mass flux. Having a strong central source of angular momentum, the accreting gas will be prevented from flowing in, and will instead accumulate near the center (Tanaka 2011). Such solutions can describe accretion disks around a compact binary (Pringle 1991), and also compact objects with strong central magnetic fields (Lynden-Bell & Pringle 1974). The latter seems similar to the physical situation imposed by the accumulated poloidal magnetic field near the innermost region of MADs, which might result in a magnetic barrier mechanism (Narayan et al. 2003; Shahamat & Abbassi 2017). Therefore, we are to apply both boundary conditions, in order to probe the effect of a magnetic barrier on the light curve of an X-ray flare.

Obtaining the appropriate Green's function, in the case of a boundary condition of zero torque, the surface density integral takes the following form:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{\Sigma }}}{{{\rm{\Sigma }}}_{0}}\,=\,{\displaystyle \int }_{0}^{\infty }{\displaystyle \frac{R}{{R}_{0}}}^{-n+(1-{b}_{1})/2}{\left(\displaystyle \frac{{R}_{0}}{{R}_{{\rm{in}}}}\right)}^{2-n}\\ \,\times \,(1-n/2){e}^{-2{\kappa }^{2}{\left(1-n/2\right)}^{2}t/{t}_{{\nu }_{{\rm{in}}}}}\\ \,\times \,[{Y}_{l}(\kappa ){J}_{l}(\kappa {x}_{0})-{Y}_{l}(\kappa {x}_{0}){J}_{l}(\kappa )]\\ \,\times \,\displaystyle \frac{[{Y}_{l}(\kappa ){J}_{l}(\kappa x)-{Y}_{l}(\kappa x){J}_{l}(\kappa )]}{[{Y}_{l}^{2}(\kappa )+{J}_{l}^{2}(\kappa )]}\kappa \,d\kappa \end{array}\end{eqnarray} \tag{ 12 }$$

and in the case of a boundary condition of zero mass flux one may achieve

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{\Sigma }}}{{{\rm{\Sigma }}}_{0}}\,=\,{\displaystyle \int }_{0}^{\infty }{\displaystyle \frac{R}{{R}_{0}}}^{-n+(1-{b}_{1})/2}{\left(\displaystyle \frac{{R}_{0}}{{R}_{{\rm{in}}}}\right)}^{2-n}\\ \,\times \,\left(1-n/2\right){e}^{-2{\kappa }^{2}{\left(1-n/2\right)}^{2}t/{t}_{{\nu }_{{\rm{in}}}}}\\ \,\times \,[{Y}_{l-1}(\kappa ){J}_{l}(\kappa {x}_{0})-{Y}_{l}(\kappa {x}_{0}){J}_{l-1}(\kappa )]\\ \,\times \,\displaystyle \frac{{Y}_{l-1}(\kappa ){J}_{l}(\kappa x)-{Y}_{l}(\kappa x){J}_{l-1}(\kappa )}{{Y}_{l-1}^{2}(\kappa )+{J}_{l-1}^{2}(\kappa )}\kappa \,d\kappa .\end{array}\end{eqnarray} \tag{ 13 }$$

To learn about new parameters and variables, and to get an insight into the adopted procedure, we refer readers to the Appendix.

What matters right now is how to fix the parameters Σ₀ and R₀.

2.1. Determination of Σ₀ and R₀

As previously stated, one can assume the ring-like clump to be a sharp concentration at radius R₀. This radius can be regarded as the same as that of gravitational instability. Therefore, we approximated Σ₀ with the expression

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}\simeq \displaystyle \frac{{M}_{\mathrm{frag}}}{2\pi {R}_{0}{l}_{{cl}}}\approx \displaystyle \frac{{M}_{\mathrm{frag}}}{2\pi {R}_{\mathrm{ins}}{l}_{{cl}}}\end{eqnarray} \tag{ 14 }$$

where l_cl denotes the clump size, which can be estimated as the local Jeans length (λ_j). First, to evaluate R₀, we need to find an expression for R_ins. This can be fulfilled by the use of the Toomre criterion (Equation (2)). Through some mathematical considerations as performed by Shahamat & Abbassi (2017), we can approximate magnetic pressure (${B}^{2}/8\pi$) with the total pressure (P). Mathematically speaking, one knows that

$$\begin{eqnarray*}&&{B}^{2}/8\pi \lt P+{B}^{2}/8\pi ,\end{eqnarray*}$$

we then have

$$\begin{eqnarray*}&&{p}_{{\rm{m}}{\rm{a}}{\rm{g}}}\approx \beta (P+{B}^{2}/8\pi )\,{\rm{w}}{\rm{i}}{\rm{t}}{\rm{h}}\,\beta \lt 1.\end{eqnarray*}$$

This provides us with the following expression:

$$\begin{eqnarray}&&{B}^{2}=8\pi \displaystyle \frac{\beta }{1-\beta }P\approx 0.3P\end{eqnarray} \tag{ 15 }$$

in which, to get the last equality, we choose β to be 0.01, which is the same value as adopted by Shahamat & Abbassi (2017). With the help of Equation (8) we reach

$$\begin{eqnarray}&&{B}^{2}=-\displaystyle \frac{0.3}{\alpha }\rho \nu R\displaystyle \frac{d{\rm{\Omega }}}{{dR}}.\end{eqnarray} \tag{ 16 }$$

Replacing these quantities, and applying relations $\dot{M}=3\pi {\rm{\Sigma }}\nu$ (valid for steady accretion) and n = 1/2 (an acceptable value for advection-dominated disks, see, e.g., Tanaka 2013), one may obtain

$$\begin{eqnarray}&&R\gt {\left(3{h}^{2}\alpha \displaystyle \frac{M}{\dot{M}}\sqrt{1.03{h}^{2}{GM}}\right)}^{2/3}={R}_{\mathrm{ins}}.\end{eqnarray} \tag{ 17 }$$

Finally, we are in need of the local magnetic Jeans length, λ_jmag, which takes the form

$$\begin{eqnarray}&&{\lambda }_{\mathrm{jmag}}={\lambda }_{j}\sqrt{1+\displaystyle \frac{{v}_{{\rm{A}}}^{2}}{{c}_{s}^{2}}}={c}_{s}\sqrt{\displaystyle \frac{\pi }{G\rho }}\sqrt{1+\displaystyle \frac{{v}_{{\rm{A}}}^{2}}{{c}_{s}^{2}}}.\end{eqnarray} \tag{ 18 }$$

Regarding all the above considerations, we will have

$$\begin{eqnarray}&&{\lambda }_{\mathrm{jmag}}=\sqrt{\displaystyle \frac{3.1{h}^{5}{\pi }^{2}\alpha M\sqrt{{{GMR}}_{\mathrm{ins}}}}{\dot{M}}}.\end{eqnarray} \tag{ 19 }$$

Hereby, we have our approximately estimated values for Σ₀ and R₀.

2.2. What about R_in and R_out?

The inner radius (R_in) differs in for the two different boundary conditions we adopted here. In the case of zero torque, we choose it to be equal to the innermost stable circular orbit, while for zero mass flux, in which the magnetic barrier can be taken into account, we considered the magnetospheric radius as the inner radius. This quantity can be defined as (Proga & Zhang 2006)

$$\begin{eqnarray}&&{R}_{m}\approx 60{\epsilon }_{-3}^{2/3}{\dot{M}}_{1}^{-2/3}{M}_{3}^{-4/3}{\phi }_{30}^{4/3}\end{eqnarray} \tag{ 20 }$$

where ϕ₃₀ ≡ ϕ/(10³⁰ cm² G). Here ϕ is the magnetic flux accumulated in the disk inner side and it has been taken to be about 10⁻², which can be supposed to be a typical value for the activity of the GRB's central engine (Xie et al. 2009; Shahamat & Abbassi 2017). Also, ₋₃ = 10³, in which is a parameter defined to be the ratio of the radial velocity of the accreted matter to its freefall velocity inside the magnetospheric radius (Shahamat & Abbassi 2017). We adopted a value of 10⁻³ for . Indeed, this parameter is less certain (for more discussion about this parameter see, e.g., Lloyd-Ronning et al. 2016; Narayan et al. 2003). ${\dot{M}}_{1}$ and M₃ denote $\dot{M}/1\,{M}_{\odot }\,{{\rm{s}}}^{-1}$ and M_BH/3 M_⊙, respectively.

On the other hand, as time passes, the accretion rate declines after the prompt phase. Metzger et al. (2008) reported that a neutrino-cooled accretion disk, in its late-time viscous evolution, experiences a decreasing mass accretion rate with a self-similar behavior $\dot{M}\propto {t}^{-4/3}$. Hence, we approximate the late-time accretion rate as

$$\begin{eqnarray}&&\dot{M}=\dot{{M}_{0}}{\left(\displaystyle \frac{t}{{t}_{0}}\right)}^{-4/3}\end{eqnarray} \tag{ 21 }$$

where $\dot{{M}_{0}}$ is the accretion rate related to a given time, t₀, specified (arbitrarily) during the late activity phase. We have fixed $\dot{{M}_{0}}$ to have the value 0.04 M_⊙ s⁻¹ at t₀ = 50 s after the prompt phase (which is supposed to have an accretion rate of about 0.1–10 M_⊙ s⁻¹, typically).

In the context of the outer boundary, we consider R_out as the radius at which the viscous evolution of the clump sets in. To estimate this radius, we take the shear, self-gravity, and tidal forces per unit mass exerted on the clump as introduced by Dall'Osso et al. (2017). We then add the Lorentz force in order to account for the large-scale effects of magnetic field.

A clump of linear size l_cl and mass M_cl will be affected by the shear force per unit mass

$$\begin{eqnarray}&&{F}_{\nu }=\displaystyle \frac{{l}_{{\rm{cl}}}\dot{M}{\rm{\Omega }}}{2\pi {M}_{{\rm{cl}}}}\end{eqnarray} \tag{ 22 }$$

in which the relation $\dot{M}=3\pi \nu {\rm{\Sigma }}$, valid for a steady-state disk, has been considered.

Given the self-gravity force per unit mass of the clump, ${F}_{{\rm{SG}}}=\tfrac{{{GM}}_{{\rm{cl}}}}{{l}_{{\rm{cl}}}^{2}}$, the tidal force due to the central object, ${F}_{T}=\tfrac{{{GM}}_{{\rm{BH}}}}{{r}^{2}}\left(\tfrac{{l}_{{\rm{cl}}}}{r}\right)$ (Dall'Osso et al. 2017), and the Lorenz force per unit mass, ${F}_{B}=\tfrac{{B}_{r}{B}_{\phi }}{4\pi {\rm{\Sigma }}}$ (Lee et al. 2000), viscous spreading of the clump will start roughly when F_ν + F_T + F_B > F_SG. It is worth noting that, making use of Equation (7) and $\dot{M}=3\pi {\rm{\Sigma }}\nu$, the Lorentz force can be rewritten in the form $-\tfrac{9{GM}\alpha h}{4{r}^{2}}$. After all, the condition for being viscously evolved reads

$$\begin{eqnarray}&&\displaystyle \frac{{{GMl}}_{{\rm{cl}}}}{{r}^{3}}+\displaystyle \frac{{l}_{{\rm{cl}}}\dot{M}}{2\pi {M}_{{\rm{cl}}}}\sqrt{\displaystyle \frac{{GM}}{{r}^{3}}}-\displaystyle \frac{9{GM}\alpha h}{4{r}^{2}}-\displaystyle \frac{{{GM}}_{{\rm{cl}}}}{{l}_{{\rm{cl}}}^{2}}\gt 0.\end{eqnarray} \tag{ 23 }$$

The solution to this inequality can provide us with an upper limit for the outer boundary R_out.

2.3. Issues on Bolometric and X-Ray Luminosities

There are some significant points to be noted. It is the relativistic jet (which can be created through the Blandford–Znajek mechanism, as the most viable scenario that can efficiently extract the rotational energy of a spinning black hole by a large-scale magnetic field threading the disk (Blandford & Znajek 1977)) that might produce the GRB's X-ray flares. To be more precise, such emissions are subject to the radiative mechanisms through which the jet power is converted into radiation. For the sake of simplicity, we assume the radiative efficiency in the jet to be constant. Therefore, the bolometric luminosity of the flare will be

$$\begin{eqnarray}&&{L}_{\mathrm{bol}}={f}_{\mathrm{rad}}{L}_{\mathrm{acc}}\end{eqnarray} \tag{ 24 }$$

where f_rad is the efficiency of conversion of energy from mass accretion into radiation, which we assume to be constant. We also adopted f_rad ≈ 0.5 in view of discussions by Fragile et al. (2012). However, in case of the zero-flux boundary condition, i.e., taking the magnetic barrier into account, and consequently using the radiation efficiency in the MAD model (Narayan et al. 2003), we adopted a value of ∼1.0 for f_rad, based on discussions by Tchekhovskoy et al. (2011) and McKinney et al. (2012).

On the other hand, we need to consider the X-ray bandpass (0.3–10 keV) in order to evaluate this model by making a comparison with observations that are performed in a finite X-ray energy band. Thus, an efficiency coefficient is required to consider the X-ray energy band instead of the bolometric light curve. Here, we follow Jin et al. (2010) and Fan & Piran (2006) by approximating the X-ray flare luminosity as a constant coefficient of bolometric luminosity, L_X = f_XL_bol, where f_X (X-ray flare efficiency) is estimated as ≈0.1.

3.1. Correlation of Model Parameters and Spectral Properties

3.2. Observations and Model Predictions

From the study by Chincarini et al. (2010) of 113 flares in the 0.3–10 keV energy band of the X-ray Telescope (XRT), and four subenergy bands, some observational characteristics of the X-ray light curves can be clarified as follows.

• (i)  
  The ratio of rise time to decay time is constant, implying that both timescales grow by the same factor, so that t_d ≈ 2t_r. Consequently, flares are self-similar in time.
• (ii)  
  The width evolves linearly with the peak time: w ≈ 0.2. These two points are the key features that strongly distinguish the flare emission from the prompt phase.

Moreover, analyzing 468 bright X-ray flares from the GRBs observed by Swift between 2005 and 2015, Yi et al. (2016) argued that the peak luminosity decreases with the peak time, following a power-law behavior ${L}_{p}\propto {t}_{p}^{-1.27}$.

In general, for both boundary conditions, our model leads to rather observationally well matched predictions, considering the mentioned shape parameters.

In the case of a zero-torque boundary condition, Figures 1 and 2 demonstrate how our model respects these three criteria. We have considered a variety of model parameters, specified within the panels, and different offset times during the early phase of flaring activity (with t_p < 1000 s), i.e., 50 s, 100 s, 200 s, and 1000 s, in order to enhance the validity of our fits.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 reflects the logarithmic behavior of X-ray luminosity over time. Four different offset times of 50 s, 100 s, 200 s, and 300 s have been considered. The best fit for peak luminosity versus peak time (the purple dotted–dashed line) indicates a power-law trend with a power index of ∼−1.44, which is interestingly close to the value discussed by Yi et al. (2016), i.e., −1.27. These light curves have been plotted for three sets of model parameters: M_BH = 3 M_⊙, h = 0.8, and α = 0.13 (green curves), M_BH = 2.5 M_⊙, h = 0.8, and α = 0.15 (orange dashed curves), and M_BH = 5 M_⊙, h = 0.6, and α = 0.22 (blue dotted curves). The inset panel shows the fit for a broader range of model parameters including those considered in Figure 2. Being of a rather similar power index, i.e., ∼−1.35, the variation of model parameters does not affect the correlation we obtained in our model, which is in good agreement with its phenomenological relation, namely ${L}_{p}\propto {t}_{p}^{-1.27}$. On the other hand, as expressed in Section 3.1, Figure 2 reveals the robustness of the shape of the light curve in our model through the considerable (acceptable) insensitivity of the asymmetry parameter (width over peak time). Furthermore, our predicted values for these two shape parameters, i.e., k ≈ 2.5 and w ≈ 0.11, are in rather good agreement with those predicted by Chincarini et al. (2010), i.e., k ≈ 2 and w ≈ 0.2.

In order to obtain results for the zero-mass-flux boundary condition, and consequently take the magnetic barrier into account, we need to probe for a set of model parameters that satisfy the condition

$$\begin{eqnarray}&&\tau =\displaystyle \frac{{t}_{\nu }}{{t}_{\mathrm{diff}}}\lt 1\end{eqnarray} \tag{ 25 }$$

where ${t}_{{\rm{diff}}}\approx H/{v}_{{\rm{A}}}$ is the diffusion timescale that estimates the magnetic field buoyancy and the time it takes to rise to the disk surface (Shahamat & Abbassi 2017). Such a limitation ensures that the magnetic barrier is likely to exist during the clump's accretion due to the accumulation of the magnetic flux inside the inner regions (Xie et al. 2009; Shahamat & Abbassi 2017). The adopted parameters in Figures 3 and 4 comply with this criterion. For instance, for M_BH = 2.5 M_⊙, h = 0.8, and α = 0.15 we found an average value of about 0.8 for τ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 illustrates the variation of X-ray luminosity with time for different sets of model parameters and for offset times of 50 s, 100 s, 200 s, and 300 s. These parameters are similar to those considered in Figure 1; to show the generality of our result, we extended the sets of parameters for the inset to include those regarded in Figure 4. Besides the fact that the power-law fit for L_p and t_p is in good agreement with phenomenological predictions, taking the magnetic barrier into account leads to an increase in the peak luminosity by an order of magnitude in comparison with similar cases for the zero-torque boundary condition. This is an expected result according to general relativistic magnetohydrodynamic simulations performed, e.g., by McKinney et al. (2012) and Tchekhovskoy et al. (2011).

Regarding Figure 4, our fits for the other two shape parameters k and w confirm a satisfactory robustness of our model in reproducing the shape of the light curve for the case of the zero-mass-flux boundary condition as well as the zero-torque boundary condition.

However, we must point out some issues on the subject of the magnetic barrier mechanism. On the one hand, Liu et al. (2014), through the study of a self-gravitating neutrino-dominated accretion disk (called NDAF), confirmed that instabilities appear in the outer parts and away from the equatorial plane. From a magnetic point of view, Shahamat & Abbassi (2017) probed for the possibility of fragmentation in a self-gravitating magnetized NDAF and also found it to be possible in the outer parts and away from the equatorial plane. Thus, we think that clumps might be thinner than what we expected from the vertical structure of the late-time advection-dominated central engine, that is, smaller values of the h parameter should be more worthy of consideration. As a consequence, clumps of smaller masses would be more likely to form. On the other hand, regarding α in its physically meaningful range, we found that the magnetic barrier is mainly plausible in the case of heaver clumps (with masses of the order of 0.4–0.8 M_⊙). Therefore, we think that such a limitation might make it lees likely for the magnetic barrier to govern the evolution of the generated clumps. Moreover, we considered ϕ to be of the order of 10⁻², which might be an overestimation for the late-time accumulated magnetic flux, since, as time passes and the accretion rate drops, both the clump's magnetic field and the accumulated magnetic flux decline. Such a drop in magnetic field, which is more likely for higher values of t_off, might weaken the possibility of the magnetic barrier even more (see, e.g., McKinney et al. 2012 and references therein).

In the context of a power-law correlation between L_p and t_p, there are some points of importance to be noticed. We need to emphasize that the time dependence we considered to set our model parameters, such as the radius at which the clump is created, cannot lead to such a power-law behavior. This is because we have used these relations to apply some physically justified constraints to our model parameters, such as $\dot{M}$ and R₀, and these relations have not been applied directly in the time evolution equation of luminosity. In particular, Equation (21) has been used only to determine the accretion rate at t_off and it has not been considered in the main formalism of the model, which provides us with the time evolution of luminosity. In other words, we only aimed to use these relations to fix some model parameters regarding the offset time at which the clump starts to accrete. Hence, their time dependences have not influenced the clump's evolution directly, so that one can think of this power-law trend as a general outcome of our model.

We reconstructed the quantitative framework established by Dall'Osso et al. (2017), through a magnetic style in which small- and large-scale effects of magnetic field have been taken into account. To this end, we added one more term corresponding to the torque exerted by the Lorentz force (see Equation (5)). To solve this evolution equation, some approximations and assumptions have been implemented in order to be able to adopt a similar semi-analytic approach to that proposed by Tanaka (2011) and followed by Dall'Osso et al. (2017). Through this method, we have provided final solutions that differ from those achieved by the previous work in some respects. First, for both cases of boundary conditions, the power-law dependence on radius changed in comparison with solutions obtained by Tanaka (2011). Second, considering h as an arbitrary parameter that is not just fixed at "1" (this affects ${t}_{{\nu }_{\mathrm{in}}}$ in the final solution) causes another difference in our result with respect to what has been achieved by Dall'Osso et al. (2017). Finally, to estimate model parameters such as clump mass, we limited ourselves to the correlations introduced between these parameters, so that our results are able to verify the framework we defined; however, in Dall'Osso et al. (2017) it was of less importance to consider such correlations in providing data. On the other hand, to investigate the possibility of a magnetic barrier occurring during each clump's viscous evolution, we compared the diffusion and viscous timescales. Moreover, two different boundary conditions have been imposed, i.e., zero torque and zero mass flux, in the hope of distinguishing the presence of the accumulated magnetic flux inside the inner regions, and consequently the magnetic barrier mechanism. We then studied the validity of our model through some key analogies to phenomenological findings. Our main conclusions can be identified as follows.

• 1.  
  The ratio of width to peak time is found to be an almost constant, although somewhat scattered, parameter with an average value ∼0.1, which is in rather good agreement with its phenomenological estimate, ∼0.2. We also came to the conclusion that a magnetic barrier, if it exists, will not change this parameter significantly; instead, it might cause our data to become less scattered (left panel in Figure 4).
• 2.  
  The insensitivity to model parameters of the skewness of the X-ray light curves (the asymmetry parameter) reveals another robust aspect of our model. This important feature is respected for both boundary conditions, with a similar value of about ∼2.4, which is in satisfying agreement with its observational prediction, i.e., ∼2.
• 3.  
  Interestingly, the power-law correlation between the maximum value of luminosity and peak time, with a value of about ∼−1.3, is found to be close to what was estimated by Yi et al. (2016) from a phenomenological point of view, i.e., ∼−1.27. Both boundary conditions adopted here respect this feature.

There are some more points worth discussing here. For one thing, Dall'Osso et al. (2017) found the bolometric luminosity to be far from satisfactory in comparison with observations. However, after conducting a spectral correction through the adoption of a time-dependent X-ray efficiency, they modified the results to be well matched with observations. Whereas, considering the large-scale effect of a magnetic field led us to obtain shape parameters in close agreement with phenomenological values in both cases of bolometric and X-ray luminosities. Note that, in our study, bolometric luminosity differs from X-ray luminosity only by a constant coefficient, L_X = f_XL_bol. Thus, the bolometric and X-ray light curves have the same shape parameters in our model.

Another point is that we assumed f_rad and f_X to be constants. Accounting for any inconstancy about these two parameters might affect the shape of the light curve. On the one hand, simulations lead us to conclude that the efficiency in converting torus mass into jet energy varies from a few per cent up to more than 100%. The black hole spin, the disk thickness, and the magnetic flux might effectively alter this parameter (see, e.g., Giacomazzo et al. 2013 and references therein). On the other hand, regarding f_X, Dall'Osso et al. (2017) argued that a time-dependent X-ray efficiency with a decreasing trend in time (which is in agreement with the "curvature effect" if the emitting region has an accelerated bulk relativistic motion (Uhm & Zhang 2015, 2016)) might affect the spectral properties of X-ray flares.

Separate from these considerations, our simple model, together with the other probable modifications, could efficiently reproduce the spectral properties of X-ray flares. However, more precise studies via simulations appear to be essential to shed more light on the dark nature of the late-time evolution of a GRB's central engine.

Following Dall'Osso et al. (2017), to find the Green function, we adopted assumed the viscosity to follow a radial power law, ν ∝ Rⁿ, and assumed a separable ansatz of the form ${\rm{\Sigma }}(R,t)={R}^{p}\sigma (R){e}^{-\lambda t}$, where p and Λ are real numbers and σ is an arbitrary function of R, so that Equation (9) can be rewritten as a Bessel differential equation (Tanaka 2011):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\partial }^{2}\sigma }{\partial {R}^{2}}+\displaystyle \frac{\partial \sigma }{R\partial R}(1.5+3/h+2(n+p))\\ +\,\sigma (R)\left[\displaystyle \frac{{\rm{\Lambda }}}{3s}{R}^{-n}+\displaystyle \frac{(n+p)(0.5+3/h+(n+p))}{{R}^{2}}\right]=0\end{array}\end{eqnarray} \tag{ 26 }$$

where s = νR⁻ⁿ is a constant. It should be mentioned that we considered the Newtonian potential with ${\rm{\Omega }}=\sqrt{{GM}/{R}^{3}}$. We also choose $p=n+\tfrac{1-{b}_{1}}{2}$, with b₁ = 1.5 + 3/h, and Λ = 3sk². All these considerations, besides comparing Equation (26) with the transformed Bessel function

$$\begin{eqnarray}\begin{array}{rcl} & & \displaystyle \frac{{\partial }^{2}y}{\partial {x}^{2}}-\displaystyle \frac{\partial y}{\partial x}\left(\displaystyle \frac{2\alpha -1}{x}\right)+y(x)\left[\displaystyle \frac{{\alpha }^{2}-{l}^{2}{\gamma }^{2}}{{x}^{2}}+{\beta }^{2}{\gamma }^{2}{x}^{2\gamma -2}\right]\\ & & \,=0,\end{array}\end{eqnarray} \tag{ 27 }$$

yield

$$\begin{eqnarray*}&&\gamma =1-n/2;\,\beta =\displaystyle \frac{k}{1-n/2};\,\alpha =2n;\,l=\displaystyle \frac{| 1-{b}_{1}| }{2(1-n/2)}.\end{eqnarray*}$$

Now, the solution to Equation (10) can be written in the following form:

$$\begin{eqnarray}&&{\sigma }_{k}(R)={R}^{-2n}[A(k){J}_{l}({ky})+B(k){Y}_{l}({ky})]\end{eqnarray} \tag{ 28 }$$

where $y=\tfrac{{R}^{1-n/2}}{1-n/2}$. It is worth noting that in case of non-integer l, Y_l must be replaced by J_−l. Integrating the above solution over all possible k-modes gives the solution

$$\begin{eqnarray}&&{\rm{\Sigma }}(R,t)={\int }_{0}^{\infty }{R}^{-n-1/4}[A(k){J}_{l}({ky})+B(k){Y}_{l}({ky})]{e}^{-3{{sk}}^{2}t}{dk}.\end{eqnarray} \tag{ 29 }$$

The mode-weighting functions A(k) and B(k) are determined by the boundary conditions and the initial surface density profile Σ(R, t = 0).

We just elaborate the zero-mass-flux boundary condition here, because the zero-torque one is implemented in a similar way. Before any further calculation to get the final solution, it is better to determine s in the power-law expression for viscosity as follows. The (magnetic) viscosity can be parameterized as $\nu =\alpha {c}_{s}^{2}/{\rm{\Omega }}$ (Shakura & Sunyaev 1973; Pringle 1991). On the other hand, invoking hydrostatic equilibrium perpendicular to the disk plane gives us $H={c}_{s}/{\rm{\Omega }}$, and consequently one can infer $\nu =\alpha {h}^{2}\sqrt{{GM}}{R}^{1/2}$, which yields $s=\alpha {h}^{2}\sqrt{{GM}}$ and n = 1/2. Furthermore, considering the continuity equation

$$\begin{eqnarray}&&R\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}+\displaystyle \frac{\partial {\rm{\Sigma }}{{Rv}}_{R}}{\partial R}=0,\end{eqnarray} \tag{ 30 }$$

together with

$$\begin{eqnarray*}&&\dot{M}=-2\pi {\rm{\Sigma }}{{Rv}}_{R},\end{eqnarray*}$$

lead us to

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}=\displaystyle \frac{1}{2\pi R}\displaystyle \frac{\partial \dot{M}}{\partial R}.\end{eqnarray} \tag{ 31 }$$

Comparing Equations (31) and (9) provides us with the following result:

$$\begin{eqnarray}&&\dot{M}\propto \left[\displaystyle \frac{\partial }{\partial R}\left({R}^{3}\nu {\rm{\Sigma }}\displaystyle \frac{d{\rm{\Omega }}}{{dr}}\right)+3\displaystyle \frac{d{\rm{\Omega }}}{{dR}}\displaystyle \frac{{\rm{\Sigma }}}{h}\nu {R}^{2}\right]{| }_{{R}_{\mathrm{in}}}=0.\end{eqnarray} \tag{ 32 }$$

By substitution of the corresponding expressions for ν, Ω, and Σ, and making use of the recurrence relations between different Bessel functions

$$\begin{eqnarray*}\begin{array}{rcl}{J}_{l-1}(x)-{J}_{l+1}(x) & = & 2{J}_{l}^{{\prime} }(x)\\ {J}_{l-1}(x)+{J}_{l+1}(x) & = & 2l\displaystyle \frac{{J}_{l}}{x}(x),\end{array}\end{eqnarray*}$$

the following relation can be achieved:

$$\begin{eqnarray}&&\displaystyle \frac{A(\kappa )}{B(\kappa )}=-\displaystyle \frac{{Y}_{l-1}(\kappa )}{{J}_{l-1}(\kappa )},\end{eqnarray} \tag{ 33 }$$

where κ = ky_in. As we proceed below, we need to obtain some more ingredients to get the final solution for Σ(R, t). First,

$$\begin{eqnarray*}&&{\rm{\Lambda }}=3{{sk}}^{2}=3\nu {R}^{-n}{k}^{2},\end{eqnarray*}$$

and, second, we know

$$\begin{eqnarray}&&{t}_{\nu }=\displaystyle \frac{2}{3}\displaystyle \frac{{R}^{2}}{\nu }=\displaystyle \frac{2}{3}\displaystyle \frac{{R}^{2-n}}{s}.\end{eqnarray} \tag{ 34 }$$

Thus, Λ can be rewritten in the form

$$\begin{eqnarray*}&&{\rm{\Lambda }}=-2{\kappa }^{2}{\left(1-n/2\right)}^{2}\displaystyle \frac{1}{{t}_{{\nu }_{\mathrm{in}}}}.\end{eqnarray*}$$

Combining all the above, the expression for surface density takes the form

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{\Sigma }} & = & {\displaystyle \int }_{0}^{\infty }{R}^{-n+(1-{b}_{1})/2}[c(\kappa ){\kappa }^{-1}][{Y}_{l-1}(\kappa ){J}_{l}(\kappa x)\\ & & -\,{Y}_{l}(\kappa x){J}_{l-1}(\kappa )]{e}^{-2{\kappa }^{2}{\left(1-n/2\right)}^{2}t/{t}_{{\nu }_{{\rm{in}}}}}\kappa \,d\kappa \end{array}\end{eqnarray*}$$

where ${t}_{{\nu }_{\mathrm{in}}}$ is the viscous timescale of the inner radius of the disk and x = y/y_in. The mode weight c(κ) might be obtainable by making use of a generalized Weber transform (Zhang et al. 2006; Tanaka 2011)

$$\begin{eqnarray*}&&\begin{array}{l}{\phi }_{l}(x)\,=\,{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{W}_{l}(k,x;a,b)}{{Q}_{l}^{2}(k,;a,b)}{{\rm{\Phi }}}_{l}(\kappa )\kappa \,d\kappa \\ {{\rm{\Phi }}}_{l}(\kappa )\,=\,{\displaystyle \int }_{1}^{\infty }{W}_{l}(\kappa ,x;a,b){\phi }_{l}(x)x\,{dx}\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{W}_{l} & = & {J}_{l}(\kappa x)[(a-{lb}){Y}_{l}(\kappa )+b\kappa {Y}_{l-1}(\kappa )]\\ & & -\,{Y}_{l}(\kappa x)[(a-{lb}){J}_{l}(\kappa )+b\kappa {J}_{l-1}(\kappa )]\end{array}\end{eqnarray*}$$

and

$$\begin{eqnarray*}\begin{array}{rcl}{Q}_{l}^{2} & = & {\left[(a-{lb}){Y}_{l}(\kappa )+b\kappa {Y}_{l-1}(\kappa )\right]}^{2}\\ & & -\,{\left[(a-{lb}){J}_{l}(\kappa )+b\kappa {J}_{l-1}(\kappa )\right]}^{2}.\end{array}\end{eqnarray*}$$

If a = 1 and b = 0, the pair is identical to the ordinary Weber transform. However, a = l and b = 1 correspond to our desired boundary condition, i.e., zero mass flux. After all, what would be the Green function? Clearly,

$$\begin{eqnarray*}\begin{array}{rcl}{R}^{n-(1-{b}_{1})/2}{\rm{\Sigma }}(R,t=0) & = & {\displaystyle \int }_{0}^{\infty }[c(\kappa ){\kappa }^{-1}][{Y}_{l-1}(\kappa ){J}_{l}(\kappa x)\\ & & -\,{Y}_{l}(\kappa x){J}_{l-1}(\kappa )]\kappa \,d\kappa .\end{array}\end{eqnarray*}$$

Finally, through applying the above considerations and running some more calculations, one may achieve the following result for c(κ):

$$\begin{eqnarray*}&&\begin{array}{c}\begin{array}{r}\begin{array}{c}c(\kappa )\\ ={\int }_{1}^{\infty }\displaystyle \frac{{R}^{n-(1-{b}_{1})/2}{\rm{\Sigma }}\left(\left.R,0)[{Y}_{l-1}(\kappa ){J}_{l}(\kappa x)-{Y}_{l}(\kappa x){J}_{l-1}(\kappa \right)\right]{xdx}}{\left[\left.{Y}_{l-1}^{2}(\kappa )+{J}_{l-1}^{2}(\kappa \right)\right]}.\end{array}\end{array}\end{array}\end{eqnarray*}$$

Thus,

$$\begin{eqnarray*}&&\begin{array}{c}{\rm{\Sigma }}(R,t)={\displaystyle \int }_{0}^{\infty }\left(\displaystyle \frac{{\displaystyle \int }_{{r}_{{\rm{in}}}}^{\infty }[{Y}_{l-1}(\kappa ){J}_{l}(\kappa {x}^{{\rm{{\prime} }}})-{Y}_{l}(\kappa {x}^{{\rm{{\prime} }}}){J}_{l-1}(\kappa )]{\rm{\Sigma }}({R}^{{\rm{{\prime} }}},0)\tfrac{{\left({R}^{{\rm{{\prime} }}}\right)}^{\left(1+{b}_{1})/2\right.}}{{R}_{{\rm{in}}}^{2-n}}{{dR}}^{{\rm{{\prime} }}}}{\left[\left.{Y}_{l-1}^{2}(\kappa )+{J}_{l-1}^{2}(\kappa \right)\right]}\right)\\ \times (1-\,n/2){R}^{-n+(1-{b}_{1})/2}[{Y}_{l-1}(\kappa ){J}_{l}(\kappa x)-{Y}_{l}(\kappa x){J}_{l-1}(\kappa )]{e}^{-2{\kappa }^{2}{\left(1-n/2\right)}^{2}t/{t}_{{\nu }_{{\rm{in}}}}}\kappa \,d\kappa .\end{array}\end{eqnarray*}$$

Considering Equation (10), the Green's function would be obtainable as

$$\begin{eqnarray}\begin{array}{rcl}G(R,{R}^{{\rm{{\prime} }}},t) & = & {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{[{Y}_{l-1}(\kappa ){J}_{l}(\kappa x{\rm{^{\prime} }})-{Y}_{l}(\kappa x{\rm{^{\prime} }}){J}_{l-1}(\kappa )]}{[{Y}_{l-1}^{2}(\kappa )+{J}_{l-1}^{2}(\kappa )]}\\ & & \times \,\displaystyle \frac{{\left({R}^{{\rm{{\prime} }}}\right)}^{(1+{b}_{1})/2}}{{R}_{{\rm{in}}}^{2-n}}(1-n/2){R}^{-n+(1-{b}_{1})/2}\\ & & \times \,[{Y}_{l-1}(\kappa ){J}_{l}(\kappa x)-{Y}_{l}(\kappa x){J}_{l-1}(\kappa )]\\ & & \times \,{e}^{-2{\kappa }^{2}{\left(1-n/2\right)}^{2}t/{t}_{{\nu }_{{\rm{in}}}}}\kappa \,d\kappa .\end{array}\end{eqnarray} \tag{ 35 }$$

Finally, regarding the initial surface density as

$$\begin{eqnarray}&&{\rm{\Sigma }}(R,t=0)={{\rm{\Sigma }}}_{0}{R}_{0}\delta (R-{R}_{0})\end{eqnarray} \tag{ 36 }$$

in which Σ₀ and R₀ are arbitrary (Tanaka 2011), we have

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{\Sigma }}}{{{\rm{\Sigma }}}_{0}}\,=\,{\displaystyle \int }_{0}^{\infty }{\displaystyle \frac{R}{{R}_{0}}}^{-n+(1-{b}_{1})/2}{\left(\displaystyle \frac{{R}_{0}}{{R}_{{\rm{in}}}}\right)}^{2-n}\\ \,\times \,(1-n/2){e}^{-2{\kappa }^{2}{\left(1-n/2\right)}^{2}t/{t}_{{\nu }_{{\rm{in}}}}}\\ \,\times \,[{Y}_{l-1}(\kappa ){J}_{l}(\kappa {x}_{0})-{Y}_{l}(\kappa {x}_{0}){J}_{l-1}(\kappa )]\\ \,\times \,\displaystyle \frac{[{Y}_{l-1}(\kappa ){J}_{l}(\kappa x)-{Y}_{l}(\kappa x){J}_{l-1}(\kappa )]}{[{Y}_{l-1}^{2}(\kappa )+{J}_{l-1}^{2}(\kappa )]}\kappa \,d\kappa \end{array}\end{eqnarray} \tag{ 37 }$$

Likewise, in the case of the zero-central-torque boundary condition, one can evaluate the Green's function as

$$\begin{eqnarray}\begin{array}{rcl}G(R,R{\rm{^{\prime} }},t) & = & {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{[{Y}_{l}(\kappa ){J}_{l}(\kappa x{\rm{^{\prime} }})-{Y}_{l}(\kappa x{\rm{^{\prime} }}){J}_{l}(\kappa )]}{[{Y}_{l}^{2}(\kappa )+{J}_{l}^{2}(\kappa )]}\\ & & \times \,\displaystyle \frac{{\left({R}^{{\rm{{\prime} }}}\right)}^{(1+{b}_{1})/2}}{{R}_{{\rm{in}}}^{2-n}}(1-n/2){R}^{-n+(1-{b}_{1})/2}\\ & & \times \,[{Y}_{l}(\kappa ){J}_{l}(\kappa x)-{Y}_{l}(\kappa x){J}_{l}(\kappa )]\\ & & \times \,{e}^{-2{\kappa }^{2}{\left(1-n/2\right)}^{2}t/{t}_{{\nu }_{{\rm{in}}}}}\kappa \,d\kappa \end{array}\end{eqnarray} \tag{ 38 }$$

and subsequently the surface density will read

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{\Sigma }}}{{{\rm{\Sigma }}}_{0}}\,=\,{\displaystyle \int }_{0}^{\infty }{\displaystyle \frac{R}{{R}_{0}}}^{-n+(1-{b}_{1})/2}{\left(\displaystyle \frac{{R}_{0}}{{R}_{{\rm{in}}}}\right)}^{2-n}\\ \,\times \,(1-n/2){e}^{-2{\kappa }^{2}{\left(1-n/2\right)}^{2}t/{t}_{{\nu }_{{\rm{in}}}}}\\ \,\times \,[{Y}_{l}(\kappa ){J}_{l}(\kappa {x}_{0})-{Y}_{l}(\kappa {x}_{0}){J}_{l}(\kappa )]\\ \,\times \,\displaystyle \frac{[{Y}_{l}(\kappa ){J}_{l}(\kappa x)-{Y}_{l}(\kappa x){J}_{l}(\kappa )]}{[{Y}_{l}^{2}(\kappa )+{J}_{l}^{2}(\kappa )]}\kappa \,d\kappa \end{array}\end{eqnarray} \tag{ 39 }$$