Additionally, some of the flares may arise in the secondary jet rather than in the primary jet. They could be recognized for their short variability time scale. However, since the secondary jet is expected to represent only a fraction of the flux of the primary jet, the fast variations would appear at a low amplitude, consistent with Valtaoja et al. [109] who found 1–10% peak-to-peak amplitude variations during the 1983/84 flaring season, with a 15.7 min period. Pihajoki et al. [89] find that variations down to 3.5 d may be associated with the primary while the variability of the secondary could be faster by at least a factor of 122, depending on the secondary spin. For a maximally spinning secondary the 15.7 min period may arise [110]. It is the period of the innermost stable orbit about the secondary black hole in the jet reference frame.

Since the secondary disk is likely to overflow its Roche lobe, we may also find interesting X-ray emission from the ring in the primary disk where the overflow stream terminates, akin to X-ray binaries. It is also worth searching for the direct accretion disk signature of the primary in the SED including a strong disk contribution in X-rays persistently detected in AGN with SMBHs of that mass.

**Figure 13.** Optical spectra of OJ 287 taken at the *Nordic Optical Telescope*. The continuum has been fitted with a low-order polynomial and subtracted and the spectra have been shifted in vertical direction for clarity. The positions of *λ*4861 H*β* and *λ*50006 [OIII] at the redshift of OJ 287 are indicated by vertical lines. The red curve below shows the positions of prominent atmospheric emission lines. Apparent deep absorption structures in the spectrum are telluric lines. The epochs are (from top to bottom): 12 February 2010, 23 October 2010, 18 November 2010, 4 January 2011, 17 January 2011, 1 October 2011, 11 February 2011, 23 April 2011, 7 May 2011, 16 November 2011, 19 February 2012, 27 December 2012, and 6 March 2013.

### **9. Are There More OJ 287-like Systems?**

The remaining lifetime of the OJ 287 black hole binary before a merger is 104 yr. Comparing this with the merger time scale of 107 yr, there is only a chance of ∼10−3 to catch the system at this late stage of evolution. Is this a problem? It depends on the parent population of OJ 287-like systems. As an order of magnitude estimate, it has been deduced that there are something like 104 binary black hole systems in the sky in the OJ 287 brightness category [111]. If this is the parent population, some ten of them would be found in the final stages of merger [112].

How could we find other OJ 287-like systems? We define an OJ 287-like binary black hole system such that it shows repeated bremsstrahlung flares due to the impacts on the accretion disk. The main requirement is that the mass ratio of the two black holes is rather large, so that the smaller black hole does not disrupt the accretion disk entirely in the process [51].

OJ 287 is a blazar with a jet pointing nearly toward us which increases its brightness many-fold, by about 3 magnitudes by our estimate above. The jets of other similar systems would likely point in some other direction. Thus it is unlikely that other OJ 287-like systems are found in lists of brightest radio sources, such as the Ohio Survey which gives the initial to the name of OJ 287.

However, if such a source is discovered through the repeated bremsstrahlung flares, it will take time before enough flares are seen to solve the orbit. In the case of OJ 287 we

are lucky that the source lies in the ecliptic zone which has been frequently photographed over 130 yrs. No such luck can be expected for the other OJ 287-like systems. It is also rather likely that these other systems contain smaller binary components by mass than OJ 287, with weaker bremsstrahlung flares. While other candidate binary black hole systems have been identified in recent years (review by [113]), none of them were explained by involving optical bremsstrahlung flares, and none of them so far has the dense polarimetry monitoring that OJ 287 has seen. Thus we may find more OJ 287-like systems, but not before full-sky monitoring programs have been operational for some time. We now briefly list GW aspects of our description for OJ 287.

### **10. Pulsar Timing Array Response to nHz GWs from OJ 287's SMBH Binary Engine**

The above listed observational campaigns in many electromagnetic observational windows clearly point to the presence of a nHz GW emitting SMBH binary in OJ 287. Therefore, it is important to develop prescriptions to model the expected pulsar timing array (PTA) response to nHz inspiral GWs from systems like OJ 287 [10]. The fact that the orbital eccentricity and primary BH spin can substantially influence the BH binary orbital evolution, as demonstrated in [33,114], prompted us to employ an improved version of GW phasing approach to model temporally evolving GW polarization states from such systems [115–117]. This is required as the PTA relevant GW-induced pre-fit pulsar timing residual, namely PTA response to inspiral GWs, induced by an isolated GW source at a given epoch *t* may be written as [118]

$$R(t) = \int\_0^t \left[ h(t' - \tau\_P) - h(t') \right] dt',\tag{9}$$

where *τP* is the geometric time delay in the solar system barycenter (SSB) frame of the underlying millisecond pulsar while *t* and *t* 0 are the SSB relevant coordinate time variables. Further, *h*(*t*) provides the effective GW strain induced by an isolated source at an epoch *t* and depends on the intrinsic parameters of the GW source like the total mass, spin values, orbital eccentricity and angles that specify its sky location, as detailed in [116]. In other words, PTA response of a pulsar to inspiral GWs from OJ 287 like system at a time *t* depends not only on the strength of the inspiral GWs at the SSB at *t* but also on the GW amplitude at location of the pulsar at time *t* − *τP*. This is another reason why we need to employ PN approximation to track the orbital evolution and the resulting GW polarization states while modeling PTA response to inspiral GWs from SMBH binaries [116].

In an ongoing effort, we are modeling temporally evolving quadrupolar order *h*(*t*) from comparable mass BH binaries inspiraling along PN-accurate eccentric orbits under the influence of dominant order spin-orbit interactions due to the primary BH spin. This is being done by adapting and extending the GW phasing efforts of Damour et al. [115], Susobhanan et al. [116]. It may be recalled that this approach usually employs Kepleriantype parametric solution to model the conservative dynamics of SMBH binaries in eccentric orbits. Thereafter, the effects of GW emission are incorporated by adapting the method of variation of constants, as detailed in [115]. In practice, we invoke the arguments of energy and angular momentum balance that lead to differential equations for orbital frequency and eccentricity with the help of orbital averaged PN-accurate far-zone energy and angular momentum fluxes [117].

In Figure 14, we plot the expected PTA response of PSR J0751 + 1807 to inpsiral nHz GWs associated with our SMBH binary central engine description for OJ 287 as detailed in [33]. The plot is based on Equation (9) where we have employed quadrupolar order *h*(*t*) expression for spinning BH binaries in eccentric orbits from Königsdörffer & Gopakumar [119] while the orbital parameters are as detailed in Dey et al. [33]. The sharp features depend on the orbital eccentricity and pulsar distance and provide unique signatures of the expected PTA response to nHz GWs from such SMBH binaries. Unfortunately, the expected timing residuals are substantially below what we can achieve with the current timing campaigns on PSR J0751 + 1807. However, it should

be possible to find and time milli-second pulsars with such exquisite timing precision during the SKA era [29]. Therefore, it should be possible to monitor persistent inspiral GWs from SMBH binaries like OJ 287 during the SKA era. This development should naturally allow us to pursue persistent multi-messenger GW astronomy with crucial inputs from the ongoing, planned and proposed multi-wavelength electromagnetic observations of OJ 287. This is because such observational campaigns should allow us to constrain various astrophysical parameters of the binary BH central engine description for the unique blazar OJ 287.

**Figure 14.** Timing residuals induced by inspiral GWs from OJ 287's SMBH binary on a typical PTA pulsar PSR J0751 + 1807. We let the pulsar distance to be 1.4 Kpc [120] while SMBH binary parameters as listed in Section 1.1 with its luminosity distance to be 1.6 Gpc. Timing of milli-second pulsars with SKA are expected to achieve similar timing residuals.

### **11. Discussion**

This article emphasizes the importance of continuing monitoring of OJ 287 at all wavelengths of the electromagnetic spectrum as well as in gravitational waves. Even though the binary orbit has now been solved [33] and confirmed [34], there are many aspects of the system which still need verification and from which we can learn interesting astrophysics.

(1) The confirmation of the 2022 optical flare: This flare does not appear in any of the simple schemes that have been proposed over the years as alternatives to the LV96 binary black hole model. Difficult as it is to observe from the ground, we should just be able to observe the tail end of this flare. Polarisation observations play an important role here since the tail end of the flare should very highly polarised [36,50,59].

(2) The detection of the signals from the 3 December 2021 disk impact: This is energetically a huge event which has at previous times been studied only by occasional photometric measurements. This time we expect a full coverage from optical polarisation to X-rays which are the expected key indicators of the impact. So far, in ongoing MOMO/Swift and ground-based monitoring observations, OJ 287 was found in a low-state in early December 2021 (Komossa et al., and Zola et al., in prep.). The timing of the impact also depends on the detailed model, which will be improved after the 2021/22 observing season.

(3) The *Hα* flares: This has been caught only once, accidentally by Sitko & Junkkarinen [65]. Later spectroscopic surveys, some of which are reported in this paper for the first time, have not been equally lucky. Now we give a rather precise prediction for the *Hα* flare which follows the 3 December 2021, disk impact.

(4) The secondary black hole: The secondary black hole has never been directly detected. One of the problems is the short separation between the two black holes in the sky, typically 10 microarcseconds, depending on its orbital phase in an eccentric orbit. It will be several years before this resolution is achieved in radio VLBI. However, an attempt was made to detect the secondary by its spectral lines, but so far we only report a null result. It is a difficult task due to the high background level of the jet component. New attempts

should be made during periods of low light level. Here we have made a prediction when such occasions may arise.

(5) Radio jet of the secondary: Even though the immediate surroundings of the secondary are hard to resolve, the current resolution in radio VLBI should be just enough to see the jet. A strong radio component was recently reported [70] at the expected position of the secondary jet [75]. It remains to be seen if this component behaves as expected. One problem is that the current direction of the primary jet is not very far from the suspected secondary jet. However, if the primary jet behaves according to our models [75], the primary jet will move away from the secondary jet, and the identification becomes more clear. Thus the continuing radio VLBI monitoring of OJ 287 is essential.

(6) Radio jet of the primary: The radio jet of OJ 287 has been mapped in radio VLBI since the 1980s. The origin of the jet and its connection to the central black hole and its accretion disk are still open problems. In the LV96 binary black hole model it is possible to calculate the direction of the spin axis of the primary as well as the normal to the innermost disk. They follow each other more or less [75], with the disk normal showing a little more short time scale wobble than the spin axis. If we make the reasonable assumption that the initial jet direction is determined either by the spin axis (spin model) or the disk normal (disk model), we obtain a unique model for the jet which wobbles predictably. When projected onto the sky, we obtain position angles at any point in time, and since the resolution at higher frequencies is better than at lower frequencies in the global VLBI, we get snapshots of the jet at different distances from the center. After solving for the unknown viewing angles, we get a full picture of the bulk flow in the jet.

(7) Optical polarisation angle: It is usually assumed that the optical emission arises in the jet closer to the black hole than the radio emission. We have applied the same model for the optical polarisation that works for the radio position angle, and find that the region of optical emission is about four times closer to the central black hole than the inner 86 GHz emission knot. It is interesting also that the optical polarisation gives the direction of the innermost jet. This is not surprising since the jet has been found to have helical magnetic field geometry [70]. When compressed along the jet axis in shock fronts, the field lies practically perpendicular to the jet axis and the polarisation is lined up with the jet. It is important to continue regular polarisation monitoring programs to learn more about the inner jet.

(8) Optical fades: The times of low light level in the jet are useful also for the detection of the host galaxy. The absolute magnitude of the host has been measured by surface photometry out to ∼30 kpc from the galaxy center, beyond which it is difficult to measure any light above the background generated by the bright point source, the AGN itself. Its brightness measured out to this distance is among the highest ever detected, and therefore it is likely that the total brightness of OJ 287 host is also among the highest. The brightness may also be measured by multicolor photometry at low light, because the color of the host is expected to be very different from the color of the jet. Thus the combined color of the source becomes more and more like galaxy colors when the source fades. This method suggests a brightness increase by one magnitude if the outer envelope of the host galaxy is included.

(9) High-energy emission and broad-band SEDs: Very dense, simultaneous X-ray–UV– optical monitoring has allowed us to identify the latest binary after-flares in 2017 and 2020, based on the brightest X-ray flares so far discovered from OJ 287 and other arguments [80]. Ongoing monitoring will reveal new outstanding X-ray–UV outbursts. Spectroscopy, SED modelling, cross-correlation analyses and structure function analyses [81,86] will allow us to probe in detail the binary and outburst physics across the electromagnetic spectrum, and will allow to search for an accretion disk contribution of the secondary at select epochs. Ongoing high-energy monitoring could also serve as trigger criterion for new VHE observations.

(10) Very high energy gamma rays: VHE gamma rays have been detected in OJ287 only once. This occurred after the two big tidal flares which followed the 2015 impact flare. Even though we do not know the mechanism which produced this event, it may useful

to monitor OJ 287 again after the tidal flares which follow the 2022 impact flare. More frequent monitoring would be also useful to see if the two phenomena are really connected.

(11) Fast variability: The shortest time scale of periodic or quasi-periodic variability that may be associated with the primary jet is 3.5 days. This has been seen in many occasions. The time scale of variability should be ∼122 times shorter in the secondary jet, due to the smaller mass of the secondary by this factor, or it could be even as low as ∼15.7 min if the secondary spins much faster than the primary. Thus the study of fast variability in radio, optical and X-rays may reveal details of the secondary jet which are otherwise hard to recover.

(12) Pulsar Timing Array: It will be important to develop accurate and efficient algorithms to model the expected PTA signals from SMBH binaries like OJ 287. Further, data analysis approaches that can efficiently look for such nHz GW sources in the PTA datasets of SKA era will need to be developed [10]. Additionally, it will be desirable to search for milli-second pulsars close to the Sky location of OJ 287 as they are expected to be more sensitive to inspiral GWs from OJ 287. Investigations should be pursued to figure out the possibility of tracking the trajectory of secondary BH during the ngEHT era.

(13) It has been argued that OJ 287 may produce a few muon neutrinos with energies more than 100 TeV in instruments like IceCube-Gen2 during ten years of Fermi flaring epochs, influenced by the modeling of the 2017 flare of TXS 0506 + 056 [121,122]. The neutrino studies will be an important part of multi-messenger astronomy with OJ 287 in the coming decades.

**Author Contributions:** M.J.V. and L.D. have been responsible for the main design of this article as well as for theoretical developments. L.D. and A.G. were in charge of the gravitational wave part, S.Z. of optical photometry, S.K. of the MOMO project and X-ray parts, T.P. and S.Z. of optical spectroscopy, R.H. of historical magnitude measurements, H.J. and A.V.B. of optical polarimetry and J.L.G. of radio mapping. All authors have read and agreed to the published version of the manuscript.

