A relativistic orbit model for temporal properties of AGN

We present a unified model for X-ray quasi-periodic oscillations (QPOs) seen in Narrow-line Seyfert 1 (NLSy1) galaxies, $\gamma$-ray and optical band QPOs that are seen in Blazars. The origin of these QPOs is attributed to the plasma motion in corona or jets of these AGN. In the case of X-ray QPOs, we applied the general relativistic precession model for the two simultaneous QPOs seen in NLSy1 1H 0707-945 and deduce orbital parameters, such the radius of the emission region, and spin parameter $a$ for a circular orbit; we obtained the Carter's constant $Q$, $a$, and the radius in the case of a spherical orbit solution. In other cases where only one X-ray QPO is seen, we localized the orbital parameters for NLSy1 galaxies REJ 1034+396, 2XMM J123103.2+110648, MS 2254.9-3712, Mrk 766, and MCG-06-30-15. By applying the lighthouse model, we found that a kinematic origin of the jet based $\gamma$-ray and optical QPOs, in a relativistic MHD framework, is possible. Based on the inbuilt Hamiltonian formulation with a power-law distribution in the orbital energy of the plasma consisting of only circular or spherical trajectories, we show that the resulting Fourier power spectral density (PSD) has a break corresponding to the energy at ISCO. Further, we derive connection formulae between the slopes in the PSD and that of the energy distribution. Overall, given the preliminary but promising results of these relativistic orbit models to match the observed QPO frequencies and PSD at diverse scales in the inner corona and the jet, it motivates us to build detailed models, including a transfer function for the energy spectrum in the corona and relativistic MHD jet models for plasma flow and its polarization properties.


