1. Introduction
The standard theory of accretion discs assumes that accreting material cools efficiently, allowing the energy generated through viscous dissipation to be largely radiated away. This radiative cooling results in disc temperature being significantly lower than the local virial temperature, enabling steady accretion under relatively cool and stable conditions (Shakura & Sunyaev Reference Shakura and Sunyaev1973; Lynden-Bell & Pringle Reference Lynden-Bell and Pringle1974). However, since the early 1990s, this classical paradigm has been fundamentally challenged by the emergence of advection-dominated accretion flows (ADAFs), now widely regarded as key models within the broader category of radiatively inefficient accretion flows (RIAFs; Narayan & Yi Reference Narayan and Yi1994; Nakamura et al. Reference Nakamura, Kusunose, Matsumoto and Kato1997; Narayan et al. Reference Narayan, Kato and Honma1997), which have reshaped our understanding of accretion physics.
ADAFs have attracted considerable attention in astrophysics due to their unique characteristics, which distinguish them from standard accretion discs. One of the most striking features of ADAFs is their significantly lower luminosity (L) compared to standard accretion discs with the same mass accretion rate (
$\dot{M}$
). Specifically, the radiative cooling efficiency (
$\eta = \frac{L}{\dot{M} c^2}$
) in ADAFs typically falls below 0.01, whereas standard accretion discs exhibit efficiencies exceeding 0.1 (Abramowicz et al. Reference Abramowicz, Chen, Granath and Lasota1996). This reduced efficiency arises from the inability of ADAFs to cool radiatively as effectively as standard accretion discs.
In ADAFs, the energy generated by viscous forces is predominantly advected inward by the accreting material rather than being radiated away. Consequently, the gas temperature rises close to the local virial temperature, resulting in several notable structural distinctions. The high thermal energy provides substantial pressure support, making ADAFs geometrically thick and optically thin. Their low density, high temperature and inefficient radiative cooling significantly reduce opacity, allowing most radiation to escape with minimal absorption or scattering.
Additionally, the angular velocity of gas in ADAFs is sub-Keplerian due to increased pressure support (Narayan & Yi Reference Narayan and Yi1994; Ghoreyshi Reference Ghoreyshi2020). This reduces the gravitational force required to maintain rotational equilibrium at Keplerian velocity. These properties distinguish ADAFs from standard accretion discs, which are typically characterised by optically thick, geometrically thin structures and Keplerian angular velocity profiles.
These optically thin ADAFs provide a compelling explanation for interpreting the low-luminosity phenomena observed across various astrophysical systems. Examples include Active Galactic Nuclei (AGNs) such as Sagittarius A* (Nakamura et al. Reference Nakamura, Kusunose, Matsumoto and Kato1997; Satapathy et al. 2023) and cataclysmic variables (Narayan & Popham Reference Narayan and Popham1993), where the accretion processes occur in environments characterised by low radiative efficiency. Observational evidence supporting ADAFs underscore the necessity of deeper investigations into the physics governing these flows, particularly their dynamical behaviour and broader implications for astrophysical models.
Early studies of ADAFs around black holes largely neglected general relativistic effects (Narayan & Yi Reference Narayan and Yi1994; Chakrabarti Reference Chakrabarti1996; Nakamura et al. Reference Nakamura, Kusunose, Matsumoto and Kato1997; Chen et al. Reference Chen, Abramowicz and Lasota1997; Narayan et al. Reference Narayan, Kato and Honma1997). A comprehensive analysis of ADAFs surrounding black holes requires the incorporation of such considerations (Popham & Gammie Reference Popham and Gammie1998 and Gammie & Popham Reference Gammie and Popham1998). Advancements in this field have since led to the development of general relativistic models for optically thin, extremely hot ADAFs surrounding black holes. One of the earliest contributions was by Abramowicz et al. (Reference Abramowicz, Chen, Granath and Lasota1996), who provided solutions describing ADAFs in the gravitational field of a Kerr black hole, although their treatment of viscosity remained non-relativistic. Subsequent studies by Gammie & Popham (Reference Gammie and Popham1998) and Takahashi (Reference Takahashi2007) employed a fully general relativistic framework to describe the shear tensor; however, these models did not incorporate the effects of heat flux. Nonetheless, given the suppression of radiative cooling in ADAFs, it is important to consider non-radiative energy transport mechanisms, such as thermal conduction, which can significantly modify the energy budget and influence the flow structure.
Developing a more realistic and comprehensive modelling for ADAFs near black holes requires incorporating the relativistic effects of both viscosity and heat flux, as their interplay is crucial for a complete theoretical understanding of ADAF physics. The influence of heat flux on the dynamics of ADAFs has long been investigated in numerous studies (Page & Thorne Reference Page and Thorne1974; Beloborodov et al. Reference Beloborodov, Abramowicz and Novikov1997; Riffert & Herold Reference Riffert and Herold1995; Manmoto Reference Manmoto2000; Compere & Oliveri Reference Compere and Oliveri2017; Potter Reference Potter2021; Faraji & Hackmann Reference Faraji and Hackmann2020). However, the majority of these works have exclusively focused on the vertical component of the heat flux, often treating it as the only non-zero component of the heat flux four-vector.
In this paper, we aim to investigate the stationary axisymmetric structure of an ADAF, taking into account all components of the viscous shear tensor and the heat flux four-vector. We will evaluate the relative importance of these components in thermal energy transfer, with particular attention to how they compare to the vertical component, which has traditionally been the primary focus of study. Additionally, we will explore the fully general relativistic formulation of the shear tensor to provide a more comprehensive understanding of ADAF dynamics.
The next sections are organised as follows: Section 2 outlines our model, including the underlying assumptions and background geometry, while Section 3 details the governing dynamical equations. Section 4 presents analytical solutions that provide suitable models for the velocity components and temperature distribution. Section 5 presents the results and their physical interpretation, and Section 6 provides a comparison with numerical studies. Finally, Section 7 summarises the key findings of the study and discusses open problems for future research.
2. The model
We are interested in the low-luminosity accretion flows formed by plasma surrounding a non-rotating black hole. Due to radiatively inefficient cooling, the viscous energy dissipated during the accretion process does not have the opportunity to escape from the disc’s surface in the form of radiation. Instead, this energy is largely advected inward with the accreting plasma, resulting in a hot, optically thin, and geometrically thick flow, commonly described by the advection-dominated accretion flow (ADAF) model.
The presence of a strong gravitational field near the black hole significantly distorts the surrounding space-time, influencing both the geometry of the flow and the underlying physical processes. To incorporate these relativistic effects in a self-consistent way, we adopt the
$3+1$
formalism of general relativity. This approach decomposes the four-dimensional curved space-time into three spatial dimensions and one time dimension, denoted by coordinates
$(x^0,x^1,x^2,x^3) = (t,r,\theta,\phi)$
.
2.1. Assumptions
The key ingredients and assumptions of our model are as follows:
-
(I) Spherical coordinates: We employ standard spherical polar coordinates for the spatial dimensions, with the black hole located at the origin. The z-axis is aligned with the rotation axis of the accretion disc.
-
(II) Axisymmetric and stationary plasma: The plasma is assumed to be axisymmetric and stationary, satisfying
$\frac{\partial}{\partial \phi}=\frac{\partial}{\partial t}=0$
. Additionally, it is symmetric about the equatorial plane. -
(III) Geometric units: We adopt geometric units where the speed of light
$c=1$
, the universal gravitational constant
$G=1$
and the mass of the central black hole
$M=1$
. Consequently, the characteristic length scale is
$m=GM/c^{2}=1$
. -
(IV) Negligible self-gravity: The self-gravity of the disc is ignorable compared to the strong gravity of the central black hole.
-
(V) Geometrically thick disc: ADAFs are generally geometrically thick, characterised by
$\frac{H}{r}\sim 1$
, where H represents the vertical thickness of the disc, and r is the radial distance from the central black hole. Furthermore, we assume no meridional velocity (
$u^{\theta}=0$
) and no meridional variation (
$\frac{\partial}{\partial \theta}=0$
) of the fluid. -
(VI) Vertical averaging: We apply the vertical averaging approximation for all physical flow variables f, which is expressed as
(1)where
\begin{equation}\int d\theta d\phi\sqrt{g}f \simeq 4\pi H_{\theta} f(\theta=\pi/2),\end{equation}
$H_{\theta}$
represents the vertical half-thickness (angular scale height) of the accretion flow near the equator and
$\sqrt{g}$
is the metric determinant for curved space-time geometries. This approximation assumes that variations in the flow along the
$\theta$
-direction are small compared to those in the radial directions, allowing the flow dynamics to be effectively captured by evaluating f at the equatorial plane (
$\theta=\frac{\pi}{2}$
). This reduction significantly simplifies the mathematical treatment of the problem while preserving the essential physics of the system (Gammie & Popham Reference Gammie and Popham1998; Takahashi Reference Takahashi2007).
-
(VII) Radial dependency: The accretion flow is treated as a one-dimensional problem, with all physical quantities being functions only of the radial distance r from the central black hole.
-
(VIII) Boyer–Lindquist coordinates: All equations are formulated in Boyer–Lindquist coordinates.
-
(IX) Sub-Keplerian rotation: The angular velocity of the accreting plasma is assumed to be sub-Keplerian, meaning that it is lower than the angular velocity of a test particle in a stable circular orbit around the black hole, where gravitational and centrifugal forces are exactly balanced. In ADAFs, however, the flow is optically thin and radiative cooling is highly inefficient, causing the viscous energy to remain trapped within the gas. This leads to a significant rise in temperature and pressure, making the disc geometrically thick and pressure-supported. As a result, less centrifugal force is required to counteract gravity, and the plasma rotates more slowly than in a thin, Keplerian disc.
-
(X) Weak magnetic fields: We assume that the magnetic field is dynamically weak, i.e. its pressure and energy density are small compared to those of the thermal plasma. Although magnetic fields could likely have played a crucial role in angular momentum transport through processes such as the magnetorotational instability (MRI), we do not include them in the present model in order to maintain a simplified description of the plasma. This assumption allows us to focus on the relativistic and thermal aspects of the flow without the added complexity of magnetic field dynamics. As a result, angular momentum transport is modelled through shear viscosity, without invoking magnetohydrodynamic (MDH) effects.
2.2. Background geometry
In the
$3+1$
formulation, the strong gravitational field of a non-rotating black hole gives rise to a curved space-time described by the Schwarzschild metric. The line element in this stationary and spherically symmetric geometry is given by
where the covariant components of the metric tensor are
$g_{\rm tt}=-(1-\frac{2}{r})$
,
$g_{\rm rr}=(1-\frac{2}{r})^{-1}$
,
$g_{\theta\theta}=r^{2}$
and
$g_{\phi\phi}=r^{2} \sin^{2}\theta$
. The Christoffel symbols
$\Gamma^{\alpha}_{\mu\nu}=\Gamma^{\alpha}_{\nu\mu}=\frac{1}{2}g^{\alpha\rho}(\partial_{\mu}g_{\nu\rho}+\partial_{\nu}g_{\rho\mu}-\partial_{\rho}g_{\mu\nu})$
associated with the Schwarzschild metric at the equatorial plane (
$\theta=\frac{\pi}{2}$
) are given by
\begin{eqnarray*}&& \Gamma ^{t}_{\rm rt}=\frac{1}{r^{2}-2r},\,\Gamma ^{r}_{\rm tt}=\frac{r-2}{r^{3}},\, \Gamma ^{r}_{\rm rr}=-\frac{1}{r^{2}-2r},\, \Gamma ^{r}_{\theta\theta}=-(r-2),\nonumber\\&& \Gamma ^{r}_{\phi\phi}=-(r-2),\, \Gamma ^{\theta}_{r\theta}=\frac{1}{r},\,\Gamma ^{\phi}_{r\phi}=\frac{1}{r}.\end{eqnarray*}
It is useful here to introduce some key kinematic quantities used in relativistic hydrodynamics: the projection tensor
$h^{\mu\nu}$
, the 4-velocity
$u^{\nu}$
, and its covariant derivatives including the fluid world line expansion
$\Theta$
, and the 4-acceleration
$a^{\mu}$
,
3. Basic dynamical equations
Perfect fluids have been particularly successful in general relativistic MHD (GRMHD) studies (Gammie et al. Reference Gammie, McKinney and Toth2003; McKinney Reference McKinney2006; McKinney Reference McKinney2004; McKinney & Gammie Reference McKinney and Gammie2004; Witzany, & Jefremov Reference Witzany and Jefremov2018). However, including dissipative effects such as viscosity and heat conduction provides a more complete view of plasma phenomena.
We are interested in investigating a non-magnetised, viscous plasma accreting onto a non-rotating black hole in the form of an ADAF taking into account the effects of heat flux.
wherein
$T^{\mu\nu}_{\rm Fluid}=\rho\eta u^{\mu}u^{\nu}+pg^{\mu\nu}$
,
$T^{\mu\nu}_{\rm Visc}=t^{\mu\nu}$
and
$T^{\mu\nu}_{\rm Heat}=q^{\mu}u^{\nu}+u^{\mu}q^{\nu}$
are the energy-momentum tensors for the perfect fluid, viscosity and heat flux, respectively. Here,
$\rho$
denotes the rest-mass density,
$\eta=\frac{\rho+p+u}{\rho}$
represents the relativistic enthalpy, u is the internal energy, p stands for the gas pressure,
$t^{\mu\nu}$
is the viscous stress-energy tensor and
$q^{\mu}$
is the heat flux 4-vector.
3.1. Heat flux four-vector and viscous stress tensor
In the study of relativistic accretion discs, the treatment of heat flux is crucial for understanding the thermodynamic behaviour and energy transport within the flow. Several models have been proposed to describe the heat transfer mechanisms, often focusing on simplified assumptions to highlight dominant effects. Among these, the vertical component of heat flux has been widely used due to its relative simplicity and the common belief that it plays the dominant role in thermal energy transport. For example, Page & Thorne (Reference Page and Thorne1974) and Faraji & Hackmann (Reference Faraji and Hackmann2020) modelled the vertical heat flux using the form
$q_z= F(r) \frac{z}{H(r)}$
, where
$H(r)=\frac{H_{\theta}}{r}$
is the disc’s relative half-thickness, and F(r) is the time-averaged radiant energy flux. Alternatively, Potter (Reference Potter2021) adopted a surface-emission-based formulation, expressing the vertical heat flux as
$q_z=\sigma T_s^4$
, where
$T_s$
is the surface temperature of the disc and
$\sigma$
is the Stefan–Boltzmann constant. In the present study, however, we aim to go beyond these simplified models by considering all components of the heat flux vector using the general relativistic formulation given by Eckart (Reference Eckart1940) and Misner et al. (Reference Misner, Thorne and Wheeler1973):
where
$\kappa$
is the thermal conductivity, T is the temperature, and
$a_{\nu}=g_{\nu\mu} a^{\mu}$
represents the covariant acceleration. In the absence of bulk viscosity, the viscous stress tensor
$t^{\mu\nu}$
is related to the shear tensor
$\sigma^{\mu\nu}$
via the dynamic viscosity coefficient
$\lambda$
as
where the shear tensor
$\sigma^{\mu\nu}$
, describing the rate of shear deformation in the flow, is defined by
3.2. Conservation equations
The motion of a viscous plasma is governed by the baryon mass conservation law (or continuity equation)
$(\rho u^{\mu})_{;\mu}=0$
, and the energy-momentum conservation law
$T^{\mu\nu}_{;\nu}=0$
. By expanding these conservation laws and projecting them along appropriate directions, we derive the fundamental dynamical equations that govern the behaviour of the accretion flow.
3.2.1. Continuity equation
The particle number conservation equation can be written as
where
$g=|Det (g_{\mu\nu})|=r^4$
. Applying the vertical averaging approximation (Equation 1) to the above expression yields
Integrating once with respect to the radial coordinate gives the vertically averaged form of the mass conservation (continuity) equation
where
$\dot{M}$
is a constant of integration identified as the mass accretion rate.
3.2.2. Mass-energy flux
The mass-energy flux is obtained by projecting the energy-momentum conservation law
$T^{\mu\nu}_{;\nu}=0$
, onto the temporal component using the time-like Killing vector
$\varepsilon^{\mu}_{t}=(1,0,0,0)$
. This yields
Applying the vertical averaging approximation to the equation above leads to
where
$\dot{E}$
represents the mass-energy flux, corresponding to the rate of increase in the black hole mass as measured at the event horizon. In this study, we adopt the assumption
$\dot{E}\approx \dot{M}=1$
, which is justified in the limit of cold and slow accretion flows at large radii (see, e.g. Gammie & Popham Reference Gammie and Popham1998; Popham & Gammie Reference Popham and Gammie1998).
3.2.3. Radial momentum conservation
The radial momentum conservation equation is derived by projecting the energy–momentum conservation law,
$T^{\mu\nu}_{;\nu}=0$
, onto the radial direction using the projection tensor
$h^r_{\mu}$
, such that
$h^r_{\mu} T^{\mu\nu}_{;\nu}=0$
. Expanding the energy–momentum tensor, we have
Applying the continuity Equation (5) to simplify the terms involving
$\rho$
, we obtain
This leads to the following expression for the radial momentum conservation
where the term
$n_{\rm HI}=-\frac{h^{r}_{\mu}}{\rho\eta}[ t^{\mu\nu}_{;\nu}+(q^{\mu}u^{\nu}+u^{\mu}q^{\nu}) _{;\nu}] $
accounts for the effect of heat inertia, as discussed in Beloborodov et al. (Reference Beloborodov, Abramowicz and Novikov1997). Thus, the radial momentum conservation equation can be compactly written as
3.2.4. Angular momentum conservation
The equation for the angular momentum conservation can be derived by projecting the energy-momentum conservation law
$T^{\mu\nu}_{;\nu}=0$
, onto the azimuthal direction using the Killing vector
$\varepsilon^{\mu}_{\phi}=(0,0,0,1)$
. This yields
Applying the vertical averaging approximation (Equation 1), we obtain the vertically averaged angular momentum conservation equation
where
$\dot{M}j$
represents the total inward flux of angular momentum and j is the specific angular momentum.
3.2.5. Energy conservation
Projecting the energy–momentum conservation law
$T^{\mu\nu}_{;\nu}=0$
, along the fluid four-velocity
$u_{\mu}$
leads to the energy conservation equation
$u_{\mu} T^{\mu\nu}_{;\nu}=0$
. After standard manipulations and simplifications, this yields
where
$q^{+}_{\rm vis}$
and
$q^{-}_{\rm rad}$
represent the viscous heating and radiative cooling rates, respectively, and are defined as
It is worth noting that, in ADAFs, the difference between viscous heating and radiative cooling represents the advected energy. Specifically, the quantity
$q_{\rm adv}=q^{+}_{\rm vis}-q^{-}_{\rm rad}$
denotes the net thermal energy carried inward by the accreting plasma, rather than being radiated away.
4. Analytical solution
By substituting the continuity Equation (6) into the mass–energy flux Equation (8) and the angular momentum conservation Equation (11), we obtain the modified forms of the mass-energy flux
and the angular momentum conservation equation
Combining these two equations, we solve for
$\rho$
and
$\eta$
in terms of the heat flux and viscous terms as
The stress-energy components
$t^{r}_{t}$
and
$t^{r}_{\phi}$
are expressed using the shear tensor components
with the shear tensor components
$\sigma_{rt}$
and
$\sigma_{r\phi}$
given by
Expanding the general relativistic heat flux expression (Equation 2), the components of
$q^{\mu}$
are obtained as
\begin{eqnarray}&&q^{t}=- \kappa\left[h^{\rm rt}\frac{dT}{dr}+T(h^{\rm tt}a_{t}+h^{\rm rt}a_{r}+h^{t\phi}a_{\phi})\right]\!,\nonumber\\ && q^{r}=- \kappa\left[h^{\rm rr}\frac{dT}{dr}+T(h^{tr}a_{t}+h^{\rm rr}a_{r}+h^{r\phi}a_{\phi})\right]\!,\nonumber\\ && q^{\theta}=0,\nonumber\\ && q^{\phi}=- \kappa\left[h^{r\phi}\frac{dT}{dr}+T(h^{t\phi}a_{t}+h^{r\phi}a_{r}+h^{\phi\phi}a_{\phi})\right]\!.\end{eqnarray}
To compute the heat flux components accurately, appropriate models for the temperature and velocity profiles must be specified.
4.1. Velocity model
The four-velocity is given by
$u^{\mu}=(\gamma,\gamma \mathbf{v})$
, where the Lorentz factor
$\gamma$
is determined from the normalization condition
$u_{\mu} u^{\mu}=-1$
. This leads to
We adopt a power-law model (
$v^{r}\propto r^{-n}$
) for the radial inflow velocity, outside the innermost stable circular orbit (ISCO), following Ozel et al. (Reference Ozel, Psaltis and Younsi2022), Satapathy et al. (Reference Satapathy and Psaltis2023). This ansatz is consistent with the scale-free nature of ADAFs and reflects the fact that far from the boundary conditions imposed at the inner and outer edges of the flow, the dynamics of a hot ADAF are dominated by advective transport of energy. In this regime, the flow remains scale-free over the large radial range. So, the inflow naturally adjusts with radius far from the boundaries in a scale-free manner, and the governing equations admit self-similar power-law behaviour for the radial velocity. Therefore, a self-similar power-law inflow of the form
$v^r= -\frac{\beta}{r^n}$
, where
$\beta$
is a constant, naturally emerges from the underlying dynamics.
To illustrate the radial behaviour of the inflow, we plot the Mach number, defined as
$\mathscr{M}\equiv \frac{u^r}{c_s}$
, where
$u^r$
is the radial component of the 4-velocity and
$c_s$
is the relativistic sound speed, defined as
$c_s=\sqrt{\frac{ p}{\eta\rho}}$
(Das & Chakrabarti Reference Das and Chakrabarti2008; Shapiro & Teukolsky Reference Shapiro and Teukolsky2024).
This definition of Mach number provides a relativistically consistent measure of the inflow velocity relative to the sound speed, enabling analysis of the transonic nature of the flow. Figure 1 demonstrates that the flow undergoes a transonic transition, remaining subsonic at large radii, passing through a sonic point, and becoming supersonic in the inner region. For standard ADAF solutions, the sonic point typically lies at
$r \sim 2$
–
$10$
m (Feng & Ke-liang Reference Feng and Ke-liang1999; Moscibrodzka et al. Reference Moscibrodzka, Das and Czerny2006; Abramowicz et al. Reference Abramowicz, Jaroszynski, Kato, Lasota, Rozanska and Sqdowski2010). Accordingly, we restrict the radial velocity power-law index to the range
$0.5\leq n\leq 0.8$
, for which the sonic point lies within this physically relevant interval.
Radial profiles of the Mach number
$\mathscr{M}$
, for different values of the radial velocity power-law index n. The golden diagonal-cross curve corresponds to
$n=1$
, the green dashed curve to
$n=0.8$
, the purple dash-dotted curve to
$n=0.6$
, the red dotted curve to
$n=0.5$
, the cyan solid curve to
$n=0.4$
, and the dark-green long-dashed curve to
$n=0.3$
. The black solid horizontal line marks the sonic condition,
$\mathscr{M}=1$
.