**Funding:** L.D. and A.G. acknowledge the support of the Department of Atomic Energy, Government of India, under project identification # RTI 4002. A.G. is grateful for the financial support and hospitality of the Pauli Center for Theoretical Studies and the University of Zurich. S.Z. acknowledge the grant NCN 2018/29/B/ST9/01793. R.H. acknowledges support by GACR grant 13-33324S and by the MSMT Mobility project 8J18AT036.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We thank Kari Nilsson for a valuable contribution in data reduction and producing one of the figures. We also appreciate the help that Vilppu Piirola has provided and continues to provide in optical polarisation observations of OJ 287. We also thank Gene Byrd and Shirin Haque for reading the manuscript and for valuable inputs.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

### *Article γ***-ray Flux and Spectral Variability of Blazar Ton 599 during Its 2021 Flare**

**Bhoomika Rajput \* and Ashwani Pandey**

Indian Institute of Astrophysics, Block II Koramangala, Bangalore 560034, India<EMAIL_ADDRESS>**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Blazars are known to emit exceptionally variable non-thermal emission over the wide range (from radio to *γ*-rays) of electromagnetic spectrum. We present here the results of our *γ*-ray flux and spectral variability study of the blazar Ton 599, which has been recently observed in the *γ*-ray flaring state. Using 0.1–300 GeV *γ*-ray data from the *Fermi* Gamma-ray Space Telescope (hereinafter *Fermi*), we generated one-day binned light curve of Ton 599 for a period of about one-year from MJD 59,093 to MJD 59,457. During this one year period, the maximum *γ*-ray flux detected was 2.24 ± 0.25 × 10−6 ph cm−2 s −1 at MJD 59,399.50. We identified three different flux states, namely, epoch A (quiescent), epoch B (pre-flare) and epoch C (main-flare). For each epoch, we calculated the *γ*-ray flux variability amplitude (F*var*) and found that the source showed largest flux variations in epoch C with F*var* ∼ 35%. We modelled the *γ*-ray spectra for each epoch and found that the Log-parabola model adequately describes the *γ*-ray spectra for all the three epochs. We estimated the size of the *γ*-ray emitting region as 1.03 × 1016 cm and determined that the origin of *γ*-ray radiation, during the main-flare, could be outside of the broad line region.

**Keywords:** galaxies; active; galaxies–quasars; individual (Ton 599)

### **1. Introduction**

Blazars are the jetted subclass of the Active Galactic Nuclei (AGNs) that are understood as AGNs with very small viewing angles to the line of sight [1,2]. Blazars emit broadband electromagnetic radiation that ranges from radio to extremely high *γ*-ray energies. The jets of blazars are highly Doppler boosted, resulting in flux variations over the entire accessible electromagnetic wavebands. Flat spectrum radio quasars (FSRQs) and BL Lacertae (BL Lacs) objects are the two subclass of blazars. The difference between these subclasses is determined by the equivalent width (EW) of the emission lines in their optical spectra, with FSRQs having EW > 5 Å and BL Lacs having EW < 5 Å [3]. A more physical criterion for distinguishing between FSRQs and BL Lacs was proposed by [4] which is based on the ratio of broad line region (BLR) luminosity (*LBLR*) to Eddington Luminosity (*LEdd*). FSRQs have the value *LBLR*/*LEdd* > 5 × 10−4 , whereas BL Lacs have the value *LBLR*/*LEdd* < 5 × 10−4 .

The broad-band spectral energy distributions (SEDs) of blazars comprise of two humps; the low energy hump and the high energy hump. The low energy hump peaks at optical/UV/X-ray region and the high energy hump peaks at MeV/GeV/TeV region [5–7]. The genesis of the low energy hump is well understood through synchrotron emission mechanism of the relativistic electrons, while the high energy hump originates through inverse Compton (IC) emission process [6]. Based on the position of the synchrotron peak frequency (*νsyn*), blazars are further classified as low synchrotron peaked (LSP; *νsyn* < 1014 Hz), intermediate synchrotron peaked (ISP; 1014 Hz < *νsyn* < 1015 Hz) and high synchrotron peaked (HSP; *νsyn* > 1015 Hz) blazars.

Due to the Doppler boosting in blazars' jets, the observed emission *Sobs* relative to the emission in the co-moving frame *Sint* is defined as e.g., [8]

$$S\_{\rm obs} = S\_{\rm int} \delta^q \tag{1}$$

Flux and Spectral Variability of Blazar Ton 599 during Its 2021 Flare. *Galaxies* **2021**, *9*, 118. https://doi.org/ 10.3390/galaxies9040118

**Citation:** Rajput, B.; Pandey, A. *γ*-ray

Academic Editor: Santanu Mondal

Received: 15 November 2021 Accepted: 10 December 2021 Published: 13 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

where *q* = 3 + *α* for a moving compact source and *q* = 2 + *α* for a stationary jet [9]. *α* is the spectral index which is defined as *fν* ∝ *ν* −*α* . *δ* is the Doppler boosting factor which is described as *δ* = 1/Γ(1 − *β*cos*θ*), where Γ is bulk Lorentz factor (Γ = 1/(1 − *β* 2 ) 1/2), *β* is the speed of jet in units of the speed of light and *θ* is the viewing angle between observer's line of sight and jet's axis. The observed time is also shortened by the effect of Doppler boosting by a factor *δ* −1 . These two consequences of the Doppler boosting increase the chances of observed variations in blazars over a wide range of wavelengths and make them the brightest objects in the extragalactic sky.

The study of flux variability is a valuable tool for determining the size of the emission zone in blazars. Blazars exhibit variability over a wide range of timescales, from few minutes to several years, across the full wavelength range i.e., from radio to *γ*-rays [10–14]. The variability in blazars can be explained by the shock-in-jet model [15]. In this scenario, the inhomogeneities in the jet flow produce relativistic shocks. These shocks travel through the jet plasma at relativistic speed and accelerate the particles. The other scenario is the magnetic reconnection, which is responsible for the rapid variations in blazar jets and has been investigated extensively in recent years [16–19]. Within the reconnection region, magnetic reconnection possibly enables small-scale ultra-relativistic flow, also called jet-injet scenario [16,19,20]. The ultra-relativistic motion of small plasmoids causes additional Doppler boosting and results in very short and bright flares.

The launch of *Fermi* Telescope in 2008 gave us an unprecedented opportunity to explore the *γ*-ray regime in blazars. Using *Fermi* data at *γ*-ray energies, flux doubling time scales have been detected as short as few minutes in both FSRQs and BL Lacs [14,21–23]. Causality considerations in these circumstances point to highly compact *γ*-ray emission zones in blazars' jets. The detection of *γ*-ray photons with energy greater than 10 GeV in FSRQs leads to the conclusion that the *γ*-ray emission region should be located at outside the cavity formed by the broad emission line (BLR) [24]. Otherwise, pair production on UV photons emitted by (BLR) clouds should severely attenuate the *γ*-ray photons. Despite several studies, the physics of *γ*-ray variability in blazars is still a captivating subject of research.

Ton 599 (4FGL J1159.5 + 2914; [25]) is an FSRQ, located at redshift z = 0.725 [26] with R.A. = 11*h*59*m*31.8*s* and Dec. = +29◦14043.800. Ton 599 is a strongly polarized and a highly optically violent variable quasar [27]. It was first detected in high energy *γ*rays by the Energetic Gamma Ray Experiment Telescope (EGRET) in the second EGRET catalogue [28] and later by the *Fermi* Large Area Telescope (*Fermi*-LAT) [29]. The correlation study between radio and *γ*-ray bands was carried out by [30] for this source to put the constraints on *γ*-ray emission region in parsec-scale jets. In 2017, for the first time, Ton 599 went through a protracted flaring condition spanning the full electromagnetic spectrum. During this flare, the maximum *γ*-ray flux observed was 1.26 × 10−6 ph cm−2 s −1 [31]. The detailed study of *γ*-ray flux variability during this flare was carried out by [31,32]. A multiwavelength study and the broad-band SED modelling of this flare was conducted by [33] using a two-component leptonic emission model. An EC mechanism with a dusty torus (DT) photon field producing seed photons was identified to be responsible for the GeV emission in this study. Ton 599 recently displayed a bright flare from July to September 2021, allowing us to investigate its *γ*-ray emission process and its consequences during this source's flaring condition.

In this paper, we describe the *γ*-ray analysis performed on the source Ton 599 utilising data collected over a one-year period from September 2020 to August 2021, with the goal of constraining the *γ*-ray emission region in blazars. The data used in this study is presented in Section 2. Section 3 describes the *γ*-ray light curve, followed by results in Section 4 and discussion in Section 5. Section 6 presents a summary of the results.

### **2.** *γ***-ray Observations and Data Reduction**

In this work, we used the *γ*-ray observations of Ton 599 taken by the *Fermi*-LAT for a period of ∼1 year, from 2020 September to 2021 August (MJD: 59,093–59,457). *Fermi*-LAT is a pair-conversion *γ*-ray telescope that can detect *γ* rays with energies ranging from 20 MeV to more than 1 TeV. It has a large field of view of about 2.4 sr and scans the full sky in every 3 h, covering 20% of the sky at any time [34]. We used the package ScienceTools v1.2.23 with the instrument response function P8R3\_SOURCE\_V31 for our analysis. We used the latest LAT Pass 8 data in the energy range 100 MeV to 300 GeV, where the photon-like events are classified as 'evclass = 128, evtype = 3'. The region of interest (ROI) is specified as a circle with a radius of 10◦ and is centred on the *γ*-ray position of the source. We used a maximum zenith angle of 90◦ to remove *γ*-ray contamination from the earth's limb. The latest isotropic model 'iso\_P8R3\_SOURCE\_V3\_v1' and the Galactic diffuse emission model 'gll\_iem\_v07' were used to analyze the data. The recommended condition '(DATA QUAL > 0)&&(LAT CONFIG= = 1)' was then used to construct the required good time intervals. Over the time period of interest, an unbinned likelihood analysis is performed to generate the 1-day binned *γ*-ray light curve of Ton 599. In the light curve, the source was considered to be detected if the test statistics (TS) > 9, which corresponds to a 3*σ* detection [35]. Our final one-day binned *γ*-ray light curve contains 256 confirmed measurements of Ton 599.

### **3.** *γ***-ray Light Curve**

The one-day binned *γ*-ray light curve of Ton 599 from 2020 September to 2021 August (MJD: 59,093–59,457) is shown in Figure 1. For the entire period of the light curve, we estimated the average *γ*-ray flux which was found to be 0.43 ± 0.36 × 10−6 ph cm−2 s −1 . Based on the average *γ*-ray flux, we visually identified three different flux states, namely, quiescent state, pre-flaring state, and the main flaring state, which were labelled as epoch A, epoch B, and epoch C, respectively. The duration of these epochs are marked by vertical lines in Figure 1 and the details of these epochs are given in Table 1. We classified epoch A as a quiescent epoch since the flux was lower than the average flux for the whole time span. During epochs B and C, the flux increased 2–3 times than the average flux, so we classified these epochs as flaring epochs. Over the whole time span of the light curve, the value of the flux was maximum during the flaring epoch C. The value of the maximum flux was determined to be 2.24 ± 0.25 × 10−6 ph cm−2 s −1 at MJD 59,399.50.

**Table 1.** The epochs considered for the *γ*-ray light curve study are listed in the table below. Here, *N*, *x*¯, and *Fvar* denote the total number of data points, the average *γ*-ray flux, and the fractional variability amplitude, respectively, for the epoch. The *γ*-ray flux is in the units of 10−6 ph cm−2 s −1 .

**Figure 1.** One-day binned *γ*-ray light curve of Ton 599 (**upper panel**) with the TS values corresponding to each flux value (**lower panel**) . The vertical black lines refer to the quiescent epoch (Epoch A), whereas vertical blue lines refer to the *γ*-ray flaring epochs (Epoch B and Epoch C). In the upper panel the horizontal cyan line represents the average flux from September 2020 to August 2021.

### **4. Results**

### *4.1. γ-ray Flux Variability*

We used the fractional variability amplitude (*Fvar*) to characterize the flux variability properties of Ton 599. The *Fvar* is commonly used to quantify the intrinsic variations in blazar light curves and is defined as, e.g., [12,36,37]

$$F\_{\text{var}} = \sqrt{\frac{S^2 - \overline{\sigma\_{err}^2}}{\mathfrak{x}^2}} \tag{2}$$

In the above equation *S* 2 is the sample variance and *σ* 2 *err* is mean square error. These are defined as

$$S^2 = \frac{1}{N-1} \sum\_{i=1}^{N} (\mathbf{x}\_i - \bar{\mathbf{x}})^2 \tag{3}$$

and

$$
\overline{\sigma\_{err}^2} = \frac{1}{N} \sum\_{i=1}^{N} \sigma\_{err,i}^2 \tag{4}
$$

Here *N* is the total number of data points in the light curve and *x*¯ is the arithmetic mean of the light curve.

The uncertainty in *Fvar* is given by the following equation

$$err(F\_{var}) = \sqrt{\left(\sqrt{\frac{1}{2N}} \frac{\overline{\sigma\_{err}^2}}{\overline{\pi}^2 F\_{var}}\right)^2 + \left(\sqrt{\frac{\overline{\sigma\_{err}^2}}{N} \frac{1}{\overline{\pi}}}\right)^2} \tag{5}$$

For each epoch, we calculated the value of *Fvar* separately by considering the average flux (*x*¯) and the data points (*N*) within the epoch only. The values of *N*, *x*¯, and *Fvar* for each epoch are listed in Table 1. The errors calculated using Equation (5) are given in Table 1. For epoch A, the error in the value of *Fvar* is relatively large which also indicates that it is a quiescent state. We considered the source to be variable in epochs if *Fvar* > 3 × *err*(*Fvar*). According to our criteria for variations, the source was not variable in epoch A (quiescent state), while significant *γ*-ray flux variations were noticed during epochs B and C with *Fvar* values of 22% and 35%, respectively.

### *4.2. γ-ray Spectral Fitting*

We modeled the *γ*-ray spectrum of Ton 599 for each epoch to investigate the inherent distribution of electrons that causes the *γ*-ray emission during the epoch. We used the power-law (PL) and the log parabola (LP) models to fit the *γ*-ray spectra of Ton 599 using the maximum likelihood analysis. The PL model is defined as follows [38]:

$$\frac{dN(E)}{dE} = N\_\circ \left(\frac{E}{E\_\circ}\right)^\Gamma \tag{6}$$

where *N*◦ is the prefactor (normalization of the energy spectrum), *E*◦ is the pivot energy (∼523.51 MeV) given in the 4FGL catalog, and Γ is the photon index, while, the LP model is given by

