Dynamics of Vortex Structures: From Planets to Black Hole Accretion Disks

: Thermo-vortices (bright spots, blobs, swirls) in cosmic fluids (planetary atmospheres, or even black hole accretion disks) are sometimes observed as clustered into quasi-symmetrical quasi-stationary groups but conceptualized in models as autonomous items. We demonstrate—using the (analytical) Sharp Boundaries Evolution Method and a generic model of a thermo-vorticial field in a rotating “thin” fluid layer in a spacetime that may be curved or flat—that these thermo-vortices may be not independent but represent interlinked parts of a single, coherent, multi-petal macro-structure. This alternative conceptualization may influence the designs of numerical models and image-reconstruction methods.


Introduction
Recently, a highly complex analysis of observations of the black hole Sagittarius A* (Sgr A*) yielded an impressive image (see Figure 1, panel A) that suggests the possible presence of three bright spots within its structure outside the shadow of a black hole (see Refs. [1][2][3][4][5][6]).Visually, this phenomenon is resemblant of the formations already observed in the atmospheres of planets (also shown in Figure 1, panels B-D).Indeed, regular in shape, long-lasting vortices and multi-vortex structures appear frequently in nature; some well-known examples are the Antarctic Polar Vortex, oceanic vortices, cyclones and anticyclones in atmospheres, cyclone storms on Jupiter poles, and analogous vortex structures in the atmospheres of outer planets.If further observational studies confirm that the bright spots in Figure 1A are not artifacts of the image-reconstruction algorithm but are indeed physical in nature (apparently associated with temperature anomalies/emissions in localized areas), then the list of known cosmic environments with multi-vortex structures would be expanded to include black hole accretion disk flows.
All vortex structures are naturally studied using frameworks and models that seem most appropriate for their native environments.Typically, however, the individual vortices are conceptualized as stand-alone items, sometimes interacting, sometimes not.
In this paper-while contemplating the image presented in Figure 1A and its visually bright spots (although the only significant real feature in the image is the shadow and its diameter, the non-uniformity of the brightness along the ring is not very significant and is subject to processing/reconstruction methods)-we analytically demonstrate that such spots may exist and represent connected "petals" of a single multi-petal thermo-vorticial macrostructure.Within such a structure, the vorticity of the "hydrodynamic" (i.e., collective) flow is interlinked with the temperature field of the medium.This vorticity, a Lagrangian characteristic of the hydrodynamic flow, is transported by the flow, thus creating spots with an elevated temperature.Earth Observatory [8].
The awareness of the existence of two interpretations of the structure and the linkage of thermo-vortices is important.The initial conceptualization-whether the thermo-vortices are treated from the start as stand-alone items or as inextricable parts of one macro-structuremay meaningfully influence the numerical models and image-reconstruction methods.
In our model, we make the following assumptions and impose the following restrictions: the flow of the disk fluid is assumed to be large-scale (the Reynolds number or its analogue is large, i.e., viscous effects are small and may be neglected); the flow is slow (relative to the sound speed); and the assumption of medium incompressibility is acceptable.To simplify the mathematical complexity of analytic solutions, our model considers the "thin"layer case: (1) it is assumed that the accretion disk is relatively thin, i.e., the characteristic size of the structure is assumed to be significantly larger than the characteristic thickness of the disk (layer) at the location at which the structure is localized; (2) the motion of the medium is assumed to be quasi-two-dimensional, i.e., the component of the collective flow velocity in the direction perpendicular to the disk plane is suppressed or significantly smaller in magnitude than the characteristic flow velocity in the disk plane (in this context, a "thick" system may also, in principle, be described using the two-dimensional treatment if a vorticity patch is broken into columnar sub-patches, as in Refs.[9,10]); (3) the temperature decreases towards the center of the disk due to the assumption that the heating of the disk mainly occurs not from the hot central body but from the periphery where the capture and destruction of captured bodies occurs due to the tidal effect.
The medium (disk layer) in our model is presumed to rotate.Vorticity is non-zero where there is velocity shear.We describe the "average" vorticity, representative of some "band" within the disk, using the symbol ω. (Obviously, this vorticity ω is not the same as the spin of the central body or the characteristic angular velocity Ω of the disk's rotation.)When Ω is "sufficiently" high, the centrifugal effects (in the rotating frame) become prominent in the leading terms, and even dominate the gravitational effect of attraction to the central body.
To avoid any confusion, for the modelers accustomed to working with trajectories for particles, let us emphasize that we work with the field, not with individual particles (or their trajectories or orbits).Perhaps what may help the reader grasp this nuance better is the reminder that the velocity of the displacement of electrons in a usual house-wire is not the same thing as the speed of propagation of the electro-magnetic field perturbation along that same wire.
Additionally, we limit our consideration to the case when the magnetic field is weak and the plasma is non-relativistic.This means that in our treatment we neglect the effect of the magnetic field on the large-scale "hydrodynamics" of the flow, i.e., in the full tensor of the energy-momentum of the macroscopic "hydrodynamic" flow combined with the electromagnetic field, we neglect the magnetic part.However, this does not mean that this "weak" magnetic field has no impact on the electromagnetic radiation (for example, on the polarization of the synchrotron radiation).Indeed, as Ref. [11] indicates, EHT images in linear polarization have identified a coherent spiral pattern around the black hole, produced from ordered magnetic fields threading the emitting plasma.
Finally, our focus is on steady-state thermo-hydrodynamic structures.This is important to keep in mind for numerical modelers who typically initiate their models at time t = 0 (and may start, for example, with the spherical Bondi accretion) and allow their models to evolve towards some steady state (hoping to minimize accumulated errors along the way).In contrast, our steady-state is derived analytically; we jump straight to the model of a flattened rotating medium (disk, layer) and examine its quasi-stationary dynamics.
The paper is organized as follows: Section 2 explains the model, Section 3 presents the results, and Section 4 concludes with the discussion.Appendices A and B explain the details of calculations; Appendix C elaborates the accretion disk model.
The derived insights may be useful for analyses of environments ranging from accretion disks (in an approximation of curved spacetime) to planetary atmospheres (in an approximation of flat spacetime): the design of numerical models and image reconstruction methods may be influenced by how the steady-state of the media is conceptualized-as multiple (separate) quasi-localized vortices or as a single (integrated) multi-petal structure.

Model
We use the quasi-2D fluid dynamics model, which describes large-scale vorticial motions of a "thin" rotating disk in a mildly curved spacetime.
Coordinates and Notations: In the model, we use the coordinate system where x ≡ x 1 and y ≡ x 2 are the axes in the horizontal plane and z points vertically up.We employ the following notations: Φ is the potential of the force field; temperature θ ≡ θ 0 + τ s + τ, where θ 0 is the (constant) baseline temperature of the disk, τ s is the spatially inhomogeneous, axially symmetrical, time-independent part of the temperature distribution, and τ = τ(t, x j ) is the temperature-a dynamical quantity related to the vortex structures in the fluid; v i signifies the components of the velocity field v; ω 3 is the z-component of vorticity [∇ × v]; and ϵ ijk is the standard alternating tensor (ϵ 123 = 1, ϵ ijk = 0 when any two indices are equal, ϵ ijk = +1 for any even number of permutations of indices from ϵ 123 , and ϵ ijk = −1 for any odd number of permutations).Note also that ϵ 3jk ∂ j τ s ∂ k Φ ≡ 0 because both τ s and Φ are axially symmetric.Furthermore, we assume that the band of the disk between r 1 and r 2 > r 1 , where the vortex structure of interest is formed, is characterized by an approximately constant angular velocity Ω = (0, 0, Ω).
"Thin" Disk: In the simplest model-in which the vortex structures are realized in a "thin" flat sheet with h ≪ L, V 2 ≪ s 2 ≪ c 2 , Re ≫ 1-the z-component of velocity, w, vanishes and may be dropped in all formulas (see Ref. [10], developing this particular model, or classical Refs.[12][13][14]).Here, h is the characteristic local thickness of the disk, which is slowly dependent on the radial coordinate; L is the characteristic spatial length over which the macroscopic flow characteristics change substantially; s is the so-called "isothermal sound speed"; V is the characteristic fluid velocity; c is the speed of light; dimensionless parameter Re is the Reynolds number (or its analogs) such that Re ≫ 1.
The System of Equations: The set of equations for the z-component of vorticity ω 3 and for variations in temperature τ becomes: where indices j, k = 1, 2. Operator D t ≡ ∂ t + v k ∂ k is the total derivative.The condition of incompressibility div v = 0 does not mean that the density of the medium is constant.This condition simply means that the density evolves according to the equation ), after some simple manipulation, the evolution equation for temperature perturbations takes the form of the third equation in Equation (1).For the equation of state, we will assume, in accordance with the assumption div v → 0, that the density essentially depends only on the temperature and not on the pressure.Thus, we write ρ ≃ ρ 0 (1 − β(T − T 0 )), where subscript 0 denotes the reference values.The coefficient of thermal expansion β = −ρ −1 (∂ρ/∂p) T is presumed to be constant in the layer.Since quantity β(T − T 0 ) is generally small (relative to 1), one may neglect the density variations and hence replace ρ with the constant value ρ 0 in all terms of the dynamic momentum equation except for the "buoyancy" term.We will presume that a collective fluid motion forms in the system, for which the relationship between the flow velocity and temperature variation is given by ( Here, tensor ϵ ij is the antisymmetric unit, second-order tensor, for which ϵ 12 = −ϵ 21 and the diagonal components are zero; and a is some parameter assuring the correct dimension for the relationship between the considered quantities.In this case, vorticity Then, notably, the presence (or absence) of the component with constant Ω 3 does not impact the core set of Equation (1).
As mentioned, the vorticial field is intertwined with the thermal field.Thus, in addition to the set of Equation ( 1), any model must include equations for the temperature field.We set the ansatz for the temperature profile as θ ≡ θ 0 + τ s + τ, where θ 0 is the (constant) baseline temperature of the disk, τ s is the spatially inhomogeneous, axially symmetrical, time-independent part of the temperature distribution, and τ = τ(t, x j ) is the temperature-a dynamical quantity related to the vortex structures in the fluid.We chose the stationary temperature deviation τ s in a form that allows for a self-consistent analytical solution.This choice does not contradict experimental observations (see, for example, Ref. [15]).Therefore, we will model (using notation r = |x| and x = (x 1 , x 2 ), and the Taylor series expansion) the background distribution of the time-independent temperature deviation τ s in the simplest form: where τ 1 is the parameter characterizing the "rapidity" of the increase in the basal temperature of the medium as the distance r from the disk rotation axis grows.(Generally speaking, the question of what shape the temperature profile takes within a black hole's accretion disk is an open one.Astronomical measurements remain challenging, despite significant progress.See, for example, Refs.[15,16].)When the leading contribution to potential Φ comes from the centrifugal effects, we can express Here, the square brackets denote the cross product of two vectors.Param- eter K 2 is determined by the expression K 2 = GM BH /r 3 = c 2 r g /(2r 3 ), where K ∼ r −3/2 is the so-called Kepler parameter (in the Newtonian approximation of gravity), which is interpreted as the angular velocity of a test particle on a circular orbit at distance r: G is the gravitational constant; c is the speed of light; M BH is the black hole "mass".Thus, when angular velocity Ω is presumed to remain constant within the band of r where the thermo-vorticial macro-structure forms, the set of the interlinked nonlinear evolution equations becomes: where, obviously, quantity v k = aϵ ik ∂ i τ is the k-th component of the flow velocity.Equations ( 4) and ( 5) have a transparent physical meaning.Their left sides describe the transport of dynamic (field) quantities: vorticity (∼∆τ) and temperature perturbation (τ).The right sides of these equations describe "sources" that generate the vortices and the temperature perturbations.In other words, the vortices are generated by the source described by the right part of Equation ( 4), which is effective (non-zero) only when there exists an inhomogeneous gravity-like force field Φ ∼ Ω 2 with which temperature perturbation τ interacts via the equation of state (when compressibility β ̸ = 0), while the quantity ∆τwhich characterizes the vortex field in the fluid-is transported by the self-induced flow (the left part of Equation ( 4)).On the other hand, in Equation ( 5), temperature perturbation τ is transported by the self-generated flow (v i ̸ = 0); temperature perturbation τ(t, x) is generated by the source (the right part of Equation ( 5)), which is non-zero only when there exists a spatial-and time-independent temperature gradient (i.e., when τ 1 ̸ = 0).The processes are interlinked because the "source" of one dynamic quantity depends on a complex combination of other dynamic quantities.Equations ( 4) and (5) show that when the stratification of temperature is absent (τ 1 = 0) and there is no disk rotation (Ω = 0, i.e., centrifugal force is zero), then Equations ( 4) and (5) degenerate into the traditional equations for vortex evolution and transport in a two-dimensional ideal fluid.Comparing Equations ( 4) and (5), we conclude that both the first and second equations describe the evolution of the same physical entity.As follows from the definition of hydrodynamical velocity, the dimension Lagrangian Quantity ("Q-Vorticity"): By combining Equations ( 4) and (5), i.e., after multiplying the second equation by βΩ 2 /aτ 1 and subtracting it from the first, we obtain: Here, parameter showing that R is a space scale factor.If the question is posed to examine the evolution of the temperature field, then the nonlinear differential equation Equation ( 6) is solved numerically, specifying the initial distribution of the temperature field, as well as the spatial boundary conditions.
On the other hand, we may define a new quantity q ≡ a(∆τ − R −2 τ), (7) which we will call "q-vorticity" because this quantity is linked to the streamfunction via definitions.This quantity is transported by the flow according to Equation (6), i.e., as a Lagrangian quantity.Thus, the integral of any function of q in the xy-plane is always conserved.
Temperature: Temperature τ is found via the Green function representation.It is presented as follows: Here, the Green function satisfies the following equation: where δ is the Dirac delta-function.The collective "hydrodynamical" velocity of the flow (orthogonal to the temperature gradient) is calculated using Equation (2).
Modeling Steady-State: So, what is the initial vorticity distribution q(x, t) and how can it be modeled?Recall that we are not talking about the evolution of the system from some state of "nothing", but are focusing on the (dynamical) steady-state to which the system eventually evolves.
Two methods are useful for modeling such a steady-state of q(x, t).
For a strongly localized vortex-see, for example, Ref. [17]-quantity q(x, t) can be parameterized by a function in the form F(x − x 0 (t)), i.e., the one with the center at coordinate x = x 0 (t), which satisfies the equation of transport ẋ0k (t) + ϵ ik ∂ i ψ = 0.In other words, in such cases, the center of the vortex moves according to ẋ0 (t) − v[x 0 (t)] = 0. (Here, the symbol "dot" signifies the derivative with respect to time.) Alternatively-as we will consider below-the distribution of q(x, t) can be described in the form of a macro-structure with "petals" (with constant q(x, t) = q 0 inside) whose moving boundaries evolve according to Equation (6).
Extended (i.e., not narrowly localized) stationary vortex structures-the ones rotating with constant angular velocity ω (≡ ω 3 )-are simpler to consider in a rotating coordinate system where the structures appear immovable, i.e., when their rotation direction is coaligned with the z-axis and the derivative with respect to time becomes ∂ t = −ωϵ ik x i ∂ k .(The procedure is laid out, for example, in Ref. [18].)Indeed, when f = f (ρ, ϕ − ωt), then calculations relying on the properties of Jacobeans produce the following: Then, Equation (5) may be rewritten as follows: which is satisfied by the ansatz The explicit expression of function F can be found from the obvious fact that temperature-driven flow perturbations must vanish when temperature perturbations vanish themselves.The suitable expression is as follows: function F(u) = −τ 1 u/ω.Then, temperature fluctuations are expressed as follows: It follows from here that parameter a = −ω/τ 1 .Combining this with Equation (4), where ∂ t = −ωϵ ik x i ∂ k is taken into account, we obtain the second evolution equation-the one which describes the rotation of a stationary vortex structure with angular velocity ω caused by the self-induced field of hydrodynamical velocity: Here, 2 , where the thermal length scale λ θ = (τ 1 β) −1/2 characterizes the background regime of this model.Obviously, an extended (multi-petal) thermo-vortex structure exists when parameter τ 1 , characterizing the "rapidity" of the change in temperature (and stream function) with the distance from the disk rotation axis of τ 1 > 0.
The "Blob" in The Method of Contour Dynamics / Sharp Boundaries Evolution Method (SBEM): The term "thermo-vorticial blob" defines, in the plane (x 1 , x 2 ), a domain D filled with constant vorticity q 0 (with the dimensions [q 0 ] = [time] −1 ) and bounded by a closed contour ∂D.
The contour can be described in the parametric form, where parameter s is the contour's arc length (starting from some initial point on the contour), such that 0 < s < L, where L is the entire length of the contour.The shape of the blob boundary is given by the following equations: where the quantities ν 1 , ν 2 , and ϕ are functions of s.Next, we introduce dimensionless variables and parameters: σ = L/l, s → Ls, α = (ω/q 0 )(l/R) 3 , η = σ/4K(m).Here, K(m) is the complete elliptic integral of the first kind.Omitting the details of the analysis, we will only point out that the functions are linked with the curvature of the contour κ via the following non-linear differential equation: Here, prime signifies a derivative with respect to dimensionless argument σs.The angular coordinate ϕ(s) of the tracing point on the contour (Figure 2) is calculated from the following obvious expression: which follows from the definition of curvature κ.
To bring the equation into a canonical form without the free term, we can make the following manipulations: make a "shift" in κ, "flip" the result "upside down", make a change in the magnitude, make one more "shift" in the result, and finally, make a change in the magnitude again.In other words, the solution has the following form: where z(λs) is the Jacobi elliptic cosinus and λ = n4K(m) for n-petal configuration: n = 2, 3, . . . .Parameter λ provides the periodicity of the solution (z(0) = z(λ)), i.e., the closing of the contour.Therefore, we change the argument, transitioning from σs to θ = λs.
The initial conditions are introduced as κ(0) = κ 0 and κ ,s (0) = 0.The solution contains six free parameters: α, η = σ/4K(m), A, B, m, ϵ.However, one parameter remains free; this fixes the values of all other parameters.We select ϵ as the controlling parameter.Parameter ϵ determines the magnitude of the initial perturbation because of κ(0) = A + B/(1 − ϵ).Obviously, for the unperturbed state, when ϵ = 0, we have Omitting the details of calculations, we found Parameter µ(n, ϵ) is fixed by the condition that the solution of equation has to satisfy the boundary conditions ϕ(0) = π/2 and ϕ(1) = 5π/2.Here, cn(ξ, m) is the Jacobi elliptic cosinus, K(m) is the complete elliptic integral of the first kind, 4K(0) = 2π, and n is the number of petals in the vortex structure.In Equation (19), parameters A, B, and others, are taken from Equation (20); dimensionless parameter s changes from s = 0 to s = 1.
The calculation, which we do not present in this section, provides the leading term in the expansion of µ = µ(n, ϵ) in a series with respect to ϵ (which does not vanish when ϵ → 0): the expression µ → n 2 − 1.

Results
To analyze the thermo-vorticial "blobs", we employed the so-called Sharp Boundaries Evolution Method (SBEM).Within the blob, q(x, t) = q 0 is constant; outside the blob, q(x, t) = 0.This method originated in the 1960s-1970s and was further advanced in subsequent years.Despite the seemingly radical simplification, the method produces a reasonably good description and estimates for large-scale vortex flows/structures ( [9,10]).Indeed, for the motion of an incompressible fluid in the xy-plane, the velocity vector v lies in the same plane, and the vorticity vector has only one non-zero component, q z .For large-scale vortex structures, it is intuitively apparent that there exist zones (which are sometimes rather spacious) where the vorticity magnitude reaches its maxima ("hilltops") when the sign of vorticity is positive and its minima ("chasms") when the sign of vorticity is negative.There are also zones where vorticity is zero or close to zero ("flats"), transitional zones between hilltops and flats ("upper slopes"), and transitional zones between flats and chasms ("lower slopes").Since typically the fluid velocity and thermodynamical quantities are analyzed in terms of the integral characteristics derived from vorticity, when the slopes are steep (i.e., the Reynolds number Re or its equivalent is enormous), these transitional zones may be treated as sharp (precipitous, vertical) moving boundaries which evolve in accordance with the laws of hydrodynamic flows.In this paper-to not complicate the calculations-we consider only the case with the hilltops, upper slopes, and flats (not including the chasms and lower slopes).
One Coherent Multi-Petal Structure: We find that, because of the intertwined nature of the thermal and vorticial fields in the fluid medium, what often looks like a group of multiple individual "circular" vortices (as depicted in Figure 1, for example) is in fact one coherent macro-structure with a complex multi-petal shape.In other words, the localized "circular" vortices are not independent; they influence each other, redistributing conservable quantities-such as vorticity q(x, t) and energy, which are linked to temperature τ(x, t)locally, within and between each other, while preserving the integrals over the entire multi-petal system.
Figure 3 illustrates the contours of thermo-vortices with three and five petals.Describing the shape of such a multi-petal structure mathematically is not a trivial exercise.We demonstrate how we characterize the shape using several key parameters: the number n of petals (counting the "fingers"), the level of non-linearity ϵ (describing whether the "fingers" are long or short; measuring the gap between the peak and trough); and parameter µ(n, ϵ) (describing whether the "fingers" are thick or thin, and sharply pointed or not).The latter is achieved by employing (instead of "smooth" sinus/cosinus functions) the Jacobi elliptic sinus/cosinus, which permits the inclusion of higher spectral harmonics into the description of the shape's "fullness".
In Figure 3, both shapes have the same peak/trough nonlinearity, i.e., they possess the same parameter ϵ = −0.28.But because they have a different number of petals n, their values of µ(n, ϵ) are different.Figure 4 illustrates how µ(n = 3, ϵ) depends on ϵ and how parameter ϵ (and therefore µ(n, ϵ)) impacts the contour shape.Once the shape of the contour of the macro-structure can be described, the properties of the system can be analyzed.

Number of Petals and Other Key Parameters:
Generally speaking, the number of petals n evolves as the system evolves because of the complex linkage between the parameters in the system.Physically, multi-petal vortex structures arise when thermal energy is introduced into the system, leading to the formation of a hot (τ 0 ̸ = 0) domain of finite (0 < l < ∞) size with non-zero vorticity (Q ̸ = 0) inside.Mathematically, when three quantities Q, l, and τ 0 are specified as the initial conditions, they "fix" three geometrical parameters: n, ϵ, and µ.Here, the conserved quantity is as follows: Recall that here, l is the characteristic size; σ ≡ L/l, where L is the contour length; parameter s is the contour's arc length.
Quantities l and τ 0 may be measured in observations.Thus, their links to other parameters in the system are informative: All coefficients are functions of ϵ, n, and µ.The dimensions of the quantities assure the correct dimensions of the written expressions for l and τ 0 .
Temperature Distribution: Determination of the contour is the first step of the analysis.Once the steady-state contour for the thermo-vorticial blob is established, the temperature distribution τ = τ(t, x j ) inside the blob may be determined.The expression for temperature distribution in terms of dimensionless variables is as follows: Figure 6 shows-for the same study-the vorticity distribution q(x, y) and the temperature profile (in the xz-plane for y = 0; obviously, temperature is distributed similarly within the planes slicing the other petals).The right panel also shows how the temperature profile varies depending on the geometrical parameter β = l/R.The greater β is, the "tighter" the thermal structure within the petal-the petal peak is higher, the sides are narrower, and also, near the center of macro-rotation, the central thermo-vortex appears.

Discussion and Conclusions
The fact that rotating cosmic fluids contain vortices has been understood for a while.Vortices have been discussed in the context of planetary atmospheres (at least since the discovery of the Jupiter's Great Red Spot), in the context of vortices and magnetic flux tubes in the accretion disks of black holes [23], galactic center hot spots [24], and so on.When the observed vortices are multiple, coherent, and long-lasting, they are also often clustered into orderly quasi-symmetric arrangements.These multiple vortices are typically interpreted as possibly interacting, but nonetheless as stand-alone items.
In our paper, we demonstrated that an additional physical interpretation of what these vortices represent may exist.What one sees in astronomical images as distinct spots may be not multiple individual quasi-circular-shaped vortices, which possibly interact and drift to form a quasi-symmetrical quasi-stationary macro-arrangement, but instead these distinct spots (thermo-vortices) may be the observable parts of one structure whose shape is a complex contour with multiple "petals".Either alternative is theoretically possible, but the latter one is typically not considered.
To bring attention to this multi-petal interpretation, we presented and analyzed a model that can be illustrative and adaptable to various environments.Specifically, in order for our model to be of interest for studies of the accretion disks of black holes, we used the metric that works for the mildly curved spacetime (naturally, this model works in flat spacetime as is the case for planets).We had to limit our considerations to the model of a rotating "thin" layer because otherwise the mathematical calculations would become unwieldy.Our focus was on the description of the steady-state, large-scale, thermo-vorticial motions, not on the evolution from the state of "nothing".We demonstrated that the steady-state may indeed look like a multi-petal structure and that the shape of individual petals may be highly nonlinear (the petals may look like "fat fingers", or be elongated or pointy).Various parameters (relative dimensions of the system, thermal gradient, initial vorticity, etc.) and their linkages determine the outcome.We also demonstrated that the petals indeed have an elevated temperature and would appear as bright spots if observed from afar.
The physical explanation for why a system would evolve towards a steady-state with N petals rather than M petals is intuitively simple.The contour of the macro-structure bounds the inner domain of higher vorticity.Naturally, the contour is closed.Therefore, if subjected to some perturbation, the contour's steady-state turns into a "standing wave" (like a standing wave on the surface of a swimming pool), with the number of crests (petals) defined by the properties of the system.The shape of the crests, however, is not sinus-wave-like but is complex (finger-like) because the system is highly nonlinear.Obviously, no steady-state is eternal, but we did not study the problem of the transition of quasi-steady-states.
Several notable multi-vortex arrangements, which are likely to be multi-petal structures, can be mentioned.As seen in Figure 1 (second row; panels C and D ), in the span of one year, the Antarctic ozone hole evolved from a quasi-circular into a two-petal structure.Based on the common knowledge of what the ozone hole represented, we can conclude with certainty that the two vortices observed in 2002 were not separate entities; they were indeed two petals of one shape-changing macro-structure.The arrangement of vortices on Jupiter's South Pole (Figure 1B) is also a candidate multi-petal macro-structure.Specialists in geophysical fluid dynamics may spot an analogy between the study of this paper and the Rossby waves (see, e.g., Ref. [25]): the interface on the thermocline η (or, equivalently, the streamfunction) corresponds to the temperature τ (which is proportional to the streamfunction due to the presumed Equation ( 2)), the pressure gradient responsible for the geostrophic flow corresponds to the gradient of the centrifugal potential Φ, and the potential vorticity corresponds to the q-vorticity.Not surprisingly, the global nature of planetary solitary waves corresponds to multi-petal structures and explains images like those in Figure 1B-D.One difference, however, is that the Rossby waves occur at the sharp boundaries between regions with different temperatures; in our case, the regions are of different vorticities.
For a strongly relativistic system-such as an accretion disk in the environment of a rapidly rotating supermassive object or a supermassive black hole, see Refs.[45][46][47]-the conditions (in covariant description) of adiabaticity and incompressibility remain the same, while the equation of fluid motion (in spacetime with a metric defined by the square of the four-interval ds 2 = m ik dq i dq k in the equatorial plane of a rotating black hole when the coordinate q 2 ≡ θ = π/2) has the following form in the leading approximation (these details of our analysis are not listed in this paper): Here, D is the operator of the total derivative with respect to the "synhronized time"; the enthalpy w is defined as w/c 2 = γρ + (1/c 2 )(ρe 1 + p), where γ is the Lorentz factor, p is the pressure, and e 1 is the non-relativistic part of the internal energy); c is the light speed; v β is the covariant β-component of the 3D velocity.Naturally, the basic elements of the model momentum in Equation (A84) are the pressure gradient (the first term), the "Coriolis force" (the second term), and the gravity and centrifugal forces, which are packed in the third term of the momentum equation.The key difference compared to the classic hydrodynamic treatment is the replacement of the classic gravity/centrifugal acceleration −∇Φ + Ω 2 r with −∇ ln √ m 00 , and the replacement of the classic averaged angular velocity Ω of the accretion disk rotation with the relativistic parameter Ω = (c/2)γu 0 curl(−g 00 g).
If the thermo-vortices are at some "distance" from the event horizon and conditions x = r/r g ≫ 1, Ωx < 1 are satisfied, then, leaving the leading terms in the expansion Here, the terms of higher orders of smallness are omitted.The first term describes the centrifugal effect, the second describes the contribution of the force of attraction in the Newtonian approximation of gravity, the third term describes the Schwarzschild model of a non-rotating black hole, the last takes into account the rotation of both the black hole and the disk.For x ≫ 1, Ωx < 1, but Ωx 3/2 > 1/2, the principal contribution provides the first term ∼ Ω 2 , i.e., the centrifugal effect, independent of the exact structure of the spacetime near the rotating, or not rotating, black hole.
To conclude, the illustrative model that we presented provides insights that are useful for a broad range of physical applications, ranging from planetary atmospheres to black hole accretion disks.Whether the observed vortices in cosmic fluids are conceptualized as stand-alone objects or inextricable parts of one macro-structure may meaningfully influence the numerical models employed in consideration of the phenomena and image reconstruction methods.
In the complex plane, the Θ(z)-function can be written in the form of the contour integral (see Refs. [19,26], for details), where the contour is traced in the counterclockwise direction: Equation (A7) is, obviously, a corollary of the Cauchy theorem.Equation (A7) can be rewritten in the following frm: With the help of two useful formulae (see Ref. [19]) the derivatives of Θ-function with respect to z * and z may be expressed in form of the contour integral: The introduction of the so-called ψ-function, defined as ψ ≡ aτ, helps simplify derivations.Indeed, vorticity q ∼ q 0 ∼ curlv is a local characteristic, which determines the magnitude of the vortex motion in the structure, whereas temperature τ depends on the characteristics of the environment (ω, τ 1 , . . .).The ψ-function is a local solution of the corresponding Green equation Here, dx = dx 1 dx 2 is a surface element of the xy-plane; function K 0 (ξ) denotes a modified Bessel function of the 0-th order.Since K 0 (ξ) is the Green function (see, for example, Ref. [48]), we can rewrite Equation (A11) as follows: The Laplace operator (with respect to the running coordinates with index 1) is obviously From here, we can obtain the following: Integrating (with respect to z * 1 ) Equation (A14) by parts and using known transformations for the Bessel functions, we find that Next, we use Equation (A10) and change the integrating order.We find the following: Since functions ψ and τ are connected by Equation ( 12), via τ = −(τ 1 /ω)ψ, we obtain, for the temperature distribution: Transitioning to dimensionless variables (s → Ls, z → lz, ẑ → l ẑ, l/R → β, L/l → σ)-note that we preserve the same symbols for efficiency-we can rewrite Equation (A17) in the following form: where Parameter τ 0 = τ 1 q 0 R 2 /ω is the parameter which determines the characteristic temperature level of the vorticial blob.The integral Equation (A19) is a function of coordinates x = (x, y) which depend on the dimensionless geometrical parameters β = l/R (the ratio of the blob size l to the "thermal" scale R) and σ = L/l (the ratio of the contour length L to the characteristic blob size l).Recall that τ 1 , with dimension [τ 1 ] = [temperature]/[length] 2 , is the parameter characterizing the "rapidity" of the increase in basal temperature τ s with distance from the rotation axis.

Appendix A.2. Value of ψ-Function for Blob Contour
To find the expression for the contour, we first transition from x → x(s).Then, from Equation (A16), we obtain the following for ψ = ψ( ẑ(s): where the variables are still dimensional.
We assume that parameter R is small in comparison with the characteristic scale of the blob.Then, the main contribution to the integral comes from a small region in the neighborhood of the tracing point ẑ(s).The contour may possess some (small) curvature (meaning | ẑ,ss | ≪ R −1 ).In the vicinity of point s (with small ν = u − s and ẑ,ss = i κ ẑ,s where κ(s) = ϕ ,s ) we can write the following: Here, all terms of order ≥ ν 3 are omitted.Equation (A20) is written (for a small curvature of ∂D) as a sum of two integrals: the first one, due to symmetrical limits, yields zero; the second-not zero: Here, both limits of integration (±∞) stretch to infinity (keep in mind the exponential nature of the decrease in the Bessel function with the increasing value of its argument).The value of the integral is πR 2 .Accordingly, we obtain the final result: This expression establishes the link between the values of ψ-function (temperature) on the contour and its local curvature κ.
The unperturbed state of Equation (A31) is the one satisfied by κ(s It may appear that one can simply solve the system in Equation (A31) numerically.Indeed, the solutions depend on one variable s and are fixed by two dimensionless parameters: parameter α = (ω/q 0 )(l/R) 3 and phase parameter σ.The local curvature of contour κ 0 at the initial point is the governing parameter.Since we are dealing with differential equations with first-order derivatives and are given initial conditions, obtaining a numerical solution is not difficult.However, we are only interested in the solution that will "close" the contour, i.e., the tracing point must return to the starting point.This means that one more requirement must be added to the system.

Appendix B. Equations Permitting Analytical Solution
Aiming to find the equation for the blob contour, we first rewrite Equation (A31) in the following form: which may be further rewritten as follows: Equation (A33) resembles the non-linear Duffing's equation: which is known (see Ref. [20][21][22], for example) to be solved analytically in terms of Jacobi elliptic functions y = cn(θ, m) with the period 4K(m).Here, K(m) is the complete elliptic integral of the first kind.The non-zero free term 8α complicates the straightforward solving of Equation (A33).
Obviously, ϕ(1) = ϕ(0).Thus, the parameters σ and α take the following values: where n = 2, 3 . . .The phase ϕ(s) in the considered approximation is given by the following expression: After this step in calculations, we find the following quantities: and derive the contour of the vortex blob in the parametric form: where ρ = µ 2 + ν 2 and ψ = arctan(ν/µ).Thus, the angle Φ = ϕ + ψ determines the polar angle.In order for there to be no overlap, it is necessary that the derivative dΦ/ds is greater than zero.Naturally, the linear approximation can only indicate the tendency.One harmonic is not enough to reveal the details of a deep contour deformation.For a more profound understanding of the problem, it is necessary to use an exact analytical solution.
On the other hand, the simple analysis that was discussed makes it possible to grasp the non-trivial fact that the parameters α and σ are not arbitrary: they are "quantized" in a certain sense.Their values are determined by integers: n = 2, 3. . . . .That is, if some localized area is filled with a quasi-homogeneous vortex with vorticity density q 0 , then an n-petal vortex structure can form in this area only if the parameters of this structure (size, intensity, background temperature gradient, angular velocity of disk rotation, etc.) satisfy the quantification conditions of Equation (A42).If these conditions are not met, no n-petal structure is formed at this location.

Appendix B.2. Strongly Nonlinear Configuration: Exact Analytical Solution
With the initial condition κ(0) = κ 0 and κ ,s (0) = 0, the basal equation is 1 If there were no linear term in curvature κ, the periodic solution of Equation (A47) would be expressed in terms of Jacobi elliptic functions cn(λs, m), or sn(λs, m), where dimensionless parameter λ provides the periodicity of the solution: z(0) = z(λ).
Let us bring the equation into canonical form.To do this, we will make the following manipulations: make a "shift" in κ, "flip" the result "upside down", make a change in

Appendix C. Model of Accretion Disk Near Black Hole
Accretion Disk: Generally speaking, different types of accretion disks are discussed in the literature: thick disks ("doughnuts"), thin disks (Shakura-Sunyaev model), slim disks, and advection-dominated accretion flows.The standard thin accretion disk model is applied for an accretion rate that is substantially smaller than the Eddington limit.This model was one of the first studies of accretion disks in the vicinity of black holes, formulated in the early 1970s by Shakura and Sunyaev.The authors articulated the basic assumptions and equations to describe classic thin disks in the Newtonian approximation of gravity.Relativistic disks were first described by Bardeen, Press, and Teukolsky, and by Novikov and Thorne.Thorne introduced the "slim disk" model, whose equations include a number of terms neglected in the thin disk model equations.For more detail, see, for example, Refs.[46,47]; these reviews cover the main aspects of the black hole accretion disk theory.
In the model presented below, the following presumptions are made about the accretion disk: (a) It is thin, i.e., its characteristic vertical thickness is much smaller than its characteristic horizontal size; (b) It is located in the equatorial plane, implying that the averaged θ-component of 4-velocity becomes suppressed; (c) It rotates with the rotation axis perpendicular to the disk plane; (d) It is heated non-homogeneously, with higher temperatures at the outer layers (due to the processes of the tidal disruption of captured matter and strong dissipation in the shear flows).
Let us assume that, due to some mechanism, a locally "twisted" region of heated matter spontaneously formed in the rotating disk.Then, in the field of the centrifugal force, within the outer layers, "hot bubbles" form (i.e., localized plasma clusters whose temperature exceeds the "average" for the disk and whose density is lower), which move towards the axis of the disk rotation.However, when these hot bubbles are also vorticial in nature (for example, due to the overall rotation of the disk), then each vortex bubble (via induction of the velocity field) "forces" other vortex bubbles to rotate around itself.Hence, the bubbles move "sideways", rather than directly towards the rotation axis of the disk.As a result, since all vortices are affected by the cumulative velocity field, they gradually self-organize into zones of stability-a symmetric thermo-vorticial macro-structure which rotates as a whole around the mutual center of symmetry.The dynamics and longevity of this structure are linked to the thermal and vorticial properties of the system and its elements.Visually, the protruding components of this structure look like bright "hot petals".For large-scale quasi-2D flows in "thin" fluid layers, with very high Reynolds numbers, the so-called "smoothing" effect is observed when the level of vorticity acquires an approximately constant value throughout the region of the vortex (see classical Ref. [26], and subsequent developments in Refs.[9,10,18,27], and references therein).
The following assumptions are also made in the model: the mass of the accretion disk (∼ 10 −5 ÷ 10 −4 M BH ) is negligible compared to the black hole "mass" M BH , the radiative cooling does not strongly affect the dynamics of the fluid motion in the disk, the electrons and ions are very weakly coupled by Coulomb interaction, and therefore, ion and electron plasma components may have different temperatures, even with T e much greater than T i (see Ref. [28]).Then, it is the electron plasma component which contributes the most to the equation of state of the accretion disk matter.Due to the large difference in the masses of electrons and protons, electrons are highly mobile and provide the quasineutrality of the plasma.Also, due to the high conductivity of the plasma, its own magnetic field can be considered "frozen-in".A range of possibilities may exist.In hot flows like those around Sgr A* or M87, electron temperature is thought to be typically lower than the ion temperature due to radiative losses to synchrotron, inverse Compton, and bremsstrahlung processes.
Model Assumptions Regarding the Thickness Of the Disk: We consider the motion of the medium relatively far from the event horizon, i.e., for hydrodynamic structures, it may be assumed that r ≫ r g , where r g is the Schwarzschild radius.In this paper, the approximation r > 3r g may be considered "far away".When we are not dealing with phenomena occurring near the event horizon (which is r ∼ r g ), then, for r > 3r g , we can restrict ourselves to the roughest approximation for the spacetime metric near a rotating (or not ) black hole.
It is assumed that the disk is geometrically thin-h ≪ r-where h is the characteristic local thickness of the disk.
When r > r g (more specifically, r > 3r g ), then the vertical component of gravity acceleration can be estimated as This is balanced by the pressure gradient.Here, G is the gravitational constant, c is the speed of light, K 2 = GM BH /r 3 = c 2 r g /(2r 3 ), and K ∼ r −3/2 is the radius dependent on the so-called Kepler parameter: the angular velocity of a test particle with a circular orbit at a distance r from M BH in the Newtonian approximation of gravity.
Gas pressure along the direction perpendicular to the disk plane (x, y) is determined by the hydrostatic equilibrium, dP = ρg z dz.Generally speaking, electromagnetic radiation and the induced magnetic field may contribute to the pressure.In the simplest case, the dominant contributor is the gas thermal pressure; then the vertical temperature distribution is isothermal, which is an acceptable approximation when the disk is optically thick and externally heated.The equation of state of gas/plasma is then P = s 2 ρ, where s is the "isothermal sound speed" (which is not a function of the transversal coordinate z).
If we further assume that z ≪ r, the equation of hydrostatic equilibrium becomes s 2 dρ = −s 2 h −2 ρzdz, the solution of which is ρ(z) = ρ 0 exp(−z 2 /2h 2 ).Here, the transversal space scale h is introduced: h = s r 3/2 /(GM BH ) 1/2 .Parameter ρ 0 represents the density at the disk mid-plane (at z = 0), and parameter h is the characteristic local thickness of the accretion disk.
In order for the accretion disk to be considered thin, it is necessary that h/r ≃ sr 1/2 /(GM BH ) 1/2 ∼ (s/c)(r/r g ) 1/2 ≪ 1.On the other hand, for great distances from the black hole (r > r g ), specific relativistic effects may be neglected or parametrized.Thus, we obtain the following natural bounds for the validity of our consideration: r g < r < r g (c/s) 2 .In other words, the presented model may be valid for the consideration of spots in accretion disks (resembling those in Figure 1A, for example) if the spots are located within these limits.
Note that we work with the field, not with individual particles (their trajectories)."Thin" Disk Temperature Distribution T(r): The problem of how physical parameters, such as temperature, are distributed within a geometrically "thin" accretion disk may be examined, of course, within the framework of the hydrodynamical approximation, introducing, among others, the concepts of turbulent flow, Reynolds tensor, turbulent viscosity, and the closure of equations.Instead, since we are aiming to obtain only a model-specific estimate, we will first (1) neglect the dependence of physical parameters on the perpendicular to the disk-plane z-coordinate (i.e., integrate/average vertically) and (2) consider the disk without radial thermal advection (i.e., the transfer of heat with matter along the radius) and without the loss of mass from the lateral surfaces of the disk.In such a disk, the angular velocity of rotation of the disk matter at each location along the radius r is approximately equal to the angular velocity of rotation of a free test particle.In other words, v r ≪ v ϕ .(Recall that key parameters that determine the structure of a geometrically thin disk are the mass of the gravitating center M BH , the internal radius of the accretion disk R in , and the accretion rate ṁ.) The next assumption is that the gas of the accretion disk spirals inward and this spiraling is very gradual, i.e., the orbits of the gas particles are almost circular.The orbital speed is Keplerian and is estimated as v ϕ = √ GM BH /r = Kr.Quantity K = GM BH /r 3 is the so-called Keplerian parameter, which is radially dependent, determining the period (T K = 2π/K) of revolution of the test particle located at distance r from the center of attraction M BH .
The categorization of a (vertically isothermal) disk with characteristic thickness h as "geometrically thin" means that h/r ≪ 1, where h 2 = 2s 2 /K 2 , s ∼ √ T is the isothermal sound speed for the gas, T is the absolute temperature, and the gas pressure is p = ρs 2 .Then, it follows from the static equilibrium for the motion in z-direction that −dp/dz = −s 2 dρ/dz = ρ(GM BH /r 2 )(z/r) = ρK 2 z.This yields ρ = ρ 0 exp(−z 2 /h 2 ), where ρ 0 is the density of the disk mid-plane.
The specific angular momentum of the test particle is r v ϕ = √ GM BH r.This means that, to flow inward (so that the radial velocity dr/dt < 0), the gas must lose its angular momentum by redistributing the angular momentum within the disk.The "inner" gas transfers its angular momentum to the "outer" gas; the gas matter flows inward.The loss of the angular momentum by the entire system (via the outward wind, for example) results in the inward flow of the remaining disk matter (with increased viscosity, the process becomes diffusive).
For the unit mass located at the radial distance r, the acting potential U = −GM BH /r produces the (inward-directed) force of attraction f = −dU/dr = −GM/r 2 .When the mass dm is transported inward over the distance dr, its potential energy E p is changed by dE p = (GM/r 2 )dm dr.The potential energy of the particle dm when moved from infinity to the location r decreases from zero to −G dm M BH /r.One half of this energy is converted into kinetic energy; the other half is radiated away (this is the consequence of the so-called "virial theorem").In fact, for any gravitationally bound system the time-averaged potential gravitational energy E p and the time-averaged kinetic energy E kin of the motion of the particles of the system satisfy the following condition: E p + 2E kin = 0. Thus, variations in the quantities are linked via the simple relationship: δE kin = −(1/2)δE p .Since the total energy of the disk is E = E kin + E p = (1/2)E p = −E kin , then δE = −δE kin .In other words, any addition of energy to the disk reduces the kinetic energy of its particles (the disk components) and, conversely, energy radiation leads to an increase in E kin , i.e., the temperature of the disk.Since one half of the variation in the gravitational energy E p goes to the kinetic energy of the gas, then the other half is radiated.The luminosity is thus L = (GM ṁ/2r 2 )dr.Divided by the radiating area, 2 × 2πr × dr (the first factor 2 appears because the disk has two sides), this expression produces the luminosity per unit area.
In the approximation where the disk radiates as a black body, its radiation power may be characterized by the effective temperature T (from σT 4 , where σ is the Stefan-Boltzmann constant).Equating the quantity L to the rate of energy loss via black-body radiation, we obtain GM ṁ/8πr 3 = σT 4 .When ṁ is presumed to be independent of the coordinate r, this expression leads to the radial temperature distribution: T = (GM ṁ/8πσr 3 ) 1/4 ∼ r −3/4 .This estimate was obtained for disks without radial advection (the transfer of heat with matter along the radius), without the loss of mass from the surfaces of the disk, without consideration of the boundary conditions at the inner edge of the disk, etc.For an extended source of mass (such that, within the disk, where R in < r < R, the rate ṁ = ṁ0 (r/R) s and parameter s > 0)-i.e., when the inward mass transfer is balanced by the mass production at the disk periphery (r ∼ R) and the mass is "devoured" by the black hole at the eventhorizon (r ∼ R in )-the radial temperature dependence is T ∼ r (s−3)/4 .This means that, depending on the value of the parameter s, the temperature curve may take various shapes.For s = 3, T is radius-independent (the disk is heated uniformly); no heat flow occurs along the r-direction.For s > 3, the disk is hotter at its periphery, where the tidal destruction of captured objects takes place.
Model of Spacetime Metrics: Obviously, equations of fluid motion in the vicinity of a black hole must be written using the concept of relativistic dynamics.In our model (see details in Refs.[29,30]), the spacetime metric is fully characterized by the black hole's mass parameter and "spin" (for an extensive discussion, see Refs.[31][32][33][34][35][36], and the bibliographies therein).Specifically, we use the Kerr metric-an exact, singular, stationary, and axially symmetric solution of the Hilbert-Einstein equations of the gravitational "field" in a vacuum.
The notations, here and presented below, are as follows: Latin indices and suffixes take the values 0, 1, 2, 3; parameter x 0 = ct, c is always the speed of light; the Greek letters take the values 1, 2, 3 and correspond to the spatial coordinates.The Galilean metric (special relativity) is characterized by a metric tensor g ik = (1, −δ αβ ).For the three-dimensional vector below, in Cartesian coordinates, there is no need to distinguish contraand covariant components.
Using the Boyer-Lindquist four-coordinates q i = (t, r, θ, ϕ)-and it is well-known that, besides the Boyer-Lindquist coordinate representation, other representations of spacetime locations exist-the square of the interval is written as ds 2 = g ik (r/r g , θ)dq i dq k , where notations are standard and the components of g ik depend only on the dimensionless combination r g /r and θ.Here, as accepted, r g = 2GM BH /c 2 is the Schwarzschild radius, c is the speed of light, G is the gravitational constant, and M BH is the "mass" of the black hole.The off-diagonal term g 03 in the metric tensor is proportional to the rate of the black hole's own rotation and to 1/r.We use the metric signature diag(+ − −−) (see, for example, Ref. [30], and references therein).To satisfy the principle of causality for moving material objects, ds 2 > 0. The four timespace coordinates q i = (t, r, θ, ϕ) provide the location of a world event from the viewpoint of a remote observer.The meaning of space coordinates r, θ, ϕ is clear once transitioned to the limit r ≫ r g , r ≫ ωr 2 g /c.When the square of the interval becomes ds 2 → c 2 dt 2 − dr 2 − r 2 (dθ 2 + sin 2 θ dϕ 2 ), i.e., at infinity, parameters r, θ, ϕ may be interpreted as the standard spherical coordinates in flat spacetime.For parameter r, strictly speaking, note that it is not the "distance", in the usual sense, from the center of the black hole.This is because, for any material object, in the spacetime defined by equation ds 2 = g ik dq i dq k , no central point r = 0 exists in the sense of a world event on a valid world-line.
Relativistic Flows Of Perfect Fluids: We take into account the effects of special and general relativity to obtain the relativistic Euler equations.To avoid misunderstandings that may arise due to the underdefinition of some concepts (for example, the definition of the signature of the metric tensor), we include many details in this section, even though an advanced reader will certainly be aware of them.
As is known (see, for example, Ref. [31]), the contravariant energy-momentum tensor of a perfect relativistic fluid is written as T ik = (e + p)u i u k − pg ik .Here, e is the internal energy of the fluid; p is the pressure.The quantity w = e + p is the heat function (enthalpy); g ik is the contravariant metric tensor.The quantity u i is the contravariant 4-velocity of the fluid flow; u i is its covariant 4-velocity.The 4-velocity vectors of the flow are normalized by the condition u i u i = 1.The covariant metric tensor (included in the definition of the interval ds via ds 2 = g ik dx i dx k ) is the tensor g ik reciprocal to the tensor g ik ; that is, g il g lk = δ i k .The usual rule of summation is always used.The metric signature is chosen as g 00 > 0. The procedure of the raising and lowering of these indices follows the standard rule.Therefore, the 4-vector components' links are as follows: u i = g ik u k and u i = g ik u k .
The relativistic internal energy e includes the rest energy of particles nmc 2 , where m is the rest mass of one particle and n is the proper number density of particles (i.e., 1/n is the volume per particle).The heat function w, normalized per mc 2 , is written as w/mc 2 → w = n + w 1 /mc 2 .Here, w 1 captures the non-relativistic part of the heat function.Normalized pressure p/mc 2 → p has the same dimension as normalized w: The explanations of why tensor T ik takes such a form can be found in Ref. [12], §133.The mixed tensor T k i is thus A more complex model of the stress-energy tensor of a viscous relativistic fluid with an energy flux may be found in Ref. [33], §22.3, or Ref. [12], §136.A method for building analytical models of relativistic accretion disks may also be found in Ref. [46], which also contains an extensive bibliography on the topic.
To better explain the idea of how to construct the equations of fluid motion, as the first step, we consider the flow in the flat spacetime.Then, if the Cartesian coordinates are used, g 00 = 1, g αβ = −δ αβ .The equations of fluid motion and the condition for the conservation of the proper number density of particles n are contained in the following equations: Reminder: this system becomes closed through the inclusion of the equation of state (for the ideal gas, Fermi gas, plasma, etc.).This is the most subtle part of the problem: one cannot just use some equation of state to "see what happens", one must consider the real astrophysical situation in the proper spacetime for the problem.The magnetic field can also be included in the consideration.Then, tensor T k i would contain an additional magnetic term.If the plasma is highly conductive, then the magnetic field may be considered "frozenin" within the fluid.Then, the additional Lorentz force appears in the dynamic equations and an additional equation for the evolution of the "frozen-in" magnetic induction H appears in the system of equations.Magneto-hydrodynamics is strictly applicable in configurations when the mean free path and mean free time of electrons and protons are much smaller than the characteristic scales and time intervals of the macroscopic motions in question.However, there exist situations when, even for systems with long free paths of current-carriers, the equations obtained from the kinetic theory are formally identical to the MHD equations.Such a situation is observed, for example, in a non-equilibrium plasma when the electron temperature considerably exceeds the ion temperature.
It is physically clear that when a conducting fluid moves in a magnetic field, electric fields are induced in it.Thus, electric currents begin to flow.The magnetic field exerts forces on these currents which, in principle, may considerably modify the flow characteristics.Conversely, the currents themself perturb and even strongly modify the magnetic field.Consequently, a complex interaction between the magnetic and fluid dynamics phenomena takes place.Thus, the fluid flow must be studied by combining the fluid dynamics equations with those of electromagnetic field equations-this is the MHD formulation.Such an approach covers a wide range of physical objects, from liquid metals in a magnetic field to cosmic plasmas.In most configurations, MHD processes are extremely complex.
The description radically simplifies if (a) all dissipative processes are neglected, i.e., no account is taken of thermal conduction and viscosity and the electrical conductivity is considered unbounded such that, in the case of a perfectly conducting fluid, the electrical field is completely screened and the magnetic field is "frozen-in" within the fluid; (b) the fluid is incompressible; (c) in the equations of conservation, the following terms are added: The MHD formulation implies that the displacement current is neglected in the Maxwell equations, i.e., c −1 |∂ t E| ≪ |curl H|; together with the "frozen-in" condition, this leads to the following condition: vl/c 2 τ ≪ 1.Here, c is the light speed, l and τ are, respectively, the characteristic space scale of the considered flow structure and of the changes in its its evolution.From the Maxwell equation for the evolution of the magnetic field H, which, in highly conductive plasma, becomes ∂ t H = curl[v, H], it follows that τ −1 ∼ vl −1 .From these two inequalities, we find that v 2 ≪ c 2 , i.e., the flow must be nonrelativistic.Because ∂ t v ∼ −(4πρ) −1 [H, curl H], i.e., ρv/τ ∼ H 2 /l, and taking into account vl/c 2 τ ≪ 1, we find that it must be that H 2 ≪ ρc 2 .Thus, if the inequality H 2 ≪ ρc 2 s ≪ ρc 2 (where c s ∼ √ T is the speed of sound) is satisfied, the hot plasma may be considered non-relativistic and the effect of the magnetic field can be taken into account as a small perturbation and can be neglected in the leading approximation.
With respect to the modeling of black hole accretion disks in general, one caveat is critically important: the present understanding of the physical conditions within the disks (for example, the equation of state of the matter) is highly uncertain; there are no observational measurements that confirm any of the existing models of the equation of state.Unfortunately, as all studies note, the sensitivity of numerical models to the uncertainties in the parameters of the equation of state is also high.Hence, until some reliable observational data appear regarding the physical conditions in the vicinity of black holes, the overzealous obsession with numerical details remains meaningless and semi-qualitative estimates will suffice.
Here, the three space components of these four equations are the relativistic generalization of Euler equations; the time component for u 0 is the consequence of the other three.Equation (A61) was obtained from Equation (A60) and the definition for T l k , as follows: Here, the second term cancels the fifth term because u k u k = 1, and the fourth term is zero because u k ∂ l u k = 0, for the same reason; from this, we can obtain Equation (A61).
The equations of relativistic fluid dynamics in the general theory of relativity are obtained from Equations (A59) and (A61) by simply replacing the ordinary derivatives with the covariant ones (see Ref. [12], §134) and keeping in mind that the 4-velocity expression is modified as well: Recall that the covariant derivatives for vectors u i and u i are respectively determined as The quantities of Γ k lm are the Christoffel symbols.Also, we recall that these quantities are expressed via the metric tensor g ik as Γ i kl = (1/2)g lm (∂ l g mk + ∂ k g ml − ∂ m g kl ).In Galilean coordinates, Christoffel symbols Γ k lm = 0; therefore, the covariant differentiation reduces to the ordinary differentiation.
Because the covariant derivative of a scalar function produces the same result as an ordinary derivative, Equations (A63) are rewritten as or, for a contravariant component of 4-velocity of flow u s , the first equation becomes When the coordinate grid and spacetime metric are specified, then the group of the last terms in the parenthesis in Equation (A65) manifests itself as a specific "force" acting on a fluid relativistic particle in the "Euler description."If, in analogy with the dynamics of relativistic particles with du i /ds, we wrote, in the left part of Equation (A65), the quantity u k ∂ k u i ≡ (dx k /ds)∂ k u i , this would have been incorrect.Properly, the left part should look as written in Equation (A65).The covariant quantity u k ∂ k u i is not the specific covariant force acting on the fluid particle; this force is the quantity It is helpful to see what will happen in the case of low velocities and in the absence of fields, i.e., when γ → 1, For the α-spatial contravariant component of 4-velocity, we obtain (with the coefficient c −2 ) the "material" derivative: To find the covariant component of the "material" acceleration a α , we use the standard rule for the lowering of indices.
Thus, we become capable of calculating the spatial "force" components, using g αs g is = δ i α , based on the following equation: which has a transparent physical meaning: the left side of the equation determines the kinematic characteristics of the process-the evolution in time and space of the 4-velocitywhich characterizes the state of the fluid particle; the right side of the equation specifies the causes of this change, namely, the pressure gradient and the presence of an external force field.
The subsequent calculations are conceptually rather transparent: for a given metric tensor g ik found from interval ds, defined via ds 2 = g ik dx i d k , we calculate the Christoffel symbols Γ k il ; everything is substituted into Equations (A66), which are solved for the selected model of the equation of state of (disk) matter for the quantities p, e, ϵ, σ.Equations (A66), which can be applied in many physically interesting situations, are quite complicated.The equations of fluid dynamics in the first approximation beyond the Newtonian (i.e., in the post-Newtonian, via expansion with respect to parameter c −2 ) were obtained in Ref. [37], and are discussed in Ref. [33].Additional information can be also found in Refs.[38][39][40][41]44].
Explicit Form of Metric Coefficients for Kerr Geometry with Rotation: For the Kerr metric (see Ref. [30] and references therein), the components of the metric tensor may be found from the following expression: If we need to transition to an uniformly rotating frame of reference, then we make the following transformation x = x 1 cos Ωt − y 1 sin Ωt, y = x 1 sin Ωt + y 1 cos Ωt, z = z 1 , where parameter Ω is the angular velocity of the rotation co-linear with the z-axis.The same notations remain: r → r, θ → θ, ϕ → ϕ + Ωt.Then, This dimensionless expression is expressed in units of length r g = 2GM BH /c 2 and time r g /c.Therefore, r → r g r, â → r g a, Ω → Ω = Ωr g /c.Thus, Equation (A68) can be written as ds 2 = m ik dq i dq k using the following: As is known, for the Kerr model, the parameter is always a 2 < 1/4.From Equation (A69), it follows that dimensionless spacial parameter x is bounded as follows: x + < x < x max (Ω, a, θ) if m 00 > 0. Obviously, setting a = 0 and Ω = 0 in Equation (A69) provides the Schwarzschild metric.Also, the limit value x + corresponds to the external event horizon.
Finally, we provide the expressions for the components of metric tensor m ik in the leading approximation for θ = π/2, a 2 < 1/4 and x ≫ 1: From these expressions, it follows that if the fluid flow is away from the event horizon (i.e., x ≫ 1 > |a| in dimensionless units) and parameter Ω satisfies conditions aΩ ≪ 1 and (approximately) Ωx 3 ≪ 1, then the principal contributions to the components of the metric tensor m ik come from the terms which depend on Ω; the stand-alone gravitational part in the metric is not the focal one.
The expressions in Equation (A70) make it possible to determine the size of the strip of a thin accretion disk, within which the approximation of the classical hydrodynamics can be determined, and where the main contribution to the flow dynamics is made not by the fine structure of the metric in the vicinity of a rotating black hole, but by the effect of the rotation of the accretion disk, characterized by parameter Ω. Parameter Ω appears in these expressions because the accretion disk mainly forms through the captured matter, which has a non-zero angular momentum; therefore, the (averaged) angular velocity of the disk-which is a material formation of finite mass and limited dimensions-is non-zero.Hence, it is logical to transition to a coordinate system rotating with such angular velocity.However, the very concept of a rotating coordinate system contains two implied caveats (see Ref. [31], § 84): (1) some material bodies must exist (or the coordinate system has no "anchors") and (2) they exist within the spacetime domain, which is necessarily bound (or the material bodies outside x max (Ω, a, θ) would rotate with velocities exceeding the speed of light, which is impossible; the limit x max (Ω, a, θ) is where component m 00 of the metric tensor turns negative).This means that a rotating coordinate system cannot extend to infinity.In other words, it is important to remember that all considerations are always made for a finite spacetime domain surrounding the axis of rotation-beyond this domain lies the "forbidden" zone.(Near the axis of rotation is the other "forbidden" zone-the ergosphere.) As the rotation speeds up (Ω → Ω max ), the "forbidden" zones (where m 00 < 0) transform from two detached zones (Figure A1A) into one merged zone (Figure A1B). between two close points can be found from dl 2 = γ αβ dx α dx β .The contravariant 3D metric tensor is γ αβ = −g αβ (see Ref. [31], §84).Also, when g α = −g 0α / ĝ00 , the contravariant component of the same vector g is determined as follows: g α = γ αβ g β = −g 0α .We can also obtain that ĝ00 = ĝ−1 00 − g α g α .Synchronized Time: The introduction of parameter z 0 defined by 0 = √ g 00 (dx 0 − g α dx α )-the so-called synchronized time-is not some kind of mathematical manipulation.It reflects the fundamental need within the framework of relativism to simultaneously measure positions at two separate points and, therefore, to measure the particle velocity.This means that the measured time interval is ĝ00 (dx 0 − g α dx α ), divided by c; the concept of the velocity of fluid particle movement is formulated using this very time interval.Equation allows us to conclude that 3D contravariant and covariant velocities take the following form: Furthermore, from Equations (A75) and (A76), the contravariant components of 4velocity are not introduced arbitrarily, but, following the procedure above, are expressed via the 3-velocity vector v as follows: Equation (A77) provides the link between the 4-velocities u i and 3-vectors of flow velocities v α .Note that, in this case, the operator u k ∂ k can be written as ≡ (γ/c)D, where D = ∂ t 0 + v β ∂ β and t 0 = z 0 /c.In fact, γ(( dz 0 + √ g 00 g α dx α √ g 00 dz 0 Thus, the operator D is nothing more than a "material" derivative with respect to time, in which, instead of differentiation with respect to the time t of a remote observer, differentiation with respect to synchronized time t 0 = z 0 /c is implied. "Specific Force" and Coriolis Effect: Let us return to Equation (A65), which we write using the following form: where, after some simple manipulations, f β becomes a remarkably symmetric expression The term f i in Equation (A78) is non-zero in a curved spacetime and manifests itself as a certain specific force acting on a moving fluid particle.
The key question that immediately arises is whether Equation (A78) provides similar expressions to the centrifugal and Coriolis forces which appear in the classical hydrodynamics.ω = curl v (in terms of components: ω λ = ϵ λµν ∂ µ v ν ) has only one non-zero component A 3 ≡ ψ, i.e., ω 3 is perpendicular to the xy-plane.
We assume that a state of dynamical equilibrium exists: We will now write p = p s + p 1 , ρ = ρ 0 + ρ 1 .Combining Equation (A85) with Equation (A84), we find that a perturbed state is described by the following: In this equation, the small term ∼ ρ 1 Dv β is omitted for similar reasons to those articulated in Refs.[12], §13, or [14]: We suppose that the space-scale of structures which are of interest to us is small in comparison with the distances over which the force-field from m 00 causes a noticeable change in density, and we can regard the fluid itself as incompressible.This means that, during this process, we can neglect the change in the density caused by the pressure change.The change in density caused by thermal expansion, ρ 1 ≃ (∂ρ/∂ T ) p (T − T 0 ), cannot be neglected, because this is the one which causes the phenomenon.
All the necessary attributes in the evolution equation Equation (A86) are presented: on the left is the kinematic quantity, which is dependent on the flow velocity v, which describes the change in the state of the fluid particle in the Euler description, and on the right are the terms characterizing the causes which drive this change-the gradient pressure and the force characteristic determined by the metric coefficient m 00 .
Next, the following well known expressions are used: The condition div v = 0 does not mean that the density of the medium is constant.This condition simply means that the density evolution takes place according to the following equation: Dρ = 0. Using ρ = ρ 0 + (∂ρ/∂ T ) p (T − T 0 ) + (∂ρ/∂p) T (p − p 0 ) ≃ ρ 0 (1 − β p (T − T 0 )), after some simple manipulation, the evolution equation for temperature perturbations takes the form of the third equation in Equation (1).For the equation of state we shall assume, in accordance with the assumption of div v → 0, that the density essentially depends only on the temperature and not on the pressure.Thus, we set ρ = ρ ( 1 − β p (T − T 0 )), where subscript 0 denotes the reference values.The coefficient of thermal expansion β p is presumed constant.Since the quantity β p (T − T 0 ) is generally small, one may neglect the density variations, and hence replace ρ with the constant value ρ 0 in all terms, except in the "buoyancy" term.
Condition div v = 0 is not merely wishful.This condition is applicable in situations where, on the one hand, the magnitude of the instantaneous current velocity v may be considered small compared to the speed of sound s, i.e., v 2 ≪ s 2 , and, on the other hand, the time τ during which the flow configuration changes meaningfully is large compared to the time which is necessary for the sound signal to travel distances of the order of the size l of the vorticial structure in the flow, i.e., τ ≫ l/s (see Ref. [12]).Then, we can assume that information about the disturbance of the medium is transmitted by an acoustic signal, as if instantly, i.e., the medium can be considered incompressible.Obviously, the requirement must also be met that the speed of sound s, although large compared to the value of the local current velocity in the medium, must be less than the speed of light c.
Condition div v = 0 means that the flow velocity v is uniquely determined by specifying the vorticity ω = curl v. Indeed, from the well known identity curl curl b = ∇div b − ∆b for any vector quantity b, it follows that when div v = 0, then ∆v = −curl ω, and, consequently, v = −curl dx 1 G(x − x 1 )ω(x) 1 .Here, G(x − x 1 ) is the Green function of the problem ∆G(x − x 1 = δ(x − x 1 with boundary conditions.The final remark is that since the velocity is determined through the integral of the product of the vorticity q and the Green's function G, even if the model vorticity distribution function is chosen to be "bumpy", the integration process-smoothing-produces a resulting function that is "smooth", i.e., the function is "well-defined" and correctly (at least qualitatively) describes the examined process.

Figure 1 .
Figure 1.First Row, Left Panel (A): Composite image of the black hole Sgr A* derived from radio (1.3 mm) data collected by the Event Horizon Telescope (EHT) Collaboration [1].First Row, Right Panel (B): Multiple cyclones on Jupiter's North Pole [7].Second Row, Left Panel (C): The Antarctic Ozone Hole in Earth's stratosphere; the vortex is quasi-circular on 24 September 2001.Second Row, Right Panel (D): The vortex evolved into a two-petal structure by 24 September 2002.Credit: NASA's Earth Observatory [8].

Figure 2 .
Figure 2. Illustration of a vortex blob.