For the azimuthal velocity, we assume a sub-Keplerian profile
$v^{\phi}=\frac{l}{r^{3/2}}$
, where the angular velocity parameter
$l\lt1$
ensures that the flow remains sub-Keplerian. Substituting these expressions into the Lorentz factor equation yields
Accordingly, the four-velocity becomes
With this form of
$u^{\mu}$
, the associated dynamical quantities, including the four-acceleration
$a^{\mu}$
, the projection tensor
$h^{\mu \nu}$
, the expansion scalar
$\Theta$
, and the shear tensor
$\sigma^{\mu\nu}$
, can now be explicitly computed. In Figure 2, the radial profiles of the non-zero components of the shear tensor are shown. Our general relativistic calculations reveal that the
$\sigma^{\rm rr}$
,
$\sigma^{\rm rt}$
, and
$\sigma^{\rm tt}$
components are larger in magnitude and more influential than the
$\sigma^{r\phi}$
component, which is traditionally regarded as the dominant contributor to disc dynamics and angular momentum transport, particularly in Newtonian frameworks. This finding indicates that a comprehensive understanding of relativistic disc dynamics requires accounting for these terms, as they can strongly influence both angular momentum redistribution and energy dissipation in the disc. Neglecting them in simplified models may lead to an incomplete, or even misleading, description of the flow dynamics.
Radial profiles of the non-zero components of the shear tensor. The cyan dash-dotted curve represents
$\sigma^{\rm rr}$
, the black dotted curve shows
$\sigma^{\rm rt}$
, the gold dashed curve corresponds to
$\sigma^{\rm tt}$
, the solid green curve indicates
$\sigma^{r\phi}$
, the red space-dotted curve denotes
$\sigma^{\phi\phi}$
, and the blue long-dashed curve depicts
$\sigma^{t\phi}$
. The constant parameters are
$n=0.5$
,
$\beta=1$
and
$l=0.8$
.