$$\frac{dN(E)}{dE} = N\_\circ \left(\frac{E}{E\_b}\right)^{-\alpha - \beta \ln\left(\frac{E}{E\_b}\right)}\tag{7}$$

In this equation *N*◦ is the normalization, *Eb* is the pivot energy (same as in PL), *α* is photon index at *Eb* and *β* is the curvature index that defines the curvature of the spectrum around the peak.

For each epoch, the model fitted *γ*-ray spectra of Ton 599 is shown in Figure 2, and the values of best fitted model parameters are given in Table 2. In Figure 2, the uncertainties are large at higher energies because of relatively low photon counts at these energies.

**Table 2.** Details of the PL and LP model fits for the three epochs. Here, *γ*-ray flux is in the units of 10−6 ph cm−2 s −1 .

To determine whether the *γ*-ray spectrum has a curvature or not, and which model (LP or PL) best describes the *γ*-ray spectrum of Ton 599, we calculated *TScurve* (curvature of the test statistics) [38], which is defined as:

$$TS\_{curve} = \mathcal{Z}(\log L\_{LP} - \log L\_{PL}) \tag{8}$$

In the above equation, *L* represents the likelihood function. The value of *TScurve* for each epoch is given in Table 2. We employed the *TScurve* > 16 threshold (i.e., 4*σ* level; [35]) to determine the presence of statistically significant curvature in the *γ*-ray spectrum. We found that the LP model best describes the *γ*-ray spectra of Ton 599 for all three epochs.

**Figure 2.** *γ*-ray spectrum for the epochs (**A**–**C**). The name of the epoch is mentioned in each plot.

### *4.3. Location of the γ-ray Emission Region*

To determine the size of the *γ*-ray emitting region, the flux doubling time scale is usually estimated, e.g., [31,39]. We calculated the flux doubling time scale for the one-day binned *γ*-ray light curve during the main-flaring epoch (epoch C) of Ton 599 as follows:

$$F(t\_2) = F(t\_1) \times 2^{\Delta t/\tau\_d} \tag{9}$$

In the above expression, *F*(*t*1) and *F*(*t*2) are the flux values at times *t*1 and *t*2 respectively, ∆*t* = *t*2 − *t*1 and *τd* represents the flux doubling time.

We found a flux doubling timescale of ∼13.2 h during the epoch C for the blazar Ton 599. Using the flux doubling time scale, we estimated the size of the *γ*-ray emitting region for Ton 599 using the following expression:

$$r \le c\tau\_d \delta / (1+z) \tag{10}$$

A gamma-ray Doppler factor of delta = 12.5 was calculated for Ton 599 by [40] using the multiwavelength data and the known radio Doppler factors with the assumption that the observed boosted emission is from the SSC model. This value is consistent with the lower limits (10.75 and 13.45) of delta obtained for the two bright flares of Ton 599 by [32]. Using *δ* = 12.5, the size of the *γ*-ray emission region is estimated to be 1.03 × 1016 cm.

According to [41], the location of the *γ*-ray emission region with respect to the central super massive black hole (SMBH) can be approximated by using the relation R = r/*φ*, where r is the size of the *γ*-ray emitting blob and *φ* is the jet opening angle.

An intrinsic opening angle of 0.58◦ (∼0.01 radian) was determined for Ton 599 by [42]. They estimated the jet opening angles for Fermi detected AGNs using 15.4 GHz Very Long

Baseline Array (VLBA) observations following two different methods: (a) by analyzing transverse jet profiles in the image plane and (b) by model fitting the data in the (u, v) plane. Using these values of *r* and *φ*, the location of the *γ*-ray emitting region for Ton 599 was found to be at 1.03 × 1018 cm from the SMBH.

### **5. Discussion**

In this section we give our interpretation of the findings of above-mentioned analyses, as well as a discussion of them.

### *5.1. γ-ray Flux Variability*

In the *γ*-ray band, blazars exhibit remarkable flux variability. For all the three epochs of the source Ton 599, we estimated the flux variability amplitude using one-day binned light curves. We found that during the main-flaring epoch (epoch C), the source showed largest variations with *Fvar* = 0.35 ± 0.03. During the pre-flaring epoch (epoch B), the source was also variable with *Fvar* = 0.22 ± 0.04 and during the quiescent epoch (epoch A), the source was not significantly variable (*Fvar* = 0.36 ± 0.12) within the error bar. During the main-flare (epoch C), the source showed a maximum flux of 2.24 ± 0.25 × 10−6 ph cm−2 s −1 at MJD 59,399.50, which is larger than the maximum flux observed during the 2017 flare [31]. We found that the photon index value is 1.99 ± 0.10 during the main-flaring epoch (epoch C), when flux is highest, harder than the 4FGL value of 2.19 ± 0.01.

The observed *γ*-ray variability in the blazars could be attributed to both intrinsic and extrinsic factors. The distribution of relativistic electrons responsible for the emission is one of the intrinsic effects. These relativistic electrons, which can be accelerated to a Lorentz factor of up to 106 , are responsible for non-thermal emission from blazar jets via synchrotron or inverse Compton emission processes. The distribution of these relativistic electrons and seed photon responsible for the inverse Compton emission process (synchrotron self Compton (SSC); [15,43] and external Compton (EC); [44,45]) are the intrinsic factors that cause short-term variability in jets. Extrinsic factors, such as the moving blob's high bulk Lorentz factor (Γ) ∼ 50, are in addition to intrinsic factors. Though the large bulk Lorentz factor is a favourable explanation in the case of BL Lacs, because the seed photons come from inside the jets and produce *γ*-rays through inverse Compton. On the other hand, the seed photons in FSRQ, emanate from outside the jets and produce *γ*-rays. However, pair-production through dense broad line region (BLR) can attenuate *γ*-ray emission. As a result, in the case of FSRQs, the large bulk Lorentz factor that causes *γ*-ray variability is not a plausible explanation [46]. The *γ*-ray flux variability in blazars jets can be explained by the magnetic reconnection (jet-in-jet) scenario, in which the mini jets generated in the jets can have a large bulk Lorentz factor and they can produce *γ*-ray flux variability on short time scale [16,19].

### *5.2. γ-ray Spectra*

For each epoch, we modelled the *γ*-ray spectra of Ton 599 using the power law and the log parabola models. We used the *TScurve* value to determine the best fitted model. Based on the *TScurve* value, we found that the LP model best describes the *γ*-ray spectra of Ton 599 during all the three epochs indicating that the GeV *γ*-ray spectrum of Ton 599 is curved.

FSRQs usually have a curved GeV *γ*-ray spectrum, e.g., [47–49]. Several theories have been proposed, in the literature, to explain the curved *γ*-ray spectra of FSRQs. The attenuation of *γ*-ray photons by pair-production (inside the BLR) could explain the curved *γ*-ray spectrum [50]. A viable scenario for defining the curve in the *γ*-ray spectrum is the Klein-Nishina effect on the inverse Compton scattering of BLR photons through relativistic electrons present in the jet [51]. However, several investigations conducted to locate the *γ*-ray emission region discovered that the observed *γ*-ray spectra in FSRQs are not caused by IC scattering of BLR photons and that the *γ*-ray emission site is located outside the BLR [19,52]. The curved *γ*-ray spectra of FSRQs could also be explained by the curved energy distribution of the electrons emitting the radiation [53,54].

### *5.3. Location of the γ-ray Emitting Region*

We calculated the minimum *γ*-ray flux doubling timescale of 13.2 h during epoch C for Ton 599. Using this minimum doubling timescale, we estimated the size of the *γ*-ray emitting region as 1.03 × 1016 cm. We have also identified the location of the *γ*-ray emitting blob to be at a distance of 1.03 × 1018 cm from the SMBH. The size of the BLR for TON 599 was estimated as ∼(2.11–2.45) × 1017 cm by [32]. So, the location of the *γ*-ray emitting blob, found in this work, is outside of the BLR. During the 2017 *γ*-ray flare of this source, using two component leptonic model, Ref. [33] found that the seed photons for the GeV emission produced by the dusty torus (DT) and the GeV emitting region is located outside of the BLR. Our results are consistent with the findings of [33].

### **6. Summary**

We investigated the *γ*-ray flux and spectral variability of the blazar Ton 599 during MJD 59,093 to MJD 59,457. For our study, we have chosen three epochs of different flux states: quiescent, pre-flaring, and main-flaring. The outcomes of the analysis of these epochs are summarised below.

**Author Contributions:** Conceptualization, B.R.; Data curation, B.R. and A.P.; Formal analysis, B.R.; Investigation, B.R.; Project administration, A.P.; Resources, B.R.; Visualization, B.R.; Writing—original draft, B.R.; Writing—review & editing, A.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data used in this work is publicly available from the Fermi-LAT data archive https://fermi.gsfc.nasa.gov/ssc/data/access/, accessed on 1 September 2021.

**Acknowledgments:** We thank the anonymous referees for their constructive comments, which helped us to improve the manuscript. This publication made use of data from the *Fermi* Gamma-ray Space Telescope, which accessed from the Fermi Science Support Center. We also acknowledge the use of High Performance Computing Facility (Nova cluster) at Indian Institute of Astrophysics.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Note**

1 https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/, accessed on 1 September 2021.

### **References**

### *Article* **The Radiative Newtonian 1** < *γ* ≤ **1.66 and the Paczy ´nski–Wiita** *γ* = **5/3 Regime of Non-Isothermal Bondi Accretion onto a Massive Black Hole with an Accretion Disc**

**Jose M. Ramírez-Velásquez 1,**∗ **, Leonardo Di G. Sigalotti 2 , Ruslan Gabbasov 3 , Jaime Klapp 4 and Ernesto Contreras 5**

**Abstract:** We investigate the non-isothermal Bondi accretion onto a supermassive black hole (SMBH) for the unexplored case when the adiabatic index is varied in the interval 1 < *γ* ≤ 1.66 and for the Paczy´nski–Wiita *γ* = 5/3 regime, including the effects of X-ray heating and radiation force due to electron scattering and spectral lines. The X-ray/central object radiation is assumed to be isotropic, while the UV emission from the accretion disc is assumed to have an angular dependence. This allows us to build streamlines in any desired angular direction. The effects of both types of radiation on the accretion dynamics is evaluated with and without the effects of spectral line driving. Under line driving (and for the studied angles), when the UV flux dominates over the X-ray heating, with a fraction of UV photons going from 80% to 95%, and *γ* varies from 1.66 to 1.1, the inflow close to the gravitational source becomes more supersonic and the volume occupied by the supersonic inflow becomes larger. This property is also seen when this fraction goes from 50% to 80%. The underestimation of the Bondi radius close to the centre increases with increasing *γ*, while the central overestimation of the accretion rates decreases with increasing *γ*, for all the six studied cases.

 Introduction**

Observations of giant elliptic galaxies provide firm evidence of the presence of supermassive black holes (SMBHs) in their centres [1,2], which accrete matter from the surroundings and liberate enormous amounts of energy that affect their environments from pc to Mpc scales (see [3] for spherically symmetric black holes with quantum corrections). In particular, the mass accretion rate onto these SMBHs, which is a key quantity to understand the galactic evolution, is usually estimated using Bondi accretion theory [4]. However, semi-analytic calculations and numerical simulations based on the Bondi accretion model do not provide information about the flow transport down to pc and sub-pc scales and the mass accretion rates (e.g., [5–8]). Since active galactic nuclei (AGNs) evolve concomitantly with their host galaxies, they affect each other. For example, good evidence for this kind of feedback has been provided by observations of AGN-starbursts [9]. On the other hand, the energy released during accretion onto a SMBH can hinder further accretion and drive the gas away, which in turn self-regulates the galaxy growth [10].

**Citation:** Ramírez-Velásquez, J.M.; Sigalotti, L.D.G.; Gabbasov, R.; Klapp, J.; Contreras, E. The Radiative Newtonian 1 < *γ* ≤ 1.66 and the Paczy ´nski–Wiita *γ* = 5/3 Regime of Non-Isothermal Bondi Accretion onto a Massive Black Hole with an Accretion Disc. *Galaxies* **2021**, *9*, 55. https://doi.org/10.3390/galaxies 9030055

Academic Editor: Dimitris M. Christodoulou

Received: 11 May 2021 Accepted: 9 August 2021 Published: 11 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

An earlier study on the structure of X-ray-irradiated accretion discs in AGNs was considered by the authors of [11]. They found that the irradiated region above and below the disc consists of a region which is supported by the radiation pressure where the UV flux is created and a warmer thin layer above this region, which is optically thin to the UV radiation. In the last years, studies of flows in AGNs with the inclusion of radiation source terms have become progressively more consistent [12–17]. All of these studies were able to compute cooling and heating functions, which are important to produce the correct opacities in these environments, and new methodologies were proposed to properly include the coupling between matter and radiation. Winds and accretion processes are of primary importance to improve our understanding of the feedback between the galaxy and its guest SMBH [18–23] and therefore on the galactic evolution [24–27]. Photoionization calculations of radiative forces due to spectral lines (i.e., [15]) have also shown the importance of a proper treatment of the coupling between matter and radiation along with the influence of the non-LTE effects (i.e., [28]).

Much effort has been devoted to study the dynamical evolution of the system in terms of the accretion rates onto a SMBH from the numerical solution of the hydrodynamics equations, where a common assumption has been to fix the boundary conditions at infinity as it is indeed required by the classical Bondi solution. However, an inconvenience with this approach is that when exploring the dynamics close to the black hole, this boundary conditions may fail to represent the finiteness of the spatial region under study. Several calculations of accretion processes with an ideal equation of state with different values of the adiabatic index (*γ*; in full symmetric and axisymmetric configurations in relativistic contexts e.g.; [29–35]) show it playing an important role in the estimation of the mass supply in cosmological simulations [21,36,37], and so does the exploration of the effects of varying the adiabatic index on the AGN accretion dynamics. In this paper we extend the analysis of [38] and present solutions for the radial Mach number and density profiles along with the critical points of transonic solutions for the radiative, radial accretion of matter in the potential well of a SMBH for adiabatic indices in the range 1.1 ≤ *γ* ≤ 1.66. It is found that variations of the adiabatic index may change the whole dynamics of the system as well as the mass supply. The methodology employed is detailed in Section 2 and the results are described in Section 3. A catalogue of pure absorption lines for 1.1 ≤ *γ* ≤ 1.66 is given in Section 4. Estimates of the Bondi radius and mass accretion rate as functions of the adiabatic index are given in Section 5. Finally, Section 6 summarizes the relevant conclusions.

### **2. Non-Isothermal Radial Bondi Accretion**

Under the effects of irradiation by X rays, the accretion flow onto the central SMBH is described by the mass and momentum conservation laws

$$\frac{d\rho}{dt} = -\rho \nabla \cdot \mathbf{v}\_{\prime} \tag{1}$$

$$\frac{d\mathbf{v}}{dt} = -\frac{1}{\rho}\nabla p + \mathbf{g} + \mathbf{F}^{\text{rad}},\tag{2}$$

where *ρ* denotes the density, **v** the velocity field, *p* the gas pressure, **g** the gravitational acceleration due to the SMBH, **F** rad = (*Fr* , *Fθ* = 0, *Fφ* = 0) the radiation force per unit mass, and *d*/*dt* = *∂*/*∂t* + **v** · ∇ the material time derivative. Using spherical coordinates and assuming that the angular components of the velocity vanish (i.e., *vθ* = *vφ* = 0), Equation (2) can be written as

$$\frac{\partial v\_r}{\partial t} + v\_r \frac{\partial v\_r}{\partial r} = -\frac{1}{\rho} \frac{\partial p}{\partial r} - \frac{\partial}{\partial r} \psi\_{\text{grav}}(r) + \frac{\partial}{\partial r} \psi\_{\text{rad}}(r), \tag{3}$$

for the radial velocity component, where *vr* = *vr*(*r*, *θ*, *φ*), *p* = *p*(*r*, *θ*, *φ*), *M*BH is the mass of the SMBH, and *ψ*rad(*r*) is a function defined below. If we further assume azimuthal symmetry (i.e., *∂*/*∂φ* = 0) and perform the analysis for a fixed angle *θ*0, such that *vr*(*r*, *θ*, *φ*) = *vr*(*r*, *θ*0) = *vr* and *p*(*r*, *θ*, *φ*) = *p*(*r*, *θ*0) = *p*, the above equation becomes

