Covariant Evolution of Gravitoelectromagnetism

The long-range gravitational terms associated with tidal forces, frame-dragging effects, and gravitational waves are described by the Weyl conformal tensor, the traceless part of the Riemann curvature that is not locally affected by the matter field. The Ricci and Bianchi identities provide a set of dynamical and kinematic equations governing the matter coupling and evolution of the electric and magnetic parts of the Weyl tensor, so-called gravitoelectric and gravitomagnetic fields. A detailed analysis of the Weyl gravitoelectromagnetic fields can be conducted using a number of algebraic and differential identities prescribed by the 1+3 covariant formalism. In this review, we consider the dynamical constraints and propagation equations of the gravitoelectric/-magnetic fields and covariantly debate their analytic properties. We discuss the special conditions under which gravitational waves can propagate, the inconsistency of a Newtonian-like model without gravitomagnetism, the nonlinear generalization to multi-fluid models with different matter species, as well as observational effects caused by the Weyl fields via the kinematic quantities. The 1+3 tetrad and 1+1+2 semi-covariant methods, which can equally be used for gravitoelectromagnetism, are briefly explained, along with their correspondence with the covariant formulations.


Introduction
In general relativity, the properties of curvatures are described by the Riemann curvature, which can be split into terms including the Ricci tensor defined by the Einstein equations [1], and the traceless part, called the Weyl tensor [2][3][4] constrained by the Ricci and Bianchi identities [5][6][7][8]. The Weyl tensor can also be split into an electric part, called the gravitoelectric field, and a magnetic part, called the gravitomagnetic field, owing to some similarities to their electromagnetic counterparts [9]. The analogies between electromagnetism and gravitational fields had been demonstrated in other works prior to the 1960s [10][11][12][13][14][15][16], which enabled the spacetime perturbations to be expressed in terms of the electric and magnetic parts of the traceless Riemann curvature with respect to timelike vectors.
The Bianchi identities, which are held by the Riemann curvature, can be decomposed into a set of constraint and propagation equations governing the dynamics of the gravitoelectric and gravitomagnetic fields in a form that rather reminds us of the Maxwell equations [7,8,17,18,18,19]. However, the Weyl gravitoelectromagnetic fields constrained by the Bianchi identities cannot completely describe a solution of the Einstein equations, so we also need the kinematic quantities that are subject to the Ricci identities [18,20,21]. Accordingly, the kinematic constraints of the gravitoelectric/-magnetic fields are provided by the Ricci identities. In this way, the Bianchi and Ricci identities serve as the main equations for analyzing the dynamics and kinematic evolution of the tidal and frame-dragging fields encoded in the Weyl tensor.
The 1 + 3 covariant formalism contains algebraic and differential identities that allow us to express exact (non-linear) solutions to the dynamical and kinematic evolution of the Weyl gravitoelectromagnetic fields. The development of the 1 + 3 covariant approach began with the work of Heckmann and Schücking [22,23], Ehlers [20,21], Kundt and Trümper [5,6], and other early works [7,8,24,25]. The first covariant analysis of perturbations in this context was carried out by Hawking [17]. An elegant covariant form of the Bianchi identities governing the gravitoelectric-/magnetic fields and matter evolution was considered by Ellis [18], which was earlier proposed by Trümper [26]. The covariant formulation was then used to study density perturbations [27], the nonlinear dynamics of comic microwave background (CMB) anisotropies [28], and other cosmological models (see the review by [29]).
The covariant approach has extensively been employed in the studies of the gravitoelectromagnetic fields, whose dynamics and kinematics are constricted by the Bianchi and Ricci identities. The covariant formulations allow the classification of cosmological models, a fluid description of the matter field, and a kinematic description of perturbations in almost-FLRW universes [19,29]. In this literature, we have a covariant study of the evolution of the gravitoelectric/-magnetic tensors in comparison to Newtonian theory [30], as well as a proof of the inconsistency of a universe with a purely gravitomagnetic field [31] and a covariant analogy between the Bianchi equations and the Maxwell equations [32]. An improved covariant method was used to demonstrate that the silent universe, where the gravitomagnetic field vanishes, is inconsistent with the exact nonlinear theory [33].
It is the aim of this paper to review the applications of the 1 + 3 covariant formalism for gravitoelectromagnetism. The paper is organized as follows. In Section 2, we describe the covariant forms of the derivatives of projected vectors and rank-2 tracefree tensors, the kinematic quantities of the fluid, and the dynamic quantities of the matter. In § 3, we covariantly express the gravitoelectromagnetic fields and other algebraic terms of the curved spacetime. In § 4, we see that the gravitoelectric/-magnetic fields, which are not locally affected by the matter field, are indeed coupled to the dynamic quantities of the matter field via the Bianchi identities. In § 5, we show how the curls and distortions of the electric and magnetic parts of the Weyl tensor characterize gravitational waves. Section 6 discusses an irrotational purely gravitoelectric dust model, where the gravitomagnetic field vanishes, as well as a purely gravitomagnetic dust model in the absence of the gravitoelectric field. Section 7 summarizes the covariant formulations for multi-fluid models, followed by a discussion of the tetrad formulation in § 8 and the 1 + 1 + 2 semi-covariant formalism in § 9. In § 10, we discuss observational effects of the Weyl fields induced by the kinematic quantities on electromagnetic waves.
Various 4-velocity vector fields are typically present in a given region of relativistic 1 + 3 spacetime and cosmological models. We choose a timelike 4-velocity field u a to be a unit vector field, i.e. u a u a = −1, and use it to perform a 1 + 3 decomposition of spacetime g ab into the time direction and the 3-dimensional space without questioning the appropriate choice of such a 4-velocity. Accordingly, the metric g ab is projected parallel and orthogonal to u a as follows (see e.g. [19][20][21]): where h ab is the projector tensor yielding the spatial metric containing 3 space quantities, and g ab and h ab have the following properties: g ab u b = u a , g ac g cb = δ a b , δ a a = 4, (2) h ab u b = 0, h a c h cb = h ab , h a a = 3.
The spatially projected alternating tensor is then defined as follows (also compare with the notations used by [39,40]): ε abc = η abcd u d , ε abc = ε [abc] , ε abc u c = 0, where η abcd is the spacetime alternating tensor defined by The above expressions are the basis mathematics for covariant irreducible decomposition of tensors and derivatives, while we have the following practical identities and contractions: e , ε ae f ε de f = 2h a d , ε ae f ε ae f = 6.
The projection of any tensor is denoted by a ⊥ symbol, i.e. ⊥ T ab···d ≡ h a c h b e · · · h d f T ce··· f , where u a ⊥ T ab···d = 0. The spatially projected vectors and projected symmetric tracefree (PSTF) rank-2 tensors are defined as (see Appendix A for higher-rank PSTF tensors): A spatially projected rank-2 tensors S ab can be split into a scalar trace, a projected vector being spatially dual to the skew part, and a PSTF part: Sh ab + ε abc S c + S ab , where S ≡ S cd h cd is the spatial trace, and S a ≡ h ab S b = 1 2 ε abc S [bc] is the projected vector dual to the skew part.
We may also define a vector product and its generalization to rank-2 tensors [41]: The covariant (spacetime) derivative ∇ a of any tensor can be split into the following time derivative and spatial derivative, respectively, Following [33], D a symbolizes the spatially projected part of the covariant derivative. 1 The Fermi derivatives, orthogonal projections of time derivatives along u a , are then denoted bẏ The spatial divergences and curls are defined as, respectively, [33] divV ≡ D a V a , (divS) a ≡ D b S ab , (13) curlV a ≡ ε abc D b V c , (curlS) ab ≡ ε cd(a D c S b) d . (14) If S ab = S (ab) , then (curlS) ab = (curlS) ab . If S ab = S [ab] , we get (S a ≡ 1 2 ε abc S [bc] ): The div curl does not typically vanish for vectors or rank-2 tensors (see [32,33,44,45]). We should remind that D c h ab = 0 = D d ε abc ,ḣ ab = 2u (aub) ,ε abc = 3u [a ε bc]du d , and u aε abc = −u a ε abc . We may also employ the produce and its generalization to define the temporal rotation relative to u a as follows The spatial distortions of vectors and rank-2 tensors are also expressed by [38] As demonstrated by [38], the covariant derivatives of scalars, vectors, and rank-2 tensors can irreducibly be decomposed into projected algebraic terms: where the various algebraic terms prescribed by the kinematic quantities, defined below, emerge from the relative motion of comoving observers. The kinematic quantities of the fluid are defined by [19][20][21] whereu a ≡ u b ∇ b u a is the acceleration vector of the fluid motion under forces, Θ ≡ D a u a is the expansion scalar, σ ab ≡ D a u b is the traceless symmetric (σ ab = σ (ab) , σ a a = 0) shear tensor describing the distortion of the matter flow, and ω ab ≡ D [a u b] is the skewsymmetric (ω ab = ω [ab] , ω a a = 0) vorticity tensor describing the rotation of the matter relative to a non-rotating frame. The vorticity vector ω a is also defined as ω a = − 1 2 ε abc ω bc , where ω a u a = 0, ω ab ω b = 0 and the magnitude ω 2 = 1 2 ω ab ω ab ≥ 0. Accordingly, we may also express the vorticity vector as ω a = − 1 2 ε abc D b u c . 2 In the frame of instantaneously comoving observers, we haveu a =u a andu a = −ḣ ab u b . The shear and vorticity products of a rank-2 PSTF tensor can also be written using the produce and its generalization, [σ, S] a = ε abc σ b d S cd , and [ω, S] ab = ε cd a S b c ω d . The vorticity vector can be split into some algebraic terms encoded by the kinematic quantities: whereω a ≡ h a bω b . The relative motion of the fluid, expansion, rotation, and local distortion are described by the kinematic quantities, while the dynamical effects of the energy and momentum are expressed by the dynamic quantities.
The total energy-momentum tensor T ab is decomposed in terms of the dynamic quantities as follows [47,48] T where ρ ≡ T ab u a u b is the energy density relative to u a , p ≡ 1 3 T ab h ab is the pressure, q a ≡ −T a b u b = −h a c T cb u b is the energy flux (q a u a = 0), and π ab ≡ T ab = T cd h c h ab h cd T cd is the anisotropic stress (π a a = 0, π ab = π (ab) , and π ab u b = 0). Taking q a = π ab = 0, we have a perfect fluid, while additionally p = 0 leads to a pressurefree matter (dust or cold dark matter). A detailed study of analytic equations governing the gravitoelectric/-magnetic fields requires the use of a number of algebraic and differential identities written in terms of the kinematic and dynamic quantities.