In contrast, the vertical component
4.2. Temperature model (Ion-dominated plasma)
In a radiatively-inefficient, hot, dilute, and weakly collisional ADAF, the coulomb energy-exchange timescale between ions and electrons exceeds the inflow timescale. As a result, the plasma thermally decouples into a two-temperature structure, in which the ions absorb most of the viscous heating and thus retain the majority of the internal energy and pressure. While the electrons cool efficiently through synchrotron, bremsstrahlung, and inverse-compton emission. Because the ions carry the dominant share of thermal energy and pressure, they control the hydrodynamic behaviour of the flow.
We therefore adopt the ion temperature
$T_i$
as the primary thermodynamic variable in our model, where weak cooling and efficient inward advection drive the ions towards virial temperature equilibrium with the gravitational potential. Electrons, while important for radiative emission, contribute negligibly to the dynamical pressure and are therefore not treated explicitly here.
Under virial balance, the gravitational binding energy per particle
$\frac{G M m_p}{r}$
approximates the thermal energy
$\frac{3}{2} k T_i$
, leading to the scaling
$T_i \approx \frac{1}{r}$
, in agreement with the virial scaling commonly adopted in analytic ADAF models (Akizuki & Fukue Reference Akizuki and Fukue2006). Accordingly, we prescribe the ion temperature as
$T_i=T_0 (\frac{r}{r_0})^{-1}$
, wherein
$T_0$
represents the ion temperature at a chosen reference radius
$r_0$
. With this prescribed ion temperature distribution, all components of the heat-flux four-vector appearing in Equation (22) become fully determined.
5. Results and physical interpretation
The adopted temperature model allows for a complete specification of the heat flux four-vector components, thereby enabling the determination of the remaining physical variables in the flow. We proceed to compute these components explicitly and analyse their relative contributions to the dynamical structure of the flow. The resulting profiles offer valuable insight into the spatial behaviour of these quantities across the disc.
In Figure 3, the radial variations of the different components of the heat flux are presented. The results indicate a clear decline in heat flux with increasing radius, reflecting a gradual reduction in thermal energy transport as the disc cools outward. This trend is consistent with the expected thermodynamic behaviour of ADAFs.
Radial profiles of the non-zero components of the heat flux four vector. The solid green curve represents
$q^r$
, the red dotted curve shows
$q^t$
, the blue dashed curve corresponds to
$q^{\phi}$
, and the gold dash-dotted curve indicates
$q^z$
. The constant parameters are
$n=0.5$
,
$\beta=1$
,
$\lambda=1.5, \kappa=1.5$
,
$l=0.8$
and
$j=1.5$
.