$$\frac{d\mathcal{H}(r)}{dr} = p\frac{d}{dr}\left(\frac{1}{\rho}\right) - \frac{\partial v\_r}{\partial t},\tag{4}$$

where

$$\mathcal{H}(r) = \frac{p}{\rho} + \frac{1}{2}v\_r^2 - \psi\_{\text{grav}}(r) - \psi\_{\text{rad}}(r) \tag{5}$$

is the Bernoulli function H(*r*), *ψ*grav = (*GM*BH)/*r* and *ψ*rad = −*C*(*r*)/*r* are the gravitational and radiation potentials, respectively. Compared to the analysis of [38], a radial dependence is now assumed for the term *C* in the definition of the radiation potential. For simplicity the accretion disc is assumed to be flat, Keplerian, geometrically thin and optically thick. The pressure is related to the density by means of an ideal equation of state

$$p = \frac{k\_B \rho T}{\mu m\_p} = p\_\infty \mathfrak{p}^\gamma,\tag{6}$$

where *kB* is the Boltzmann constant, *T* is the gas temperature, *µ* is the mean molecular weight, and *mp* is the proton mass. Here *γ* is varied in the range 1 ≤ *γ* < 5/3 and *ρ*˜ ≡ *ρ*/*ρ*∞, where *p*∞ and *ρ*∞ are the values of the pressure and density at infinity, respectively. The sound speed is defined according to

$$c\_s^2 = \frac{\gamma p}{\rho}.\tag{7}$$

The radiative contributions in Equations (2) and (4) are here implemented using the strategies used by [13,39–41] (hereafter P07). Although the radiation field from both the disc and the central black hole are modelled as in P07 by adding a radiation force in the momentum equation, there are important differences with P07 that are worth commenting. For instance, in P07 the analysis is based on time-dependent, axisymmetric simulations with *γ* = 5/3, while our analysis consists of steady-state calculations along radial streamlines for fixed *θ* and varying values of the adiabatic index in the interval 1 . *γ* < 5/3. These calculations have been designed to serve as initial conditions for fully three-dimensional simulations of accretion discs [42]. As in P07, the disc is assumed to emit only UV light (i.e., *f*disc = *f*UV) and the central SMBH to emit only X-rays (i.e., *f*? = *f*X = 1 − *f*disc) so that the disc and central luminosities are given by *L*disc = *f*disc*L* and *L*? = *f*?*L*, respectively, where *L* is the total accretion luminosity. The radial component of the radiation force is approximated according to the relation (e.g., P07)

$$F\_{r\_{\rm rad}}(r, \theta\_0) = \frac{\sigma\_T L}{4\pi r^2 c m\_p} [f\_\star + 2 \cos \theta\_0 f\_{\rm disk} (1 + M(r, t))] \,\tag{8}$$

where *c* is the speed of light in vacuum, *σT*/*mp* is the mass scattering coefficient for free electrons, *θ*0 is some fixed and constant polar angle measured from the rotational axis of the disk, and *t* is the optical depth

$$t = \frac{(\sigma\_T / m\_p) \rho v\_{\rm th}}{|dv\_r / dr|},\tag{9}$$

where *v*th is the thermal velocity (set using a temperature of 25,000 K), *dvr*/*dr* is the velocity gradient along the radial direction, and *M*(*r*, *t*) is the so-called force multiplier [43], defined by

$$M(r,t) = kt^{-\alpha} \left[ \frac{(1+\tau\_{\text{max}})^{(1-\alpha)} - 1}{\tau\_{\text{max}}^{(1-\alpha)}} \right],\tag{10}$$

where *α* = 0.6 is the ratio of optically thick to optically thin lines. According to [44], the parameter *k* is given by the minimun between

$$k = 0.03 + 0.385 \exp\left(-1.4 \xi^{0.6}\right),\tag{11}$$

and

$$\log k = \begin{cases} -0.383, & \text{for } \log T \le 4, \\ -0.630 \log T + 2.138, & \text{for } 4 < \log T \le 4.75, \\ -0.870 \log T + 17.528, & \text{for } \log T > 4.75, \end{cases} \tag{12}$$

which is based on detailed photoionization calculations performed using the XSTAR code (P07), while *τ*max = *tη*max with

$$\log \eta\_{\text{max}} = \begin{cases} 6.9 \exp\left(0.16 \xi^{0.4}\right), & \text{if } \log \xi \le 0.5, \\ 9.1 \exp\left(-7.96 \times 10^{-3} \xi\right), & \text{if } \log \xi > 0.5. \end{cases} \tag{13}$$

where *ξ* is the photoionization parameter. The parameter *η*max defined by Equation (13) is used to determined the maximum force multiplier, i.e., *M*max = *k*(1 − *α*)*η α* max. The local X-ray flux

$$\mathcal{F}\_X = \frac{L\_\star}{4\pi r^2} \exp(-\tau\_X),\tag{14}$$

is used to estimate the photoionization parameter *ξ* = 4*π*F*X*/*n*, where

$$
\pi\_{\mathbf{X}} = \int\_0^r \kappa\_{\mathbf{X}} \rho dr,\tag{15}
$$

is the X-ray optical depth, evaluated between the centre (*r* = 0) and a radius *r* in the accreting flow, and *n* = *ρ*/(*µmp*) is the number density of the accreting gas. The absorption coefficient *κX* is set to 0.4 g−1 cm2 for all values of *ξ*.

The *γ* = 5/3 regime of non-isothermal Bondi accretion is modelled using the Paczy ´nski– Wiita (PW) potential [45,46]

$$
\psi\_{\rm grav} = \frac{GM\_{\rm BH}}{r - R\_s}'\tag{16}
$$

where *Rs* = 2*GM*BH/*c* 2 is the gravitational radius of the black hole. While this pseudo-Newtonian potential does not obey the Poisson equation, it has become a standard tool because it accurately models the general relativistic accretion of matter onto a nonrotating SMBH. In particular, this form of the potential reproduces the radii of a marginally stable Keplerian orbit (*r* = 3*Rs*) and of a marginally bound orbit (*r* = 2*Rs*) as predicted by Einstein's gravity in the Schwarzschild metric [46].

Under the assumption of steady-state motion, Equations (4) and (5) can be combined to produce after a few algebraic steps the integral equation

$$\int \frac{d}{dr} \left( \frac{\mathcal{M}^2 \tilde{c}\_s^2}{2} + \frac{\tilde{\rho}^{(\gamma - 1)}}{\gamma - 1} - \frac{1}{\mathfrak{x}} + \frac{I\_{\text{tot}}^{\text{rad}}|\_{\theta = \theta\_0}(\mathfrak{x})}{\mathfrak{x}} \right) dr = 0,\tag{17}$$

where

$$\infty \equiv \frac{r}{r\_B}, \quad \mathcal{E}\_\mathbf{s} \equiv \frac{c\_\mathbf{s}}{c\_\infty} = \not{\mathfrak{p}}^{(\gamma - 1)/2}, \quad \mathcal{M} \equiv \frac{v\_r}{c\_\mathbf{s}} \,\tag{18}$$

and

$$l\_{\rm tot}^{\rm rad}|\_{\theta = \theta\_0}(\boldsymbol{\pi}) = l\_{\rm Edd}^{\rm rad} f\_{\rm rad}|\_{\theta = \theta\_0}(\boldsymbol{\pi}). \tag{19}$$

Here *C*|*θ*=*θ*0 (*x*) = *GM*BH*l* rad tot |*θ*=*θ*0 (*x*),

$$r\_B = \frac{GM\_{\rm BH}}{c\_{\infty}^2} \,\prime \tag{20}$$

is the Bondi radius, *c* 2 ∞ = *γp*∞/*ρ*∞, M is the Mach number, *l* rad Edd = *L*/*L*Edd with *L*Edd = 4*πcGM*BH*mp*/*σT* being the Eddington luminosity, *σT* = 6.6524 × 10−25 cm2 is the Thomson cross section, and *f*rad|*θ*=*θ*0 (*x*) is the radiative force parameter given by the relation

$$f\_{\rm rad}|\_{\theta=\theta\_0}(\mathbf{x}) = f\_\star + 2\cos\theta\_0 f\_{\rm disk} [1 + M(\mathbf{x}, t)].\tag{21}$$

Although the force multiplier is an explicit function of *x* and *t*, for simplicity hereafter we omit the *x*-dependence and write it as *M*(*t* = *τ*).

Equation (17) can be integrated to give

$$\rho^{(\gamma - 1)} \left( \frac{\mathcal{M}^2}{2} + \frac{1}{\gamma - 1} \right) = \frac{1}{\mathfrak{x}} - \frac{l\_{\text{tot}}^{\text{rad}}|\_{\theta = \theta\_0}(\mathfrak{x})}{\mathfrak{x}} + \frac{1}{\gamma - 1}. \tag{22}$$

This equation describes the non-isothermal (*γ* > 1), steady-state accretion along a radial streamline for a fixed angle *θ* = *θ*0 when the effects of radiation emission due to electron scattering and spectral discrete lines with appropriate boundary conditions at infinity and ionization changes due to temperature corrections are considered. A sketch showing the geometry of the system is displayed in Figure 1, where a radial streamline at a meridional angle *θ*0 is shown. In terms of the above assumptions and the normalized parameters (18), Equation (1) becomes

$$
\alpha^2 \mathcal{M} \mathfrak{gl}^{(\gamma+1)/2} = \lambda\_\prime \tag{23}
$$

where *λ* is the accretion parameter given by

$$
\lambda = \frac{\dot{M}\_B}{(4\pi f\_{\text{solid}})r\_B^2 \rho\_{\infty} c\_{\infty}},
\tag{24}
$$

which determines the accretion rate for given boundary conditions and mass of the SMBH. In this expression, 4*π f*solid = R sin *θdθdφ* is the solid angle covered by the streamline at the polar angle *θ*0 (see Figure 1). When *f*solid = 1 (full solid angle), Equation (24) reduces to the classical accretion parameter. However, for a *θ*0-dependent force *f*solid << 1. This dependence is important only for the final calculation of *M*˙ *B* = *λ*(4*π f*solid)*r* 2 *B ρ*∞*c*∞, which is fixed once the value of *θ*0 is chosen. Replacing *ρ*˜ from Equation (23) into Equation (22), the radiative-radial Bondi problem reduces to solving the equation

$$\mathbf{g}(\mathcal{M}) = \Lambda f(\mathbf{x}), \quad \text{with} \quad \Lambda = \lambda^{2(1-\gamma)/(\gamma+1)}, \tag{25}$$

where *χ* rad tot |*θ*=*θ*0 = 1 − *l* rad tot |*θ*=*θ*0 , *λ* = *χ* rad tot | 2 *θ*=*θ*0 *λ*cr, and

$$\lg(\mathcal{M}) \quad = \; \mathcal{M}^{2(1-\gamma)/(\gamma+1)} \Big( \frac{\mathcal{M}^2}{2} + \frac{1}{\gamma-1} \Big)\_{\cdot} \tag{26}$$

$$f(\mathbf{x}) = \mathbf{x}^{4(\gamma - 1)/(\gamma + 1)} \left( \frac{\chi\_{\text{tot}}^{\text{rad}}|\_{\theta = \theta\_0}}{\mathbf{x}} + \frac{1}{\gamma - 1} \right),\tag{27}$$

$$
\lambda\_{\rm cr} = \frac{1}{4} \left( \frac{2}{5 - 3\gamma} \right)^{(5 - 3\gamma)/[2(\gamma - 1)]} \quad . \tag{28}
$$

**Figure 1.** Model geometry. (**a**) Integration is performed along a radial streamline for a fixed meridional angle *θ*0. (**b**) The full angular dependence can be obtained by integrating along different radial streamlines.

For this case Equations (17)–(27) remain the same with the only change being the substitution of *x* by *x* − *xs* , where *xs* = *Rs*/*rB*. Here we take *Rs*/*rB* = 1.75 × 10−5 for the calculations of Figure 2 and 2.18 × 10−6 for the rest of the calculations. The radiative Bondi-like problem then reduces to solving Equations (25)–(28) for the Mach number, M, as a function of radius. On the other hand, it must be noticed that the PW potential, in addition to reproducing the innermost stable circular orbit and the marginally bound orbit, also accounts for a non-vanishing critical radius for *γ* = 5/3, in contrast to the purely Newtonian radial accretion, for which

$$\mathbf{x}\_{\text{crit}} = \mathbf{x}\_{\text{s}} + \sqrt{\frac{2\chi\_{\text{total}}^{\text{rad}}\mathbf{x}\_{\text{s}}}{\mathfrak{Z}}},\tag{29}$$

✡

where now the factor *χ* rad total (>0) appears, which is not present in [38,47].

**Figure 2.** (**Left**) Angular dependence of the UV emission for models *M*1,2,3 and *M*? 1,2,3. The horizontal black solid line at *χ* rad tot = 0 mark the transition from type 5 to types 1 and 2 solutions. (**Right**) Mach number as a function of radius for model *M*? 1 and *γ* = 1 (analytical solution; orange curves), *γ* = 1.1 (numerical solution from [38] with *ψ*rad = −*C*/*r*; gray curves, and this work with *C* = *C*(*r*); green curves), *γ* = 3/2 (blue curves), *γ* = 1.55 (black curves), *γ* = 1.66 (magenta curves) and *γ* = 5/3 (non-radiative (NR), with *xs* = 2.5 × 10−3 ; red curves). The thick solid curves represent supersonic inflow (*x* < *xs*) and outflow (*x* > *xs*) solutions, while the dashed lines represent inflow and ourflow subsonic solutions.