Introduction
Active galactic nuclei (AGN), at the center of most galaxies, are known to be powered by black holes of masses M • = 10 5 − 10 9 M [1 -3]. These systems are believed to be the scaled-up version of black hole X-ray binaries (BHXRB), possessing the same physical process of accretion [4]. The riveting evidence of this conjecture is the similarity between the X-ray variability in AGN and BHXRB [4,5]. However, a complete understanding of the physical processes of accretion and the jet emission in AGN requires the variability analysis in various wavelength bands, from optical to γ ray.
X-ray Power spectral density (PSD) shape: The X-ray variability is a key diagnostic for understanding the physical processes in the innermost regions of the accretion flow. The similarity in the behavior of X-ray variability in AGN and BHXRB is an important aspect of the AGN-BHXRB connection. The PSD of both BHXRB and AGN are known to show red noise, which decreases steeply at high frequencies (small timescales) following a power law, P(ν) ∝ ν −α , where α ∼ 2 typically [28][29][30]. At lower frequencies, below a characteristic frequency, called the break frequency (ν b ), the PSD flattens (α ∼ 1) [28][29][30]. Such a characteristic PSD shape is well defined by a bending power-law model [28] and found in various types of AGN with ν b ranging from ∼10 −6 -10 −4 Hz [29,31]. This break frequency is expected to approximately scale as an inverse of the black hole mass. However, the bending power-law shape of the PSD shape in BHXRB is known to be associated only with the soft spectral state [32]. Hence, the understanding of such a characteristic shape of the PSD is fundamental for probing the inner regions close to the black hole.
X-ray QPOs: The QPOs detected so far in the X-ray light curves of AGN are seen to be mostly associated with the Narrow-Line Seyfert 1 (NLSy1) galaxies, which are identified by the narrow width of their broad Hβ emission line with FWHM< 2000 kms −1 , strong FeII lines, and weak forbidden lines [33,34]. NLSy1 galaxies are also known to show rapid X-ray variability and near Eddington accretion rates [35]. The first detection of a significant QPO in an X-ray light curve was reported in RE J1034+396 with timescale ∼ 3730 s using the XMM-Newton data [8]. Another significant QPO at ∼ 3.8 hour timescale was reported in an Ultrasoft Active Galactic Nucleus Candidate 2XMM J123103.2+110648 [9]. A QPO with ∼ 2 h timescale was detected in MS 2254.9-3712 [10]. Later, 1H 0707-495 also showed the detection of a significant QPO at ∼ 3800s and another at timescale ∼ 8265 s with relatively low significance in the X-ray light curve [36,37]. A highly significant QPO of timescale ∼ 6450 s was reported in NLSy1 Mrk 766 [38], while another QPO (but not simultaneous) with a period of ∼ 4200 s was also reported [39], making these two signals be in a ∼ 3:2 resonance. Another significant X-ray QPO was reported in NLSy1 MCG-06-30-15 of timescale ∼ 3600 s [40]. Very recently, the detection of two QPOs was reported at timescales ∼ 8064.5 s and ∼ 14706 s in ESO 113-G010 [41]. All these X-ray QPOs, discussed above, were found in XMM-Newton data (0.3-10 keV). The connection between QPOs and 3:2 twin peaks in BHXRB, ultra-luminous X-ray sources (ULXs), and AGN was shown in [42] as a universal scaling of these frequencies with the black mass and spin.
Optical and γ ray QPOs: The optical and γ ray QPOs are also known to be discovered in a few BL Lacertae objects, also known as BL Lac. These objects are a class of AGN characterized by their large polarization, high variability, and weak emission lines [43,44]. These objects are interpreted as systems with a relativistic jet pointing directly towards the line of sight of the observer; hence, the jet emission dominates in these systems, and the discovered QPOs are thought to be associated with the jets. The γ ray QPOs are majorly discovered using the FERMI-LAT observations (100MeV-300GeV). A γ ray QPO of timescale T ∼ 630 days was reported in PKS2155-304 [11], where this timescale was found to be twice the optical period originally proposed by [45]. Later, the presence of both these QPO timescales was confirmed [12]. A QPO of timescale 2.18 ± 0.08 years was discovered in γ ray light curve of PG 1553+113 [13], where correlated oscillations were found in the radio and optical fluxes. Later, an optical QPO of similar timescale, ∼810 days, was confirmed in PG 1553+113 [14]. Another γ ray QPO with a timescale of a few months, T∼280 days, was reported in PKS 0537-441, where an optical QPO of timescale ∼T/2 was also discovered [15]. A pair of optical and γ ray QPO was also reported in BL Lac, having similar timescales of ∼ 680 days [14,16]. Another QPO of period 34.5 days is observed in the γ-ray light curve of blazar PKS 2247-131 [46]. Recently, an optical QPO, having a temporal period of 44 days, was detected in the Kepler light curve of an NLSy1 galaxy KIC 9650712, which may or may not be a jet-based QPO [21]. Another interesting case is OJ 287, which is a quasar with a quasi-periodic optical outburst emission cycle of 12 years. This prominent outburst is explained by a black hole binary model, where a secondary black hole interacts with the accretion disk of a much more massive primary black hole [22][23][24][25][26][27]. A comprehensive analysis of PSDs of 11 blazars was carried out recently, using the Fermi-LAT gamma-ray 10-years-long light curves, where a QPO in PKS 2155-304 was confirmed [47].
In this paper, we present a model that unifies the origin of multiwavelength QPOs originating in the disk and the jet and also probes the genesis of the PSD shape of the X-ray light curve due to a corona. We study the association of X-ray QPOs discovered in NLSy1 type AGN with the relativistic circular and spherical orbits around a Kerr black hole using the generalized relativistic precession model (GRPM) [48][49][50]. We also motivate that the non-equatorial trajectories (for example, spherical orbits), which are the consequence of axisymmetry of the Kerr spacetime, are also the viable solutions to the QPO frequencies using the GRPM [50]. We also apply a relativistic jet model [51,52] to study the optical and γ ray QPO timescales in BL Lacertae objects. This model describes the simultaneous QPOs with 1:2 or 3:2 frequency ratio as the harmonics obtained in the Fourier series of the Doppler factor of the pulse profile from a blob rotating along with the jet. The Doppler factor includes the relativistic effects, such as Doppler boost, relativistic aberration, gravitational redshift, and light bending. We also present a model to describe the typical bending power-law profile of the PSD observed in AGN. Assuming the bending power-law profile of the PSD shape, we find the intrinsic profile of the energy distribution of the particles orbiting in circular and spherical trajectories in the corona around a Kerr black hole, which results in a distribution in the fundamental frequency space. The X-ray flux gets modulated at this fundamental frequency, which is a function of distance from the black hole, and consequently also a function of E. A unified picture of these models of multiwavelength QPOs and PSD shape is shown in Figure 1, where r M is the marginally bound spherical orbit (MBSO) radius, r I is the innermost stable circular orbit (ISCO) radius, and r X is the outer disk radius. Figure 1. The figure shows a unified picture of the models for X-ray, optical, and γ ray QPOs and the origin of X-ray power spectral density (PSD) shape in AGN. The X-ray QPOs observed in NLSy1 galaxies are associated with the fundamental frequencies of the equatorial orbits in the accretion disk sandwiched by a corona region, which we call the outer corona (OC) region, r I < r < r X ; the inner corona (IC) region, r M < r < r I , is associated with the fundamental frequencies of the spherical orbits around a Kerr black hole. The optical and γ ray QPOs in Blazars are shown as the harmonics of the timescale of a blob of matter moving along the jet. The shape of the PSD is studied using the fundamental frequency of matter which is governed by the radial effective potential, V e f f (E, L, a, Q), providing the gravitational background responsible for the geodesic motion, in IC and OC regions to derive the energy distribution of the orbiting matter, N(E), which is directly related to the observed intensity, I(ν), where ν is the temporal frequency.
The structure of this paper is as follows: In Section 2, we discuss the generalized relativistic precession model (GRPM) for the X-ray QPOs [48][49][50]. In Section 2.1, we present the method for the error estimation in the parameters calculated for the case of AGN with two simultaneous X-ray QPOs, 1H 0707-495. In Section 2.2, we discuss the association of X-ray QPO frequencies with the equatorial circular orbits using the GRPM, whereas we study the spherical orbits as the origin of X-ray QPOs using the GRPM in Section 2.3. In Section 3, we apply a basic jet model [51,52] to study the timescales of optical and γ ray QPOs in Blazars. We then study the genesis of the bending power-law shape of the PSD in AGN and derive the intrinsic energy distribution of the orbiting particles in Section 4. We summarize our results in Section 5 and draw conclusions in Section 6. A glossary of symbols is provided in Table 1.