A comparison among these components highlights the dominant role of the radial heat flux in the overall energy transport. In contrast, the vertical component,
$q^z$
, defined through
$q^z=q^r \cos \theta$
with
$\cos \theta= \frac{H_{\theta}}{r}$
, which was previously assumed to contribute (Page & Thorne Reference Page and Thorne1974; Faraji & Hackmann Reference Faraji and Hackmann2020; Potter Reference Potter2021), appears to have a minimal influence throughout the disc. Additionally, the azimuthal component,
$q^{\phi}$
, is comparatively weak, suggesting that azimuthal transport of thermal energy plays only a minor role in the disc’s energy redistribution.
In addition to the spatial components of the heat flux responsible for transporting thermal energy across the disc, there exists a temporal component
$q^t$
, whose role – although not spatial in nature – is more significant than that of the azimuthal
$q^{\phi}$
and vertical
$q^z$
components. While
$q^t$
does not contribute to the directional flow of heat, it represents the local energy density of the heat flux perceived by a local observer in the comoving frame and is a purely relativistic effect. Our results show that
$q^t$
consistently lies between the dominant radial component
$q^r$
and the weaker components
$q^z$
and
$q^{\phi}$
in magnitude.
This anisotropy in the heat flux components underscores the importance of considering direction-dependent conduction, particularly in relativistic regimes where space-time curvature and flow geometry can strongly influence transport processes. Overall, these results emphasise that in the context of ADAFs, accurate modelling of the radial component of heat flux is essential for capturing the dominant mechanisms of energy transfer. Neglecting or oversimplifying this component may lead to significant inaccuracies in predicting the thermal structure and dynamics of the accretion flow. Furthermore, the non-negligible contribution of the temporal component of heat flux emphasises the necessity of incorporating the full relativistic form of the heat flux vector to accurately capture the energy dynamics in ADAFs.
Figures 4 and 5 present the radial profiles of the system’s physical variables. All variables exhibit a declining trend with increasing radius, consistent with the expected outward cooling behaviour of the disc. In the inner regions, where the disc is hotter and more thermodynamically active, the influence of thermal conductivity becomes more pronounced (Figure 4).
Radial variations of the physical variables – rest-mass density
$\rho$
, relativistic enthalpy
$\eta$
, vertical half-thickness
$H_{\Theta}$
, pressure P, and sound speed
$c_{s}$
– for different values of thermal conductivity
$\kappa$
. The lines correspond to: solid blue for
$\kappa=2$
, red dotted for
$\kappa=1$
, and green dashed for
$\kappa=0$
. The constant parameters are
$n=0.5$
,
$\beta=1$
,
$\lambda=2$
,
$l=0.8$
, and
$j=1.5$
.