✑

The angular variation of the UV emission, *χ* rad tot , for models with and without the force multiplier is shown in the top panel of Figure 2. Model *M*1 has *f*? = 0.5 and *f*disc = 0.5, model *M*2 has *f*? = 0.2 and *f*disc = 0.8, while model *M*3 has *f*? = 0.05 and *f*disc = 0.95, with *M*(*τ*) = 0 so that both types of radiation can be evaluated in detail. Models *M*? 1 , *M*? 2 , and *M*? 3 differ from the preceding models in that *M*(*τ*) 6= 0. For different fractions of *f*?, *f*disc, and position *x*, we estimate the critical angles for which a transition from type 5 to type 1 and 2 solutions is allowed1 . This plot, which corresponds to *ρ*˜ = 0.001, *x* = 3, and *γ* = 1.55, is an example of the complexity of the system. In particular, *M*? 2 (green line with empty squares) and *M*? 3 (red line with full triangles) have physical transonic solutions for *θ* > *π*/8, while model *M*? 1 (blue line with empty circles) does not admit solutions for which *χ* rad tot < 0, and so this flow can be seen in any angular direction. The transition for *M*? 2 occurs at a lower angle than for *M*? 3 , meaning that we should find more collimated streamlines for the *M*? 2 flows than for the *M*? 3 flows. It is interesting to see that model *M*1, for which only electron scattering is present (*M*(*τ*) = 0) and *f*? = *f*disc = 0.5, behaves like model *M*? 1 with no transition to the unphysical region where *χ* rad tot < 0. As the fraction of disc photons increases (i.e., from *f*disc = 0.8 to 0.95), critical angles for which unphysical type 5 solutions occur (for *θ* phase *M*2,3 < *π*/8) are found even when only electron scattering is dominating.

In both sets of models *M*1,2,3 and *M*? 1,2,3, the radiation acceleration depends on the ionization parameter defined as

$$\mathfrak{F} = \frac{(f\_{\star}L)\exp\left(-\tau\_{\mathbf{X}}\right)}{n r^{2}},\tag{30}$$

where *L* = 7.45 × 1045 erg s−1 , which is appropriate for a 108*M* SMBH accreting at an efficiency of 8%. The gas density is *n*(*x*) = (*ρ*∞*ρ*˜(*x*))/(*µmp*) cm−3 , with *ρ*∞ = 10−20 gr cm−3 , *µ* = 0.7, and the optical depth *τX* is integrated using Equation (15) with *κX* = 0.4 g −1 cm2 . The physical parameters for all models computed in this study are listed in Table 1 (a full analysis of the angular dependence is beyond the scope of the present work, and will be presented in a forthcoming paper). In all cases we choose *θ*0 = *θ* phase acc . The fractions *f*? and *f*disc in Table 1 are exactly the same used by P07, which comply with observational results from [49,50].

**Table 1.** Parameters of the accretion models.

(a) *M*(*t*) 6= 0 (yes) and *M*(*t*) = 0 (no). (b) The values of *γ* = {1.1, 1.2, 1.3, 13/9, 3/2, 1.53, 1.55, 1.63, 1.65, 1.66}. (c) The angle *θ* phase acc is given in terms of *θ* phase, i.e., *θ* phase acc = 1.1*θ* phase. All models *M*? ∗ and *M*∗ shared the same angle.