Algebraic Identities of the Curvature
In this section, we discuss the algebraic properties of the Riemann curvature. The gravitational fields are locally described by a set of algebraic relations between the Ricci tensor and the energy-momentum tensor, the so-called Einstein equations. However, the tidal forces and frame-dragging effects, which are not directly influenced by the matter field, are described by the traceless part of the Riemann curvature, called the Weyl tensor.
The Einstein field equations [1] describe how the Ricci curvature is affected by the nearby matter field, which, in the absence of a cosmological constant, take the following form where R ab ≡ R c acb is the Ricci tensor, R abcd is the Riemann curvature, and T ab is the energymomentum tensor of the matter field (T ≡ T c c ). Taking the definition (26) of T ab , a successive contraction of the Einstein field equations leads to the following relations: Additionally, the trace of (27) implies R = −T, where R = R a a and T = T a a = −ρ + 3p. 2 Note that the vorticity vector and tensor here have an opposite sign compared to definitions often found in some papers (e.g. [46]). Following [18,20,21], the sign convention is similar to the Newtonian definition The Riemann curvature tensor R abcd is decomposed into the Ricci tensor terms and the Weyl conformal tensor C abcd as follows [2][3][4]: where the Weyl tensor has the following properties The part of the Riemann curvature that is not directly affected by gravitational sources is represented by the Weyl curvature tensor, which describes tidal forces, frame-dragging effects, and gravitational waves. The Weyl tensor is irreducibly split into the gravitoelectric and gravitomagnetic fields [9,49] (see also [17,18,50]): The gravitoelectric tensor represents a covariant Lagrangian description of tidal forces, while frame-dragging effects are provided by the gravitomagnetic tensor. Both of them support the propagation of gravitational waves. These fields are spacelike and traceless symmetric: The gravitoelectric/-magnetic tensors are irreducible, and each has five independent components. They are, in principle, physically measurable in the frame of a comoving observer in such a way that they are encoded in the Weyl tensor: In this way, they enable gravitational action at a distance (tidal forces, frame-dragging, and waves) and affect the motion of matter via the geodesic deviation equation [11,[51][52][53]. Along with the Ricci tensor R ab constrained locally by the matter via the Einstein field equations, the gravitoelectric and gravitomagnetic fields completely describe the algebraic properties of the Riemann curvature. Accordingly, we may decompose the Riemann curvature into its perfect/imperfect matter terms, and gravitoelectric/-magnetic parts, i.e. R ab cd = R ab P cd + R ab I cd + R ab E cd + R ab H cd , as follows [29] R ab P cd = where R ab P cd (and R ab I cd ) is the perfect (and imperfect) matter term, and R ab E cd (and R ab H cd ) is the gravitoelectric (and gravitomagnetic) curvature part.

Spatial Curvature
We define the spatial Riemann curvature as follows [54] (see also [55][56][57][58][59][60] for alternative definitions): where k ab = D b u a is the relative flow tensor between two neighboring observers as defined by (23). In irrotational spacetime (ω a = 0), R ab cd represents the 3-curvature of the space orthogonal to u a with the typical Riemann curvature symmetries. In the presence of an irrotational matter fluid, the tangent planes of fundamental observers closely fit together to create spacelike hypersurfaces orthogonal to their worldlines [61]. They are typical of the u a -congruence and simultaneously represent the hypersurfaces for all comoving observers. However, according to Frobenius' theorem, a rotational spacetime does not have such an integrable hypersurface [62,63], so the observers' planes cannot smoothly go along with each other.
Assuming non-vanishing vorticity, the spatial curvature has the following algebraic properties: In the absence of vorticity, we see that R abcd = R cdab , where the spatial Riemann curvature have the same symmetries of its 4-dimensional counterpart. From the spatial Riemann curvature, one may also define the corresponding spatial Ricci tensor [54]: with the following property [54]: as well as the spatial Ricci scalar describing the local scalar of the space orthogonal to u a [54]: Substituting them into the Einstein field equations yields the following generalized Friedmann [64,65] equation: which depicts how the matter field is associated with the spatial curvature. In the absence of vorticity (ω ab = 0), R is the Ricci scalar of the hypersurfaces orthogonal to the fluid lines.
Considering the decomposition of the Riemann curvature, Eqs. (36)- (39), one could also rewrite the spatial Riemann curvature in terms of the gravitoelectric tensor, along with kinematic and dynamic quantities [61], Contracting (47) on the first and third indices, we derive the spatial Ricci tensor expressed by the so-called Gauss-Codacci equation [61]: An additional contraction of the above equation results in the generalized Friedmann equation (46).

Commutation Laws
Using the spatial derivative definition (11), one can arrive at the commutators of the spatial derivatives for scalars, vectors, and rank-2 tensors. The key point to be considered is that for any scalar function f we have, Applying the Leibniz rule of ∇ a to the last bracket and using Eq. (23), one obtains the following purely relativistic expression [54] The above scalar commutator corresponds to the relativistic behavior of rotating spacetimes in general relativity. Similarly, for any vector field V a orthogonal to u a (V a u a = 0), we derive the following vector commutator, Furthermore, for any rank-2 tensor field S ab orthogonal to u a (S ab u a = S ab u b = 0), we obtain The above fully nonlinear commutators hold at all perturbative levels. In the absence of rotation (ω ab = 0), R abcd corresponds to the Riemann curvature of 3-D hypersurfaces orthogonal to the u a -congruence. However, time derivatives do not typically commute with their spacelike counterparts. In particular, for any scalars, we have the following expression at all perturbative levels [28] which is the key equation for solving the dynamical evolution of spatial gradients that are covariantly associated with inhomogeneities [27]. Similarly, for any vector field V a and tensor field S ab orthogonal to u a to first order, we get: Contracting the above equations leads to their divergences to first order: In an almost-FLRW background, the orthogonally projected gradient and time derivative of the first-order vector V a and spacelike tensorṠ ab possess the linear commutation laws aD aVb = (aD a V b ) · and aD aṠbc = (aD a S bc ) · . From Eq. (50), for any scalar, it also follows which depicts the relation between vorticity and non-integrability, implying that in the case of non-zero vorticity there is no constant-time 3-surfaces everywhere orthogonal to u a , since the instantaneous rest spaces cannot smoothly mesh each other.

Identities in Irrotational Dust Spacetimes
One of the most commonly used fluids in cosmology is the dust model (p = q a = π ab = 0), while the irrotational dust model (ω ab = 0) also appears to be useful for the late universe (see [33] for relevant calculations). Here we summarize some properties of irrotational dust spacetimes.
In the case of dust spacetimes, the Ricci tensor is R ab = 1 2 ρ(u a u b + h ab ), so the Riemann curvature is written as The dab S ce , together with (59) in irrotational spacetime (∇ b u a = 1 3 Θh ab + σ ab andu a = 0 = ω ab ) lead to the following key identities [33]: Eq. (65) generalizes the Ricci identities for the commutation of spatial derivatives of rank-2 tensors. There are the further important identities [33]: where R ab is the spatial Ricci tensor defined by (48), and D 2 = D a D a is the spatial Laplacian operator. Linearized forms (σ ab = E ab = H ab = 0) of (68)- (70) for the spatial and time derivative of the curl, as well as the curl of the curl can be written, respectively, For a dust congruence, we imposeu a ≡ D a Φ = 0, so Eq. (50) leads to D [a D b] Φ = −ω abΦ = 0. Therefore, for the fundamental observers seeing an isotropic radiation field in a dust spacetime, the spacetime is locally FLRW (EGS theorem [25], see also [66][67][68][69]). 3

Analytic Identities for the Weyl Tensor
In this section, we describe the analytic identities that provide dynamical and kinematic constraints for the gravitoelectric/-magnetic parts of the Weyl tensor. The Ricci and Bianchi identities are the main equations that govern the gravitoelectric and gravitomagnetic fields, while Einstein's equations correspond to the algebraic relation between the Ricci curvature and the matter field. We see that the decomposition of Bianchi identities facilitated by the covariant formalism leads to the dynamical and evolutionary equations for the gravitoelectric/-magneticc tensors.
We consider the Bianchi identities as field equations for the free gravitational field, As shown by Kundt and Trümper [5,6], the Bianchi identities are equivalent to the following form (see also [18,19,50,53]) Substituting the dynamic quantities (26) into the Bianchi identities and decomposing them according to the formalism introduced in the previous sections, the following equations are obtained [7,8,18,19,29,46]: The decomposition of R abcd into R ab and C abcd combined with the once-contracted Bianchi identities yield two constraints and two propagation equations for the gravitoelectric and gravitomagneticc tensors, which are analogous to the Maxwell equations in an expanding spacetime (see [17,19]). The twice contracted Bianchi identities describe the conservation of the total energy momentum tensor, which are decomposed according to the covariant formalism as follows [19,29]: In a perfect-fluid model, they are reduced to the typical energy-density and momentumdensity conservation laws:ρ Eq. (81), which is the conservation of total energy-momentum law, describes how the matter field determines the geometry, i.e., the motion of the matter. Considering the Einstein field equations defined by G ab ≡ R ab − 1 2 g ab R = T ab − Λg ab , the twice-contracted Bianchi identities (81) leads to ∇ b G ab = 0 (∇ b Λ = 0 if Λ is constant in space and time) that describes the dynamical evolution of the spacetime and the matter in it.
The kinematic equations are provided by the Ricci identities of the vector field u a , i.e.

Perfect-fluid Models
Using the 1 + 3 covariant notations introduced in Sec. 2, Eqs. (77)-(80) for a perfectfluid matter can simply be rewritten as follows: Based on Eq. (93), the spatial gradient of the energy density covariantly emerges as a source for the gravitoelectric field, while the shear and vorticity products of the gravitomagnetic field are associated with the inhomogeneity in the general-relativistic fluid model. This implies that the gravitoelectric field represents a generalization of the Newtonian tidal force [18]. In Eq. (94), the gravitomagnetic field originates from the angular momentum density, i.e., (ρ + p)ω a , as well as the shear and vorticity products of the gravitoelectric field.
In a perfect-fluid model, the Ricci equations (87)-(92) can also be simplified using the covariant conventions: According to Eq. (99), the gravitoelectric tensorial field can directly induce the acceleration distortion and the time derivative of the shear. It can be seen in Eq. (102) that the vorticity distortion and the shear curl can directly originate from the gravitomagnetic tensorial field (see [32] for more discussions). Although other equations do not involve the gravitoelectric and gravitomagnetic fields, they are important constraints for the kinematic quantities. The Raychaudhuri equation (97) is the only equation that includes the dynamic quantities, which can be employed to describe the gravitational attractive force of the nearby matter (see e.g. [18][19][20][21]).

Nonperturbative Shearless Spacetimes
We now consider small perturbations of the fluid motion and the Weyl gravitoelectric/magnetic fields in a way that neglects products of small quantities and performs derivatives relative to the undisturbed metric. In a nonperturbative shearless (σ ab = 0) perfectfluid model, we avoid perturbations that are merely associated with coordinate transfor-mation and have no physical significance, so Eq. (93)- (99) can be rewritten in first order as follows (also compare with [71]): It can be seen that the evolution of the gravitoelectric/-magnetic fields or their curls do not lead to perturbations of the expansion, vortecity, or density according to the first-order dynamical equations in nonperturbative shearless spacetimes (see also [17]).

Irrotational Dust Spacetimes
Irrotational dust spacetimes have widely been used to study the late universe, as well as for the evolution of density perturbations [36] and gravitational waves [33]. Here we rewrite the Bianchi and Ricci equations for an irrotational (ω ab = 0 =u a ) dust (p = 0) model using the covariant formulations, and similar to the nonperturbative shearless model we also perform derivatives with respect to the undisturbed metric, so we have Irrotational dust models with realistic inhomogeneities accommodate both the gravitoelectric and gravitomagnetic fields. In the linearized theory, where the model is almost close to an FLRW spacetime, gravitational waves in irrotational dust spacetime are typically characterized by transverse traceless tensor modes. We should note that an FLRW spacetime is covariantly described by D a ρ = 0 = D a Θ and σ ab = E ab = H ab = 0, whereas the density and expansion gradients and the gravitoelectric/-magnetic tensorial fields are first order of smallness in almost-FLRW spacetimes such as a linearized irrotational model. In the linearized theory, imposing E ab = 0 leads to the vanishing of anisotropy and inhomogeneity, resulting in an FLRW spacetime, where H ab = 0, so a linearized purely gravitomagnetic irrotational dust spacetime deos not exist [31]. Covariant calculations can be conducted to show that the constraint equations (109) and (113) are preserved under the evolution in irrotational dust spacetimes. Let us assign C A = 0 (where A = 1, . . . , 4) to the following constraints: The time derivatives of the above constraints (C A ) generate a set of equations that can be denoted byĊ A = F A (C A ), where the terms F A do not contain any time derivatives, since they are removed by using the propagation equations and relevant identities [31,33]: Assuming an initial spatial surface {t = t 0 }, where t is a proper time along the worldlines, is satisfied by the constraints C A , i.e. C A | t 0 = 0, the evolution of the constraints (Ċ A ) implies that the constraints should be satisfied for all time, as C A = 0 is a solution for the initial condition. The constraint set is linear, so the solution should be unique. Thus, the constraints C A are preserved under evolution. For example, if we provide σ ab (t 0 ) and D a ρ(t 0 ), C 1 a produces D a Θ(t 0 ), C 2 ab generates H ab (t 0 ), and C 3 a presents D b E ab (t 0 ). A consistency condition is expected to be imposed by C 4 a , though this is not the case, since we have Let us consider that D a Θ is determined by C 1 a , H ab by C 2 ab , and D a ρ by C 3 a , so the constraint equations are consistent with each other owing to the presence of C 4 a . Therefore, if we have a solution to the constraints on {t = t 0 }, it is consistent and evolves consistently. Eq. (123) implies that no new vectorial constraint emerges from the divergence of the tensorial constraint C 2 ab , so the tensorial constraint is characteristically transverse traceless. Supposing the gravitomagnetic tensorial field that is divergence-free (a condition for gravitational waves), the constraint (118) leads to The covariant formulations can be employed to show that the condition (124) is satisfied under the evolution without considering further conditions. In particular, the condition (124) does not impose H ab = 0, so it establishes spacetimes with D b H ab = 0 = H ab .

Imperfect-fluid Models
We now covariantly write the exact (nonlinear) Bianchi equations in imperfect-fluid models. Substituting the total energy-momentum tensor (26) into the Bianchi equations (77)- (80) and applying the 1 + 3 covariant notations yield: It can be seen in Eqs. (125) and (126) that the shear and vorticity coupled with the energy flux and gravitomagnetic field act like sources for the gravitoelectric field, whereas the shear and vorticity coupled with the anisotropic stress and gravitoelectric field appear as sources for the gravitomagnetic field. The long-range gravitational fields, associated with tidal forces, frame-dragging effects, and gravitational waves, especially making tensorial contributions to CMB anisotropies [28], are ruled by Eqs. (125)- (128). The implications of inhomogeneous coupling terms in the nonlinear Bianchi equations are beyond the scope of this review and deserve further discussion (see e.g. [32]). In irrotational spacetimes (ω a = 0) under homogeneous and isotropic conditions (σ ab = π ab = 0), Eqs. (125)-(128) simply reduce to: where H = 1 3 Θ is the Hubble parameter. We see that the gradient of the energy density and the energy flux constrained by the Hubble expansion rate act as sources for the gravitoelectric divergence, while the curl of the energy flux is a source for the gravitomagnetic divergence. Additionally, Eqs (131) and (132) present wave solutions if we take the curl of curl(H ab ) and the time derivative of curl(E ab ), and applying the linearized identities (73) and (74) to them, resulting inḦ ab − D 2 H ab . Similarly, the covariant calculations of curlcurl(E ab ) and curl(Ḣ ab ) in the linearized forms give a wave solution for the gravitoelectric tensor, i.e.,Ë ab − D 2 E ab .

Linearization Scheme
The exact dynamics of the gravitoelectric/-magnetic tensorial fields with the full matter including flux and anisotropic terms are covariantly governed by the nonlinear Bianchi equations (125)- (128). We may also linearize these equations around any chosen background such as an FLRW metric. In particular, in the FLRW background, all the inhomoge-neous and anisotropic terms (q a and π ab ) vanish, and only the first order of the quantities appear in the linearized Bianchi equations.
Linearization of Eqs. (125)- (128) leads to the following system of constraint and propagation equations: Furthermore, the Ricci identities (89) and (92) according to the linearization scheme are written as follows: Once again, the gravitoelectric and gravitomagnetic fields originate from the energy density gradient and the angular momentum density in Eqs. (133) and (134), respectively. Moreover, we see that the anisotropic stress gradient and the energy flux constrained by the expansion scalar are additional sources for (divE) a , while the energy flux curl acts as a source for (divH) a . The anisotropic stress π ab still remains in the propagation equations (135) and (136) after linearization. We also have the energy flux distortion in the curl(H ab )equation. The full energy-momentum tensor that includes the flux and stress allows us to investigate inhomogeneities and anisotropies in cosmological models. In particular, the anisotropic stress can be used to analyze the CMB (see e.g. [28,72,73]).

Purely Tensorial Perturbations under Linearization
We now consider the so-called purely tensorial perturbations [74] under a linearization regime, which are characterized by vanishing vector and scalar variables, i.e. f = D a f = D a D b f = 0 and V a = D b V a = 0 for any scalars f and vectors V a . Imposing only tensor perturbations, the linearized equations (133)-(138) shall be: The time derivative and curl of equations (140) and (141), as well as the time derivative of (142a), together with the conditions D b E ab = 0 and D b H ab = 0, lead to [74] (see [75] for detailed calculation): where c s is the adiabatic sound speed of the fluid flow, and w is the effective barotropic index [74] (c 2 s =ṗ/ρ and w = p/ρ in FLRW spacetimes [37]). Solving Eqs. (143)-(145) yields the evolution of tensor perturbations. Note that Eq. (144) is effectively third order due to the presence of a term with σ ab coupled with the energy density, which makes it impossible to write a closed wave equation for E ab . However, taking π ab = 0 for consistency, Eqs. (143) and (144 ) should result in wave solutions since the shear also maintains a wave solution provided by Eq. (145).

Gravitational Waves
Gauge-invariant transverse tensor perturbations in homogeneous and isotropic FLRW models are known as gravitational waves. The wave solutions of the gravitoelectric/magnetic tensorial fields were first introduced by Hawking [17] in a model, where perturbations of the Weyl tensor do not originate from rotational or density perturbations, and were later covariantly extended by Ellis and Bruni [27] in the gauge-invariant perturbation approach.
Considering the divergence equations (103), to have a wave solution, we require the condition: Taking a time derivative from the curlH ab equation (104) and a spatial curl from the curlE ab equation (105), after simplification, we have [76,77] In empty, non-expanding (Θ = 0) space, it reduces to the typical form seen in the linearized theory, i.e. ∇ 2 E ab ≡ ∇ a ∇ a E ab = 0. From Eq. (147), it follows that that a necessary covariant condition for the propagation of gravitational waves should be: Therefore, the vanishing divergences and non-vanishing curls of the gravitoelectric and gravitomagnetic tensors are the necessary covariant conditions for the existence of gravitational waves. The gravitomagnetic field obviously has an essential role in supporting gravitational waves. Based on the similarity with electromagnetism, where neither the electric nor magnetic fields independently provide a full description of electromagnetic waves, it was argued that both the gravitoelectric and gravitomagnetic fields are simultaneously necessary for supporting tensor perturbations [78]. Indeed, as seen in Eq. (147), gravitational waves are exactly characterized by the spatial curls of the electric and magnetic parts of the Weyl tensor.
The nonlinear evolution of E ab and H ab , with a perfect-fluid source are determined from Eqs. (95) and (96). The only differences between these two equations are the signs of the right sides, besides the energy density and pressure coupled with the shear in the curlH ab equation. Once they are linearized about an FLRW background, they reduce to Eqs. (104) and (105). On taking the time-derivative of the former and substituting the curl of the latter, the curlE ab and curlH ab terms support the propagation of gravitational waves. However, Eqs. (95) and (96) do not close up in general, owing to the presence of the 1 2 (ρ + p)σ ab term in the curlH ab equation. It is necessary to add the shear evolution equation to the curlH ab equation to obtain a third-order equation for E ab , once its time derivative is calculated (see e.g. [77]).
In the linearization scheme, purely tensorial perturbations are obtained by requiring the density perturbations and rotational perturbations vanish to first order. The condition that the terms on the right hand side vanish from the constraint (103) is similar to the transverse condition of tensorial perturbations in the metric approach. We remind that the Weyl tensor is the traceless part of the Riemann tensor, so both E ab and H ab are traceless, again similar to the tensorial perturbations of the Bardeen formalism [79]. 4 As the condition (146) is required for gravitational waves, one might expect that both E ab and H ab possess the SO(2) electric-magnetic symmetry (see the review by [81]), so their propagation equations should be invariant under the transformation E ab ⇔ H ab . However, this is not generally true; indeed, it is valid when the equation of state of the background satisfies special conditions. In particular, the evolution equations for E ab and H ab are not of the same order, the E ab wave equation is of third order, whereas H ab has a second-order wave equation. However, both of them reduce to second order if we suppose (ρ + p) = 0, i.e. in vacuum (also de Sitter spacetime), as well as for the equation of state p = − 1 3 ρ. Thus, under those special conditions, both E ab and H ab provide the second-order wave solutions [77], which also have the SO(2) electric-magnetic invariance [81].

Distortions of the gravitoelectromagnetic fields
The distortion of the curvature is a part of its derivative that is dissociated from the Bianchi identities and is not locally connected to the matter. The divergences of the locally free fields E ab and H ab are pointwise connected to the dynamic quantities of the matter, such as the energy density gradient. The curl of E ab (or H ab ) are pointwise associated with the time derivative of H ab (or E ab ) and the matter terms. Additionally, the distortion of E ab (or H ab ), similar to the curl, covariantly characterizes a locally free part of the spacetime gradient of the gravitational field. It was shown by [38] that the spatial Laplacian (D 2 ≡ D a D a ), which is essential for describing a wave solution, can be produced by the distortion. Thus, the non-vanishing distortions of E ab and H ab , along with the non-zero curls, is another condition for the propagation of gravitational waves. 4 In the Bardeen formalism [79], the equation of motion isḦ T = 0 that is derived for the first-order gauge-invariant amplitude of the tensor perturbation H (2) T in the absence of anisotropic stress perturbations, where ℓ is the cosmological scale factor, K is the spatial curvature of the FLRW background, and k is the physical wavenumber with K = 0 [80]. The metric tensor perturbation is H Considering the decomposition (21) of the covariant distortion of rank-2 tensors, the covariant and spatial distortions of E ab are split into various terms enclosed by the kinematic quantities [38]: whereÊ cab is the symmetric, traceless gravitoelectric distortion defined byÊ cab =Ê (cab) ≡ D c E ab , andĖ ab ≡ h (a c h b) dĖ ab − 1 3 h cdĖ cd h ab is the time derivative of the spatially projected gravitoelectric tensor. 5 Similarly, ∇ c H ab can also be decomposed. Comparison of Eq. (150) with Eq. (23) indicates that the gravitoelectric and gravitomagnetic distortions,Ê cab andĤ cab , which are the divergence-less and curl-less spatial variation of the Weyl tensor, could be analogous to the shear σ ab of the fluid flow. We may also apply the decomposition (149)-(150) to any PSTF tensors such as the shear σ ab : Although the divergence and curl of σ ab are previously seen in the Ricci identities (100) and (102), respectively, its distortion does not appear in the Ricci and Bianchi equations. If we take a spatial distortion from (99), it is noticed that the evolution ofσ cab is controlled byÊ cab . The Laplacian equation, which plays a key role in the wave solution, also emerges from the distortion. Let us consider the gravitomagnetic distortion: which has the following divergence equation: Now if we use the commutation (65), the last term on the right of Eq. (153) can be transferred into a divergence term and curvature correction terms [38]: where N [H] ab is a non-linear term defined by [38], In an almost-FLRW model, these covariant terms in N [H] ab are very small [27], so this non-linear term could be negligible. As the non-linear term N [H] ab vanishes in the linearization mode, while the divergence is also zero for gravitational waves, we get If the gravitomagnetic distortionĤ cab vanishes, a spatial Helmholtz equation is satisfied by H ab , which as argued by [38] excludes general gravitational waves (see also [17]). A similar result follows if we start with zero distortion of E ab . Hence, a necessary condition for gravitational waves is non-vanishing distortions of the gravitoelectric and gravitomagnetic fields, i.e., Therefore, the non-zero curls and distortions are essential for the propagation of gravitational waves. Additionally, the distortions of the gravitoelectric and gravitomagnetic fields can also broaden our understanding of cosmological models, which are typically characterized by the divergences and curls of E ab and H ab (see e.g. [18,33,[82][83][84]).

Purely Gravitoelectric/-magnetic Models
The gravitoelectric field E ab is a general-relativistic generalization of Newtonian tidal forces, while there is no analogy for the gravitomagnetic field (H ab ) in Newtonian theory [18,30]. However, a model with H ab = 0 cannot support gravitational waves as shown in dust spacetimes [76,77], so called silent universe owing to the absence of propagating signals [85][86][87]. The silent universe (H ab = 0) is sometimes called Newtonian-like due to the presence of a purely Newtonian counterpart, i.e. the gravitoelectric field. The Newtonianlike model with only E ab corresponds to the general-relativistic generalization of Newtonian theory. However, a consequence of vanishing H ab in post-Newtonian models is that the locally free characteristics of E ab cannot restore in the Newtonian limit, implying that it is generally incorrect to suppose H ab = 0. The purely gravitoelectric model (H ab = 0) has limited applications in cosmological models, particularly in studies of gravitational instabilities [83,88,89], so it is essential to have the gravitomagnetic field in physically realistic inhomogeneous models of the late universe. Models with H ab = 0 cannot be Newtonianlike, but they consistently satisfy the locally free fields in general [30].
A purely gravitomagnetic model with E ab = 0 = H ab is sometimes called anti-Newtonian [31] due to the presence of only a field with no Newtonian counterpart. A purely gravitomagnetic model (E ab = 0) was first studied in [24] in which either the shear or the vorticity are shown to be non-vanishing. We note that E ab = 0 = H ab is associated with an FLRW model (see § 4.1.2). A linearized irrotational dust model with E ab = 0 was also found to be exactly FLRW, resulting in H ab = 0 [31], so linearized anti-Newtonian irrotational dust models are inconsistent and cannot generally exist.

Purely Gravitoelectric Dust Spacetimes
Here we discuss an irrotational dust Newtonian-like model (H ab = 0), so-called silent [87], where E ab acts as Newtonian tidal forces [30,88]. We know that the gravitomagnetic tensor H ab has no Newtonian analogue. The silent models are known to be generally inconsistent and not likely to surpass the spatially homogeneous spaces [30,31].
Consider the evolution of the constraints (Ċ A ) in irrotational dust spacetimes, namely Eqs. (119)- (123). Although placing H ab = 0 modifies the constraints C A , they are still consistent, since theĊ A -equations are still valid. Eq. (111) offers an extra constraint for the Newtonian-like model: which should be satisfied, along with the propagation. The spatial divergence and the time derivative of C 5 ab are calculated by applying the algebraic identities (67)-(69) as follows [31]: where H ab is defined by [31] H The consistency of the evolution of the constraints in the Newtonian-like model is satisfied if H ab = 0. Eq. (161) implies the condition H ab = 0 and its evolution is equally fulfilled by linearization around an FLRW background, which is also characterized by D a ρ = D a Θ = 0 and σ ab = E ab = H ab = 0 [27]. The purely gravitoelectric models generally have inconsistent solutions that are contingent on linearization instabilities in such a way that their linearized solutions are not constrained by any consistent solutions of the exact (nonlinear) equations. According to [83], it is unlikely that silent solutions could be a physically realistic model for the late universe or gravitational instabilities. Consistent, realistic models require the gravitomagnetic field, as independently confirmed by [88]. As shown by [31], the Newtonian-like silent models have a few limited applications, so general relativity does not straightforwardly correspond to Newtonian theory.

Newtonian Limit
A Newtonian-like universe with only the Poisson equation for gravitational potential is impractical without some extensions to Newtonian theory as pointed out by Heckmann and Schücking [22,23,90]. Newtonian theory should be extended in a way to satisfy the essence of the non-locality, so that local physical laws cannot be detached from instantaneous boundary conditions at infinity. The singular limit of general relativity, where gravitational interaction is at an infinite speed, reduces to the Newtonian theory. Obviously, this Newtonian limit of general relativity excludes the gravitomagnetic field.
Here we discuss a covariant Newtonian approximation of general relativity, known as the Heckmann-Schücking limit, having boundary conditions for spatially homogeneous spacetimes. Under certain boundary conditions, the gravitoelectric tensor can covariantly be approximated in the Newtonian limit by where E (N) ab is the gravitoelectric tensor equivalent to tidal forces in the Newtonian limit [30,88], and Φ is the Newtonian potential. 6 In irrotational spacetimes, Eq. (98) reduces to curlu a = 0, which implies the presence of an acceleration scalar potentialΦ defined by [20,21,91]: It can be seen that the acceleration vector is covariantly proportional to the spatial gradient of a scalar. According to the equivalence principle in general relativity, it is indeed the Newtonian potential (Φ = Φ). Eq. (163) can also be solved for the velocity field u a defined by the normalization of the stationary Killing vector ξ a = ξu a . Accordingly, Killing's equations result in Θ = σ ab = 0 [20,21,92], so Eq. (23) reduces to D b u a = − 1 3 ε abc ω c . As Φ is invariant under ξ a , we should haveΦ = 0 of the Newtonian limit, and Eq. (163) satisfies curlu a = 0. Thus, as covariantly described by Eq. (163), the acceleration vector in irrotational spacetimes strictly corresponds to the gradient of the acceleration scalar potential Φ.
The Newtonian limit of E ab in the Heckmann-Schücking approach satisfies the condition, lim c→∞ E ab = E ab (t)| ∞ , in a way that the gravitoelectric field (162) prescribed as a function of time arbitrarily propagates to any point with infinite speed (c → ∞), although is constrained by Eq.(103), resulting in the Poisson equation of Newtonian gravity: However, we have no analogue of H ab in the Newtonian limit [18], as shown by strictly mapping general relativity onto Newtonian theory in the singular limit [93]. In Newtonian theory, we do not haveΦ , as well asĖ ab seen in the propagation (104). We may also approximate the gravitoelectric field to first order as E ab = E ab is the Newtonian term given by Eq. (162), and E (PN) ab is the first post-Newtonian term. In the Newtonian limit, the gravitoelectric field E ab simply reduces to E (N) ab , which is tidal forces.

Quasi-Newtonian Gravitational Attraction
The Raychaudhuri equation (97) [34,70], which is the basic equation for gravitational attraction [18][19][20][21] and the singularity theorem [19,29], can also be employed to provide a simpler interpretation of the Newtonian limit [18,94]: where A a ≡u a = D a Φ is the acceleration vector in the Newtonian limit. In Eq. (165), the (ρ + 3p) acts as the gravitational source. In the static case, where the expansion of the timelike congruence vanishes (Θ = 0), and the Raychaudhuri equation represents the general-relativistic generalization of the Poisson equation of gravity. In a quasi-Newtonian model [45], where a congruence of worldlines is irrotational and shearless (ω a = 0 = σ ab ), it reduces to the well-known static gravitational attraction: The above equation covariantly generalizates the Poisson equation of Newtonian gravity, correlating the accelerationu a with the active gravitational mass (ρ + 3p). In quasi-Newtonian dust spacetimes (p = 0), Eqs. (164) and (166) imply D au a +u au a = D 2 Φ, which is the Laplace equation of the Newtonian potential. In quasi-Newtonian dust models, the Raychaudhuri equation determines the acceleration divergence. However, in expanding spacetimes, it includes the evolution of the expansion, along with the divergence of the acceleration.

Purely Gravitomagnetic Dust Spacetimes
We here consider a case of purely gravitomagnetic dust models, known as anti-Newtonian universes [31]. Substituting E ab = 0 and curlE ab = 0 into Eqs. (111) and (112) yields the following two constraints: The former equation determines the evolution of the shear σ ab , whereas the latter equation is associated with the propagation of H ab without influencing the other quantity evolution. It is seen that the the matter evolution is entirely detached from the gravitomagnetic fields, i.e. no coupling between the Weyl tensorial fields and the matter source in the geodesic deviation equation [95]. The evolution of the dynamic and kinematic quantities is controlled by Eqs. (114) and (167), which is consequently supplied to the propagation of the gravitomagnetic field via Eq. (168).
In the case of irrotational dust anti-Newtonian spacetimes ( § 4.1.2), the constraint (118) implies D b H ab = 0, which still maintains Eq. (123) with E ab = 0. However, the propagation is no longer provided by Eq. (110) with vanishing the gravitoelectric terms, which turns into a new constraint analogous to the Newtonian-like equation (158): The spatial divergence and the time derivative of C 6 ab are calculated as follows [31]: where J a and E ab are defined by [31] J a ≡ΘD a ρ − 3ρD a Θ − The evolution of the constraints is consistent in the anti-Newtonian model if the two conditions J a = 0 and E ab = 0 are satisfied. Thus, the consistency of the constraints on an initial surface is covariantly fulfilled by the following condition, We see that an algebraic relation between the spatial gradients of the energy density and the expansion scalar appears as a main integrability condition for the consistency. There is at least one special situation where this condition is identically satisfied (174). Nevertheless, only specific initial conditions {ρ, Θ, σ ab } on the initial surface can be incorporated into the covariant condition (174), which lead to inconsistencies of the evolution equations (114) and (167) in general cases. Thus, anti-Newtonian spacetimes are generally inconsistent. Linearization about an FLRW background under the consistency condition J a = 0 implies that, unlike the Newtonian-like model, linearized integrabilities in the anti-Newtonian case are not insignificant [27]. Indeed, all the linearized anti-Newtonian so-lutions appear to be inconsistent, since the linearized integrability conditions are satisfied only with H ab = 0. As we already impose E ab = 0 in the anti-Newtonian case, H ab = 0 shall make an FLRW spacetime. Linearization about an FLRW space under the condition E ab = 0 implies Θσ ab = 0 that requires either Θ = 0 or σ ab = 0. Moreover, the linearized form of the covariant condition J a = 0 is 3ρD a Θ − ΘD a ρ = 0 that holds if one imposes Θ = 0. If we additionally have σ ab = 0, there is an FLRW model, which is is not anti-Newtonian. We also see that vanishing the expansion scalar results in ρ = 0 according to the linearized form of Eq.(114a), subsequently σ ab = 0 based on Eq.(114b). Hence, we must rule out any possibilities of linearized anti-Newtonian models (see [31] for full discussion).

Anti-Newtonian Limit
In the Newtonian limit, the interaction propagates to any point at an infinite speed, and the gravitoelectric field to first order is roughly written as E ab = E  (2) electric-magnetic duality transformation E ab → H ab (see e.g. [81])), the gravitomagnetic tensor may also be written similar to Eq. (162) in the anti-Newtonian limit as follows: where Ψ is the anti-Newtonian potential. The above equation is based on the gravitoelectric/magnetic duality invariance (see [32,81,96]). The electric-magnetic duality has important implications in quantum field theory.
Considering (100) and (101) in an anti-Newtonian irrotational shear-free model (Θ = σ ab = 0), along with the velocity field u a defined by the normalization of the stationary Killing vector ξ a = ξu a , the constrain (100) then yields curlω a = −2[u, ω] a [20,21]. It can be seen that the vorticity curl corresponds to the product of vorticity and acceleration, which implies the presence of a vorticity scalar potentialΨ defined by The vorticity vector is then covariantly characterized by the vorticity scalar. Suppose the equivalence between the anti-Newtonian potential and the vorticity scalar (Ψ =Ψ), by analogy with the equivalence principle (Φ =Φ), the ani-Newtonian limit (175) of H ab constrained by Eq.(103b), results in the Helmholtz equation [41]: The Bianchi equations relate the gravitomagnetic tensor to the source (ρ + p)ω a , in another word, the angular momentum density [18]. The gravitomagnetic tensor H ab , which has no Newtonian analogue, is associated with either gravitational radiation [78] or cosmic inflation [97]. In § 5, we saw that the spatial curls of both the gravitoelectric and gravitomagnetic fields characterize gravitational waves [38]. In post-Newtonian models, the non-locality -in a way that local physics cannot be decoupled from boundary conditions at infinity -of Newtonian tidal forces induced by the gravitoelectric field cannot be restored without the gravitomagnetic field.

Multi-fluid Models
There are many situations in cosmology where it is necessary to describe a system containing multiple fluids of various matter species rather than a single fluid, which is made possible by studies of multi-component systems (see the review by [98]). There are many examples of multi-fluid systems such as CMB inhomogeneities, including radiation, baryonic matter and neutrinos [72], and incorporating peculiar velocities into nonlinear gravitational collapse [99]. In either multi-component systems or analyses of peculiar velocities, it is required to include the velocity tilt between the matter species and the fundamental observers (see [37,73,[98][99][100][101]).
In a multi-fluid system, the 4-velocity u a of the fundamental observers and the 4velocities u (i) a of i-th species satisfy u a u a = u (i) a u a (i) = −1, so the projector tensors orthogonal to u a and u (i) a are defined, respectively, as follows: The velocity u a is related to u where v (i) a is the peculiar velocity of the i-th species with respect to u a , and a v a (i) ) −1/2 is the Lorentz-boost factor, so γ (i) = 1 in non-relativistic peculiar motion. The boost relation is a reaction to the hyperbolic angle of tilt β (i) [100] between u a and u (i) a in a way that satisfies which turns Eq. (179) to [100] u (i) It can easily be seen that v (i) = tanh β (i) , so v (i) ≃ β (i) for small tilt angles (β (i) ≪ 1) in non-relativistic peculiar motion. The dynamic quantities in a multi-fluid system should be defined to include all dynamically significant species: where ρ (i) , p (i) , q (i) (a , and π (i) ab are the energy density, pressure, energy flux, and anisotropic stress measured in the frames of i-th species with u (i) a defined by Eq. (179). As an example, a multi-fluid system may include the electromagnetic radiation (i = R), baryonic matter (i = B) modeled by a perfect fluid, cold dark matter (i = C) described by a dust model over the era of CMB anisotropies, neutrinos (i = N), and dark energy modeled by a cosmological constant (i = V) (see [28] for further details).
The inverse form of Eq. (179) can be expressed as (wherev The above equations directly lead to the following nonlinear energy-momentum tensors of i-th species measured in the fundamental u a -frame: where the nonlinear dynamic quantities of i-th species are (v 2 ≡ v a v a ) [28,102]: These nonlinear dynamic quantities generalize the well-known linearized results (see [37,103]). The total dynamic quantities are simply obtained using To linear order, the dynamic quantities measured in the i-th frame are roughly equal to those in the fundamental frame, expect for the energy flux q a that needs a velocity correction:ρ In the nonlinear situation, as seen in Eqs. (186)-(189), the dynamic quantities observed in the i-th frame and the fundamental frame are not generally identical.

Multi-Perfect Fluids
Here we consider a multi-component system consisting of i-th perfect-fluid species defined by energy density ρ (i) and pressure p (i) moving along the timelike 4-velocity u (i) a . The energy-momentum tensor of each i-th fluid with respect to u (i) a is given by where h (i) ab defined by Eq. (178). With respect to the fundamental frame u a , Eq. (192) turns into an imperfect fluid as follows where the dynamic quantitiesρ (i) ,p (i) ,q ab are given by [37] In non-relativistic peculiar motion (β (i) ≪ 1), quadratic terms in v (i) are negligible, resulting in γ (i) ≃ 1, so Eqs. (194)-(197) reduce tô The energy-momentum tensors of i-th species for a combination of interacting and noncomoving perfect fluids are therefore written as where v

Tilted Frames
Let us assume u a andũ a to be two timelike, (u a u a =ũ aũ a = −1) 4-velocities of two observers O andÕ satisfying the following projection tensors (h a b u b = 0 =h a bũ b ): which represent the spatial parts of the local rest frames of the O andÕ observers. For each hypersurface-orthogonal, the projection tensor is a spatial metric in the surface. In other words, the projectors are the metrics in the subspaces of the tangent space orthogonal to the corresponding 4-velocities.
Suppose that there is a relation between u a andũ a determined by the hyperbolic angle of tilt β [100]: also constrained by the direction of tilt, either given by the directionc a of the O observer motion (i.e. projection of u a ) in theÕ local rest frame: or specified by the c a direction of theÕ motion (i.e. projection ofũ a ) in the O local rest frame: h a bũ b = sinh βc a ⇒ c a u a = 0, c a c a = 1.
We thus have: In the case, where u a andũ a are orthogonal to the homogenous surface everywhere (u a = u a ), the hyperbolic angle vanishes (β = 0) [104][105][106]. However, in the case of tilted models [100], c a andc a are uniquely expressed by (202), (203) and (204) in a way that u a is tilted with respect to the homogenous surface (u a =ũ a ), so β > 0. It is also useful to consider γ ≡ cosh β rather than β, where γ is the contraction factor for the relativistic peculiar velocity of the fluid relative to the homogeneous surface [100]: We note that sinh β = (γ 2 − 1) 1/2 . The 3-velocity v a = vc a is related to γ by [100] From Eq. (204), it follows some algebraic relations: c a u a = sinh β = −c aũ a , c ac a = cosh β.
There is a spacelike difference vector d a that holds the following relations [100]: The projection tensorh ab is also related to h ab as follows: [37]: Linearization of the above relations under the condition β ≪ 1, where the motion ofÕ with respect to O is non-relativistic, yields: The change between the frames O andÕ with a small relative velocity is referred to as a first-order change in β. For a set of timelike worldlines tangent to the 4-velocity u a , normalized such that u a u a = −1, we consider the matter field moving withũ a [99,107]: where u a is the 4-velocity with respect to observers comoving with the matter (i.e. v a u a = 0), and v a is the peculiar velocity of the matter relative to u a . In an FLRW background, as u a andũ a are chosen parallel to the canonical time direction, the peculiar velocity vanishes (v a v a ≪ 1), meaning that γ ≃ 1 to ensureũ aũ a ≃ −1. For two independent frames u a and u a , the 4-velocity fieldũ a correspond to an timelike direction and an associated projection tensor given by Eq. (200). As v a is not orthogonal toũ a , even for γ = 1, Eq. (217) implies thatũ a v a = v 2 = 0, and then Eq. (200) ensures thath ab v b = v a + v 2ũ a = v a .

Tilted Frame Transformations
For tilted O andÕ frames described by 4-velocity u a andũ a , respectively, the kinematic and dynamic quantities observed by O undergo some transformations as measured byÕ. Here, we summarize the exact forms of these transformations.
Transformation of u a to the comoving frameÕ performed by Eq. (217), have the following algebraic identities: The above relations, together with the decomposition (20), and ∇ a γ = γ 3 v b ∇ a v b , result in the following nonlinear transformation of the kinematic quantities [102] Θ =γΘ + γ(D a v a +u a v a ) + γ 3 Moreover, the dynamic quantities are nonlinearly transformed as follows [102] (also compare with [108]) Eqs. (222)-(229) to linear order reduce to: Furthermore, the transformation of the Weyl tensor relative to the comoving frameÕ yields: The kinematic quantities in the above equation are measured in the frame ofÕ, whose projection tensor and velocity are defined by (200) and (217) respectively, so the gravitoelectric and gravitomagnetic tensors are transformed to the frameÕ as follows: To linear order, we haveẼ ab ≈ E ab andH ab ≈ H ab . In a multi-fluid system, the same transformations can be employed to analyze the gravitoelectric/-magnetic fields of multiple matter species relative to the fundamental frame (see [102] for details). In multi-component models, it is essential to consider the 4-velocities of various species with respect to the fundamental frame. Each matter component has a distinct 4-velocity, which needs to be analyzed in a fundamental frame. Each of the 4-velocities leads to slightly different covariant formulations of the gravitoelectromagnetic fields. These different variations may also be regarded as partial gauge-fixings that can be solved using a covariant method. However, in an FLRW spacetime, any differences among the 4-velocities of different matter species vanish, as covariantly demonstrated by the consistent and gauge-invariant linearization around an FLRW background [27,42].

1 + 3 Tetrad Formalism
In addition to the 1 + 3 covariant formalism, there is the so-called 1 + 3 tetrad approach that also allows us to express geometric interpretations in terms of the projected vectors and PSTF tensors coupled with the dynamic and kinematic quantities. In particular, it is also useful to consider the tetrad formulations along with the covariant formalism that represents covariant relations of quantities with geometrical and/or physical meanings. In the covariant formalism, we do not have a full set of equations that demonstrates how the metric and connection are linked to each other. A tetrad description allows us to formulate a complete set of equations linking the metric and connection that are equivalent to the covariant equations. The development of the tetrad formulations was started by Pirani [52], and followed by Newman and Penrose [109], Ellis [39], Stewart and Ellis [110], and MacCallum [111,112] (see reviews by [62,[113][114][115][116][117][118]). This approach has been extensively employed by a number of authors [29,46,82,[119][120][121][122][123]. In this section, we briefly introduce the tetrad formalism and summarize its key identities, which are useful for analyzing the gravitoelectric/-magnetic tensorial fields. Rewriting the Bianchi equations using the tetrad approach puts the constraint and evolution equations of gravitoelectromagnetism on another identical footing.
To establish the tetrad method, a vector basis is chosen for an orthonormal tetrad {e a } with the timelike vector e 0 as follows: A local coordinate system is built by x i and a tetrad by {e a } in a way that [39] e a = e a i ∂/∂x i ≡ e a i ∂ i , where e a i are the functions (det |e a i | = 0) including components of the tetrad vectors e a relative to the basis ∂/∂x i and directional derivatives of the coordinate functions x i as e a i = ∂ a (x i ) [39] . The local tetrad transformation, e a ′ = Λ a a ′ e a , where Λ a a ′ is a position-dependent Lorentz matrix, along with the coordinate transformation, results in changes of the functions e a i [39]. The components e a i of ∂/∂x i are typically referred to as the basis {e a } specified by [110] As the tetrad is orthonormal, the tetrad components of the spacetime metric are written by where η ab = diag(−1, 1, 1, 1) is the Minkowski metric, so the basis unit vectors e a and e b are orthogonal to each other. As the metric tensor g ab are numerically equivalent to g ab , tetrad indices can be raised and lowered by useing g ab = η ab , and vice versa.
In the tetrad method, a set of four linearly independent 1-forms {ω a } may be chosen at each point of the spacetime manifold in such a way that the line elements can locally be written by ds 2 = η ab ω a ω b , so the vector fields {e a } are dual to the 1-form fields {ω a } such that ω a , e b = δ a b [46]. Suppose the components of the metric tensor expressed by δ b a , from Eq. (239) it follows that δ b a = e a i e i b = e a · e b , so e a i and e j b are inverse matrices [110]. Thus, we may write δ i j = e j b e b i that implies g ij = e i a e j b g ab .
Accordingly, Eqs. (238) and (240) are inverse to Eqs. (237) and (239), respectively. We may obtain components of any vectors V a or tensors T ab cd with respect to the basis ∂/∂x i or e a (see e.g. [29,48,124,125]): T ab cd =e a i e a j e c k e d l T ij kl .
The Ricci rotation coefficients are typically regarded as tetrad components of the Christoffel symbols [39] or connection components for the tetrad [29]: where Γ c ab calculates the c-component of the covariant derivative in the e b -direction of the basic vector e a .
Using the rotation coefficients, the covariant derivatives of any vectors V a or tensors T bc can be expressed in terms of tetrad components as follows [48,124]: But, for any scalars f , we shall calculate the derivative of f in the e a -direction by We know that e a (g bc ) = ∇ a g bc = 0, so the rotation coefficients are symmetric on the two indices Γ (ab)c = Γ abc + Γ bac = 0, where they are only raised and lowered.

Commutators
The Lie derivative of e b relative to e a is expressed by the [e a , e b ] commutator calculated by Following [39], the [e a , e b ] commutator is connected to a basic vector e c via commutation functions γ c ab : where the commutation functions γ c ab satisfy Eq. (247) leads to the inverse of Eq. (250), the so-called the Christoffel relation: It can be seen that the rotation coefficients Γ abc are linearly connected to the commutation functions γ abc , and vice versa. Suppose that the timelike direction of the orthonormal frame {e a } is aligned with the preferred timelike reference congruence e 0 = u (u a = δ a 0 , u a = −δ 0 a ), the commutation functions with one or two indices set to zero can be written in term of the frame components of the kinematic quantities (22) and (23) of the timelike congruence, so we have [126,127] where Ω a is the local angular velocity in the rest-frame of an observer with 4-velocity u a , andė a ≡ u b ∇ b e a is a Fermi-propagation axes relative to {e a }. For the Fermi-derivatives e a ·ė b ≡ e i a u j ∇ j e bi , we get e a ·ė b = −e b ·ė a . The spatial functions γ a bc are split into a symmetric tensor n ab = n (ab) and an antisymmetric term specified by a vector a a : 7 where n ab and a b are defined by where ε abc is the spatial permutation tensor with ε 123 = 1 = ε 123 . From Eq. (251), the rotation coefficients have the following components [39,104,123]: The first two equations contain the kinematic quantities, whereas the latter two equations include the rotation rate of the spatial frame {e a } relative to a Fermi-propagation basis and the quantities (a b and n bc ) determining the nine rotation coefficients, respectively. Considering the commutators (249) and the variables introduced in Eq. (253), we have [e a , e b ] =2ε abc ω c e 0 + 2a [a δ c b] + ε abd n dc e c .
Under the conditions, where a spatial frame vector e a is hypersurface orthogonal (see e.g. [62]), we get [122]

Tetrad Curvature
Let us consider u a as the basis vector e b (i.e. e b · u a = δ a b ). Applying Eq. (245) to the Ricci identities (86) results in the tetrad Riemann curvature given by Contracting the tetrad Riemann curvature leads to the tetrad Ricci curvature: The antisymmetry properties of the Riemann curvature (R a Substituting the tetrad forms (262) and (263) into Eq. (31), we derive the Einstein equations, along with the 16 Jacobi identities and constraint and evolution equations for E ab and H ab in terms of the tetrad derivatives (e a ), the tetrad variables (Ω a , a a , and n ab ), and the kinematic and dynamical quantities (see [29,46,82,122,123] for full details).

Correspondence between Covariant and Tetrad Formulations
Here we summarize how the 1 + 3 covariant notations correspond to the 1+3 tetrad analogues. To obtain the equivalent tetrad formulations, we can use the following conversion rules [46,128,129]: Applying the above rules to the Bianchi equations (93)-(96) with a perfect-fluid matter yields The above equations can be expressed in a symmetric normal hyperbolic form that determines their hyperbolic characteristics [123].
While the 1+3 covariant approach split spacetime into time and space using a timelike 4-velocity vector u a (u a u a = −1), the 1 + 1 + 2 semi-covariant formalism additionally decomposes space into one dimension and 2D surface with the help of a unit spacelike vector n a (i.e. n a n a = 1 and u a n a = 0). The projection tensor h ab given by Eq. (1) can be split into the spacelike vector n a and a new projection tensor N ab as follows: N ab ≡ h ab − n a n b = g ab + u a u b − n a n b . (280) The tensor N ab , which projects vectors orthogonal to n a and u a onto a 2-surface referred to as the sheets [132], has the following properties: The alternating Levi-Civita 2-tensor ε ab , is determined from the volume element of the observers' rest-spaces ε ab ≡ ε abc n c = u d η dabc n c , which possesses the following identities and contractions: The projection tensor N ab can be employed to irreducibly decompose any 3-vector V a into a scalar part V of the vector parallel to n a , and a 2-vector V a sitting on the sheet orthogonal to n a : where V and V a can be specified by Similarly, a PSTF tensor W ab can be split into a scalar W, a vector W a , and a 2-tensor W ab as follows: where W, W a , and W ab can be obtained using W =n a n b W ab , The curly brackets are utilized to denote the projection of a vector with N ab and the transverse-traceless (TT) part of a tensor. This approach allows us to decompose any object into scalars, 2-vectors laying in the sheet, and TT 2-tensors again in the sheet. Using n a and N ab , one can also define the following two new derivatives of any tensors:T where the hat (ˆ) denotes the derivative along n a in the surfaces orthogonal to u a , and δ denotes a projected derivative on the sheet, with projection on every free index, which is analogous to the spatial derivative D in the 1 + 3 covariant formalism The above derivatives can be employed to provide the following algebraic identities [133]:Ṅ where we haveṅ The spatial derivatives D of any scalars, projected 2-vectors, and TT 2-tensors, their spatial derivatives can be split into the 1 + 1 + 2 notations: By analogy with Eqs. (22) and (23), the spatial derivative D of n a orthogonal to u a is split into the following irreducible kinematic quantities: D b n a =δ b n a +n a n b , wheren a = n b D b n a is the acceleration of the sheet moving along n a , φ = δ a n a is the expansion of the sheet, ζ ab = δ {a n b} represents the distortion of the sheet calculated by the shear of n a , and ξ = − 1 2 ε ab δ a n b is the twist of the sheet described by the rotation of n a . These kinematic quantities of n a in the 1 + 1 + 2 approach can be treated in the same way as the kinematical variables of u a in the 1 + 3 method.
To formulate the 1 + 1 + 2 approach, we should also decompose the kinematical quantities defined by Eqs. (22) and (23), the Weyl fields, and the energy-momentum tensor. The acceleration is provided by (298), but for the vorticity and shear, we have the following decomposition: σ ab =Σ(n a n b − 1 2 N ab ) + 2Σ (a n b) + Σ ab .
The gravitoelectric tensor is splits into E , E a , and E ab , similarly the gravitomagnetic tensor, as follows E ab =E (n a n b − 1 2 N ab ) + 2E (a n b) + E ab , H ab =H(n a n b − 1 2 N ab ) + 2H (a n b) + H ab .
The shear scalar σ is calculated by Finally, the energy flux q a and the anisotropic stress π ab of the energy-momentum tensor (26) are split into π ab =Π(n a n b − 1 2 N ab ) + 2Π (a n b) + Π ab .
We see how the 1 + 3 variables are split to formulate their 1 + 1 + 2 counterparts, which produce new variables associated with the irreducible parts of ∇ a n b . Accordingly, the 1 + 1 + 2 Bianchi and Ricci equations can be constructed for the timelike vector u a and the spacelike vector n a in terms of all the quantities introduced above (see e.g. [131][132][133]).

Correspondence between Covariant and Semi-covariant Formulations
The 1 + 3 covariant formulations can simply be mapped into their 1 + 1 + 2 semicovariant analogues using the following rules [133]: (curlW) ab →ε c{aŴb} c + ε cd δ c W d(a n b) The above relations can be put into the covariant formulas to get their equivalences in the semi-covariant formalism.
The above linear equations describe the evolution of the Weyl fields in a perfect-fluid model according to in the 1 + 1 + 2 approach (see [133] for a generic matter field). We see that the 1 + 1 + 2 formalism offers an alternative description of the gravitoelectric/magnetic fields in terms of quantities with clear physical or geometrical meanings, namely scalar, 2-vector, and TT 2-tensors. This instrument has been used for a number of studies, such as electromagnetic perturbations on locally rotationally symmetric spacetimes [132], electromagnetic radiation made by gravitational waves interacting with a strong magnetic field of vibrating massive compact objects [131], and perturbations in rotationally symmetric spacetimes [133].

Observational Effects
Astrophysical observations are performed by detecting electromagnetic waves (light) and gravitational waves that have traveled through the past lightcone. We know from where the light originates and how it reaches us. Although astronomical observations allow us to collect information about the universe, there is a strict limit on the information due to constraints imposed by the past lightcone [138]. In addition to astrophysics, observations also enable us to collect information about the kinematic properties of large-scale structures such as the CMB anisotropies, where are described by global properties of u a e.g. the shear (see [40,128,129,139]). The propagation of electromagnetic waves and gravitational waves is affected by the spacetime metric, in other words the kinematic quantities, so the description of it in terms of the kinematic quantities will help us to better interpret the observed effects of the gravitoelectric/-magnetic field on small scales near compact objects such as around black holes, as well as large-scale structures (e.g. CMB).
Covariant handling of observations was first implemented by Kristian and Sachs [140], MacCallum and Ellis [105]. They tried to link different distance measurements to redshift based on the source's intrinsic and observed brightness. The expansion of the universe results in an energy loss of the electromagnetic waves reaching us. The energy loss of each photon is described by the redshift z, i.e. ν obsv = (1 + z) −1 ν emit , where ν emit is the emitted frequency and ν obsv is the observed frequency. Cosmological redshift is one of the important measurable quantities associated with the stretching effect of space expansion in cosmology.
The redshift of electromagnetic waves can be investigated in a more careful way via solving the Maxwell equations, ∇ a F ab = J b (where F ab = 2∇ [a A b] ) on a pseudo-Riemannian manifold for current and current-free assumptions, respectively, Suppose that the light wavelengths are small relative to the curvature, light is then taken to be tangent to null surfaces of the constant phase φ, so travels on null geodesics. If the velocity of a photon is specified by a null (geodesic) vector k a ≡ ∇ a φ, we then have ∇ [a k b] = 0, and ∇ a φ∇ a φ = k a k a = 0, k a ∇ a k b = 0.
The former is an eikonal equation [141]. The angle between the vector k a and the traveling velocity u a of an observer is associated with the relative angular frequency ω = 2πν = −u a k a of electromagnetic waves measured by the observer [142,143]. A photon traveling between two points is then redshifted based on the relative frequency difference between the two points, i.e., 1 + z = (u a k a ) emit / u b k b obsv . This redshift effect can be explicitly described by the change rate of frequency along the photon path: As k a is a null vector, it may be expressed as a linear mixture of a parallel term and a term orthogonal to u a : k a = −u b k b (u a + e a ), e a e a = 1, u a e a = 0, where e a is a normalized spacelike vector orthogonal to u a that represents the propagation direction of the electromagnetic waves with respect to u a . Substituting the kinematic quantities into the above equation yields k a ∇ a ω = −ω 2 1 3 Θ + e au a + e a e b σ ab − ε abc e a e b ω c , which results in Θ + e au a + e a e b σ ab − ε abc e a e b ω c δl where δl 2 = h ab dx a dx b = (u a k a )dv 2 = −ωdv 2 . The relative direction of electromagnetic waves is affected as their frequency changes, which can be explained by the kinematic quantities. The expansion Θ reduces the energy and frequency of photons and increases the wavelength, which causes a red shift. However, as seen in Eq. (327), the changes in frequency caused by the shear σ ab , vorticity ω a , and accelerationu a are direction-dependent. Depending on the direction of the incoming photon, relative to the direction of vorticity and acceleration, ω a andu a can make either an increase or a decrease in the frequency. Sim-ilarly, the shear σ ab also modifies the light frequency according to the difference between its stress direction and the incoming light direction [20,21,66,105,140]. Furthermore, the Ricci equations (89) and (92) imply that the gravitoelectric and gravitomagnetic fields, which are governed by the Bianchi equations (77)- (80), induce the shear and vorticity. In this way, the gravitoelectric/-magnetic fields, along with the kinematical quantities, indirectly have some observational consequences [140]. The effects of the Weyl conformal tensor occur in the geodesic deviation equation, where the distortions affect bundles of null geodesics, leading to strong gravitational lensing, multiple images and Einstein rings, as well as weak gravitational lensing and image distortions (see e.g. [144,145]). 8 These effects have profound implications for cosmology, where observations allow us to determine the large-scale structure of the universe in the past null cone [138].

Conclusions
In this paper, we have reviewed the 1 + 3 covariant formalism [18,[20][21][22][23] (in addition to alternative methods in § 8 and § 9), which can be employed to covariantly describe the locally free fields of the (traceless) Weyl conformal tensor in general relativity. In this approach, we have the projected vectors and PSTF tensors instead of the spacetime metric, together with the kinematic quantities of the fluid flow, and the dynamic quantities of the matter field. This formalism leads to the definitions of spatial divergences, curls, and distortions of the projected vectors and PSTF tensors as introduced in § 2. Moreover, we have the kinematic quantities that describe the relative motion of comoving observers [47,48], and the dynamic quantities, including the energy density and pressure in a perfect fluid, in addition to the energy flux and anisotropic pressure in an imperfect fluid. These quantities are contingent upon either a single-type of matter or a multiple-type of matter species [28,37,73,100,102] (as seen in § 7), such as implications of inhomogeneities for CMB [72] and exact (nonlinear) solutions of gravitational collapse [99].
As seen in Section 4, the Bianchi identities are split into constraint and propagation equations, which covariantly characterize the roles played by the gravitoelectric and gravitomagnetic fields. These equations includes the spatial divergences, spatial curls, and time derivatives of the gravitoelectric and gravitomagnetic tensors [18,19,29]. The gravitoelectric divergence equation can simply be interpreted as tidal forces in Newtonian theory (see § 6). The Newtonian analogy occurs due to the presence of variables in the gravitoelectric dynamics that have obvious Newtonian counterparts. In this way, the gravitoelectric field that originates from the energy density gradient acts as a general-relativistic generalization of Newtonian tidal forces [30]. In § 6.1, a Newtonian-like model was built from the gravitoelectric field using the Heckmann-Schücking approach [22,23,90] under the assumption of instantaneous gravitational interaction that produces the Poisson equation of Newtonian theory [30,88]. However, in particular, the gravitomagnetic field, which orginates from the angular momentum density [18], has no analogy in Newtonian theory. In § 6, a purely gravitoelectric (Newtonian-like) model and a purely gravitomagnetic (anti-Newtonian [31]) model were shown to be generally inconsistent without their counterpart field.
The decomposition of the Bianchi identities also generates the gravitoelectric evolution (Ė ab ), which has no analogy in Newtonian gravity, as well as the gravitomagnetic evolution (Ḣ ab ). The absence ofĖ ab from Newtonian theory is obvious since the instantaneous interaction prevents any possibility for wave solutions. We also notice a close analogy between the evolution equations of E ab and H ab and the Maxwell equations of E a and H a in expanding spacetimes [32]. As discussed in Section 5, wave solutions of E ab and H ab appear under some certain conditions, where their curls [33] and distortions [38] support gravitational waves linearly characterized in free space byË ab − D 2 E ab = 0 and H ab − D 2 H ab = 0 via the propagation equations [17,27,76,77].
The geodesic deviation equation can be used to perceive how gravitoelectric and gravitomagnetic fields affect astrophysical and cosmological observations of electromagnetic waves, such as strong and weak lensings [140]. Although the effects were covariantly explained in terms of the kinematic quantities in § 10, gravitoelectric /-magnetic fields can also induce shear and vorticity, which subsequently modify the frequency and relative direction of light. The expansion of the universe reduces the photon energy (frequency), resulting in cosmological redshift [20,21,105]. Nevertheless, the observational effects on light caused by either shear or vorticity are direction-dependent. The light frequency can be changed based on its direction relative to the shear or vorticity direction.
Following [165,166], the symmetric tracefree (STF) part of a tensor T (A ℓ ) can be obtained by (i) constructing the symmetric part: and (ii) removing all traces: B ℓn g (a 1 a 2 · · · g a 2n−1 a 2n S b 1 b 1 ···b n b n a 2n+1 ···a ℓ ), where [ℓ/2] denotes the largest integer less than or equal to ℓ/2, and B ℓn is calculated by Similarly, the STF part of a spatial tensor is determined using the spatial metric h ab instead of the spacetime metric g ab [165]. Following [128,129], the PSTF part of a tensor is denoted by which can recursively be built using a vector basis e a (see [167]; e a e a = 1, e a u a = 0): T A ℓ = e A 1×ℓ = e a 1 e a 2 e a 3 · · · e a ℓ−1 e a ℓ .
Similarly, a (ℓ + 2)-rank PSTF tensor can be generated using a ℓ-rank PSTF tensor T A ℓ and a rank-2 PSTF tensor S ab [128]: h eg h f h S gh T e f (A ℓ−2 h a ℓ−1 a ℓ h ab) .
Contraction of h (A 2×ℓ ) also leads to Following [164], the integrals of the multipoles over the direction e a are [28,128]: where dΩ = d 2 e is the volume element covered by two independent de a . From [128], we also have: where h A 1×ℓ B 1×ℓ ≡ h a 1 b 1 · · · h a ℓ b ℓ , and ∆ ℓ is given by Therefore, we have From (A26), we find (compare to Eq. A20), Consider the relation between the directions e a andẽ a at x i specified by [167] e aẽ a = cos(β) ≡ X.