As the thermal conductivity parameter
$\kappa$
increases, the relativistic enthalpy
$\eta$
, vertical half-thickness
$H_{\theta}$
, pressure P, and sound speed
$c_s$
are elevated, indicating enhanced thermal transport and increased internal energy. Conversely, the rest-mass density
$\rho$
decreases in the inner regions, suggesting that higher thermal conduction leads to a more expanded and less dense medium.
This implies that thermal conduction plays a crucial role mainly in the inner, high-temperature regions of ADAFs, where it reduces steep temperature and pressure gradients by redistributing heat more efficiently. As a result, it enhances thermal pressure support, helping to maintain a thicker vertical structure in these regions. In summary, Figure 4 demonstrates that thermal conduction can substantially influence both energy transport and the structural configuration of the ADAFs, particularly in the inner regions where thermal gradients are steep and advection dominates the energy flow. On the other hand, at larger radii – where the disc becomes cooler and advection less significant – the effect of the thermal conductivity diminishes, and the profiles of the physical variables tend to converge regardless of the value of
$\kappa$
.
In contrast to thermal conductivity, which primarily influences the inner regions, viscosity affects the dynamical structure throughout the entire disc. An increase in the dynamic viscosity coefficient
$\lambda$
, which governs the strength of viscous effects in the system, leads to an increase in both the rest-mass density
$\rho$
and pressure P. However, the vertical half-thickness
$H_\theta$
decreases, while the relativistic enthalpy
$\eta$
and sound speed
$c_s$
are only weakly influenced (Figure 5). This behaviour suggests that enhanced viscosity tends to compress the disc vertically while increasing its pressure and density. It is noteworthy that viscosity not only facilitates angular momentum transport, enabling accretion to proceed, but also contributes to energy generation through dissipation. The observed decrease in the vertical half-thickness with increasing viscosity points to the essential characteristic of ADAFs that the energy generated by viscous dissipation is not efficiently radiated away but is instead advected inward with the flow. This generated heat does not significantly contribute to vertical expansion, resulting in a denser and more vertically compressed disc structure.
Radial variations of the physical variables – rest-mass density
$\rho$
, relativistic enthalpy
$\eta$
, vertical half-thickness
$H_{\Theta}$
, pressure P, and sound speed
$c_{s}$
– for different values of the dynamic viscosity coefficient
$\lambda$
. The lines correspond to: solid blue for
$\lambda=2$
, red dotted for
$\lambda=3$
and green dashed for
$\lambda=4$
. The constant parameters are
$n=0.5$
,
$\beta=1$
,
$\kappa=2$
,
$l=0.8$
, and
$j=1.5$
.