The bottom panel of Figure 2 shows the dependence of the Mach number of the inflow and outflow transonic solutions with the normalized radial distance from the black hole for model *M*? 1 and various characteristic values of the adiabatic index, namely *γ* = 1 (full analytical solution; orange curves), *γ* = 1.1 (numerical solution from [38] with *ψ*rad = −*C*/*r*; gray curves, and this paper with *ψ*rad = −*C*(*r*)/*r*; green curves), *γ* = 3/2 (blue curves), *γ* = 1.55 (black curves), *γ* = 1.66 (magenta curves), and *γ* = 5/3 (red curves). Note that the latter curves correspond to the PW model with no radiation, i.e., with *f*? = *f*disc = 0. The intersection of curves of the same colour define the position of the sonic point (*xs*). The thick solid curves for *x* < *xs* represent supersonic inflow solutions, while those at larger radii (i.e., *x* > *xs*) represent supersonic outflow solutions. The dashed curves at radii *x* > *xs* are subsonic inflow solutions and those for *x* < *xs* are subsonic outflow solutions. The Mach number profiles are shown up to radial distances of ≈10−3.1*rB* from the centre. It is evident from this figure that isothermal (*γ* = 1) inflow occurs at slightly higher values of Macc(*r*) than its non-isothermal counterparts. As the adiabatic index is increased the accretion rate occurs at relatively smaller values of the Mach number everywhere and closer to the central source so that the region of supersonic accretion involves progressively smaller central volumes as *γ* increases. In contrast, the isothermal outflow takes place at much lower Mach numbers compared to the non-isothermal cases. These differences are clearly expected to influence the rates for non-isothermal accretion. As the adiabatic index is increased above *γ* = 3/2 the outflow close to *x* = 10−3 becomes supersonic only when *γ* = 1.66 and keeps subsonic for all values <1.66. In the interval 0.01 . *x* . 0.1, the outflow proceeds supersonically for all non-isothermal cases (*γ* > 1) with a small dependence of the Mach number on *γ*, while towards larger *x*(&0.2) the outflow becomes more supersonic as the adiabatic index is increased.

### **3. Results**

The Mach number as a function of distance from the source is shown in the top panel of Figure 3 for model *M*? 1 and varied *γ*. The solid curves depict the inflow solutions, while the dotted lines with arrows correspond to the outflow solutions. In all cases, the complete subsonic and supersonic solutions are shown. The position of the sonic point, *x*crit, is defined by the intersection between the inflow and the outflow solutions. When matching the subsonic portion of the solution at *x* ≤ *x*crit with the supersonic one at *x* > *x*crit, we then obtain the "Bondi-outflow" solution. The difference with the Parker solution is the location of the integration limits in the Bernoulli equation. In the Bondi problem they are located at infinity, while in the Parker problem they are located at the base of the wind. For all models, the Bondi outflows are always depicted by dotted lines with small arrows pointing out of the system. On the other hand, when matching the subsonic part of the solution at *x* > *x*crit with the supersonic one at *x* ≤ *x*crit, we then obtain the "Bondi-inflow" solution. Increased values of *γ* produce slower inflows and faster outflows everywhere, while the position of the sonic point is shifted towards smaller radii. The left panel of Figure 4 shows a close comparison of the Mach number profiles between the Newtonian *γ* = 1.66 case with spectral line forces (gray curves) and the *γ* = 5/3 case with *M*(*τ*) = 0 (orange curves). For *M*(*τ*) 6= 0 the inflow and outflow Mach number profiles essentially overlap those for the Newtonian *γ* = 1.66 case with *M*(*τ*) 6= 0 regardless of the radiation field. In spite of differences in the value of *M*(*τ*), the inflow and outflow profiles for both cases look almost identical.

Figure 3 (right panel) displays the Mach number as a function of radius for model *M*? 2 and varied *γ*. Compared to the left panel of Figure 3, we see that increasing the fraction of UV emission (from *f*? = 0.5 to 0.2) increases the inflow Mach number everywhere, while the position of the sonic point is shifted towards larger radii. Consequently, for given *γ*, when the ionizing flux decreases faster inflows are produced everywhere. In particular, at small radii *x* . 0.01 the inflow becomes more supersonic. A similar trend is observed for model *M*? 3 with the weakest X-ray heating (*f*? = 0.05) and the strongest UV emission (*f*disc = 0.95), where the sonic point is moved farther away compared to model *M*? 2 and therefore the supersonic inflow close to the centre occupies a larger volume. A comparison of the right panel of Figure 4 with that of left panel shows that as the UV heating dominates over the X-ray luminosity, the (*γ* = 5/3) inflow becomes faster than the *γ* = 1.66 − *M*? 1 inflow. It is interesting to note that the increase of the UV emission from *f*disc = 0.8 to *f*disc = 0.95 has only mild effects on the inflow and outflow Mach number profiles as well. The Mach number profiles for models *M*1 and *M*2 with no force multiplier (i.e., *M*(*τ*) = 0) and varied *γ* are now displayed in Figure 3 (right), respectively. The profiles for model *M*3 are very similar to those shown in Figure 3 for model *M*2. For model *M*1, where the central X-ray heating is the strongest and no line driving is present, the solution does

not admit type 5 flows regardless of the incident angle. As shown in Figure 3 (bottom) the critical points occur far from the central source between *x* ≈ 0.2 (for *γ* = 1.1) and *x* ≈ 0.007 (for *γ* = 1.66). These distances are relatively similar to those found for model *M*? 1 with force multiplier. For *θ*0 = 3 8 *π* and close to the SMBH (*x* ≈ 10−3 ), the Mach number ranges from ≈1 (for *γ* = 1.66) to ≈8 (for *γ* = 1.1). When the fraction of UV emission is increased to *f*disc = 0.8, the position of the sonic point is shifted outwards so that the flow becomes supersonic in a larger central volume as can be seen by comparing bottom with top. However, model *M*2 is characterized by smaller supersonic volumes than model *M*3. Similarly to the cases with force multiplier, as the UV emission is increased, the flow close to the central black hole becomes more supersonic. Moreover, for *γ* = 1.66 models *M*2 and *M*3 both exhibit outflows that are slower at large radii (*x* > 1) than their Bondi-inflow counterparts close to the central source, while the inverse is seen to occur for larger *γ*.

**Figure 3.** (**Top**) Radial Mach number profiles of the inflow (solid lines) and outflow (dotted lines with arrows) solutions for model (**a**) *M*? 1 with *M*(*τ*) 6= 0 and (**b**) *M*1 with *M*(*τ*) = 0, *f*? = 0.5, *f*disc = 0.5, and varied *γ* in the interval 1.1 ≤ *γ* ≤ 1.66. (**Bottom**) (**c**) *M*? 2 with *M*(*τ*) 6= 0 and (**d**) *M*2 with *M*(*τ*) = 0, *f*? = 0.2, *f*disc = 0.8.

**Figure 4.** Comparison of the Mach number profiles between the Newtonian (*γ* = 1.66) case with spectral line forces and the (*γ* = 5/3) case with *M*(*τ*) = 0.

For all models the density profiles are similar and almost independent of the adiabatic index. In all cases the density increases for *x* . 0.1 from ≈(1−7) to values as high as ≈103 when the UV emission rises from *f*disc = 0.5 to *f*disc = 0.95 for models with *M*(*τ*) 6= 0. In contrast, when *M*(*τ*) = 0 the highest central densities are seen to occur for model *M*1 with *ρ*˜ ≈ 5 × 102 at *x* = 10−3 , while models *M*2 and *M*3 with stronger UV emission and weaker central X-ray heating reach central densities of ≈1000 at the same radial distance from the SMBH, respectively. Additionally, in these cases the density profiles show a very little dependence on the adiabatic index. It is clear from the present results that the flow dynamics are sensitive to the radiative process that dominates the environment.

### **4. Catalogue of Pure Absorption Spectral Line Shapes for 1.1** ≤ *γ* ≤ **1.66**

As resolution improves with the emergence of more powerful telescopes, it becomes mandatory to predict the shape of absorption spectral lines as seen by a distant observer. This can be done by imagining a radially falling atom onto the SMBH according to the predictions of models *M*? 1,2,3, which would absorb photons emitted by the inner region at 10 in the rest-frame. Then we use the Doppler-shift formula [51]:

$$w = \frac{w\_0}{\gamma L (1 - \frac{\mathbf{v}(x)}{c} \cos \phi)},\tag{31}$$

where *w* is the angular frequency of an emitted photon with rest-frame frequency *w*0 as measured by the observer, *γL* is the Lorentz factor, v(*x*) is the velocity of the absorbing atom, and *φ* is the angle between the streamline and the line-of-sight towards the observer, which we set to *φ* = 0. For falling particles we set v(*x*) = −M(*x*)*c*s, where *c*s is given by Equation (7)2 , for values of ≈360, 420, 430, 440, and 420 km s−1 . A simple description of the absorption spectrum is given by

$$F\_{\lambda}(\mathbf{x}) = \exp\left\{-A(\mathbf{x})\exp\left[-\frac{\left(\lambda(\mathbf{x}) - \lambda\_0\right)^2}{\sigma^2}\right]\right\},\tag{32}$$

where *ν*0*λ*0 = *w*0/(2*π*)*λ*0 = *c*, *σ* is given by 0.1%*λ*(*x*), *c* is the speed of light, and *λ*0 is the photon wavelength in the rest frame. The intensity of the absorption is modelled by setting *A*(*x*) = *ρ*(*x*)/*ρ*max, where *ρ*(*x*) = *ρ*˜(*x*)*ρ*∞ and *ρ*max = *ρ*˜(*x* = 3 × 10−4 )*ρ*∞. As the particle is getting closer and closer to the black hole, the density increases and the shape of the spectral line tends to have a deeper deep shifted towards the red. Figure 5 depicts the predicted accretion absorption line shape for models *M*? 1 , *M*? 2 , and *M*? 3 with varied values of *γ*. From this figure, it is clear that the most asymmetrical lines are for model *M*? 3 corresponding to *f*? = 0.05, for which the high-energy radiation is weaker making the particle to attain higher velocities and therefore more asymmetry. Telescopes with higher resolutions will be needed to resolve the asymmetry of these lines. As the UV emission dominates over the X-ray heating the absorption lines become broader, while they become narrower with increasing *γ* regardless of the radiation field. For *γ* = 1.1 the symmetry is enough to distinguish which source of radiation is dominant with a resolution of ≈2500 km s−1 . However, as *γ* increases this property is lost as changes in the fraction of UV heating from *f*disc = 0.5 to 0.95 cannot be easily distinguished from the shape of the spectral absorption lines.

**Figure 5.** Absorption line shapes as predicted for models *M*? 1 , *M*? 2 , and *M*? 3 with varied adiabatic indices in the range 1.1 ≤ *γ* ≤ 1.66. The solid red lines depict the profiles of a typical absorption line for an atom falling off onto a black hole.

✠

✠

### **5. Estimated Bondi Radius and Mass Accretion with Respect to Conditions at Infinity**

We are now in the position to quantify the ratio between the estimated (*re*) and the true (*rB*) values of the Bondi radius and that between the estimated (*M*˙ *e*) and the true (*M*˙ rad) accretion rates when the boundary conditions are not at infinity. These ratios are given as functions of the radial distance from the SMBH by the following relations

$$\frac{r\_{\varepsilon;\{\text{acc,out}\}}(\mathbf{x})}{r\_B} = \left(\frac{\mathbf{x}^2 \mathcal{M}\_{\{\text{acc,out}\}}}{\lambda}\right)^{2(\gamma - 1)/(\gamma + 1)}\,,\tag{33}$$

✠

$$\frac{\dot{M}\_{\varepsilon;\{\text{acc,out}\}}(\mathbf{x})}{\dot{M}\_{\text{rad}}} = \frac{1}{\chi\_{\text{tot}}^{\text{rad}}|\_{\theta=\theta\_0}^2} \left[ \frac{r\_{\varepsilon;\{\text{acc,out}\}}(\mathbf{x})}{r\_B} \right]^{-(5-3\gamma)/[2(\gamma-1)]},\tag{34}$$

respectively, where *χ* rad tot (*θ*, *x*) is the angle- and space-dependent radiative factor. Figure 6 (left) shows the estimated Bondi radii for model *M*? 1 for adiabatic indices between *γ* = 1.1 and 1.66. At small radial distances (*x* = 8 × 10−3 ), model *M*? 1 with *γ* = 1.1 has an estimated Bondi radius of *re*/*rB* ≈ 0.50. As the adiabatic index is increased, the ratio *re*/*rB* decreases reaching values slightly less than 0.5 for *γ* = 1.1 and as low as ≈0.03 for *γ* = 1.66 at *x* = 10−2 . The same trends are also seen for models *M*? 2 (see Figure 6 (right)) and *M*? 3 as *γ* is increased, except that now the estimated Bondi radii achieve slightly lower values as the radiation field is dominated by the UV emission from the accretion disc. The profiles for model *M*? 3 are very similar to those displayed in Figure 6 with the ratios *re*/*rB* achieving slightly smaller values close to the centre. For all models regardless of the adiabatic index the ratio *re*/*rB* → 1 at large radii (*x* = 3).

Figure 7 shows the estimated accretion rates for models *M*? 1 , and *M*? 2 , for varied values of *γ*. When the UV emission is increased from *f*disc = 0.8 to 0.95 for model *M*? 3 the trends are almost the same as those in Figure 7 (right) for model *M*? 2 . In all cases, the radiative effects lead to an overestimation of the accretion rates for *x* . 0.1 (model *M*? 1 ) and *x* . 0.12 (models *M*? 2 and *M*? 3 ). The differences between *M*˙ *e*(*x*) and *M*˙ rad grow rapidly as we get closer to the central source. However, the rapid growth of the ratio *M*˙ *e*(*x*)/*M*˙ rad towards

the centre becomes progressively less steep as *γ* is increased, with *M*˙ *e*(*x*)/*M*˙ rad becoming almost flat for *γ* = 1.66. For instance, model *M*? 1 has an accretion rate at *x* = 10−2 for *γ* = 1.1 that is about 100 times larger than for *γ* = 1.66. A similar difference between both values of *γ* happens when comparing models *M*? 2 and *M*? 3 , where the UV emission from the disc dominates over the incident X-ray flux. The observed large overestimations in the estimated accretion rates make sense because as the closer the accreting matter is to the source of gravitation, the less accurate are the determination of *ρ* and *c*s from their values at infinity. However, the authors of [38] found that the level of overestimation of the accretion rates compared to the classical Bondi problem is slightly attenuated when the effects of line driving are ignored (model *M*1 with *γ* = 1.1). The same is true for higher values of *γ*.

**Figure 6.** (**Left**) Estimated Bondi radii for model *M*? 1 with *M*(*τ*) 6= 0. (**Right** *M*? 2 with *M*(*τ*) 6= 0), *f*? = 0.5(0.2) and *f*disc = 0.5(0.8) for varied values of the adiabatic index in the range 1.1 ≤ *γ* ≤ 1.66. The vertical solid lines mark the position of the sonic point and the symbols on them identify to which inflow solution they correspond to.

**Figure 7.** (**Left**) Estimated accretion rates for model *M*? 1 with *M*(*τ*) 6= 0 (**Right** *M*? 2 with *M*(*τ*) 6= 0), *f*? = 0.5 and *f*disc = 0.5 for varied values of the adiabatic index in the range 1.1 ≤ *γ* ≤ 1.66. The vertical solid lines mark the position of the sonic point and the symbols on them identify to which inflow solution they correspond to.

### **6. Conclusions**

In this paper, we have investigated the non-isothermal Bondi accretion onto a supermassive black hole (SMBH) when the adiabatic index is varied in the interval 1.1 ≤ *γ* ≤ 1.66 by including the effects of X-ray heating and UV emission. The radiation field is modelled using the prescriptions introduced by [13] (P07), where a central SMBH is considered as a source of X-ray heating and a surrounding optically thick, geometrically thin, standard accretion disc as a source of UV radiation. The radiative, non-isothermal radial accretion is solved semi-analytically using a method similar to that developed by the authors of [41]. The X-ray heating from the SMBH is assumed to be isotropic, while the UV emission from the accretion disc is assumed to have an angular dependence. The effects of both types of radiation are studied for different flux fractions (*f*? and *f*disc) and an incident angle

*θ*0 = *θ* phase acc in the presence and absence of the effects of spectral line driving. The *γ* = 5/3 regime is also explored using the Paczy´nski–Wiita (PW) potential for the case when the effects of line driving are ignored.

The relevant conclusions are summarized as follows:

**Author Contributions:** Conceptualization, J.M.R.-V. and L.D.G.S.; methodology, J.M.R.-V., L.D.G.S. and E.C.; software, J.M.R.-V.; formal analysis, J.M.R.-V., L.D.G.S. and J.K.; validation, J.M.R.-V., L.D.G.S., J.K., R.G. and E.C.; writing—review and editing, J.M.R.-V., L.D.G.S., J.K., R.G. and E.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been partially supported by Yachay Tech under the project for the use of the supercomputer QUINDE I and by UAM-A through internal funds. One of us (J. K.) acknowledges financial support from CONACYT under grant 283151.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The code to reproduce the average h*ξ*i calculations underlying this article are available in Zenodo at doi:10.5281/zenodo.4075323. For the new varying *ξ*, it will be available under request to the corresponding author.

**Acknowledgments:** IMPETUS is a collaboration project between the ABACUS-Centro de Matemática Aplicada y Cómputo de Alto Rendimiento of Cinvestav-IPN, the School of Physical Sciences and Nanotechnology of the Yachay Tech University, and the Área de Física de Procesos Irreversibles of the Departamento de Ciencias Básicas of the Universidad Autónoma Metropolitana–Azcapotzalco (UAM-A) aimed at the SPH modelling of astrophysical flows.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Notes**

### **References**

### *Article* **Geodesic Circular Orbits Sharing the Same Orbital Frequencies in the Black String Spacetime**

**Sanjar Shaymatov 1,2,3,4,5,\* and Farruh Atamurotov 3,6**

**Abstract:** We consider isofrequency pairing of geodesic orbits that share the same three orbital frequencies associated with Ω*r*ˆ , Ω*ϕ*ˆ , and Ω*ω*ˆ in a particular region of parameter space around black string spacetime geometry. We study the effect of a compact extra spatial dimension on the isofrequency pairing of geodesic orbits and show that such orbits would occur in the allowed region when particles move along the black string. We find that the presence of the compact extra dimension leads to an increase in the number of the isofrequency pairing of geodesic orbits. Interestingly we also find that isofrequency pairing of geodesic orbits in the region of parameter space cannot be realized beyond the critical value J*cr* ≈ 0.096 of the conserved quantity of the motion arising from the compact extra spatial dimension.

**Keywords:** isofrequency; geodesic orbits; black string

### **1. Introduction**

In general relativity, black holes are regarded as one of the most interesting and fascinating objects due to their geometric and remarkable gravitational properties. However, the point to be noted is that their existence had been predicted as a potential tool in explaining observational phenomena until the discovery of the gravitational waves detected by LIGO and Virgo scientific collaborations [1,2] and as well as the first direct shadow of supermassive black hole at the center of galaxy 87 by the Event Horizon Telescope (EHT) collaboration [3,4]. These observations opened a qualitatively new stage to reveal black hole's unknown properties and test remarkable nature of the background geometry around black hole's horizon irrespective of the fact that there are still fundamental problems that general relativity faces, i.e., the occurrence of singularity, spacetime quantization, etc. In this framework, the motion of test particles in the strong gravitational field regime has been a productive field of study for several years [5–19]. On the other hand, there is an extensive body of work devoted to understand the nature of radiative inspirals around black holes as a source of gravitational waves and binary systems [20–28].

It is worth noting that experiments and observations always allow to test the nature of the geometry through alternative theories of modified gravity in the strong field regime. The spactime geometry that can depart from the spherically symmetric Schwarzschild black hole was considered by Grunau and Khamesra [29]. This spacetime metric describes a five dimensional black string due to the compact extra dimension added to the well known Schwarzschild metric. Thus, it is one of reasons why this solution is interesting object to

**Citation:** Shaymatov, S.; Atamurotov, F. Geodesic Circular Orbits Sharing the Same Orbital Frequencies in the Black String Spacetime . *Galaxies* **2021**, *9*, 40. https://doi.org/10.3390/ galaxies9020040

Received: 21 April 2021 Accepted: 25 May 2021 Published: 28 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

investigate its properties. The motion of test particles in the different gravity string models has been investigated thoroughly by several authors [30–35]. Furthermore, the geodesic motion of test particles in the black string spacetime for both rotating and non-rotating cases has been studied in detail [29]. The geodesic motion of colliding particles was also studied in the vicinity of black string spacetime [36].

The theoretical investigation of isofrequency pairing of geodesic circular orbits has been considered in the vicinity of Schwarzschild and Kerr black holes, i.e., the region of the parameter space where such geodesic orbits occur has been discussed by providing an intuitive explanation of their occurrence. The above investigation was addressed by taking an external asymptotically uniform magnetic field in the Schwarzschild spacetime [25]. It turns out that the surface of allowed region where isofrequency pairing orbits occur decreases as a consequence of the presence of the external magnetic field. It was later extended to the spinning particles in the Schwarzschild-de Sitter spacetime [26]. In the present paper, we consider isofrequency pairing of geodesic orbits around black string spacetime geometry, as shown by the line element in [29]. We then study the effect of the compact extra spatial dimension *ω* on the isofrequency pairing of geodesic spiral orbits sharing the same three orbital frequencies.

In this paper, we consider a new frequency Ω*ω*ˆ arisen from the compact extra spatial dimension associated with the motion in the extra direction along the black string. Generally, the bound orbits are confined to the interior of a compact special spiral given by *hz*1(*e* = 0) ≤ *hz* ≤ *hz*2(*e* = 1); frequencies, in particular, Ω*r*ˆ and Ω ˆ*θ* are considered as "libration"-type frequencies associated with the radial and longitudinal periods, while Ω*ϕ*ˆ is a "rotation"-type frequency. In the case of black string, new frequency Ω*ω*ˆ can be defined by a "libration-rotation" type frequency which involves two kinds of type frequencies, as described above. For the sake of clarity, we will further focus on the motion of test particles at the equatorial plane (i.e., *θ* = *π*/2), and hence we omit Ω ˆ*θ* as it loses its meaning. Then bound geodesic orbits around black string spacetime are respectively characterized by three orbital frequencies, i.e., Ω*r*ˆ , Ω*ϕ*ˆ , and Ω*ω*ˆ associated with the periodic motions. Thus, these orbits are considered as triperiodic orbits around the black string spacetime. In doing so we predict that an infinite number of pairs of such orbits may occur in a region of parameter space around the black string spacetime where physically distinct orbits possessing E, L, and J share the same orbital frequencies associated with Ω*ϕ*ˆ , Ω*r*ˆ , and Ω*ω*ˆ , thereby referring to the "isofrequency" bound geodesic orbits synchronized in phase *ϕ* while passing periastra at the same time. Isofrequency pairing of geodesic orbits occuring in the vicinity of black string spacetime may in principle play an important role in an astrophysical scenario in understanding the nature of radiative inspirals around black holes.

The paper is organised as follows: In Section 2 we briefly describe the metric for black string and provide the geodesic equations of motion. We further study the effect of the compact extra dimension on the isofrequency pairing of geodesic spiral orbits in Section 3. We end up with conclusion in Section 4. Throughout the paper we use a system of geometric units in which *G* = *c* = 1.

### **2. The Metric and the Geodesic Equations of Motion**

The metric describing five dimensional black string spacetime with the compact extra spatial dimension *ω* added to the Schwarzschild metric in the Boyer-Lindquist coordinates (*t*,*r*, *θ*, *ϕ*, *ω*) is given by [29]

$$\begin{split} ds^2 &= -\left(1 - \frac{2M}{r}\right)dt^2 + \left(1 - \frac{2M}{r}\right)^{-1}dr^2 \\ &\quad + r^2(d\theta^2 + \sin^2 d\varphi^2) + d\omega^2 \end{split} \tag{1}$$

with *M* being the total mass of black string and *ω* corresponding to the compact extra spatial dimension.

The radial motion of geodesic test particles in the equatorial plane of black string spacetime satisfies the following equation

$$
\dot{\tau}^2 = \mathcal{E}^2 - V\_{eff}(r; \mathcal{L}, \mathcal{J}) \, , \tag{2}
$$

where

$$N\_{eff}(r; \mathcal{L}, \mathcal{I}) = \left(1 - \frac{2M}{r}\right) \left(1 + \mathcal{J}^2 + \frac{\mathcal{L}^2}{r^2}\right),\tag{3}$$