Relativistic Circular and Spherical Orbits as Solutions to X-Ray QPOs
The X-ray emission from NLSy1 galaxies is believed to have originated from the inner region of the accretion disk in the context of the unification model of AGN [3]. We apply the (G)RPM (RPM: [48,49]; GRPM: [50]) to study the X-ray QPOs discovered in a few cases of NLSy1 type of AGN and one Type-2 AGN candidate; see Table 2. The GRPM associates fundamental frequencies of the relativistic particle orbits in the accretion disk, close to a rotating black hole, with the QPO frequencies. Using this model, we estimate the parameters: the spin of the black hole, a, and radius of an equatorial circular orbit, r, in Kerr spacetime, where QPOs originate. We also implement the GRPM [50] to associate the frequencies of relativistic spherical orbits (non-equatorial) with the QPO frequencies to calculate the corresponding parameters (r s , a, Q), where r s is the radius of a spherical orbit and Q is the Carter's constant [53], which is the fourth integral of motion in the Kerr geometry, and defined as where p θ is the conjugate momentum in θ coordinate, L z is the z-component of particle's angular momentum, and E is its energy per unit rest mass. For the astrophysically relevant bound orbits, Q > 0 is a valid condition for which θ obeys 0 < θ − < θ < θ + < π, where θ − + θ + = π/2, so that the orbit is symmetric with respect to the equatorial plane [53,54]. In the equatorial plane, when θ = π/2, p θ vanishes and cos θ = 0; hence from Equation (1) we obtain Q = 0 for the equatorial orbits. The GRPM associates two simultaneous high-frequency QPOs (HFQPOs) observed in BHXRB [6,7,[48][49][50] with the fundamental frequencies: the azimuthal frequency ν φ and the periastron precession frequency 1 , ν pp = ν φ − ν r , where ν r is the frequency of radial oscillation. There are a few cases of BHXRB where a third low-frequency QPO (LFQPO) is also detected simultaneously to HFQPOs [60,61]; in such cases, the LFQPO is associated with the nodal precession frequency 2 , ν np = ν φ − ν θ , where ν θ is the frequency of vertical oscillation. See Figure 2 for an illustration of the precession frequencies.
In the GRPM, these frequencies are associated with the non-equatorial bound orbits (Q = 0), whereas only equatorial orbits (Q = 0) were studied in the RPM. There is no such known case in AGN where three QPOs are detected simultaneously; however, the X-ray QPO detected in Type 2 AGN 2XMM J123103.2+110648 (see Table 2) was suggested as a LFQPO because of its large rms value (25-50%) [9], which is the typical characterstic of LFQPOs observed in BHXRB [6].
Therefore, for AGN having a single QPO detection, we associate the QPO frequency with ν φ , except for 2XMM J123103.2+110648 where we also analyze the ν np frequency. For the cases of AGN with two simultaneous QPO detections, we use ν φ and ν pp frequencies in the GRPM. In Table 2, we have summarized the cases of AGN with either one or two simultaneous QPO detections in X-rays. 1 pp stands for the periastron precession. . Ω pp represents the periastron precession and Ω np represents the nodal precession frequency. The particle starts from the initial point A, and follows an eccentric and non-equatorial trajectory before completing one (a) radial, or (b) vertical oscillation to reach point B, where it sweeps an extra ∆φ azimuthal angle during one (a) radial, or (b) vertical oscillation because the azimuthal motion is faster than the radial or vertical motion. Image courtesy: [50].

Method for the Error Estimation
Here, we describe a generic procedure [50] which we have used to estimate errors in the orbital parameters for NLSy1 AGN with two simultaneous X-ray QPOs, 1H 0707-495; see Table 2.

1.
We assume that the frequencies, ν 1 and ν 2 , of QPOs are Gaussian distributed with their mean values at the centroid of observed QPO frequencies, ν 10 and ν 20 (with ν 10 > ν 20 ). The joint probability density distribution of these frequencies is given by where P i (ν i ) represents the Gaussian distribution of ith QPO frequency, given by where σ i is the observed standard dispersion (error) of the ith QPO.

2.
We find the Jacobian, J , of the transformation from frequency to orbital parameter space using the formulae of fundamental frequencies, which is given by ( where x 1 and x 2 represent the orbital parameters. For the equatorial circular trajectories (Q = 0), we have {x 1 , x 2 }={r, a}; whereas for the spherical trajectories (Q = 0), we have {x 1 , x 2 }={r s , a}.
The Jacobian is completely expressible in an analytic form and can be easily evaluated from Equation (3), and using the frequency formulae. We utilize Equation (7) for circular orbits in Section 2.2, and Equation (8) for spherical orbits in Section 2.3, to evaluate J (Equation (3)), where ν 1 = ν φ and ν 2 = ν pp according to the RPM and GRPM. 3.
Next, we write the probability density distribution in the parameter space given by where [x] represents the set of parameters {x 1 , x 2 } and J is given by Equation (3); and {ν 1 , ν 2 } are substituted in terms of parameters using the analytic formulae, Equation (7) for the circular orbits and Equation (8) for the spherical orbits.

4.
We calculate the exact solutions for parameters by solving ν φ = ν 10 and ν pp = ν 20 using Equation (7) for circular trajectories {r 0 , a 0 }, and Equation (8) for spherical trajectories {r s0 , a 0 } for fixed Q. We fix M • to the previously known values. We find 1σ errors in the parameters by taking an appropriate parameter volume around the exact solution, and generate sets of parameter combinations with resolution ∆x j in this volume. The chosen parameter range, exact solutions, and corresponding resolutions are summarized in Tables 3 and 4. We then calculate the probability density using Equation (4), for all the generated parameter combinations and normalize the probability density by the normalization factor where k varies from 1 to the number of total parameter combinations taken in the parameter volume; [x] k is the kth combination of the parameters in the parameter volume. Hence, the normalized probability density is given by The normalization of the probability density in the parameter space, discussed above, is done because only a sub-volume in the parameter space is astrophysically allowed for bound orbits, which is discussed below. 5.
The allowed parameter combinations for the bound orbits is governed by the condition given by [54] where e is the eccentricity and µ is the inverse latus-rectum of the general non-equatorial trajectory.
We have e = 0 for spherical orbits; hence, we ensure that the parameters (r s = 1/µ, a, Q) for spherical trajectories follow the above bound orbit condition. If any parameter combination does not obey the bound orbit condition, then P ([x]) is taken to be zero at that point in the parameter volume.

6.
For the circular orbit case, there are two parameters to estimate {r, a} using two QPO frequencies.
For the case of spherical orbits, there are three unknown parameters {r s , a, Q}; hence, we first take Q ={1, 4, 8, 12} for the spherical trajectory solutions, where the extrema of θ coordinate deviates away from the equatorial plane with an increase in Q. For each fixed value of Q, we find the normalized probability density distribution in the parameter space {x 1 , x 2 } = {r s , a} using Equation (5b). Later, using the calculated spin values and their errors for each fixed Q, we estimate the distribution of spin and the most probable spin. Using this distribution and the most probable value of the spin, we then determine the probability distribution in the {x 1 , Next, we integrate the normalized probability density, P ([x]), Equation (5b), in one dimension to obtain the profile in the other dimension. Thus, we finally obtain the one dimensional distributions {P 1 (r), P 1 (a)} for circular orbits, and {P 1 (r s ), P 1 (a)} for spherical orbits. 8.
Finally, we fit the normalized probability density profiles in each of the parameter dimensions to find the corresponding mean values and quoted errors are obtained such that it contains a probability of 68.2% about the peak value of the probability density. The results of these fit are given in Tables 3 and 4.

Circular Orbits
In this section, we use the GRPM (Q = 0) for QPOs to estimate the (r, a) parameters of the circular orbits using their fundamental frequencies, which are given by [61][62][63] where {ν φ ,ν r ,ν θ } are the dimensionless frequencies and M • is mass of the black hole. We use the dimensionless parameters: r is scaled by R g = GM • /c 2 and a ≡ J/ GM 2 • /c , where J is the angular momentum of the black hole. We use the convention a > 0 for the prograde and a < 0 for the retrograde orbits in this article.
We discuss our results below: 1.
We have computed the contours of ν φ (r, a), using Equation (7a), for the QPO frequencies (given in Table 2) of RE J1034+396 (blue), MS 2254.9-3712 (red), and MCG-06-30-15 (magenta), shown in the (r, a) plane in Figure 3a. The masses of these black holes were assumed from the previous estimations (see Table 2). We see that the QPO emission originates from a very narrow region of the accretion disk, where r ∼ (9.4 − 9.9) for RE J1034+396, r ∼ (10.4 − 11.4) for MS 2254.9-3712, and r ∼ 14.2 for MCG-06-30-15 even though a ranges from 0 to 1. This implies that the QPO emission region is very close to the black hole, and this emission region remains very narrow and nearly independent of the spin of the black hole.

2.
For the case of Mrk 766, two QPO frequencies were detected (see Table 2), but at different epochs. We have shown ν φ (r, a) contours for both these frequencies in Figure 3a, where ν 1 = 2.38 × 10 −4 Hz (orange) and ν 2 = 1.55 × 10 −4 Hz (green). The mass of the black hole was fixed to M • = 4.3 × 10 6 M [58]. The QPO origin range is r ∼ 10 for ν 1 and r ∼ (12.6−14) for ν 2 , which is again found to be in a narrow range and very close to the black hole. Although these QPOs were not detected simultaneously, we tried to estimate a simultaneous solution for (r, a) by equating ν φ = ν 1 and ν pp = ν 2 as per GRPM. We show them as curves in the (r, a) plane in Figure 3b, and we see that these contours do not cross each other, implying that there is no simultaneous solution for (r, a).

3.
For the Type-2 AGN 2XMM J123103.2+110648, the detected QPO (see Table 2) was suggested as an LFQPO type because of its large rms value [9]. If this QPO frequency is equated to the high-frequency component, ν φ (r, a), of the GRPM, we found that r ∼ 200, which is far from the black hole to emit X-rays. Hence, the GRPM predicts that this should be an LFQPO. We show the contours of the LFQPO component of the GRPM, ν np (r, a), in the (r, a) plane for the QPO frequency of 2XMM J123103.2+110648 in Figure 3c, where we fixed M • = 10 5 M [56]. We see that the emission region for this LFQPO is r ∼ (6 − 20), for the whole range of a. Hence, the detected QPO of 2XMM J123103.2+110648 is an LFQPO that originated very close to the black hole. 4.
For the case having two simultaneous X-ray QPOs, 1H 0707-495 (see Table 2), we first solve the equations {ν φ (r, a) = ν 10 , ν pp (r, a) = ν 20 } (using Equations (7a) and (7b)), assuming M • = 5.2 × 10 6 M [36], as per GRPM to estimate the exact solution for (r, a), which is found to be (r 0 = 8.214, a 0 = 0.0662). We then apply the method, described in Section 2.1, to estimate the errors in the parameters (r, a) implied due to the errors of the QPO frequencies. The range of (r, a) and corresponding resolutions used for our simulations are summarized in Table 3. Finally, we generate the probability density profiles in each parameter dimension {P 1 (r), P 1 (a)}, shown in Figure 4, where we have also shown the probability contours in the parameter space. The results of the model fits to the probability density profiles are summarized in Table 3. The errors in the parameters are quoted with respect to the exact solution (r 0 , a 0 ), whereas the simulated {P 1 (r), P 1 (a)} profiles peak at (r = 8.092, a = 0.038), which slightly differs from the exact solution. Hence, our analysis assuming the circular orbit frequencies as the origin of QPOs, using the GRPM, in NLSy1 1H 0707-495, suggests that it harbors a slowly rotating black hole (a ∼ 0.0662) at the center, and that the X-ray QPOs originate in the inner region of the accretion disk and very close to the black hole (r ∼ 8.214).

Spherical Orbits
In this section, we apply the GRPM for simultaneous QPOs of 1H 0707-495 to estimate the (r s , a, Q) parameters of the spherical orbits using their fundamental frequencies, which are given by [50,63] where ∆ = r 2 s + a 2 − 2r s and z ± are given by where L z is the z-component of particle's angular momentum and E is its energy per unit rest mass, which can be explicitly expressed as the functions of {r s , a, Q} parameters (see Equation (16) in [54]). The definitions of the Elliptic integrals are [64] F ϕ, We discuss our results below:

1.
We explore the parameter space (r s , a, Q) for the spherical orbits. Since there are two input QPO frequencies, we first vary the Q value to find various solutions of {r s , a} by solving equations {ν φ = ν 1 , ν pp = ν 2 } as per GRPM. Q = 13 is at the limit of astrophysically allowed bound orbits, Equation (6); Q < 13 in the case of 1H 0707-495. The Q = 13 orbit is an unstable orbit very close to the separation of bound and unbound (called a separatrix orbit), and such an unstable orbit is not relevant to our study; hence, we fix our parameter exploration between Q= 1 and 12. In Figure 5, we have shown these solutions in the (Q, a) and (Q, r s ) planes.

2.
Next, we fix Q = {1, 4, 8, 12} and find the errors in the {r s , a} parameters using the method described in Section 2.1. The range of {r s , a}, resolution taken in the simulations, along with the exact solutions and their errors obtained by fitting P 1 (r s ) and P 1 (a) are summarized in Table 4.

3.
The ranges of {a, r s , Q}, shown in Table 4 and Figure 5, span the complete parameter volume for QPO frequencies of 1H 0707-495. As the spin of the black hole does not change in the timescale of a few months or years, we need to find the most probable value of spin. We first find the variance of P 1 (a) with respect to the exact solution of a for each Q, given in Table 4, which is given by where P 1i (a) is the probability density ditribution in a parameter space for each value of Q. We have summarized the values of σ a for each Q in Table 4. We then minimize the likelihood function to obtain the most probable value of the spin given by We find the peak value to be a p = 0.139 for 1H 0707-495, and corresponding solution of {r s , Q} for the QPO frequencies is {r sp = 8.246, Q p = 9.814}.

4.
Next, we obtain the χ 2 a distribution function of a given by A plot of χ 2 a /χ 2 p is shown in Figure 6a, where χ 2 p = χ 2 a a p . We obtain the 2σ errors with respect to a p by normalizing the χ 2 a (a) function and obtain 0.139 0.183 −0.139 , where the region of 95% probability is indicated by the vertical dashed line in Figure 6a. We also show the range of r s and Q in Figure 6b,c within the 2σ region of a, as seen in Figure 6a, where the parameter ranges are r s = (8.214 − 8.323) and Q = (0.0001 − 12.264).

5.
Hence, we conclude that the spherical orbits, close to the black hole in the region, r s = (8.214 − 8.323) with Q values between (0.0001 − 12.264), are possible sources of the QPO frequencies observed in 1H 0707-495, while the most probable spin value to be a p = 0.139 0.183 −0.139 with 2σ confidence.

Relativistic Jet Model for the Optical and γ Ray QPOs
In a simple kinematic approach inspired by the lighthouse model [51,52,65], the basic periodicity is set by where m 6 = M • / 10 6 M and m 8 = M • / 10 8 M and r F is the radius of the footpoint of the magnetic field anchored in the equatorial plane. An important radius is the light cylinder radius, which given in geometrical units is r L = r 3/2 F + a. The plasma is expected to relativistically follow the field lines upto the light cylinder rigidly beyond which the field lines would be bend. A reasonable estimate of the cylindrical radius of the plasma motion is expected to be typically r 0 ≡ χr L where χ = (0.1 − 10). Taking an angular momentum conservation beyond the Alfven radius, r A = x A r L , where x A < 1, will lead to r 2 0 Ω 2 0 = x 2 A r 2 L Ω 2 F , setting an observed periodicity of T 0 = (χ/x A ) T F . The value of (χ/x A ) depends on details of the relativistic MHD models and x A is determined by the relativistic Bernoulli equation, but a range of (χ/x A ) = 1 − 20 is reasonable [52]. This is illustrated by estimating T F (r F ) for the range of r F = (30 − 80) (see Table 5); we see that the observed is in the range of (1 − 20).
This agreement motivates the study of the plasma motion in the background of relativistic MHD models, and its comparison with fits to the light curves in the future. Another clue of the jet physics will come from polarization models, as evidenced by the promising but simplistic cylindrical relativistic polarization signatures of the EVPA, DOP, and Doppler boosted flux profiles, as predicted by [51]; this will be an additional and useful tool to extract jet properties by doing detailed fits to polarization observations. There is an oscillatory behavior seen in both γ-ray and optical light curves [11][12][13][14][15][16]66,67] that supports the above trend. There is also evidence of the radio structure that is supported by the basic model of [52] as observed by [68,69]. Table 5. A list of statistically significant QPOs detected in the γ ray and optical energy bands in BL Lacertae type of AGN, along with their redshifts and black hole masses. The theoretical timescales are calculated, using Equation (12), such that the lower and upper limit correspond to r F = 30 and r F = 80 respectively. The lower-case letters (a to m) are links to references given in the last column.

Relativistic Orbit Model (ROM) and PSD Shape
The X-ray timing analysis of NLSy1 galaxies has been proven to be an essential tool for probing the emission region and the underlying mechanism of the variability process of the X-ray flux in these sources. The shape of the power spectral density is found to have a shape which is well fit by a bending power-law model given by [28] where P 0 is the normalization constant, and α l , α h are the PSD slopes below and above the break frequency, ν b . The power density spectrum shows that the low-frequency power spectrum is significantly flatter (α l ∼ 1) than the high-frequency power spectrum (α h > 2). The break frequencies were found to be near ν b ∼ 6.7 × 10 −6 Hz for PKS 0558-504 [29] and ν b ∼ 8 × 10 −4 Hz for NGC 4051 [28].
Here, we present a plausible relativistic orbit model to generate such a power density spectrum. As argued before, the non-equatorial orbits, such as spherical orbits, are the natural consequence of the axisymmetry of the Kerr space-time [53,54]. We assume that inside a spherical corona region of relativistic electrons (the inner corona, IC, r M < r < r I ), existing inside the radius of innermost stable circular orbit (ISCO) (see Figure 1), the particles are in non-equatorial orbits. The thin accretion disk spans the region outside ISCO, where the fluid motion is confined to the equatorial plane. We also assume that an outer corona region (OC, r I < r < r X ) of relativistic particles envelopes this accretion disk, lying almost in the equatorial plane (see Figure 1). The energy per unit rest mass of these relativistic particles, E, orbiting in the equatorial circular trajectories, is given by [62] E (r, a) = r 2 − 2r + a √ r We see that E increases with r outside ISCO, and it decreases with r inside ISCO, where it has minima at the ISCO radius; see Figure 7a. The stable circular orbits exist outside the ISCO radius, whereas the unstable circular orbits are found inside the ISCO radius. The mechanical energy per unit rest mass of the relativistic plasma, E, orbiting in the spherical trajectories, is given by [54] where E increases with r s outside the innermost stable spherical orbit (ISSO), and it decreases with r s inside ISSO, where it has minima at the ISSO radius; see Figure 7b. The stable spherical orbits, for a fixed Q, exist outside the ISSO radius, whereas the unstable spherical orbits are found inside the ISSO radius. A comparison of the ISCO and ISSO radii is shown in the (r, a) plane in Figure 8, where we see that these radii converge to r = 6 for the Schwarzschild black hole (a = 0), as the spherical orbits or ISSOs are possible only outside a Kerr black hole (a = 0) because of the axisymmetry of the Kerr space-time. Additionally, as the value of Q increases, the ISSO radii move further away from the black hole as compared to the ISCO radius. This implies that one will always find the unstable circular and the unstable spherical orbits inside the ISCO radius.

The ROM
The underlying assumptions of our model are:

1.
We associate the temporal frequency, ν, in the observed power spectral density with the fundamental azimuthal frequency of the particles orbiting in the circular orbits in the accretion disk outside ISCO, r I , and both circular and spherical trajectories between r I and marginally bound spherical orbit (MBSO) radius, r M . These frequencies are functions of the orbital radius, r or r s , (Equations (7a) and (8a)), and hence they are also fundamentally related to the mechanical energy of the orbit through Equations (14) and (15).

2.
We assume a prior distribution of the energy of particles (or electrons) given by a power-law , IC : radial range r M < r < r I , (16a) , OC : radial range r I < r < r X .
where N (E) represents the number of particles having energy E, α 1 and α 2 are the power-law indices inside and outside r I respectively, E I is the particle energy at r I , and A is the normalization constant. The energy distribution, N (E) (Equation (16)), is constructed so that it is continous at r I . Assuming that the total number of particles are N 0 (however, the PSD solution is independent of this), we have the normalization condition given by where the first and the second terms contribute for the regions inside and outside r I respectively, and E X corresponds to the energy of the particles at the outer radius of the equatorial circular accretion disk, r X . Subsequently, we obtain We redefine A (a) such that A (a, α 1 , where 3.
We assume that the break frequency of the PSD corresponds to the temporal frequency at the ISCO radius.

4.
We also assume that the particle distrbution in the temporal frequency space, F (ν), directly translates to the observed intensity for a given temporal frequency, so that the power density is given by Next, we derive the distribution of the temporal frequency, F (ν), as follows: We obtain (dN (E) /dE) from Equation (16), and numerically obtain E (ν) to derive (dE (ν) /dν), using Equations (7a) and (14) for circular and Equations (8a) and (15) for spherical orbits. In Figure 9, we have shown E(ν φ ), where it is clear that the behaviour of E(ν φ ) changes inside and outside r I . Outside r I , Figure 9a, the radius of the circular orbits goes from r I to r X ≈ 50 (1 keV≈ GM • m e r X k B ), whereas inside r I , Figure 9b, the radius of the circular orbits varies between r I and marginally bound circular orbit (MBCO), r M , which corresponds to E = 1. Next, we obtain the temporal frequency distribution given by whereν I (a) is frequency at r I andν X (a) is frequency at r X (where the energy is E X defined in Equation (17c)). N I (α 1 , α 2 , a, Q) and N X (α 1 , α 2 , a, Q) correspond to the number of particles having frequency at r I and r X respectively, where we have scaled the functions F 1 , F 2 , N I , and N X by N 0 . The expressions for the functions Φ 1 α 1 ,ν , a, Q and Φ 2 α 2 ,ν , a are given by We scale the distribution functions, Equation (19), by N I for simplicity, which yields where and We employ the condition that f 1 (α 1 , α 2 ,ν, a, Q) = f 2 (α 1 , α 2 ,ν, a, Q) at the frequency corresponding to {r I ,ν I (a)}, which gives or N I (α 1 , α 2 , a, Q) = −W (a, α 1 , α 2 ) C 2XI (α 2 , a) + N X (α 1 , α 2 , a, Q) .
Now, we describe the procedure to obtain the model parameters α 1 and α 2 using observations: 1.
If β 1 is the average slope of the observed PSD after the break frequency,ν >ν b , given by where ∆ represents the difference of values at the end points defined by MBSO and ISCO in our model: the end point of the PSD forν >ν b (whereν b =ν I (a)) is at the MBSO radius (ν M (a)), so that where where N I (α 1 , α 2 , a, Q) can be substituted using Equation (24), which yields Hence, for a given combination of {a, Q}, we obtain a relation, given by Equation (25f), where {α 1 , α 2 } are unknowns.

3.
We compute the slopes {α 1 , α 2 } by the above mentioned criteria for different combinations of (a, Q), which are shown in Table 6. We find that α 1 ranges between ∼[2.3 − 4] and α 2 is in the range ∼[3.7 − 8.9], indicating that a power-law model for the intrinsic mechanical energy of the orbiting matter describes the shape of the observed PSD reasonably well. Additionally, if we reverse the analysis to estimate {β 1 , β 2 } by fixing {α 1 = 2.5, α 2 = 3.5} for (a = 0.5, Q = 2), we find {β 1 = −1.97, β 2 = −0.77} which are in good agreement with observations. We also show contours of α 1 and α 2 in the (Q, a) plane in Figure 10, where the values of α 1 and α 2 increase with a. We also see that contours are independent of Q for small a, which is expected because the non-equatorial orbits do not exist in Schwarzschild spacetime.

4.
The examples of PSD profile obtained in the scaled frequency space,ν, are shown in Figure 11. We see that the PSD profiles for given parameter combinations in Table 6 show good fits to the expected bending power-law model, Equation (13). The PSD represents a general power spectrum obtained independent of the mass of the black hole; hence, it applies to the stellar-mass black holes also. This validates the ROM as a plausible model for PSD observed in black holes.   Table 6. The red curve shows the bending power-law model fit, given by Equation (13), where the fitting parameters are shown in Table 6. The vertical black dashed curve corresponds to the ISCO (break) frequency.

Summary
The results are summarized below: • In Section 1, we summarized the observations for X-ray QPOs, which are traditionally associated with the accretion disk and corona, γ-ray QPOs normally attributed to a jet, and the X-ray PSD usually connected with the inner and outer corona (see Figure 1).

•
In Section 2, we motivated the creation of (G)RPM models for X-ray QPOs and extracted the spins and radii for the sources, listed in Table 2, based on the model given in [48][49][50]. The GRPM model confirms that the detected QPO in Type-2 AGN 2XMM J123103.2+110648 is an LFQPO, as it was also suggested by [9]. In a statistical analysis, we were able to determine these parameters and their errors for 1H 0707-495, the case of two simultaneous QPOs, based on the observed QPO frequencies and their errors. The results are presented in Table 3 for circular orbits and in Table 4 for spherical orbits. We found non-planar orbits, with Q ∼ (1 − 12), which are very close to a Kerr black hole, that (r s ∼ (8.2 − 8.3); a ∼ 0.14) are the possible solutions for QPO frequencies of 1H 0707-495. • Next, in Section 3, we applied the relativistic kinematic jet model to check its validity by comparing the basic frequency with the observed QPO periods in BL Lac objects, given in Table 5. The ratio T 0 /T F is typically in the range 1 − 20, which is reasonable, given the range of footpoint radii of the field lines and typical location of the Alfvén point up to which the field line is rigid [52]. It motivates detailed relativistic MHD models along with polarization profile predictions (as given in [51]) to compare with observations. • In Section 4, we built a relativistic orbit model consisting of circular and spherical orbits that have a power-law distribution, and its mechanical energy is split into two parts (above and below the energy at ISCO). This formulation leads to unique results relating to the PSD slopes (before (ν <ν b ) and after (ν >ν b ) the break) with those of the energy spectrum for the given spin and mass of the black hole ( Figures 10 and 11). We plan to test this model against observations to extract {a, M • }.

Discussion and Conclusions
We add the following points of discussion of our results and conclusion: 1.
The periastron and nodal precession of the particle orbits is an intrinsic phenomenon in Kerr geometry, which is a consequence of strong gravity and axisymmetry of the spacetime. We propose in the GRPM [50] that the precession frequencies of matter blobs orbiting in these trajectories, very close to the Kerr black hole, modulate the X-ray flux, from the thin accretion disk where the flow is hot. The origin of these non-equatorial orbits of blobs in a slim torus region having a single radius is motivated in [50], where a model of fluid flow in the general relativistic thin accretion disk [73] is studied. In this study, we suggest that the edge region, defined in [73], is a launchpad for plasma instabilities, where blobs orbit with fundamental frequencies of the geodesics near the edge and in the geodesic region (defined in [73]), in which Hamiltonian dynamics is applicable. We also show in the GRPM that these geodesics span a torus region, which overlaps with the edge and geodesic region of [73]. 2.
The QPOs in NLSy1s are usually observed when L/L Edd is very high; for example, L/L Edd ∼ 10 in the case of RE J1034+396 [8] implies a high accretion rate, but the association of L/L Edd with the QPO frequencies is not clear. Moreover, even if one assumes that the accretion process in AGN and BHXRB is the same and that both show similar characteristic Q shape in the hardness-intensity diagram [6], over a timescale, T, this would be 10 5 −10 6 times more than BHXRB timescales, as T ∝ M • .

3.
Our relativistic orbit model (ROM) is built on the formulation of the intrinsic mechanical energy distribution of the plasma in motion, where three frequencies ν X < ν I < ν M correspond to the low-frequency end, break frequency, and the high-frequency end of the PSD. However, there is a noise component to be added at higher frequencies of the PSD to obtain a more realistic PSD shape to the intrinsic energy distribution related to the frequencies of the unstable orbits inside MBSO. A more generalized approach will be to incorporate frequencies of the more general eccentric and non-planar orbits (e = 0, Q = 0) contributing to the PSD shape. This is planned as future work. 4.
The fundamental frequencies of the spherical geodesics in the Kerr geometry seem to explain the PSD in the Inner Corona (IC) region, where P (ν I < ν < ν M ); whereas the frequencies of the Outer Corona (OC) region are associated with the circular orbits, where P (ν X < ν < ν I ). The results of this toy statistical model, ROM, seem promising. A detailed physical model is required to predict the power law indices in the energy spectrum. Furthermore, including a more ellaborate transfer function taking into account the GR effects like light bending and Doppler boosting, is in order for further study.

5.
The paradigm of the ROM can be tested against observations by extracting {M • , a} from observed {ν X , ν I , ν M , β 1 , β 2 }, and by exploring the parameter space {α 1 , α 2 } which is the basis of the PSD for the ROM model. In the future, we plan to apply and test this model against several observed PSD of various AGN sources. 6.
The total power of a PSD having a power-law profile is given by where X = ν/ν c , τ is the power-law index, and ν c is the upper frequency cut-off of the PSD. On the other hand, from the Wiener-Khinchin theorem, P T = F 2 var ∝ σ 2 − σ 2 N , where F var gives a measure of the time signal variance above the noise and σ 2 N is the variance in the noise measurable from observations. This gives the relation between the measured quantity and ν c as, F var ∝ ν 0.5 c , where the cutoff ν c provides a measure of the spin and mass of the black hole if the disk cuts off at the ISCO or MBSO radius; this implies ν 0.5 c ∝ M 0.5 • . Using a more complicated PSD distribution expected from the ROM and using F var , we can give better estimates for ν c and hence extract {M • , a}, using F var , and study statistical trends from a sample of sources with known {M • , a}.