Figure 5. Long description
The image contains five line graphs showing radial variations of physical variables around black holes. Each graph compares different values of the dynamic viscosity coefficient. Panel A: The line graph shows the rest-mass density (ρ) against the radial distance (r/m). The solid blue line represents a specific value, the red dotted line represents another value, and the green dashed line represents a third value. Panel B: The line graph shows the relativistic enthalpy (h) against the radial distance (r/m). The lines follow the same color scheme as Panel A. Panel C: The line graph shows the vertical half-thickness (H_θ) against the radial distance (r/m). The lines follow the same color scheme as Panel A. Panel D: The line graph shows the pressure (P) against the radial distance (r/m). The lines follow the same color scheme as Panel A. Panel E: The line graph shows the sound speed (c_s) against the radial distance (r/m). The lines follow the same color scheme as Panel A. The constant parameters for all graphs are specified but not visually detailed in the image.
Overall, Figures 4 and 5 highlight the complementary roles of viscosity and conduction in shaping the structure of ADAFs. While viscosity acts globally by transporting angular momentum and generating internal energy throughout the disc, thermal conduction operates more locally, particularly in the hot inner regions, by redistributing this energy and regulating steep thermal gradients. Therefore, a fully general relativistic, self-consistent treatment of both mechanisms is essential for accurately modelling the thermal and structural behaviour of ADAFs.
The radial profiles of the radiative cooling rate
$q^{-}_{\rm rad}$
, viscous heating rate
$q^{+}_{\rm vis}$
, and advective energy transport rate
$q_{\rm adv}$
are shown in Figures 6 and 7. All three rates, consistently diminish with radial distance, indicating a progressive reduction in thermal activity as the disc becomes cooler and less dense at larger radii. The plots quantitatively confirm the inefficiency of radiative cooling compared to the viscously generated heat. This imbalance clearly highlights the advection-dominated nature of the flow, where a significant portion of the viscously generated energy is advected inward rather than radiated away.
Radial variations of radiation cooling rate
$q_{\rm rad}$
, viscous heating rate
$q_{\rm vis}$
and advection rate
$q_{\rm adv}$
from left to right, illustrating the effects of the radial velocity parameter n. The curves correspond to: solid blue for
$n=0.5$
, red dotted for
$n=0.6$
, and green dashed for
$n=0.8$
. The constant parameters are
$\beta=1$
,
$\kappa=2$
,
$\lambda=2$
,
$l=0.8$
and
$j=1.5$
.