with E = *E*/*m*, L = *L*/*m*, and J = *J*/*m* being the conserved constants per unit mass and related to the particle's specific energy, angular momentum, and conserved quantity arising from *ω*, respectively.

It is well known that bound orbits exist in the case when L > 2 √ 3*M* with 2 √ 2/3 < E < 1 for the Schwarzschild black hole. Whereas, in the case of black string, bound orbits only exist when L > 2 √ 3(1 + J 2 ) 1/2*M* with 2 √ 2(1 + J 2 ) 1/2/3 < E < (1 + J 2 ) 1/2 is satisfied, and the motion of particles is then restricted by the turning points being labeled as the periastron *r*p and the apastron *r*a, respectively.

We now consider the innermost stable circular orbit (ISCO) for test particles orbiting around black string spacetime. To find radii of the ISCO one needs to solve the following equation

$$W\_{eff}^{\prime\prime}(r; \mathcal{L}, \mathcal{J}) = 0 \,. \tag{4}$$

In our case *rISCO* can be determined implicitly from the condition

6L 2 (*r* − 4*M*) − 4(1 + J 2 )*Mr*2 = 0 , (5)

with L obtained on solving *V* 0 *e f f*(*r*; L, J ) = 0. Thus, we have the ISCO radius *rISCO* = 6*M* for black string spacetime. It is worth noting that the ISCO is the same in contrast for the Schwarzschild black hole.

The standard condition *Ve f f*(*rp*,L, J ) = *Ve f f*(*ra*,L, J ) = E 2 leads to the specific energy E and angular momentum L of the particle, given in terms of *p* and *e*. Here, *rp* = *Mp*/(1 + *e*) and *ra* = *Mp*/(1 − *e*) respectively refer to the periastron and the apastron [25]. In between these two turning points lies bound orbits with corresponding specific energy and angular momentum. Note that parameters *p* and *e* respectively measure the size of the orbit and its degree of noncircularity [37]. Then we have

$$\mathcal{L}^2 = \frac{(p - 2 - 2e)(p - 2 + 2e)}{p(p - 3 - e^2)}(1 + \mathcal{J}^2),$$

$$\mathcal{L}^2 = \frac{p^2 M^2}{p - 3 - e^2}(1 + \mathcal{J}^2). \tag{6}$$

Following Cutler et al. [38], we turn to the integration of the geodesic equations in terms of *t*(*r*), *ϕ*(*r*), and *ω*(*r*) coordinates

$$t(r) = \mathcal{E} \int \frac{r dr}{(r - 2M) \left[\mathcal{E}^2 - V\_{eff}(r; \mathcal{L}, \mathcal{J})\right]^{1/2}} \tag{7}$$

$$\varphi(r) = \mathcal{L} \int \frac{dr}{r^2 \left[\mathcal{E}^2 - V\_{eff}(r; \mathcal{L}, \mathcal{I})\right]^{1/2}},\tag{8}$$

and

$$
\omega(r) = \mathcal{I} \int \frac{dr}{\left[\mathcal{E}^2 - V\_{eff}(r; \mathcal{L}, \mathcal{I})\right]^{1/2}}.\tag{9}
$$

Using Equations (6)–(9) and substituting *r*(*χ*) = *Mp*/(1 + *e* cos *χ*) into Equation (3), we then have the following expressions together with *χ* related to the coordinates *t*, *ϕ*, and *ω* as

$$\frac{dt}{d\chi} = \frac{Mp^2[(p-2)^2 - 4e^2]^{1/2}(p - 6 - 2e\cos\chi)^{-1/2}}{(p - 2 - 2e\cos\chi)(1 + e\cos\chi)^2},\tag{10}$$

$$\frac{d\varphi}{d\chi} = \frac{p^{1/2}}{(p - 6 - 2e \cos \chi)^{1/2}} \, \, \, \, \tag{11}$$

$$\frac{d\omega}{d\chi} = \frac{\mathcal{I}M p^{3/2} (p - 3 - \varepsilon^2)^{1/2}}{(1 + \mathcal{I}^2)^{1/2} (p - 6 - 2\varepsilon \cos \chi)^{1/2} (1 + e \cos \chi)^2} \,. \tag{12}$$

The radial period and frequency are defined by [24,25,38]

$$T^{\circ} = \int\_{0}^{2\pi} \frac{dt}{d\chi} \, d\chi \, \, \quad \Omega^{\circ} = \frac{2\pi}{T^{\circ}} \,. \tag{13}$$

The azimuthal frequency of the orbit will then have the form as

$$
\Omega^{\Phi} = \frac{1}{T^{\sharp}} \int\_{0}^{T^{\sharp}} \frac{d\varphi}{dt} dt = \frac{\Delta \varphi}{T^{\sharp}} \,, \tag{14}
$$

with azimuthal phase ∆*ϕ*, which is given by

$$\begin{split} \Delta \boldsymbol{\varrho} &= \int\_{0}^{2\pi} \frac{d\boldsymbol{\varrho}}{d\boldsymbol{\chi}} d\boldsymbol{\chi} = \\ &= \int\_{0}^{2\pi} \frac{\sqrt{p}}{\sqrt{p-6-2e\cos\chi}} d\chi = 4\sqrt{\frac{p}{\varepsilon}} K \left( -\frac{4e}{\varepsilon} \right), \end{split} \tag{15}$$

where *e* = *p* − *ps*(*e*) = *p* − 6 − 2*e* and *K*(*x*) = R *π*/2 0 *dθ*(1 − *x* sin2 *θ*) −1/2 is the complete elliptic integral of the first kind.

As a consequence of the presence of the compact extra dimension *ω*, the new frequency can be defined by

$$
\Omega^{\hat{\omega}} = \frac{1}{T^{\hat{\mathbb{P}}}} \int\_0^{T^{\hat{\mathbb{P}}}} \frac{d\omega}{dt} dt = \frac{\Delta \omega}{T^{\hat{\mathbb{P}}}} \, \, \, \, \tag{16}
$$

where ∆*ω* is the new extra phase accumulated over time interval *T r*ˆ .

### **3. Isofrequency Pairing Orbits in Black String Spectime Geometry**

Now we investigate the isofrequency and bound orbits of the particles in the vicinity of the black string spacetime. We consider orbits lying extremely close to the separatrix *e* in order to estimate the divergent quantities of the azimuthal and the new extra phases, the radial period, and their ratio in Equations (14) and (16). Considering the near-separatrix analytic expansions we obtain the corresponding expression for ∆*ω* as

$$
\Delta\omega = \frac{\mathcal{I}\mathcal{M}(3+2\varepsilon-\varepsilon^2)^{1/2}}{(1+\mathcal{I}^2)^{1/2}(1+2\varepsilon)}\sqrt{\frac{(6+2\varepsilon)^3}{e}}
$$

$$
\times \left[1+\mathcal{O}\left(\frac{\varepsilon}{4e}\right)\right]\ln\left(\frac{64e}{\varepsilon}\right). \tag{17}
$$

The previous investigations suggest that the occurrence of isofrequency pairing of geodesic orbits is strongly related to the presence of boundary regions referred as *separatrix* and *circular-orbit duals* (COD), and each and every circular orbit in the open range between the separatrix and singular curve and between the singular curve and COD as well has a dual isofrequency pairing. Note that the existence of the *separatrix* plays an important role for occurrence of these dual orbits in this isofrequency scenario. Considering the singular curve and the circular orbit duals in the particular region of the gravitational object is particularly important, and the COD can also play an important role due to the fact that as a boundary region allowing to keep pairs of isofrequency orbits that exist in this region [24–26].

Let us then study the effect of the compact extra dimension on the occurrence of the isofrequency pairing of geodesic orbits in black string spacetime. Based on the above discussions and Equations (13)–(15) we keep the result for which the range of frequencies of any pairs of orbits on Ω*ϕ*ˆ = (*M*/*r* 3 *b* ) 1/2 in the black string is given by

$$
\boldsymbol{\Omega}\_{\varepsilon=0}^{\phi} = 0.062 < \boldsymbol{\Omega}^{\phi} < \boldsymbol{\Omega}\_{\varepsilon=1}^{\phi} = 0.125 \,\text{.}\tag{18}
$$

where we have defined Ω˜ *ω*ˆ = *M*Ω*ϕ*ˆ as a dimensionless parameter.

Further, we determine the values of radial and azimuthal frequencies numerically and tabulate their values in Table 1.

**Table 1.** Values of Ω*r*ˆ and Ω*ϕ*ˆ are tabulated for different values of eccentricity *e* for isofrequency pairing of geodesic spiral orbits around the black string.

In Table 2, we show numerical values of the frequency Ω˜ *ω*ˆ of the bound orbits around the black string. As can be seen from Table 2 the value of the frequency Ω˜ *ω*ˆ arisen from the compact extra spatial dimension is increasing with an increase in the value of the forth conserved quantity J , but the height of the horizontal direction along black string is decreasing for given values of eccentricity *e*. However, numerical analysis leads to the fact that when the value of the fourth conserved quantity attains J*cr* ≈ 0.096 geodesic spiral orbits are not allowed all through to occur around the black string. Thus, the pairs of geodesic circular spiral orbits of particles with the same radial Ω˜ *r*ˆ , azimuthal Ω˜ *ϕ*ˆ , and new Ω˜ *ω*ˆ frequencies occur in the particular region of the gravitational object under this critical value irrespective of the fact that these orbits are physically distinct with different values of conserved quantities (i.e., E, L, and J ). Then the pairs of these spiral orbits of particles with three frequencies associated with (Ω˜ *r*ˆ , Ω˜ *ϕ*ˆ , Ω˜ *ω*ˆ ) would move and oscillates in a cylindrical form along the black string.

It turns out that as a consequence of the decrease in the value of J the above mentioned cylindrical region becomes larger around the black string. It could then play an important role in allowing the particles to move along the black string, thus increasing the amount of isofrequency pairs of geodesic spiral orbits in the region between the *separatrix*, *circular-orbit duals* and *horizontal length hz*, i.e., in the particular region of the black string. Hence, one can keep in mind that there would exist an infinite number of pairs of circular spiral orbits having the same three orbital frequencies around the black string regardless of physically distinct obits with three conserved quantities. As a consequence of the presence of compact extra spatial dimension, the occurrence of isofrequency pairing of spiral orbits would be of particularly importance in understanding and explaining the behavior of black string spacetime and the nature of the radiative inspirals as a source of gravitational waves or binary systems.

**Table 2.** The values of Ω*ω*ˆ are tabulated for different values of the fourth conserved parameter J for isofrequency pairing of geodesic spiral orbits around the black string. The point to be noted here is that *hz* corresponds to the horizontal length and will be valid only for the non vanishing compact extra dimension, i.e., *ω* 6= 0.

### **4. Conclusions**

We considered isofrequency pairing of geodesic spiral orbits in the black string spacetime, which recovers the Schwarzschild one in the case of vanishing *ω*. Thus, in an astrophysical scenario, it is worth investigating these geodesic orbits in the vicinity of black string as it allows to test the effects arising from the five dimensions. Our analysis suggests that the compact extra dimension *ω* leads to an increase in the number of pairs of circular spiral orbits. These pairs of orbits having the same three orbital frequencies (i.e., Ω*r*ˆ , Ω*ϕ*ˆ , and Ω*ω*ˆ ) occur in the particular region and oscillate in the cylindrical form around the black string in spite of the fact that they are physically distinct orbits possessing different conserved quantities (i.e., E, L, and J ). Furthermore, we showed that pairs of geodesic spiral orbits are not allowed all through to occur in the black string vicinity in the case when the fourth conserved quantity attains J ≈ 0.096.

The obtained results suggest that an infinite number of pairs of such orbits may occur in the particular region around the black string and it would be particularly important in explaining not only the possibility of occurrence of pairs of geodesic spiral orbits around the black holes, but also the nature of the black string spacetime. These theoretical results and discussions would be useful to the possible interpretation of astrophysical observations, and they could help to provide details of the validity of alternative models to black holes and also explain the nature of the radiative inspirals since isofrequency paring of geodesic spiral orbits can give the same information as such inspiral objects.

In a future work, we would be considering the possible extensions of recent results to the case of rotating black string and studying the effects of the rotation and a compact extra dimension on the specific isofrequency pairing of geodesic orbits. This would be of primary astrophysical importance.

**Author Contributions:** Conceptualization, S.S. and F.A.; investigation, S.S. and F.A.; methodology, S.S. and F.A.; validation, S.S. and F.A.; formal analysis, S.S. Both authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing not applicable—no new data generated.

**Acknowledgments:** S.S. acknowledges the support of the Uzbekistan Ministry for Innovative Development. F.A. acknowledges the support of INHA University in Tashkent.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

### *Article* **Study of Accretion Flow Dynamics of V404 Cygni during Its 2015 Outburst**

**Arghajit Jana 1,2,\* , Jie-Rou Shang 3 , Dipak Debnath 2 , Sandip K. Chakrabarti 2 , Debjit Chatterjee 2,4 and Hsiang-Kuang Chang 5**

	- Indian Institute of Astrophysics, Koramangala, Bangalore 560034, India

**Abstract:** The 2015 Outburst of V404 Cygni is an unusual one with several X-ray and radio flares and rapid variation in the spectral and timing properties. The outburst occurred after 26 years of inactivity of the black hole. We study the accretion flow properties of the source during its initial phase of the outburst using *Swift*/XRT and *Swift*/BAT data in the energy range of 0.5–150 keV. We have done spectral analysis with the two component advective flow (TCAF) model fits file. Several flow parameters such as two types of accretion rates (Keplerian disk and sub-Keplerian halo), shock parameters (location and compression ratio) are extracted to understand the accretion flow dynamics. We calculated equipartition magnetic field *Beq* for the outburst and found that the highest *Beq* ∼ 900 Gauss. Power density spectra (PDS) showed no break, which indicates no or very less contribution of the Keplerian disk component, which is also seen from the result of the spectral analysis. No signature of prominent quasi-periodic oscillations (QPOs) is observed in the PDS. This is due to the non-satisfaction of the condition for the resonance shock oscillation as we observed mismatch between the cooling timescale and infall timescale of the post-shock matter.

**Keywords:** X-rays: binaries—stars individual: (V404 Cygni)—stars:black holes—accretion; accretion disks—shock waves—radiation:dynamics

### **1. Introduction**

Transient black hole candidates (BHCs) have two phases in their lives: the quiescence phase and the outbursting phase. They spend most of their lifetimes in the quiescence phase. A sudden rise in viscosity leads to an outburst when the X-ray intensity rises by a factor of thousands or more that of the quiescence phase. Matter from the companion star accretes to the central black hole, and in this process, gravitational potential energy is converted to heat and radiation. A black hole (BH) X-ray spectrum generally consists of two components: a multi-color blackbody and a hard power-law tail. The multi-color blackbody part is believed to originate from a Shakura-Sunyaev type standard thin disk [1]. The power-law tail is believed to originate from a Compton corona [2,3]. In the two-component advective flow (TCAF) solution, the CENBOL or CENtrifugal pressure supported BOundary Layer [4,5] replaces the Compton corona used in other models, such as disk-corona model [6,7] or evaporated disk in ADAF [8]. In this paper, we used the TCAF solution to study the accretion flow dynamics of V404 Cygni during its first outburst in 2015 after a long quiescent of ∼26 years.

V404 Cygni is one of the most studied black hole X-ray binary systems. It is also known as GS 2023+338. It was first identified as an optical nova in 1938. In 1956, another nova outburst was reported in this system [9]. In 1989, V404 Cygni went through another outburst. The 1989 outburst was discovered by the all sky monitor onboard *Ginga*. It is

**Citation:** Jana, A.; Shang, J.-R.; Debnath, D.; Chakrabarti, S.K.; Chatterjee, D.; Chang, H.-K. Study of Accretion Flow Dynamics of V404 Cygni during Its 2015 Outburst. *Galaxies* **2021**, *9*, 39. https://doi.org/ 10.3390/galaxies9020039

Academic Editors: Yosuke Mizuno and Emilio Elizalde

Received: 10 April 2021 Accepted: 21 May 2021 Published: 25 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

located at RA = 306◦ .01 and Dec = 33◦ .86. After spending ∼26 years in the quiescent state, V404 Cygni went through a short but violent outburst on 15 June 2015. In December 2015, another short activity was observed [10]. The binary system V404 Cygni harbors a black hole of mass 9–12 *M* at the center with a K-III type companion of mass ∼1 *M* [11–13]. The inclination angle of the binary system is ∼67◦ [12,13]. The orbital period of the system is 6.5 days [11]. V404 Cygni is located at a distance of 2.39 kpc, measured by parallax method [14]. V404 Cygni has a high spinning black hole with spin parameters *a* ∗ > 0.92 [15].

The 2015 outburst of V404 Cygni was discovered on 15 June simultaneously by *Swift*/BAT [10] and *MAXI*/GSC [16]. This outburst of the source was extensively observed in multi-wavelength bands, such as in radio [17], optical [18] and X-ray [19,20]. *INTE-GRAL* observation reported multiple X-ray flares during the outburst [19]. Several radio flares were also observed [17]. Extreme variable and rotating jet is observed during the outburst [21–23]. The source showed rapid changes in the spectral properties in a very short time [24]. *INTEGRAL* observation detected *e* − − *e* + pair annihilation on 20 June 2015 [25]. *FERMI*/LAT detected high energy *γ*-ray jet in the source on 26 June 2015 [26]. King et al. [27] reported detection of emission lines with Chandra-HETG, indicating strong disc wind emission. Motta et al. [24] argued that a slim disk is formed at the inner accretion flow due to the super-Eddington accretion in the source.

In this *paper*, we study the timing and the spectral properties of V404 Cygni with combined *Swift*/XRT and *Swift*/BAT data in the broad energy range of 0.5–150 keV during the initial phase of the 2015 outburst. We have done spectral analysis with the TCAF model-based *fits* file to extract physical flow parameters. The nature of these model-fitted accretion flow parameters allowed us to investigate the physical reasons behind the origin of the several flares and their variability and turbulent features. We have also calculated the equipartition magnetic field and cooling time during the outburst of the source.

The paper is organized in the following way. In Section 2, we briefly discuss the disk structure prescribed by TCAF and the way flow parameters decide on the spectral shape. In Section 3, we discuss the observations and the data analysis procedure. In Section 4, we present the results of our analysis. In Section 5, we make a discussion based on our result, and finally, in Section 6, we summarize our findings.

### **2. TCAF Solution**

TCAF configuration is based on the solution of a set of equations which govern viscous, transonic flows around a black hole [28]. In the TCAF solution, accreting flow consists of two components: (*i*) high viscous, high angular momentum, optically thick and geometrically thin Keplerian disk flow ( . *m*d), (*ii*) low viscous, low angular momentum, optically thin sub-Keplerian halo flow ( . *m*h). The Keplerian disk accretes on the equatorial plane and is immersed within the sub-Keplerian flow. Centrifugal force rises as flow moves close to the black hole and becomes comparable with the gravitational force. Low viscous halo matter is fast-moving and forms an axisymmetric shock at the centrifugal barrier [29]. The post-shock region i.e., CENBOL 'puffed-up', and becomes 'hot'. The Keplerian disk truncates at the shock location. The multi-color black body part of the observed spectrum that is formed by the soft photons is emitted from the Keplerian disk. A fraction of these soft photons are intercepted by the CENBOL. Depending on the temperature and size of the CENBOL, soft photons become hard photons through repeated inverse-Comptonization scattering at the CENBOL. Conversely, a fraction of Comptonized photons is reprocessed at the Keplerian disk and produces a reflection hump. Thus in the TCAF solution, the reflection component is self-consistently incorporated. However, a Gaussian line may be required to add if an iron emission line is present. CENBOL is also considered to be the base of the jets or outflows [30]. Toroidal magnetic flux tubes are responsible for the collimation of jet [31,32]. Oscillation of CENBOL can be triggered when the cooling and heating times inside CENBOL are similar, and the emerging photons produce the quasi-periodic oscillations (QPOs) [33–35] hereafter C15.

Transient BHCs generally show different spectral states during their outbursts. In TCAF, these observed spectral states are controlled by the flow parameters (two types of accretion rates and two shock parameters). A typical outbursting BHC generally goes through spectral state transitions to form a hysteresis loop as follows: hard state (HS) → hardintermediate state (HIMS) → soft-intermediate state (SIMS) → soft state (SS) → softintermediate state (SIMS) → hard-intermediate state (HIMS) → hard state (HS) [36,37]. In the upper panel of Figure 1 (adopted from [38]), we show a cartoon diagram of the above four spectral states under the TCAF solution. In the lower panel, typical spectra of each spectral state correspond to the diagrams are shown. In the cartoon diagrams, brown, light green, dark green, and grey regions represent Keplerian disk, sub-Keplerian halo, CENBOL and jet, respectively. It is to be noted that the cartoon diagram and corresponding spectra are only for illustrative purposes. TCAF model fits file assumes the shock location as the truncation radius of the Keplerian disk ; however, it may not occur all the time. For example, in the HS, a disk may not be formed at all, or in the SS, the shock may not be formed with the negligible sub-Keplerian flow. However, these scenarios are the example of extreme cases. In general, both Keplerian disk and sub-Keplerian halo matters are present in the accretion flow, and the TCAF model fits file covers the most probable parameter space of the system. quasi-periodic oscillations (QPOs) [33–35] hereafter C15. Transient BHCs generally show different spectral states during their outbursts. In TCAF, these observed spectral states are controlled by the flow parameters (two types of accretion rates and two shock parameters). A typical outbursting BHC generally goes through spectral state transitions to form a hysteresis loop as follows: hard state (HS) → hardintermediate state (HIMS) → soft-intermediate state (SIMS) → soft state (SS) → softintermediate state (SIMS) → hard-intermediate state (HIMS) → hard state (HS) [36,37] . In the upper panel of Figure 1 (adopted from [38]), we show a cartoon diagram of the above four spectral states under the TCAF solution. In the lower panel, typical spectra of each spectral state correspond to the diagrams are shown. In the cartoon diagrams, brown, light green, dark green , and grey regions represent Keplerian disk, sub-Keplerian halo, CENBOL and jet, respectively. It is to be noted that the cartoon diagram and corresponding spectra are only for illustrative purposes . TCAF model fits file assumes the shock location as the truncation radius of the Keplerian disk ; however, it may not occur all the time. For example, in the HS, a disk may not be formed at all, or in the SS, the shock may not be formed with the negligible sub-Keplerian flow. However, these scenarios are the example of extreme cases. In general, both Keplerian disk and sub-Keplerian halo matters are present in the accretion flow, and the TCAF model fits file covers the most probable parameter space of the system.

and heating times inside CENBOL are similar, and the emerging photons produce the

*Galaxies* **2021**, *1*, 0 3 of 21

**Figure 1.** In upper panel, cartoon diagrams of four commonly observed spectral states under TCAF solution are shown (adopted from [38]). Brown, light green, dark green and grey region represent Keplerian disk, sub-Keplerian halo, CENBOL and jet, respectively. An outbursting black hole evolves as HS → HIMS → SIMS → SS → SIMS → HIMS → HS. Generally both rising and declining HS and HIMS show monotonically evolving type-C (low frequency) QPOs, and SIMS shows sporadic type-B or type-A QPO. In soft sate these QPOs are absent. Here, EQPO means evolving QPO; SQPO means sporadic QPO; NQPO means no QPO. In the bottom panel, theoretical spectra corresponding to the top paneled spectral states are shown. These spectra are generated using five input parameters (*M*BH in *M*; *m*˙*d* in *M*˙ Edd; *m*˙*h* in *M*˙ Edd ; *X*s in *rs*, *R*), whose values are marked inset. The cartoon diagram and corresponding spectra are only for illustrative purpose. **Figure 1.** In upper panel, cartoon diagrams of four commonly observed spectral states under TCAF solution are shown (adopted from [38]). Brown, light green, dark green and grey region represent Keplerian disk, sub-Keplerian halo, CENBOL and jet, respectively. An outbursting black hole evolves as HS → HIMS → SIMS → SS → SIMS → HIMS → HS. Generally both rising and declining HS and HIMS show monotonically evolving type-C (low frequency) QPOs, and SIMS shows sporadic type-B or type-A QPO. In soft sate these QPOs are absent. Here, EQPO means evolving QPO; SQPO means sporadic QPO; NQPO means no QPO. In the bottom panel, theoretical spectra corresponding to the top paneled spectral states are shown. These spectra are generated using five input parameters (*M*BH in *M*; . *md* in . *M*Edd; . *mh* in . *M*Edd ; *X*s in *rs*, *R*), whose values are marked inset. The cartoon diagram and corresponding spectra are only for illustrative purpose.

When an outburst is triggered due to lower viscosity and angular momentum, the sub-Keplerian matter moves with free-fall velocity, whereas the Keplerian flow moves in viscous time. Thus, the sub-Keplerian flow dominates the accretion process since it moves faster than the flow in the Keplerian disk. The Keplerian disk is truncated very far away by a large CENBOL. A strong shock (higher *R*) is formed at hundreds of Schwarzschild radius (*r*s) away from the BH. Thus, it is difficult to cool the CENBOL by the Keplerian component. When an outburst is triggered due to lower viscosity and angular momentum, the sub-Keplerian matter moves with free-fall velocity, whereas the Keplerian flow moves in viscous time. Thus, the sub-Keplerian flow dominates the accretion process since it moves faster than the flow in the Keplerian disk. The Keplerian disk is truncated very far away by a large CENBOL. A strong shock (higher *R*) is formed at hundreds of Schwarzschild radius (*r*s) away from the BH. Thus, it is difficult to cool the CENBOL by the Keplerian component. Hard X-ray flux dominates, and a hard state is observed. Compact jet is launched in this

state from the CENBOL [30]. Evolving type-C QPO is produced in this state due to the resonance oscillation of the shock [33].

The source enters the HIMS after HS. The Keplerian disk accretion rate continues to rise and becomes comparable with the sub-Keplerian halo accretion rate. As a result, accretion rate ratio (ARR = . *m*h/ . *m*d) decreases. Due to the rise of the Keplerian disk accretion rate, CENBOL shrinks and becomes cool as the shock moves farther inward. Shock strength decreases as the Compton cooling reduces the post-shock thermal pressure. Further, in this intermediate state, in the presence of jets, mass outflow rate to inflow rate ratio becomes maximum. In this state, evolving type-C QPOs may also be observed.

In the SIMS, the Keplerian rate keeps on increasing, although the sub-Keplerian flow rate started to decrease. This is because more and more sub-Keplerian flow becomes Keplerian by viscous transport. The shock becomes weak in this state. The shock further moves in, and the CENBOL becomes small. The soft X-ray flux increases, and the hard X-ray flux decreases in this state due to rapid rise in . *m*d and slow decrease in . *m*h. Generally, type-A or type-B QPOs are observed sporadically in this state due to weak oscillation of the CENBOL (type-B) or due to oscillation of shock-less centrifugal barrier (type-A) [34,35]. The time (day) difference between . *m*h and . *m*d peaks gives us rough estimation of the viscous time scale of the source [39,40]. In the SS, the Keplerian disk dominates and completely cools down the CENBOL. The soft X-ray flux dominates over the hard X-ray flux. No shock is formed. As a result, the jet is completely quenched in this state [30,41]. No QPO is produced in the soft state.

The spectral states evolve in the reverse order during the declining phase of the outburst as flow parameters show opposite nature. Starting from the SS to the SIMS transition day, both the Keplerian disk accretion rate and the sub-Keplerian halo accretion rate decreases, although the Keplerian disk rate decreases faster. As a result, ARR increases. As in the SIMS of the rising phase, one may see sporadic type-B or A QPOs in the declining SIMS. In the declining HIMS and HS, evolving type-C QPOs could be seen. Similar to the rising phase, one could observe a compact jet in the HIMS and HS in the declining phase.

TCAF solution is implemented in XSPEC [42] as an additive table model to analyze the spectral properties around the black holes [43,44]. TCAF model has four input flow parameters ( . *m*d, . *m*h, *X*s, *R*) other than mass of the black hole (*M*BH) and normalization (*N*). In the recent years, accretion flow dynamics around several black holes are studied quite successfully using TCAF model [45–48]. Frequencies of the dominating QPOs are predicted from the TCAF model fitted shock parameters [49]. Masses of the black holes are estimated quite successfully from spectral analysis with the TCAF model [49–51]. Jet contribution in the X-rays are also calculated using TCAF solution [52–54].

While fitting with the TCAF model, one must take care of degeneracy. Like any other spectral model, TCAF model is degenerate between different parameters. Thus, one must be very careful while fitting with TCAF model. One must always remember the physical condition of the system. For example, a single spectrum can be fitted with two different set of parameters: ( . *m*d = 0.1 . *M*Edd, . *m*h = 0.5 . *M*Edd, *M*BH = 10 *M*, *X*s = 100 *r*s, *R* = 3.5) and ( . *m*d = 0.2 . *M*Edd, . *m*h = 0.8 . *M*Edd, *M*BH = 10 *M*, *X*s = 10 *r*s, *R* = 4). Clearly, the second set of parameters is unphysical (the Xs can not be low while the the shock is strong and mh is higher than md). Also, one needs to find the parameters with the steppar command to avoid the degeneracy.

### **3. Observation and Data Analysis**

We analyzed *Swift* data for 19 observations between 2015 June 15 (MJD 57,188.77) and 26 June (MJD 57,199.52). We studied the source in 0.5–150 keV energy band with combined XRT and BAT data for five observations (MJD 57,191.01, 57,194.54, 57,197.21, 57,197.33, and 57,198.02). The 15–150 keV BAT data was used for four observations (MJD 57,188.77, 57,193.56, 57,198.15 and 57,199.52) when only BAT observations were available. For the rest of the ten observations, we studied 0.5–10 keV using XRT data.