Same as Figure 6, but showing the effects of the angular velocity parameter l. The curves correspond to: solid blue for
$l=0.4$
, red dotted for
$l=0.6$
, and green dashed for
$l=.8$
. The constant parameters are
$n=0.5$
,
$\beta=1$
,
$\kappa=2$
,
$\lambda=2$
and
$j=1.5$
.

Figure 7. Long description
Panel A: A line graph shows the radial energy advection rate (q_rad) as a function of radius (r) in meters. The x-axis represents the radius (r) in meters, and the y-axis represents the radial energy advection rate (q_rad). The graph includes three curves: a solid blue line, a red dotted line, and a green dashed line. The curves show a steep decline in q_rad as the radius increases, with the solid blue line showing the highest values, followed by the red dotted line and the green dashed line. Panel B: A line graph shows the viscous energy dissipation rate (q_vis) as a function of radius (r) in meters. The x-axis represents the radius (r) in meters, and the y-axis represents the viscous energy dissipation rate (q_vis) on a logarithmic scale. The graph includes three curves: a solid blue line, a red dotted line, and a green dashed line. The curves show a linear decline in q_vis as the radius increases, with the solid blue line showing the highest values, followed by the red dotted line and the green dashed line. Panel C: A line graph shows the energy advection rate (q_adv) as a function of radius (r) in meters. The x-axis represents the radius (r) in meters, and the y-axis represents the energy advection rate (q_adv) on a logarithmic scale. The graph includes three curves: a solid blue line, a red dotted line, and a green dashed line. The curves show a linear decline in q_adv as the radius increases, with the solid blue line showing the highest values, followed by the red dotted line and the green dashed line.
The effect of radial inflow velocity on the thermal rates is examined in Figure 6 through the radial velocity parameter n, while Figure 7 explores the influence of angular velocity via the angular momentum parameter l. Both a decrease in n, which indicates a faster radial inflow (Figure 6), and an increase in l, which corresponds to a more rapidly rotating disc (Figure 7), lead to an enhancement of all thermal rates, including radiative cooling, viscous heating, and advective transport. This enhancement arises from the stronger and more accelerated inward motion of the accreting material, which intensifies compression, shear, and energy transport within the flow.
6. Comparison with numerical studies
We compare our analytical results with a representative set of numerical studies that span the principal numerical approaches used to investigate hot relativistic accretion flows, including ideal-GRMHD simulations (Foucart et al. Reference Foucart2017; Dhruv et al. Reference Dhruv, Prather, Wong and Gammie2025) and extended relativistic fluid simulations that incorporate anisotropic heat conduction (Chandra et al. Reference Chandra, Foucart and Gammie2017).
While our analytical model computes the relativistic heat flux four-vector from temperature gradients in a stationary background flow, Chandra et al. (Reference Chandra, Foucart and Gammie2017) consider a weakly collisional, magnetised plasma in which the heat flux is anisotropic and constrained to flow strictly parallel to magnetic field lines. This anisotropy fundamentally alters the stability properties of the flow and allows heat transport to actively drive buoyant motions through instabilities such as the magneto-thermal instability and the heat flux driven buoyancy instability. In that framework, heat flux contributes directly to the energy-momentum tensor and plays an active role in the dynamics.
By contrast, the GRMHD simulations of Foucart et al. (Reference Foucart2017) and Dhruv et al. (Reference Dhruv, Prather, Wong and Gammie2025) adopt the ideal fluid approximation and do not include an explicit conductive heat flux. In these simulations, thermal energy is transported primarily by advection and turbulent dissipation, while angular momentum transport arises self-consistently from magnetic stresses generated by MRI-driven turbulence. The global structure and variability of the accretion flow are therefore controlled mainly by magnetic fields and turbulence rather than by conductive transport.
Our analytical results are complementary to existing numerical studies. Rather than modelling the full time-dependent evolution of magnetised accretion flows, our approach isolates the impact of relativistic heat conduction on the stationary thermodynamic structure of the disc. The solutions obtained here can therefore be regarded as equilibrium background states, which capture the mean thermodynamic structure established by relativistic heat conduction. Time-dependent numerical simulations may then be viewed as describing fluctuations and instabilities that develop around this background due to magnetic stresses, turbulence, and anisotropic transport effects.
7. Conclusion
In this paper, we have developed an axisymmetric, stationary, one-dimensional model of advection-dominated accretion flows (ADAFs) around Schwarzschild black holes by incorporating a fully general relativistic formulation of both the viscous shear tensor and the complete heat flux four-vector. In contrast to previous studies, which either neglected a fully relativistic treatment or considered only the vertical component of the heat flux, our model accounts for a fully general relativistic framework and also all components of the heat flux vector. However, our results indicate that the vertical component is, in fact, relatively insignificant compared to the radial and temporal components. This comprehensive approach enables a more accurate representation of anisotropic heat transport and its interaction with viscous dissipation in relativistic accretion flows.
Our analysis demonstrates that thermal conduction significantly influences the dynamical structure of ADAFs, especially in the hotter inner regions of the disc where thermal gradients are steep and advection is strong. As the thermal conductivity parameter
$\kappa$
increases, the relativistic enthalpy, pressure, sound speed, and vertical thickness all rise, while the density decreases. These changes reveal that conduction enhances internal energy and thermal pressure support, leading to a more expanded, vertically extended, and less dense structure. At larger radii, where the disc becomes cooler and less dynamically active, the effect of thermal conduction diminishes, and the profiles of physical quantities converge for different values of
$\kappa$
. Unlike thermal conduction, which predominantly influences the inner disc regions, viscosity impacts the dynamical structure throughout the entire disc. Enhancing viscosity leads to increased rest-mass density and pressure, accompanied by a decrease in vertical half-thickness, resulting in a denser and more vertically compressed disc. The results for conduction and viscosity highlight their complementary roles in shaping disc structure: viscosity governs angular momentum transport and energy generation on a global scale, whereas conduction redistributes thermal energy locally, particularly in the hotter inner regions. A comprehensive and fully relativistic treatment of both mechanisms is therefore essential for accurately capturing the thermal, structural, and dynamical properties of accretion flows around black holes.
The energy transport analysis confirms the advection-dominated nature of the flow, showing that radiative cooling is consistently weaker than viscous heating, and a substantial portion of the generated energy is transported inward by advection. Moreover, both faster radial inflow (smaller n) and higher angular momentum (larger l) amplify all thermal processes, namely cooling, heating, and advection, due to enhanced compression and shear in the flow.
A detailed comparison of the heat flux components suggests a strong anisotropy. The radial component
$q^r$
dominates the thermal transport throughout the disc. In contrast, the vertical component
$q^z$
, often emphasised in earlier works, is found to be negligible. Notably, the temporal component
$q^t$
, though non-spatial, contributes more significantly than both
$q^z$
and
$q^\phi$
, indicating its essential role as a relativistic measure of local heat energy density in the comoving frame.
These results emphasise the importance of including all components of the heat flux four-vector and shear stress tensor in relativistic ADAF models. Our analysis reveals that
$\sigma^{\rm rr}$
,
$\sigma^{\rm rt}$
, and
$\sigma^{\rm tt}$
can exceed
$\sigma^{r\phi}$
in both magnitude and dynamical influence, challenging the long-standing assumption of
$\sigma^{r\phi}$
dominance in disc evolution and angular momentum transport. Simplified treatments that neglect spatial anisotropy or relativistic effects may overlook key aspects of energy transport and misrepresent the thermal structure of the flow.
Despite these advances, several questions remain open that are expected to be addressed in the continuation of this work in future research. Incorporating MHD effects into a fully relativistic conduction framework remains a challenging but necessary next step for a more complete understanding of ADAFs. Magnetic fields are believed to play a central role in angular momentum transport and energy dissipation, and their interaction with anisotropic conduction could significantly alter the dynamics and thermal structure of the flow. Moreover, extending the model to Kerr geometry would allow for the investigation of black hole spin effects and relativistic frame-dragging on heat transport and disc behaviour. Addressing these issues will be crucial for building more realistic models and for connecting theoretical predictions with observations of accreting black hole systems.
Our analytical results can naturally be interpreted within the context of existing numerical studies of relativistic accretion flows. By isolating relativistic heat conduction, our analysis clarifies its role in setting the stationary thermodynamic structure of hot accretion flows. The solutions presented here may therefore be regarded as physically motivated equilibrium background states, providing a well-defined reference configuration for time-dependent numerical simulations in which magnetic stresses, turbulence, and heat flux driven instabilities are subsequently introduced.
Acknowledgements
We thank the anonymous referee for constructive comments that helped improve the clarity and presentation of this work.
Data availability statement
The data supporting the findings of this study were generated using Maple software. All relevant data are provided within the paper and its cited references, and no additional datasets were generated or analysed for this work.

M
n=1
n=0.8
n=0.6
n=0.5
n=0.4
n=0.3
M=1
σrr
σrt
σtt
σrϕ
σϕϕ
σtϕ
n=0.5
β=1
l=0.8
qr
qt
qϕ
qz
n=0.5
β=1
λ=1.5,κ=1.5
l=0.8
j=1.5
ρ
η
HΘ
cs
κ
κ=2
κ=1
κ=0
n=0.5
β=1
λ=2
l=0.8
j=1.5
ρ
η
HΘ
cs
λ
λ=2
λ=3
λ=4
n=0.5
β=1
κ=2
l=0.8
j=1.5
qrad
qvis
qadv
n=0.5
n=0.6
n=0.8
β=1
κ=2
λ=2
l=0.8
j=1.5
l=0.4
l=0.6
l=.8
n=0.5
β=1
κ=2
λ=2
j=1.5