1 Introduction

Perelman defined the W-entropy [1,2,3] as a functional with a non-decreasing Lyapunov-type property from which Hamilton’s equations [4,5,6] for Ricci flows can be derived following the variational procedure. The approach was elaborated upon for the geometric evolution of three dimensional (3-d) Riemannian metrics. There were obtained a number of fundamental results in geometric analysis and topology. Such directions in modern mathematics became famous after the elaborated methods allowed one to prove the Poincaré and Thorston conjectures. In this paper we show that using nonholonomic double \(3+1\) and \(2+2\) splitting in general relativity (GR) the geometric and statistical thermodynamics methods considered in Perelman’s work can be developed for theories of generalized relativistic geometric flows. We consider how such constructions can be applied in modern cosmology and astrophysics.

There are different ways for generalizing models of 3-d Ricci flow evolution for 4-d spacetimes with pseudo-Euclidean signature. For instance, there were formulated theories of stochastic/diffusion and kinetic processes with local anisotropy, fractional geometric evolution etc [7,8,9]. It is possible to construct thermo field models of Ricci flow evolution on imaginary time \(\varsigma =-\mathrm{{it}}(0\le \varsigma \le 1/\kappa T,\) where \(\kappa \) is Boltzmann’s constant and T is the temperature). In such a case, the pseudo-Riemannian spacetime is transformed into a Riemannian configuration space-like in thermal and/or finite temperature quantum field theory (see [10, 11] and the references therein). Here we recall that Perelman treated \(\tau =\varsigma ^{-1}\) as a temperature parameter and derived his W-entropy by analogy to formulas for the entropy in statistical mechanics.Footnote 1 In his work, it was not specified what type of underlying microstates and their energy should be taken in order to explain the geometric flows corresponding to certain thermodynamical and gravity models.

The (non-relativistic) Ricci flow evolution equations postulated heuristically by Hamilton can be written in the form

$$\begin{aligned} \frac{\partial g_{\grave{\imath }\grave{j}}}{\partial \tau }=-2\ R_{\grave{ \imath }\grave{j}}. \end{aligned}$$
(1)

In these formulas, \(\tau \) is an evolution real parameter and the local coordinates \(u^{\grave{\imath }}\) with indices \(\grave{\imath },\grave{j} =1,2,3 \) are defined on a real 3-d Riemannian manifold. We can consider that Eq. (1) describe a nonlinear diffusion process for geometric flow evolution of 3-d Riemannian metrics. For small deformations of a 3-d Euclidean metric \(g_{\grave{\imath }\grave{j}}\approx \delta _{\grave{\imath } \grave{j}}+\) \(h_{\grave{\imath }\grave{j}},\) with \(\delta _{\grave{\imath } \grave{j}}=\mathrm{{diag}}[1,1,1]\) and \(h_{\grave{\imath }\grave{j}}|\ll 1\), the Ricci tensor approximates the Laplace operator \(\Delta =\frac{\partial ^{2}}{ (\partial u^{1})^{2}}+\frac{\partial ^{2}}{(\partial u^{2})^{2}}+\frac{ \partial ^{2}}{(\partial u^{3})^{2}}\). We obtain a linear diffusion equation, \(R_{\grave{\imath }\grave{j}}\sim \Delta h_{\grave{\imath }\grave{j}}.\) In modified and normalized form, equations of type (1) can be proven following a corresponding variational calculation for Perelman’s W- and F-functionals. Using the W-entropy, analogous statistical mechanics and thermodynamics was formulated. The respective thermodynamic values (mean energy, entropy and fluctuation dispersion) can be considered as certain physical characteristics of flow evolution of Riemannian metrics. Summaries of the most important mathematical results and methods can be found in [12,13,14].

Geometric flow evolution models of pseudo-Riemannian metrics have not been formulated and studied in modern physical mathematics. Such ideas have not been developed and do not exist among gravitational and related relativistic thermodynamics/diffusion/kinetic theories. In quantum field theory some examples were considered in relativistic form of low dimensional geometric flow equations of type (1). That was even before mathematicians formulated in rigorous form respective directions in geometric analysis and topology which are related to the Ricci flow theory. Friedan published during 1980–1985 a series of works on nonlinear sigma models, \(\sigma \)-models, in two \(+\) epsilon dimensions; see [15,16,17]. Certain topological properties were studied of the \(\beta \) function and solutions of the fixed-point equation (the latter called the Ricci soliton equation) and further developments on renormalization of the \( O (N)\)-invariant nonlinear \(\sigma \)-models in the low-temperature regime dominated by small fluctuations around ordered states [18].

Generalized Perelman’s functionals were studied for various models of non-Riemannian geometries and (modified) gravity theories; see [19, 21, 22] and the references therein. The corresponding generalized Ricci soliton equations describe modified Einstein equations for certain classes of modified gravity theories, MGTs. Such models possess an important decoupling property of the fundamental evolution/dynamical equations and can be integrated in general form. We can study new classes of exact solutions and search for application in modern gravity and cosmology. Nevertheless, we cannot argue that certain analogs of generalized/modified Hamilton’s equations would describe in a self-consistent manner certain evolution processes if additional assumptions on geometric and physical properties of GR flows are not considered. We cannot treat any formal relativistic modification of Perelman’s functionals defined only in terms of the Levi-Civita connection for pseudo-Riemannian metrics as an entropy functional. For the geometric flow evolution of 4-d metrics with Lorentz signature and non-stationary solutions in GR, it is not possible to formulate a statistical thermodynamic interpretation like in the case of 3-d Riemannian ones.

The main goal of this paper is to study how the concept of W-entropy can be generalized in order to characterize 3-d hypersurface gravitational thermodynamic configurations and their GR evolution determined by exact solutions in Einstein gravity. The approach involves two other less established (general) relativistic theories: the relativistic statistical thermodynamics and the nonlinear diffusion theory on curved spacetimes and in gravity. Historically, the first relativistic generalizations of thermodynamics due to Planck [23] and Einstein [24] were subject to certain criticisms and modifications more than half century later [25,26,27,28,29,30,31]. Various ideas and constructions exclude each other and various debates continue even today [32,33,34,35,36,37,38]. For such models, different covariant relativistic thermodynamical and/or statistical thermodynamical values were postulated, with respective transformation laws under local Lorentz transforms. There is an explicit dependence on the fourth time-like coordinate and the main issue is how to define the concept of temperature and provide a physical interpretation. We plan to elaborate on a general thermodynamic treatment of relativistic Ricci flows using methods of relativistic kinetics and nonlinear diffusion theory in our further work. In this paper, we consider a generalization of Perelman’s W-thermodynamic model to flows of entropy and effective energy in the framework of “most simple” relativistic fluids theory and hydrodynamics with stability and causality.

The relativistic 4-d geometric flows depend, in general, on a time-like coordinate being described also by evolution on a temperature-like parameter. We argue that an appropriate redefinition of effective thermodynamic variables for the Ricci flow theory allows us to compute the entropy of gravitational fields and elaborate upon an effective statistical mechanics and thermodynamical formalism both for cosmology and black holes. We shall use the \(3+1\) decomposition formalism (see [39] and the references therein) and specify the conditions when a 3-d geometric evolution is “driven” in relativistic form by solutions of the Einstein equations in GR. Vacuum stationary solutions with Killing symmetries and horizons consist of a special class of gravitational configurations when the thermodynamic models are constructed for the entropy determined by the horizon gravity for a corresponding \(2+1+1\) splitting.

We compare our relativistic geometric flow approach to a recently proposed thermodynamical theory of gravitational fields with a measure of gravitational entropy was proposed in Ref. [40] (the so-called CET model). Those constructions are based on the square-root of the Bel–Robinson tensor when the measure is non-negative in contrast to other proposals. There were analyzed some examples, for instance, for the Schwarzschild black hole and Friedmann–Lemaître–Robertson–Worker, FLRW, cosmology. The black hole thermodynamics was derived as a particular case for stationary spacetimes. We study how the CET model can be obtained by corresponding nonholonomic parameterizations from the statistical mechanics and thermodynamic models based on the concept of W-entropy (see further developments in [19,20,21,22] and the references therein). There is not a unique way to generalize standard black hole thermodynamical constructions in order to include cosmological solutions. A geometric self-consistent variant is to address such problems using relativistic models with generalized W-entropy. In this work, we prove that there are such nonholonomic double \(3+1\) and \(2+2\) splitting when a relativistic generalization of statistical thermodynamics of geometric flows is possible for general classes of cosmological and other type solutions.

A definition of gravitational entropy which would be compatible with time depending and structure formation cosmological processes needs to be valid for very general classes of solutions with non-stationary and/or non-vacuum spacetimes. In order to prove that our approach really provides such a possibility, we shall apply the so-called anholonomic frame deformation method, AFDM (see a recent review in [41] and the references therein) for constructing generic off-diagonal exact solutions in various models of gravity theories and geometric flows with commutative and noncommutative variables [19, 21, 42]. This method involves \(2+2\) nonholonomic fibrations and a geometric techniques which allows us to integrate systems of partial differential equations (PDEs) with functional and parametric dependencies on generating and integration functions and constants depending, in general, on all spacetime coordinates and with various types of Killing and non-Killing symmetries.

The article is organized as follows: in Sect. 2, we provide an introduction into the geometry of double nonholonomic \(3+1\) and \(2+2\) fibrations of Lorentz manifolds. The main geometric and physical objects and the Einstein equations are written in nonholonomic variables adapted to a general double splitting.

Section 3 is devoted to the theory of relativistic Ricci flows with nonholonomic constraints for \(3+1\) splitting and auxiliary connections completely defined by the metric and nonlinear connection structures. The geometric evolution of the Levi-Civita configurations (with zero nonholonomically induced torsion) is considered as a special case defined by nonholonomic constraints. There are considered generalizations of Perelman’s functionals on Lorentz manifolds and derived the equations for the geometric relativistic evolution.

Section 4 is a summary of the anholonomic frame deformation method (AFDM) of constructing generic off-diagonal solutions in GR. Examples are considered of inhomogeneous and locally anisotropic cosmological solutions and black hole/ellipsoid solutions with self-consistent off-diagonal deformations and solitonic interactions.

In Sect. 5, a model of statistical thermodynamics is elaborated for gravitational 3-d hypersurface configurations and with relativistic evolution on (associated) Einstein manifolds. We show how to compute Perelman’s thermodynamical values for explicit examples of generic off-diagonal solutions in GR. It is analyzed how the relativistic W-entropy thermodynamics can be parameterized in order to model the CET thermodynamics and standard black hole physics. A general relativistic thermodynamic model for geometric flows of exact solutions in GR is elaborated following ideas and methods of relativistic hydrodynamics. Final remarks and conclusions are presented in Sect. 6.

2 Double \(2+2\) and \(3+1\) nonholonomic fibrations in GR

We provide an introduction into the geometry double \(2+2\) and \(3+1\) fibrations of Lorentz spacetime manifolds defined as follows in this section. The \(2+2\) splitting with nonholonomic deformations of the local frame and linear connection structures allowed one to decouple the Einstein equations in general form and to construct exact solutions with generic off-diagonal metrics depending on all space coordinates. Such nonholonomic \(2+2\) methods were applied for deformation and A-brane quantization of theories of gravity [19, 41, 42]. The \(3+1\) decomposition was introduced in GR with the aim to elaborate canonical approaches, for instance, to perturbative quantum gravity, relativistic thermodynamics etc. (for review of results see [39]). Working with nonholonomic distributions adapted to a conventional double \(3+1\) and \(2+2\) splitting, we can formulate an unified geometric approach to general relativistic Ricci flow and thermodynamical theories and apply the geometric methods for generating exact solutions of fundamental evolution and/or gravitational and matter fields equations.

2.1 Nonholonomic \(2+2\) splitting and nonlinear connections in GR

We consider a spacetime \(\, _{1}^{3}V\) in GR as a 4-d real smooth, i.e. \( \mathcal {C}^{\infty }\), Lorentz manifold \((V,\mathbf {g})\) determined by a pseudo-Riemannian metric \(\mathbf {g}\) of signature \((+,+,+,-).\) Footnote 2 It is assumed that it is possible to divide continuously over V each light cone of the metric \(\mathbf {g}\) into past and future paths. This means that \(\, _{1}^{3}V\) is time orientable. The tangent bundle of \(V_{1}^{3}\) is defined as the union of all tangent spaces \( T_{u}(\, _{1}^{3}V)\) is considered for all points \(u\in \, _{1}^{3}V,\) i.e. \( T(\, _{1}^{3}V):=\bigcup _{u}T_{u}(\, _{1}^{3}V).\) The typical fibers of \( T(\, _{1}^{3}V)\) are Minkowski spaces with pseudo-Euclidean signature. The dual bundle (equivalently, the cotangent bundle) of \(T(\, _{1}^{3}V)\) is denoted \(T^{*}(\, _{1}^{3}V).\) We shall write TV for the space of smooth vector fields and, respectively, \(T^{*}V\) for the space of 1-forms on a real 4-d spacetime manifold V, omitting left labels on signature if that will not result in ambiguities. In rigorous mathematical form, we can use the axiomatic approach to GR beginning in 1996 due to Ehlers et al. (the so-called EPS axioms) [43]; see further developments in [44,45,46]. In order to formulate in standard form the gravitational field equations in GR, the Einstein equations, it is used the unique metric compatible and torsionless Levi-Civita (LC) connection \(\nabla .\) The theory of GR was reformulated in various type of tetradic, spinor etc. variables and corresponding connection structures introduced with various purposes, for instance, to study solutions of the Einstein–Yang–Mills–Higgs–Dirac (EYMHD) systems [47, 48]. In so-called nonholonomic variables with an auxiliary canonical distinguished connection structure \(\widehat{\mathbf {D}}\) (see below the definition given in Eq. (9)), the EYMHD equations can be decoupled and integrated in very general forms with solutions depending, in general, on all spacetime coordinates. Fortunately, the EPS axiomatic can be generalized for spacetimes enabled with nonholonomic distributions and fibrations (for GR and various MGTs, see [49,50,51]). This provides a rigorous mathematical background for elaborating theories of general relativistic geometric evolution with nonholonomic variables on Lorentz manifolds, which is the main subject for study in this work.

A nonholonomic \(2+2\) splitting of \(_{1}^{3}V\) is determined by a nonholonomic distribution into local 2-d horizontal, h, and 2-d vertical, 2, subspaces (the local subspaces are of different signature). This defines a nonlinear connection (N-connection) structure,

$$\begin{aligned} \mathbf {N}:T\mathbf {V}=h\mathbf {V}\oplus v\mathbf {V,} \end{aligned}$$
(2)

where \(\mathbf {\oplus }\) is the Whitney sum and \(h\mathbf {V}\) and \(v\mathbf {V }\) are conventional horizontal, h, and vertical, v, subspaces. In our approach, boldface symbols are used in order to emphasize that certain spaces and/or geometric objects are for spaces endowed with N-connection structure. For simplicity, we shall write \(\mathbf {V}=(\, _{1}^{3}V,\mathbf {N} )\) for a Lorentz manifold with a hv-decomposition (2). This is an example of a nonholonomic manifold consisting, in our case, from a pseudo-Riemannian manifold and a nonholonomic (equivalently, anholonomic, or non-integrable) distribution \(\mathbf {N.}\) Footnote 3 In local form, a N-connection is stated by a set of coefficients \( N_{i}^{a}(u)\) when \(\mathbf {N}=N_{i}^{a}\mathrm{d}x^{i}\otimes \partial _{a},\) where \( \partial _{a}=\partial /\partial y^{a}.\)

There are structures of N-adapted local bases, \(\mathbf {e}_{\nu }=(\mathbf {e }_{i},e_{a}),\) and co-bases, \(\mathbf {e}^{\mu }=(e^{i},\mathbf {e}^{a}),\) when

$$\begin{aligned} \mathbf {e}_{i}= & {} \partial /\partial x^{i}-\ N_{i}^{a}(u)\partial /\partial y^{a},\ e_{a}=\partial _{a}, \end{aligned}$$
(3)
$$\begin{aligned} e^{i}= & {} \mathrm{d}x^{i},\quad \mathbf {e}^{a}=\mathrm{d}y^{a}+\ N_{i}^{a}(u)\mathrm{d}x^{i}. \end{aligned}$$
(4)

In general, N-adapted frames are nonholonomic because a frame basis \( \mathbf {e}_{\nu }=(\mathbf {e}_{i},e_{a})\) satisfies the relations

$$\begin{aligned}{}[\mathbf {e}_{\alpha },\mathbf {e}_{\beta }]=\mathbf {e}_{\alpha } \mathbf {e}_{\beta }-\mathbf {e}_{\beta }\mathbf {e}_{\alpha }=W_{\alpha \beta }^{\gamma }\mathbf {e}_{\gamma }, \end{aligned}$$
(5)

with nontrivial anholonomy coefficients \(W_{ia}^{b}=\partial _{a}N_{i}^{b}, W_{ji}^{a}=\Omega _{ij}^{a}=\mathbf {e}_{j}\left( N_{i}^{a}\right) -\mathbf {e}_{i}(N_{j}^{a}).\) We obtain holonomic (integrable) configurations if and only if \(W_{\alpha \beta }^{\gamma }=0.\) Footnote 4

On any nonholonomic spacetime \(\mathbf {V,}\) we can consider covariant derivatives determined by affine (linear) connections which are adapted to the N-connection splitting. A distinguished connection, d-connection, is a linear connection \(\mathbf {D}=(h\mathbf {D},v\mathbf {D} )\) which preserves under parallel transport the splitting (2). In general, a linear connection D is not adapted to a prescribed hv-decomposition, i.e. it is not a d-connection (we do not use a boldface symbol for non-N-adapted connections).

For any d-connection \(\mathbf {D}\) and using any d-vectors \(\mathbf {X,Y} \in T\mathbf {V,}\) we can define and compute in standard form the tensors of the d-torsion, \( \mathbf {T,}\) the nonmetricity, \(\mathbf {Q},\) and the d-curvature, \(\mathbf {R },\)

$$\begin{aligned}&\mathbf {T}(\mathbf {X,Y}):=\mathbf {D}_{\mathbf {X}}\mathbf {Y}-\mathbf {D}_{ \mathbf {Y}}\mathbf {X}-[\mathbf {X,Y}],\mathbf {Q}(\mathbf {X}):=\mathbf {D}_{ \mathbf {X}}\mathbf {g,} \\&\mathbf {R}(\mathbf {X,Y}):=\mathbf {D}_{\mathbf {X}}\mathbf {D}_{\mathbf {Y}}- \mathbf {D}_{\mathbf {Y}}\mathbf {D}_{\mathbf {X}}-\mathbf {D}_{\mathbf {[X,Y]}}. \end{aligned}$$

Any d-connection \(\mathbf {D}\) acts as an operator of covariant derivative, \( \mathbf {D}_{\mathbf {X}}\mathbf {Y}\), for a d-vector \(\mathbf {Y}\) in the direction of a d-vector \(\mathbf {X}.\) We omit boldface symbols and consider similar formulas for a linear connection which is not a N-connection.

We can compute in N-adapted form (with respect to (3) and (4)) the coefficients of any d-connection \(\mathbf {D}=\{{\varvec{\Gamma } }_{\ \alpha \beta }^{\gamma }=(L_{jk}^{i},L_{bk}^{a},C_{jc}^{i},C_{bc}^{a})\}.\) The N-adapted coefficients of torsion, nonmetricity, and curvature d-tensors are, respectively, labeled using h- and v-indices,

$$\begin{aligned}&\mathcal {T} =\{\mathbf {T}_{\ \alpha \beta }^{\gamma }=( T_{\ jk}^{i},T_{\ ja}^{i},T_{\ ji}^{a},T_{\ bi}^{a},T_{\ bc}^{a}) \},\quad \mathcal {Q}=\{ \mathbf {Q}_{\ \alpha \beta }^{\gamma }\}, \nonumber \\&\mathcal {R} =\{ {\mathbf {R}}_{\ \beta \gamma \delta }^{\alpha } =( R_{\ hjk}^{i}\mathbf {,}R_{\ bjk}^{a}\mathbf {,}R_{\ hja}^{i} \mathbf {,}R_{\ bja}^{c},R_{\ hba}^{i},R_{\ bea}^{c}) \}.\nonumber \\ \end{aligned}$$
(6)

The coefficient formulas for such values can be obtained using \({\varvec{\Gamma }}_{\ \alpha \beta }^{\gamma }\) e determined for the hv-components of \(\mathbf {D}_{\mathbf {e}_{\alpha }}\mathbf {e}_{\beta }:= \mathbf {D} _{\alpha }\mathbf {e}_{\beta }\) using \(\mathbf {X}=\mathbf {e}_{\alpha }\) and \( \mathbf {Y}=\mathbf {e}_{\beta }.\)

Any metric tensor \(\mathbf {g}\) on \(\mathbf {V}\) can be written as a d-tensor (d-metric), \(\mathbf {g}=(h\mathbf {g},v\mathbf {g}),\) i.e.

$$\begin{aligned} \mathbf {g}=g_{\alpha }(u)\mathbf {e}^{\alpha }\otimes \mathbf {e}^{\beta }=g_{i}(x)\mathrm{d}x^{i}\otimes \mathrm{d}x^{i}+g_{a}(x,y)\mathbf {e}^{a}\otimes \mathbf {e} ^{a}, \quad \quad \end{aligned}$$
(7)

for a N-adapted \(\mathbf {e}^{\mu }=(e^{i},\mathbf {e}^{a})\) (4). With respect to a dual local coordinate basis \(\mathrm{d}u^{\alpha },\) the same metric field is expressed

$$\begin{aligned}&\mathbf {g}=\underline{g}_{\alpha \beta }\mathrm{d}u^{\alpha }\otimes \mathrm{d}u^{\beta },\quad \nonumber \\&\quad \text{ where } \underline{g}_{\alpha \beta }=\left[ \begin{array}{cc} g_{ij}+N_{i}^{a}N_{j}^{b}g_{ab} &{} N_{j}^{e}g_{ae} \\ N_{i}^{e}g_{be} &{} g_{ab} \end{array} \right] . \end{aligned}$$
(8)

Using frame transforms (in general, not N-adapted), we can transform any metric into a d-metric (7) or in an off-diagonal form with N-coefficients.

For any metric field \(\mathbf {g}\) on a \(\mathbf {V}=(\, _{1}^{3}V,\mathbf {N}),\) there are two ‘preferred’ linear connection structures. The first one is the well-known Levi-Civita connection, \(\nabla ,\) and the canonical d-connection, \(\widehat{\mathbf {D}}.\) Such geometric objects are defined following the respective geometric conditions:

$$\begin{aligned}&\mathbf {g}\rightarrow \nonumber \\&\quad \left\{ \begin{array}{ll} \nabla : \nabla \mathbf {g}=0;\ ^{\nabla }\mathcal {T}=0, &{} \text{ the } \text{ Levi-Civita } \text{ connection; } \\ \widehat{\mathbf {D}}: \widehat{\mathbf {D}}\ \mathbf {g}=0;\ h\widehat{ \mathcal {T}}=0,\ v\widehat{\mathcal {T}}=0, &{} \text{ the } \text{ canonical } \text{ d-connection. } \end{array} \right. \end{aligned}$$
(9)

In these formulas, \(h\widehat{\mathcal {T}}\) and \(\ v\widehat{\mathcal {T}}\) are the respective torsions on conventional h- and v-subspaces. We note that there are non-zero torsion components, \(hv\widehat{\mathcal {T}},\) with non-zero mixed indices with respect to a N-adapted basis (3) and/or (4). Nevertheless, this torsion field \(\widehat{\mathcal {T}}\) is completely defined by the metric field following the parameterization (8) with \((h\mathbf {g},v\mathbf {g;N}).\)

All geometric constructions on \(\mathbf {V,}\) can be performed equivalently using \(\nabla \) and/or \(\widehat{\mathbf {D}}\) and are related via the canonical distorting relation

$$\begin{aligned} \widehat{\mathbf {D}}\mathbf {[g,N]}=\nabla [\mathbf {g}]+\widehat{ \mathbf {Z}}[\widehat{\mathcal {T}}(\mathbf {g,N)}]. \end{aligned}$$
(10)

By squared brackets [...],  we state a functional dependence when both linear connections \(\nabla \) and \(\widehat{\mathbf {D}}\) and the distorting tensor \(\widehat{\mathbf {Z}}\) are uniquely determined by the data \((\mathbf { g,N)}\) and an algebraic combination of the coefficients of the torsion \(\widehat{ \mathcal {T}}(\mathbf {g,N).}\)

The Ricci tensors of \(\widehat{\mathbf {D}}\) and \(\nabla \) are defined and computed in the standard way and denoted, respectively, by \(\ \widehat{ \mathcal {R}}ic=\{\widehat{\mathbf {R}}_{\ \beta \gamma }:=\widehat{\mathbf {R}} _{\ \alpha \beta \gamma }^{\gamma }\}\) and \(Ric=\{R_{\ \beta \gamma }:=R_{\ \alpha \beta \gamma }^{\gamma }\}.\) The N-adapted coefficients for \( \widehat{\mathbf {D}}\) and the corresponding torsion, \(\widehat{\mathbf {T}}_{\ \alpha \beta }^{\gamma },\) Ricci d-tensor, \(\widehat{\mathbf {R}}_{\ \beta \gamma },\) and Einstein d-tensor, \(\widehat{\mathbf {E}}_{\ \beta \gamma },\) are computed in [41, 42]. Any (pseudo) Riemannian geometry can be equivalently described by both geometric data \(( \mathbf {g,}\varvec{\nabla } ) \) and \((\mathbf {g,N,}\widehat{\mathbf {D}}),\) when the canonical distortion relations \(\widehat{\mathcal {R}}=\ ^{\nabla }\mathcal {R+}\ ^{\nabla }\mathcal {Z}\) and \(\widehat{\mathcal {R}}ic=Ric+\widehat{\mathcal {Z}} ic,\) with respective distortion d-tensors \(\ ^{\nabla }\mathcal {Z}\) and \( \widehat{\mathcal {Z}}ic,\) are computed for the canonical distortion relations \(\widehat{\mathbf {D}}=\nabla +\widehat{\mathbf {Z}}.\)

By the unique distortion relations (computed by introducing (10) into (6) and re-grouping the terms with \(\nabla \) and \(\widehat{\mathbf {D} }),\) we can relate, for instance, \(\widehat{\mathbf {R}}_{\ \beta \gamma }\) to \(R_{\ \beta \gamma },\)

$$\begin{aligned} \widehat{\mathbf {R}}_{\ \beta \gamma }\mathbf {[g,N]}=R_{\ \beta \gamma } \mathbf {[g,N]}+\widehat{\mathbf {Z}}_{\ \beta \gamma }\mathbf {[g,N].} \end{aligned}$$

We note that the Ricci d-tensor \(\widehat{\mathcal {R}}ic\) is not symmetric, \( \widehat{\mathbf {R}}_{\alpha \beta }\ne \widehat{\mathbf {R}}_{\beta \alpha },\) being characterized by four subsets of h-v N-adapted coefficients,

$$\begin{aligned}&\widehat{\mathbf {R}}_{\alpha \beta }=\{\widehat{R}_{ij}:=\widehat{R}_{\ ijk}^{k},\ \widehat{R}_{ia}:=-\widehat{R}_{\ ika}^{k},\nonumber \\&\quad \widehat{R}_{ai}:= \widehat{R}_{\ aib}^{b},\ \widehat{R}_{ab}:=\widehat{R}_{\ abc}^{c}\}. \end{aligned}$$
(11)

It is possible to compute the scalar of canonical d-curvature, \(\ \widehat{ \mathbf {R}}:=\mathbf {g}^{\alpha \beta }\widehat{\mathbf {R}}_{\alpha \beta }=g^{ij}\widehat{R}_{ij}+g^{ab}\widehat{R}_{ab}.\) This geometric object is different from the LC-scalar curvature, \(\ R:=\mathbf {g}^{\alpha \beta }R_{\alpha \beta }.\)

The Einstein equations in GR are written in standard form,

$$\begin{aligned} R_{\alpha \beta }-\frac{1}{2}g_{\alpha \beta }R=\varkappa \ ^{m}T_{\alpha \beta }, \end{aligned}$$
(12)

using the Ricci tensor \(R_{\alpha \beta }\) and the scalar R; they are taken for the Levi-Civita connection \(\nabla \) of \(g_{\alpha \beta }.\) In these formulas, \(\ ^{m}T_{\alpha \beta }\) is the energy-momentum tensor of matter fields \(\ ^{A}\varphi \) determined by a general Lagrangian \(\ ^{m}\mathcal {L} (\mathbf {g,}\nabla ,\ \ ^{A}\varphi )\) and \(\varkappa \) is the gravitational coupling constant for GR.Footnote 5 The gravitational field equations in GR can be rewritten equivalently using the canonical d-connection [41, 42],

$$\begin{aligned} \widehat{\mathbf {R}}_{\alpha \beta }= & {} \widehat{\varvec{\Upsilon }}_{\alpha \beta }, \end{aligned}$$
(13)
$$\begin{aligned} \widehat{\mathbf {T}}_{\ \alpha \beta }^{\gamma }= & {} 0. \end{aligned}$$
(14)

In these formulas, the effective matter fields source \(\varvec{\Upsilon } _{\mu \nu }\) is constructed via a N-adapted variational calculus with respect to (4) for \(\ ^{m}\mathcal {L}(\mathbf {g,}\widehat{\mathbf {D }},\ \ ^{A}\varphi )\) in such a form that

$$\begin{aligned} \widehat{\varvec{\Upsilon }}_{\mu \nu }=\varkappa \left( \ ^{m}\widehat{\mathbf {T}} _{\mu \nu }-\frac{1}{2}\mathbf {g}_{\mu \nu }\ ^{m}\widehat{\mathbf {T}} \right) \!\rightarrow \! \varkappa \left( \ ^{m}T_{\mu \nu }-\frac{1}{2}\mathbf {g}_{\mu \nu }\ ^{m}T\right) \end{aligned}$$

for [coefficients of \(\widehat{\mathbf {D}}]\) \(\rightarrow \) [coefficients of \(\nabla \)] even, in general, \(\widehat{\mathbf {D}}\ne \nabla .\) In these formulas, \(\ ^{m}\widehat{\mathbf {T}}=\mathbf {g}^{\mu \nu }\ ^{m}\widehat{ \mathbf {T}}_{\mu \nu }\) for

$$\begin{aligned} \ ^{m}\widehat{\mathbf {T}}_{\alpha \beta }:=-\frac{2}{\sqrt{|\mathbf {g}_{\mu \nu }|}}\frac{\delta \left( \sqrt{|\mathbf {g}_{\mu \nu }|}\ \ ^{m}\mathcal {L}\right) }{ \delta \mathbf {g}^{\alpha \beta }}. \end{aligned}$$
(15)

The canonical d-connection \(\widehat{\mathbf {D}}\) has a very important role to play in our approach. With respect to N-adapted frames of reference, it allows one to decouple in general form the gravitational and matter field equations in the form (13) with (15).Footnote 6 We can integrate nonholonomic deformations of the Einstein equations in very general form and construct exact solutions parameterized by generic off-diagonal metrics depending on all spacetime coordinates via the respective classes of generating functions and integration functions and constants. Having constructed certain general classes of solutions, we can impose at the end the LC-conditions (14) and extract LC-configurations \(\widehat{\mathbf {D}}_{\mid \widehat{\mathcal {T}}=0}=\nabla .\) This allows one, for instance, to construct new classes of generic off-diagonal solutions of (12) in GR and various MGTs. We note that to find nontrivial off-diagonal solutions is important to impose the LC-conditions (14) after a class of solutions of (13) for \(\widehat{\mathbf {D}}\) are constructed in general form. If we work only with \(\nabla ,\) we are not able to decouple the Einstein equations in general form.

2.2 \(3+1\) Decompositions adapted to nonholonomic \(2+2\) splitting

We foliate a 4-d Lorentzian nonholonomic manifold \(\mathbf {V}=(\, _{1}^{3}V, \mathbf {N})\) enabled with a pseudo-Riemannian metric \(\mathbf {g} =\{g_{\alpha \beta }\}\) of signature \((+++-)\) into a family of non-intersecting space-like 3-d hypersurfaces \(\Xi _{t}\) parameterized by a scalar field, i.e. the “time function”, \(t(u^{\alpha }),\) as described as follows. Such a \(3+1\) spacetime decomposition is necessary for elaborating various thermodynamic and flow models when a conventional splitting into time- and space-like coordinates is important for definition of physical important values (like entropy, effective energy etc.) and fundamental geometric evolution equations. We have to generalize the well-known geometric \(3+1\) formalism [39] to the case of spacetimes enabled with the nontrivial N-connection structure [41, 42].

A hypersurface \(\Xi \subset \mathbf {V}\) is considered as an one-to-one image of a 3-d manifold \({\,}_{\shortmid }\Xi .\) This image is given by an embedding \(\Xi =\zeta (\, _{\shortmid }\Xi )\) constructed as an homeomorphism with both continuous maps \(\zeta \) and \(\zeta ^{-1}.\) This guarantees that \( \Xi \) does not intersect itself. We shall use a left “up” or “low” label by a vertical bar “\(\, _{\shortmid }\ \)” in order to emphasize that a manifold is 3-d, or certain geometric objects refer to 3-d manifolds/hypersurfaces. Locally, a hypersurface is considered as the set of points for which a scalar field t on \(\mathbf {V}\) is constant, for instance, i.e. \(t(p)=0,\forall \) \(p\in \Xi .\) We assume that t spans \(\mathbb {R}\) and \( \Xi \) is a connected submanifold of \(\mathbf {V}\) with topology \(\mathbb {R} ^{3}.\) The local coordinates for a \(3+1\) splitting are labeled \(u^{\alpha }=(x^{\grave{\imath }},t),\) where the indices \(\alpha ,\beta ,\ldots =1,2,3,4\) and \( \grave{\imath },\grave{j},\ldots =1,2,3.\) In brief, we shall write \(u=(\breve{u} ,t).\) The mapping \(\zeta \) “carries along” curves/vectors in \(\, _{\shortmid }\Xi \) to curves/vectors in \(\mathbf {V,}\) for \(\zeta :(x^{\grave{\imath } })\longrightarrow (x^{\grave{\imath }},0).\) This defines respective local bases \(\partial _{\grave{\imath }}:=\partial /\partial x^{\grave{\imath }}\in T {(\, _{\shortmid }\varvec{\Xi } )}\) and \(\partial _{\alpha }:=\partial /\partial u^{\alpha }\in T\mathbf {V.}\) Correspondingly, the coefficients of 3-vectors and 4-vectors are expressed \(\, _{\shortmid }\mathbf {a}=a^{\grave{ \imath }}\partial _{\grave{\imath }}\) and \(\,\mathbf {a}=a^{\alpha }\partial _{\alpha }\) (for convenience, we shall use also capital letters, for instance, \(\, _{\shortmid }\mathbf {A}=A^{\grave{\imath }}\partial _{ \grave{\imath }}\) and \(\mathbf {\ A}=A^{\alpha }\partial _{\alpha }).\) For dual forms to vectors, 1-forms, we use the respective dual bases \(\mathbf {d}x^{ \grave{\imath }}\in T^{*}{(\, _{\shortmid }\varvec{\Xi } )}\) and \( \mathrm{d}u^{\alpha }\in T^{*}\mathbf {V.}\) We shall write for 1-forms \(\, _{\shortmid }{\tilde{\mathbf {A}}}=A_{\grave{\imath }}\mathbf {d}x^{\grave{\imath } }\) and \({\tilde{\mathbf {A}}}=A_{\alpha }\mathrm{d}u^{\alpha }\) and omit the left/up label by a tilde \({\sim } \) (writing \(\, _{\shortmid }\mathbf {A}\) and \( \mathbf {A)}\) if that will not result in ambiguities.

Using the push-forward mapping,

$$\begin{aligned} \begin{array}{cccc} \zeta _{*}:&T_{u}\, _{\shortmid }\Xi \longrightarrow T_{p} \mathbf {V;}&\,&\, _{\shortmid }\mathbf {v}=(v^{\grave{\imath } })\longrightarrow \zeta _{*}\, _{\shortmid }\mathbf {v}=(v^{ \grave{\imath }},0), \end{array} \end{aligned}$$

we can transport geometric objects from \(\, _{\shortmid }\Xi \) to \(\Xi ,\) and inversely. In dual form, the pull-back mapping acts as

$$\begin{aligned} \begin{array}{ccccccc} \zeta ^{*}: &{} T_{p}^{*}\mathbf {V} &{} \longrightarrow &{} T_{u}^{*} \, _{\shortmid }\Xi , &{} &{} &{} \\ &{} \, _{\shortmid }{\tilde{\mathbf {A}}} &{} \longrightarrow &{} \zeta ^{*}\, _{\shortmid }\tilde{\mathbf {A}:} &{} T_{u}\, _{\shortmid }\Xi \longrightarrow \mathbb {R}; &{} &{} \, _{\shortmid }\mathbf {A}\longrightarrow \langle \, _{\shortmid }{\tilde{\mathbf {A}}},\zeta _{*}\, _{\shortmid }\mathbf {A}\rangle , \end{array} \end{aligned}$$

for \(\langle \cdots \rangle \) denoting the scalar product and \(T_{p}^{*}\mathbf {V} \ni \) \( {\tilde{\mathbf {A}}}=(A_{\grave{\imath }},A_{4})\longrightarrow \zeta ^{*} {\tilde{\mathbf {A}}}=(A_{i})\in T_{u}^{*}\, _{\shortmid }\Xi .\) In this work, we identify \(\, _{\shortmid }\Xi \) and \(\Xi =\zeta (\, _{\shortmid }\Xi )\) and write simply a d-vector \(\, \mathbf {v}\) instead of \(\zeta _{*}(\, _{\shortmid }\mathbf {v).}\) For holonomic configurations, the same maps and objects are labeled in non-boldface form.

2.2.1 Induced N-adapted 3-d hypersurface metrics

We define the first fundamental form (the induced 3-metric) on \(\Xi :\)

$$\begin{aligned} \begin{array}{ccc} \mathbf {q}:=\zeta ^{*}\ \mathbf {g}, &{} &{} \text{ i.e. } \mathbf {q}_{\grave{ \imath }\grave{j}}:=\mathbf {g}_{\grave{\imath }\grave{j}} \\ \forall (\ ^{1}\mathbf {a,}\ ^{2}\mathbf {a})\in T_{u}\,\Xi \times T_{u}\,\Xi , &{} &{} \ ^{1}{{\mathbf {a}}{\cdot } }\ ^{2}\mathbf {a=g}(\ ^{1} \mathbf {a,}\ ^{2}\mathbf {a})=\mathbf {q}(\ ^{1}\mathbf {a,}\ ^{2}\mathbf {a}). \end{array} \end{aligned}$$

The 4-d metric \(\mathbf {g}\) is constrained to be a solution of the Einstein equations (13) written in nonholonomic variables. The hypersurfaces are classified following the types of induced 3-metric (for a nontrivial N-connection, this can be represented as an induced d-metric):

$$\begin{aligned} \left\{ \! \begin{array}{ll} \text{ space-like } \Xi , &{} \mathbf {q} \; \text{ is } \text{ positive } \text{ definite } \text{ with } \text{ signature } (+,+,+); \\ \text{ time-like } \Xi , &{} \mathbf {q} \; \text{ is } \text{ Lorentzian } \text{ with } \text{ signature } (+,+,-); \\ \text{ null } \Xi , &{} \mathbf {q} \; \text{ is } \text{ degenerate } \text{ with } \text{ signature } (+,+,0). \end{array} \right. \end{aligned}$$

In this article, we shall work with continuous sets of space-like hypersurfaces \(\Xi _{t},t\in \mathbb {R},\) covering some finite, or infinite, regions on \(\mathbf {V}\). For simplicity, we use only space-like hypersurfaces \(\Xi \) (which can be closed and compact if necessary) endowed with Riemannian 3-metric \(\mathbf {q}\) if other conditions will be not stated.

We introduce the concept of a unit normal d-vector, \(\mathbf {n}\), to a \(\Xi \) which is constructed following such a procedure. Such d-vectors can be used for various models of geometric flow evolution and thermodynamic models. Let us consider a scalar field \(t(u^{\alpha })\) on an open region \(U\subset \mathbf {V}\) such as the level surface is identified to \(\Xi .\) We construct in N-adapted form the gradient 1-form \(\mathbf {d}t\) and its dual d-vector \( \overrightarrow{\mathbf {e}}t=\{\mathbf {e}^{\mu }t=\mathbf {g}^{\mu \nu } \mathbf {e}_{\nu }t=\mathbf {g}^{\mu \nu }(\mathbf {d}t)_{\nu }\};\) see the operators (3) and (4). For any d-vector \(\mathbf {v}\) which is tangent to \(\Xi ,\) the conditions \(\langle \mathbf {d}t,\mathbf {v}\rangle =0\) and \( \overrightarrow{\mathbf {e}}t\) allow one to define the unique direction normal to a not null \(\Xi .\) Normalizing such a d-vector, we define

$$\begin{aligned}&\mathbf {n:=\pm }\overrightarrow{\mathbf {e}}t/\sqrt{|\overrightarrow{\mathbf {e }}t|},\nonumber \\&\quad \text{ where } \left\{ \begin{array}{ll} \mathbf {n}{\cdot } \mathbf {n}=-1, &{} \text{ for } \text{ space-like } \Xi ; \\ \mathbf {n}{\cdot } \mathbf {n}=1, &{} \text{ for } \text{ time-like } \Xi . \end{array} \right. \end{aligned}$$
(16)

The unit normal vector to supersurfaces, \(\mathbf {n}_{\alpha }\propto \partial _{\alpha }t,\) when \(\partial _{\alpha }:=\partial /\partial u^{\alpha },\) can be constructed for a future-directed time-like vector field. We can use t as a parameter for a congruence of curves \(\chi (t)\subset \, _{1}^{3}V\) intersecting \(\Sigma _{t},\) when the vector \(\mathbf {t}^{\alpha }:=\mathrm{d}u^{\alpha }/\mathrm{d}t\) is tangent to the curves and \(\mathbf {t}^{\alpha }\partial _{\alpha }t=1.\) For any system of coordinates \(u^{\alpha }=u^{\alpha }(x^{\grave{\imath }},t),\) there are defined the vector \(\mathbf {t }^{\alpha }:=(\partial u^{\alpha }/\partial t)_{x^{\grave{\imath }}}\) and the (tangent) vectors \(e_{~\grave{\imath }}^{\alpha }:=(\partial u^{\alpha }/\partial x^{\grave{\imath }})\); and the Lie derivative along \(\mathbf {t} ^{\alpha }\) results in \(\pounds _{t}e_{~^{\grave{\imath }}}^{\alpha }=0.\)

Any 2+2 splitting \(\mathbf {N}:T\mathbf {V}=h\mathbf {V}\oplus v\mathbf {V,}\) (2) induces a N-connection structure on \(\Xi ,\) \(\, _{\shortmid } \mathbf {N}:T\, _{\shortmid }\Xi =h\, _{\shortmid }\Xi \oplus v\, _{\shortmid }\Xi \mathbf {.}\) As a result, any induced 3-metric tensor \(\mathbf {q}\) can be written in N-adapted frames as a d-tensor (d-metric) in the form

$$\begin{aligned} \mathbf {q}= & {} (h\mathbf {q},v\mathbf {q})\nonumber \\= & {} q_{\grave{\imath }}(u)\mathbf {e}^{ \grave{\imath }}\otimes \mathbf {e}^{\grave{\imath }}\nonumber \\= & {} q_{i}(x^{k})\mathrm{d}x^{i}\otimes \mathrm{d}x^{i}+q_{3}(x^{k},y^{3})\, _{\shortmid }\mathbf {e}^{3}\otimes \, _{\shortmid } \mathbf {e}^{3}, \end{aligned}$$
(17)
$$\begin{aligned} \, _{\shortmid }\mathbf {e}^{3}= & {} \mathrm{d}u^{3}+\, _{\shortmid }N_{i}^{3}(u)\mathrm{d}x^{i}, \end{aligned}$$
(18)

where \(\, _{\shortmid }N_{i}^{3}(u)\) can be identified with \(N_{i}^{3}(u)\) choosing common frame and coordinate systems for \(\Xi \subset \mathbf {V}.\) We can naturally embed such a metric into a d-metric (7) reparameterized in a form adapted both to the \(2+2\) and the \(3+1\) nonholonomic splitting,

$$\begin{aligned} \begin{aligned} \mathbf {g}&=(h\mathbf {g},v\mathbf {g})=\breve{g}_{\grave{\imath }\grave{j}} \mathbf {e}^{\grave{\imath }}\otimes \mathbf {e}^{\grave{j}}+g_{4}\mathbf {e} ^{4}\otimes \mathbf {e}^{4}\\&=q_{\grave{\imath }}(u)\mathbf {e}^{\grave{\imath } }\otimes \mathbf {e}^{\grave{\imath }}-\breve{N}^{2}\mathbf {e}^{4}\otimes \mathbf {e}^{4}, \\ \mathbf {e}^{3}&=\, _{\shortmid }\mathbf {e}^{3}=\mathrm{d}u^{3}+\, _{\shortmid }N_{i}^{3}(u)\mathrm{d}x^{i},\\ \mathbf {e}^{4}&=\delta t=\mathrm{d}t+\, N_{i}^{4}(u)\mathrm{d}x^{i}. \end{aligned} \end{aligned}$$
(19)

Let us explain this construction. The lapse function \(\breve{N} (u)>0\) is defined as a positive scalar field which ensues from the fact that the unit d-vector \(\mathbf {n}\) is a unit one; see (16). We use an “inverse hat” in order to distinguish such a symbol N is used traditionally in the literature on GR [39] but, in another turn, the symbol \(N_{i}^{a}\) is used traditionally for the N-connection. We write

$$\begin{aligned} \mathbf {n}:=-\breve{N}\overrightarrow{\mathbf {e}}t\quad \text{ and/or } \quad \underline{ \mathbf {n}}:=-\breve{N}\mathbf {d}t, \end{aligned}$$

with \(\breve{N}:=1/\sqrt{|\overrightarrow{\mathbf {e}}t\cdot \overrightarrow{ \mathbf {e}}t|}=1/\sqrt{|<\mathbf {d}t\cdot \overrightarrow{\mathbf {e}}t>|}.\) For geometric constructions, it is convenient to use also the normal evolution d-vector

$$\begin{aligned} \mathbf {m}:=\breve{N}\mathbf {n} \end{aligned}$$
(20)

subject to the condition \(\mathbf {m}\cdot \mathbf {m}=-\breve{N}^{2}.\) This is justified by the property of the Lie N-adapted derivative that \(\mathcal { L}_{\mathbf {m}}\mathbf {a}\in (h\Xi \oplus v\Xi )_{t},\forall \mathbf {a}\in T\Xi _{t}.\) For a \(3+1\) spacetime splitting, we consider also the shift functions (a 3-vector \(\breve{N}^{\grave{\imath }}(u),\) or a d-vector \({\breve{\mathbf {N}}}^{\grave{\imath }}(u)\)). It is useful to define the unit normal \(\breve{n}^{\alpha }\) to the hypersurfaces when \(\breve{n} _{\alpha }=-\breve{N}\partial _{\alpha }t\) and \(\breve{n}_{\alpha }e_{~^{ \breve{\imath }}}^{\alpha }=0.\) In N-adapted form, we can consider that \( \breve{n}_{\alpha }\) is a normalized version of \(\mathbf {n}\) used in (20). This allows us to consider the decompositions

$$\begin{aligned} \mathbf {t}^{\alpha }=\breve{N}^{\grave{\imath }}e_{~^{\grave{\imath } }}^{\alpha }+\breve{N}\breve{n}^{\alpha } \end{aligned}$$

and

$$\begin{aligned} \mathrm{d}u^{\alpha }=e_{~^{\grave{\imath }}}^{\alpha }\mathrm{d}x^{\grave{\imath }}+\mathbf {t} ^{\alpha }\mathrm{d}t=(\mathrm{d}x^{\grave{\imath }}+\breve{N}^{\grave{\imath }}\mathrm{d}t)e_{~^{\grave{ \imath }}}^{\alpha }+(\breve{N}\mathrm{d}t)\breve{n}^{\alpha }. \end{aligned}$$

Using the quadratic line element \(\mathrm{d}s^{2}=g_{\alpha \beta }\mathrm{d}u^{\alpha }\mathrm{d}u^{\beta }\) of a metric tensor \(\mathbf {g,}\) we can choose such frame transforms when \(\breve{g}_{\grave{\imath }\grave{j}}=q_{\grave{\imath }\grave{ j}}=\) \(g_{\alpha \beta }e_{~\grave{\imath }}^{\alpha }e_{~\grave{j}}^{\beta }\) is the induced metric on \(\Xi _{t}.\) For the determinants of 4-d and 3-d metrics parameterized in the above-mentioned form, we compute \(\sqrt{|g|}=\breve{ N}\sqrt{|\breve{g}|}=\breve{N}\sqrt{|q|}.\) Using the coordinates \((x^{\grave{ \imath }},t),\) the time partial derivatives are \(\pounds _{t}q=\partial _{t}q=q^{*}\) and the spatial derivatives are \(q_{,^{\grave{\imath } }}:=e_{~^{^{\grave{\imath }}}}^{\alpha }\) \(q_{,\alpha }.\)

2.2.2 Induced 3-d hypersurface ‘preferred’ linear connections

There are two induced linear connections completely determined by an induced 3-d hypersurface metric \(\mathbf {q,}\)

$$\begin{aligned}&\mathbf {q}\rightarrow \nonumber \\&\quad \left\{ \begin{array}{ll} \, _{\shortmid }\nabla : \quad \, _{\shortmid }\nabla \mathbf {q}=0;\ \, _{\shortmid }^{\nabla }\mathcal {T}=0, &{} \text{ LC-connection } \\ \, _{\shortmid }\widehat{\mathbf {D}}: \quad \, _{\shortmid }\widehat{\mathbf {D}} \ \mathbf {q}=0;\ h\, _{\shortmid }\widehat{\mathcal {T}}=0,\ v\, _{\shortmid } \widehat{\mathcal {T}}=0. &{} \text{ canonical } \text{ d-connection. } \end{array} \right. \!\!\nonumber \\ \end{aligned}$$
(21)

Such formulas are induced from the 4-d similar ones; see (9). Both linear connections, \(\, _{\shortmid }\nabla \) and \(\, _{\shortmid } \widehat{\mathbf {D}},\) are involved in a distortion relation,

$$\begin{aligned} \, _{\shortmid }\widehat{\mathbf {D}}\mathbf {[q,\, _{\shortmid }N]}=\, _{\shortmid }\nabla [\mathbf {q}]+\, _{\shortmid }\widehat{\mathbf {Z}}[\, _{\shortmid }\widehat{ \mathcal {T}}(q\mathbf {,\, _{\shortmid }N)}], \end{aligned}$$
(22)

induced by (10).

There are two classes of trivial or nontrivial intrinsic torsions, nonmetricity and curvature fields for any data \((\Xi ,\mathbf {q,\, _{\shortmid }N}),\) defined by corresponding hypersurface linear connections when for any \(\mathbf {a,b}\in T\Xi \)

$$\begin{aligned} \, _{\shortmid }T(\mathbf {a,b}):= & {} \, _{\shortmid }\nabla _{\mathbf {a}} \mathbf {Y}-\, _{\shortmid }\nabla _{\mathbf {Y}}\mathbf {a}-[\mathbf {a,Y}]=0, \, _{\shortmid }Q(\mathbf {a}):=\, _{\shortmid }\nabla _{\mathbf {a}} \mathbf {q=0,} \\ \, _{\shortmid }R(\mathbf {a,b}):= & {} \, _{\shortmid }\nabla _{\mathbf {a}}\, _{\shortmid }\nabla _{\mathbf {b}}-\, _{\shortmid }\nabla _{\mathbf {b}}\, _{\shortmid }\nabla _{\mathbf {a}}-\, _{\shortmid }\nabla _{\mathbf {[a,b]}}, \end{aligned}$$

and

$$\begin{aligned} \, _{\shortmid }\widehat{\mathbf {T}}(\mathbf {a,b}):= & {} \, _{\shortmid }\mathbf {D}_{\mathbf {a}}\mathbf {b}-\, _{\shortmid }\mathbf {D}_{\mathbf {Y}} \mathbf {a}-[\mathbf {a,b}]=0,\nonumber \\ _{\shortmid }\mathbf {Q}(\mathbf {a}):= & {} \, _{\shortmid }\mathbf {D}_{\mathbf {a}}\mathbf {q=0,} \\ \, _{\shortmid }\widehat{\mathbf {R}}(\mathbf {a,b}):= & {} \, _{\shortmid }\widehat{\mathbf {D}}_{\mathbf {a}}\, _{\shortmid }\widehat{\mathbf {D}}_{ \mathbf {b}}-\, _{\shortmid }\widehat{\mathbf {D}}_{\mathbf {b}}\, _{\shortmid } \widehat{\mathbf {D}}_{\mathbf {a}}-\, _{\shortmid }\widehat{\mathbf {D}}_{ \mathbf {[a,b]}}. \end{aligned}$$

We can compute the N-adapted coefficient formulas for nonholonomically induced torsion structure \(\mathbf {\, _{\shortmid }}\widehat{\mathbf { T}}=\{\, _{\shortmid }\widehat{\mathbf {T}}_{\ \grave{j}\grave{k}}^{\grave{ \imath }}\},\) determined by \(\, _{\shortmid }\mathbf {D,}\) and for the Riemannian tensors \( \, _{\shortmid }R={\{ _{\shortmid }}R_{\, \grave{j}\grave{k}\grave{l}}^{\grave{\imath }}\}\) and \( {\, _{\shortmid }}\widehat{\mathbf {R}}={\{\, _{\shortmid }}\widehat{\mathbf {R}}_{\ \grave{j}\grave{k}\grave{l}}^{\grave{ \imath }}\},\) determined respectively, by \(\, _{\shortmid }\nabla \) and \(\, _{\shortmid }\mathbf {D.}\) Using 3-d variants of coefficient formulas, we can compute the N-adapted coefficients of the Ricci d-tensor, \(\, _{\shortmid }\widehat{\mathbf {R}}_{\ \grave{j}\grave{k}},\) and the Einstein d-tensor, \(\, _{\shortmid }\widehat{\mathbf {E}}_{\ \grave{j}\grave{k}}.\) Contracting indices, we obtain the Gaussian curvature, \(\, _{\shortmid }R=q^{ \grave{j}\grave{k}}\, _{\shortmid }R_{\grave{j}\grave{k}}\), and the Gaussian canonical curvature, \(\, _{\shortmid }^{s}R=q^{\grave{j}\grave{k}}\, _{\shortmid }\widehat{\mathbf {R}}_{\grave{j}\grave{k}},\) of \((\Xi ,\mathbf { q,\, _{\shortmid }N}).\) Such geometric objects do not depend on the type of embedding of the nonholonomic manifold \((\Xi ,\mathbf {q,\, _{\shortmid }N})\) in \((\mathbf {V},\mathbf {g,N}).\)

There are other types of curvatures which describe the (non) holonomic bending of \(\Xi ,\) i.e. dependent on the embedding. Such geometric entities are considered for any type of \(3+1\) splitting and constructed using N-adapted Weingarten maps (shape operator). In our approach, these are N-adapted endomorphisms,

$$\begin{aligned}&\chi : T_{u}\Xi \rightarrow T_{u}\Xi ; \quad \text{ and } \quad \widehat{\chi } : \quad T\, _{\shortmid }\Xi \rightarrow T\, _{\shortmid }\Xi =h\, _{\shortmid }\Xi \oplus v\, _{\shortmid }\Xi \\&\qquad \quad \;\;\mathbf {a} \rightarrow \, _{\shortmid }\nabla _{\mathbf {a}} \mathbf {n;} \qquad \qquad \qquad \qquad \quad \mathbf {a} \rightarrow \, _{\shortmid } \widehat{\mathbf {D}}_{\mathbf {a}}\mathbf {n}=h\, _{\shortmid } \widehat{\mathbf {D}}_{\mathbf {a}}\mathbf {n}\oplus v\, _{\shortmid } \widehat{\mathbf {D}}_{\mathbf {a}}\mathbf {n.} \end{aligned}$$

Such maps are self-adjoint with respect to the induced 3-metric \(\mathbf {q, }\) i.e. for any \(\mathbf {a,b}\in T_{p}\Xi \times T_{p}\Xi ,\)

$$\begin{aligned} {\mathbf {a}}{\cdot }\chi \mathbf {(b)}=\chi ({\mathbf {a}}){\cdot } {\mathbf {b}}\quad \text{ and } \quad {\mathbf {a} }{\cdot }\widehat{\chi }\mathbf {(b)=}\widehat{\chi }(\mathbf {a){\cdot } b,} \end{aligned}$$

where a dot means the scalar product with respect to \(\mathbf {q.}\) This property allows one to define two second fundamental forms (i.e. corresponding extrinsic curvature d-tensors) of hypersurface \(\Xi \)

$$\begin{aligned} \begin{array}{cccc} \, _{\shortmid }K: &{} T_{p}\Xi \times T_{p}\Xi &{} \longrightarrow &{} \mathbb {R} \text{ LC-configurations }, \\ &{} ( \mathbf {a,b}) &{} \longrightarrow &{} -\mathbf {a}{\cdot } \chi \mathbf {(b);} \\ \, _{\shortmid }\widehat{\mathbf {K}}: &{} (h\, _{\shortmid }\Xi \oplus v\, _{\shortmid }\Xi )\times (h\, _{\shortmid }\Xi \oplus v\, _{\shortmid }\Xi ) &{} \longrightarrow &{} \mathbb {R\oplus R} \text{ canonical } \text{ configurations }, \\ &{} ( \mathbf {a,b}) &{} \longrightarrow &{} -\mathbf {a}{\cdot }\widehat{ \chi }\mathbf {(b).} \end{array} \end{aligned}$$

In explicit operator form, \(\, _{\shortmid }K( \mathbf {a,b}) =- \mathbf {a}\,{\cdot } \, _{\shortmid }\nabla _{\mathbf {b}}\mathbf {n}\) and \(\, _{\shortmid }\widehat{\mathbf {K}}( \mathbf {a,b}) =-\mathbf {a}\,{\cdot } \, _{\shortmid }\widehat{\mathbf {D}}_{\mathbf {b}}{\mathbf {n}},\) which allows us to compute in coefficient form \(\, _{\shortmid }K_{\grave{\imath },\grave{j} }\) and \(\, _{\shortmid }\widehat{\mathbf {K}}_{\grave{\imath },\grave{j}}.\) Any pseudo-Riemannian geometry can be written equivalently in terms of \(( \mathbf {q,}\, _{\shortmid }\widehat{\mathbf {K}}) ,\) or \(( \mathbf { q,}\, _{\shortmid }K) ,\) for any data \((\Xi ,\mathbf {q,\, _{\shortmid }N, }\breve{N}^{\grave{\imath }},\breve{N}).\) Intuitively, we can work on 3-d space-like hypersurfaces as in Riemannian geometry. Unfortunately, such \(3+1\) splitting nonholonomic variables are not convenient for decoupling the gravitational field equations in general form.

2.3 Important formulas on space-like N-adapted hypersurfaces

Let us consider 3-surfaces \(\Xi _{t}\) enabled with Riemannian d-metrics when the normal d-vector \(\mathbf {n}\) is time-like. We follow a 4-d point of view treating d-tensor fields defined on any \(\Xi \) as they are defined for \((\mathbf {V},\mathbf {g,N}).\) This avoids the obligation to introduce special frame/coordinate systems and complicated notations depending on double fibration parameterizations etc.

We consider \(\mathrm{Vect}(\mathbf {n})\) as the 1-d subspace of \(T_{p} \mathbf {V}\) generated by the d-vector \(\mathbf {n}\) using \(\underline{\mathbf {n }}=\{\mathbf {n}_{\alpha }\}\) for its dual 1-form. In a point \(p\in \Xi ,\) the spaces of all spacetime vectors can be decomposed as \(\ T_{p}\mathbf {V=} T_{p}\Xi \oplus \mathrm{Vect}(\mathbf {n}).\) For \(\mathbf {n}{\cdot } \mathbf {n}=-1\) and any \( \mathbf {v}\in T_{p}\Xi ,\) we have \(\overrightarrow{\mathbf {q}}(\mathbf {n})=0\) and \(\overrightarrow{\mathbf {q}}(\mathbf {v})=\mathbf {v}.\) We construct the orthogonal projector onto \(\Xi \) as the operator \(\overrightarrow{\mathbf {q}} \) following the rule

$$\begin{aligned} \overrightarrow{\mathbf {q}}:\, T_{p}\mathbf {V}\rightarrow T_{p}\Xi ;\quad \mathbf {v}\rightarrow \mathbf {v+(n}{\cdot } \mathbf {v)n}. \end{aligned}$$

Similarly, we can consider a mapping \(\overrightarrow{\mathbf {q}}_{\mathbf {V} }^{*}:T_{p}^{*}\Xi \rightarrow T_{p}^{*}\mathbf {V}\) setting for any linear form \(\underline{\mathbf {z}}\in T_{p}^{*}\Xi \)

$$\begin{aligned} \overrightarrow{\mathbf {q}}_{\mathbf {V}}^{*}:\, T_{p}\mathbf {V} \rightarrow \mathbb {R};\quad \mathbf {v}\rightarrow \underline{\mathbf {z}}( \overrightarrow{\mathbf {q}}(\mathbf {v})). \end{aligned}$$

The maps can be extended to bilinear forms. Taking the induced 3-metric \( \mathbf {q}\) on \(\Xi ,\) we can consider

$$\begin{aligned} \begin{aligned} \mathbf {q}&:=\overrightarrow{\mathbf {q}}_{\mathbf {V}}^{*}\mathbf {g},\qquad \quad \; \text{ i.e. } \quad \mathbf {q=g+}\underline{\mathbf {n}}\otimes \underline{ \mathbf {n}}, \\ q_{\ \beta }^{\alpha }&=\delta _{\ \beta }^{\alpha }+\mathbf {n}^{\alpha } \mathbf {n}_{\beta }\quad \text{ and } \quad \mathbf {q}_{\alpha \beta }=\mathbf {g} _{\alpha \beta }+\mathbf {n}_{\alpha }\mathbf {n}_{\beta }. \end{aligned} \end{aligned}$$
(23)

We get the same results as for the d-metric \(\mathbf {q}\) if the two arguments of \(\mathbf {q}({\cdot } ,{\cdot } )\) are tangent d-vectors to \(\Xi .\) Such an operator gives zero if a d-vector is orthogonal to \(\Xi ,\) i.e. parallel to \(\mathbf {n.}\)

The nonholonomicmatter energy density is defined

$$\begin{aligned} \ ^{m}\widehat{E}:=\ ^{m}\widehat{\mathbf {T}}(\mathbf {n,n}),\quad \text{ i.e. } \quad ^{m}\widehat{E}=\ ^{m}\widehat{\mathbf {T}}_{\alpha \beta }\mathbf {n}^{\alpha }\mathbf {n}^{\beta }. \end{aligned}$$

Such values are written “without hat” (\(\ ^{m}\widehat{E}\rightarrow \ ^{m}E \) and \(\ ^{m}\widehat{\mathbf {T}}_{\alpha \beta }\rightarrow \ ^{m} \mathbf {T}_{\alpha \beta }\)), if \(\widehat{\mathbf {D}}\rightarrow \nabla \) in \(\ ^{m}\widehat{\mathbf {T}};\) see (15).

Similarly, the nonholonomic matter momentum density is

$$\begin{aligned} \widehat{\mathbf {p}}:=-\ ^{m}\widehat{\mathbf {T}}(\mathbf {n,}\overrightarrow{ \mathbf {q}}(...)),\quad \text{ i.e. } \quad \langle \widehat{\mathbf {p}},\mathbf {a}\rangle =-\ ^{m} \widehat{\mathbf {T}}(\mathbf {n,}\overrightarrow{\mathbf {q}}(\mathbf {a})), \end{aligned}$$

or \(\widehat{\mathbf {p}}_{\alpha }=-\ ^{m}\widehat{\mathbf {T}}_{\mu \nu } \mathbf {n}^{\mu }q_{\ \alpha }^{\nu },\) for any \(\mathbf {a}\in (h\Xi \oplus v\Xi ).\)

The nonholonomic stress d-tensor is the bilinear form

$$\begin{aligned}&\widehat{\mathbf {S}}:=\overrightarrow{\mathbf {q}}^{*}\ ^{m}\widehat{ \mathbf {T}} \text{ and } \text{ in } \text{ N-adapted } \text{ components } \\&\quad \widehat{\mathbf {S}}_{\alpha \beta }=\ ^{m}\widehat{\mathbf {T}}_{\mu \nu }q_{\ \alpha }^{\mu }q_{\ \beta }^{\nu }. \end{aligned}$$

We can consider also the trace of this field, \(\widehat{S}:=\mathbf {q}^{ \grave{\imath }\grave{j}}\widehat{\mathbf {S}}_{\grave{\imath }\grave{j}}= \mathbf {g}^{\alpha \beta }\widehat{\mathbf {S}}_{\alpha \beta }.\)

Both \(\widehat{\mathbf {p}}\) and \(\widehat{\mathbf {S}}\) are d-tensor fields tangent to \(\Xi _{t}.\) The data \((\widehat{E},\widehat{\mathbf {p}},\widehat{ \mathbf {S}})\) allow one to reconstruct

$$\begin{aligned} \ ^{m}\widehat{\mathbf {T}}=\widehat{\mathbf {S}}+\underline{\mathbf {n}} \otimes \widehat{\mathbf {p}}+\widehat{\mathbf {p}}\otimes \underline{\mathbf {n }}+\ ^{m}E\underline{\mathbf {n}}\otimes \underline{\mathbf {n}}, \end{aligned}$$

when \(\widehat{T}=\widehat{S}-\widehat{E}.\)

2.3.1 Einstein equations in nonholonomic variables for double \(2+2\) and \(3+1\) splitting

Summarizing the above formulas, the nonholonomic version of Einstein equations (13) with matter sources of type (15) can be written in the form

$$\begin{aligned} \begin{array}{lll} \begin{array}{l} \text{4-d } \text{ indices } \\ \\ \text{3-d } \text{ indices } \\ \end{array} &{} &{} \begin{array}{l} \mathcal {L}_{\mathbf {m}}\widehat{\mathbf {K}}_{\alpha \beta }:=-\, _{\shortmid }\widehat{\mathbf {D}}_{\alpha }\, _{\shortmid }\widehat{\mathbf {D}}_{\beta }N+N\{\, _{\shortmid }\widehat{\mathbf {R}}_{\alpha \beta }+\, _{\shortmid } \widehat{K}\widehat{\mathbf {K}}_{\alpha \beta }-2\widehat{\mathbf {K}} _{\alpha \mu }\widehat{\mathbf {K}}_{\ \beta }^{\mu }+4\pi [(\widehat{S }-\widehat{E})\mathbf {q}_{\alpha \beta }-2\widehat{\mathbf {S}}_{\alpha \beta }]\}, \\ \\ \mathcal {L}_{\mathbf {m}}\widehat{\mathbf {K}}_{\grave{\imath }\grave{j}}:=-\, _{\shortmid }\widehat{\mathbf {D}}_{\grave{\imath }}\, _{\shortmid }\widehat{ \mathbf {D}}_{\grave{j}}N+N\{\, _{\shortmid }\widehat{\mathbf {R}}_{\grave{ \imath }\grave{j}}+\, _{\shortmid }\widehat{K}\widehat{\mathbf {K}}_{\grave{ \imath }\grave{j}}-2\widehat{\mathbf {K}}_{\grave{\imath }\grave{k}}\widehat{ \mathbf {K}}_{\ \grave{j}}^{\grave{k}}+4\pi [(\widehat{S}-\widehat{E}) \mathbf {q}_{\grave{\imath }\grave{j}}-2\widehat{\mathbf {S}}_{\grave{\imath } \grave{j}}]\}; \\ \end{array}\\ \text{ Hamilt. } \text{ constr. } &{} &{} \, _{\shortmid }\widehat{R}+(\, _{\shortmid } \widehat{K})^{2}-\, _{\shortmid }\widehat{\mathbf {K}}_{\grave{\imath }\grave{j} }\, _{\shortmid }\widehat{\mathbf {K}}^{\grave{\imath }\grave{j}}=16\pi \widehat{E}, \\ \text{ momentum } \text{ constr. } &{} &{} \, _{\shortmid }\widehat{\mathbf {D}}_{\grave{j} }\, _{\shortmid }\widehat{\mathbf {K}}_{\ \grave{\imath }}^{\grave{j}}-\, _{\shortmid }\widehat{\mathbf {D}}_{\grave{\imath }}\, _{\shortmid }\widehat{K} =8\pi \widehat{\mathbf {p}}_{\grave{\imath }}. \end{array} \end{aligned}$$
(24)

Such systems of nonlinear partial differential equations (PDEs) are useful for stating the Cauchy problem, defining energy and momentum type values and elaborating methods of canonical and/or loop quantization. To decouple and solve such systems of equations in general analytic forms is a very difficult technical task even if such constructions are used in numeric analysis. Nevertheless, we can compute all values in (24) using any class of solutions found for (13); see Sect. 4.

2.3.2 Weyl’s tensor \(3+1\) projections adapted to nonholonomic \(2+2\) splitting

Taking any time-like unit vector \(\mathbf {v}^{\alpha }\), we can define a projection tensor, \(\mathbf {h}_{\alpha \beta }=\mathbf {g}_{\alpha \beta }+ \mathbf {v}_{\alpha }\mathbf {v}_{\beta }.\) A general d-vector \(\mathbf {v} _{\alpha }\) is necessary for studying relativistic thermodynamic and hydrodynamic models; see Refs. [7, 9, 25, 31, 33,34,35, 37, 39]. For constructions in this work, we can consider \(\mathbf {v}_{\alpha }=\mathbf {n} _{\alpha }\) and use a 3-d hypersurface metric \(\mathbf {q}_{\alpha \beta }= \mathbf {g}_{\alpha \beta }+\mathbf {n}_{\alpha }\mathbf {n}_{\beta }.\)

A \((3+1)+(2+2)\) covariant description of gravitational field is possible by splitting, respectively, into irreducible parts such that

$$\begin{aligned} \widehat{\mathbf {D}}_{\beta }\mathbf {v}_{\alpha }=-\mathbf {v}_{\beta } \mathbf {v}^{\gamma }\widehat{\mathbf {D}}_{\gamma }\mathbf {v}_{\alpha }+\frac{ 1}{3}\widehat{\mathbf {\Theta }}\mathbf {h}_{\alpha \beta }+\widehat{\sigma } _{\alpha \beta }+\widehat{\omega }_{\alpha \beta }, \end{aligned}$$

where \(\mathbf {v}^{\gamma }\widehat{\mathbf {D}}_{\gamma }\mathbf {v}^{\alpha } \) is the acceleration vector, \(\widehat{\mathbf {\Theta }}:=\mathbf {h} ^{\alpha \beta }\widehat{\mathbf {D}}_{\beta }\mathbf {v}_{\alpha }\) is the expansion scalar, \(\widehat{\sigma }_{\alpha \beta }=[\mathbf {h}_{~(\alpha }^{\gamma }\mathbf {h}_{\beta )}^{~\delta }-\frac{1}{3}\mathbf {h}_{\alpha \beta }\mathbf {h}^{\gamma \delta }]\widehat{\mathbf {D}}_{\gamma }\mathbf {v} _{\delta }\) is the shear tensor and \(\widehat{\omega }_{\alpha \beta }= \mathbf {h}_{~[\alpha }^{\gamma }\mathbf {h}_{\beta ]}^{~\delta }\widehat{ \mathbf {D}}_{\gamma }\mathbf {v}_{\delta }\) is the vorticity tensor (in the last two formulas (...) and [...] mean, respectively, symmetrization and anti-symmetrization of indices). We shall use the following decompositions of the Weyl tensor:

$$\begin{aligned} \widehat{\mathbf {C}}_{\ \alpha \beta }^{\tau \quad \gamma }:=\widehat{ \mathbf {R}}_{\ \alpha \beta }^{\tau \quad \gamma }+2\widehat{\mathbf {R}} _{[\alpha }^{\quad [\tau }\delta _{\beta ]}^{\gamma ]}+\frac{1}{3} \widehat{\mathbf {R}}\delta _{[\alpha }^{\gamma }\delta _{\beta ]}^{\tau } \end{aligned}$$
(25)

into, respectively, electric and magnetic like parts (similar to the Maxwell theory),

$$\begin{aligned} \widehat{\mathbf {E}}_{\alpha \beta }=\widehat{\mathbf {C}}_{\alpha \beta \gamma \delta }\mathbf {v}^{\gamma }\mathbf {v}^{\delta }\quad \text{ and } \quad \widehat{ \mathbf {H}}_{\alpha \beta }=\frac{1}{2}\eta _{\alpha \gamma \delta }\widehat{ \mathbf {C}}_{\ \quad \alpha \beta }^{\gamma \delta }\mathbf {v}^{\gamma } \mathbf {v}^{\varepsilon }. \end{aligned}$$
(26)

In these formulas, \(\eta _{\alpha \beta \gamma }=\eta _{\alpha \beta \gamma \delta }\mathbf {v}^{\delta }\) are for the spatial alternating tensor \(\eta _{\alpha \beta \gamma \delta }=\eta _{[\alpha \beta \gamma \delta ]}\) with \(\eta _{1234}=\sqrt{|g_{\alpha \beta }|}.\) For any space-like unit vectors \(\mathbf {x}^{\alpha },\mathbf {y}^{\alpha },\mathbf {z}^{\alpha }\) that together with \(\mathbf {v}^{\alpha }\) form an orthonormal basis, we can introduce null tetrads, thus:

$$\begin{aligned} \mathbf {m}^{\alpha }:= & {} \frac{1}{\sqrt{2}}(\mathbf {x}^{\alpha }-i\mathbf {y} ^{\alpha }),\mathbf {l}^{\alpha }:\nonumber \\= & {} \frac{1}{\sqrt{2}}(\mathbf {v}^{\alpha }- \mathbf {z}^{\alpha })\quad \text{ and } \quad \mathbf {k}^{\alpha }:=\frac{1}{\sqrt{2}}( \mathbf {v}^{\alpha }+\mathbf {z}^{\alpha }), \end{aligned}$$

for \(i^{2}=-1,\) and we express the metric \(\mathbf {g}_{\alpha \beta }=2\mathbf {m }_{(\alpha }\mathbf {m}_{\beta )}-2\mathbf {k}_{(\alpha }\mathbf {l}_{\beta )}.\) We can transform a double fibration, for instance, into a \(3+1\) fibration if we substitute \(\widehat{\mathbf {D}}\rightarrow \nabla ,\) omit hats, and change boldface symbols into similar non-boldface ones and work with arbitrary bases instead of N-adapted ones.

Equations (25) and (26) are important for constructing the thermodynamical values for gravitational fields following the CET model [40].

3 Geometric evolution of Einstein gravitational fields

Nonholonomic Ricci flows for theories with distinguished connections were considered in [20, 21]. To study general relativistic geometric flows and elaborate on thermodynamical models we use a N-adapted \(3+1\) decomposition for the canonical d-connection, \(\widehat{ \mathbf {D}}=(\, _{\shortmid }\widehat{\mathbf {D}},\ ^{t}\widehat{D})\) and d-metric \(\mathbf {g}:=(\mathbf {q,}\breve{N})\) of a 4-d spacetime \(\mathbf {V }\). In this section we formulate the general relativistic geometric flow theory in nonholonomic variables with double splitting. The fundamental functionals and Ricci flow evolution equations are constructed on 4-d Lorentz manifolds determined by exact solutions in GR.

3.1 Distortion relations on induced linear connections

On closed 3-d space-like hypersurfaces, geometric flow and gravitational field theories can be formulated in two equivalent forms using the connections \(\, _{\shortmid }\nabla \) and/or \(\, _{\shortmid }\widehat{\mathbf { D}}.\) The evolution of such connections and N-adapted frames is determined by the evolution of the hypersurface metric \(\mathbf {q}\). In N-adapted variables, we can introduce the canonical Laplacian d-operator, \(\, _{\shortmid }\widehat{\Delta }:=\) \(\, _{\shortmid }\widehat{\mathbf {D}}\) \(\, _{\shortmid }\widehat{\mathbf {D}}\) and define the canonical distortion tensor \( \, _{\shortmid }\widehat{\mathbf {Z}}\). The distortions of the Ricci d-tensor and the corresponding Ricci scalar are computed by introducing \(\, _{\shortmid }\nabla =\, _{\shortmid }\widehat{\mathbf {D}}-\, _{\shortmid } \widehat{\mathbf {Z}}\) (22) in

$$\begin{aligned} \, _{\shortmid }\widehat{\Delta }=\, _{\shortmid }\widehat{ \mathbf {D}}_{\alpha }\ \, _{\shortmid }\widehat{\mathbf {D}}^{\alpha }=\, _{\shortmid }\Delta +\, _{\shortmid }^{Z}\widehat{\Delta }. \end{aligned}$$

We obtain

$$\begin{aligned} \, _{\shortmid }\Delta= & {} \, _{\shortmid }\nabla _{\grave{\imath }}\ \, _{\shortmid }\nabla ^{\grave{\imath }} =\, _{\shortmid }\nabla _{\alpha }\ \, _{\shortmid }\nabla ^{\alpha }, \nonumber \\ \, _{\shortmid }^{Z}\widehat{\Delta }= & {} \, _{\shortmid } \widehat{\mathbf {Z}}_{\grave{\imath }}\ \, _{\shortmid }\widehat{ \mathbf {Z}}^{\grave{\imath }}-[\, _{\shortmid }\widehat{\mathbf {D}}_{\grave{ \imath }}(\, _{\shortmid }\widehat{\mathbf {Z}}^{\grave{\imath }})+ \, _{\shortmid }\widehat{\mathbf {Z}}_{\grave{\imath }}(\, _{\shortmid } \widehat{\mathbf {D}}^{\grave{\imath }})]\nonumber \\= & {} \, _{\shortmid }\widehat{ \mathbf {Z}}_{\alpha }\, _{\shortmid }\widehat{\mathbf {Z}}^{\alpha }-[\, _{\shortmid }\widehat{\mathbf {D}}_{\alpha }(\, _{\shortmid } \widehat{\mathbf {Z}}^{\alpha })+\, _{\shortmid }\widehat{\mathbf {Z}} _{\alpha }(\, _{\shortmid }\widehat{\mathbf {D}}^{\alpha })]; \nonumber \\ \, _{\shortmid }\widehat{\mathbf {R}}_{\grave{\imath }\grave{j}}= & {} \, _{\shortmid }R_{\grave{\imath }\grave{j}}-\, _{\shortmid } \widehat{\mathbf {Z}}ic_{\grave{\imath }\grave{j}},\, _{\shortmid } \widehat{\mathbf {R}}_{\ \beta \gamma }=\, _{\shortmid }R_{\ \beta \gamma }-\, _{\shortmid }\widehat{\mathbf {Z}}ic_{\beta \gamma },\ \nonumber \\ \, _{\shortmid }\widehat{R}= & {} \, _{\shortmid }R-\mathbf {g}^{\beta \gamma }\, _{\shortmid }\widehat{\mathbf {Z}}ic_{\beta \gamma }=\, _{\shortmid }R-\mathbf {q}^{\grave{\imath }\grave{j}}\, _{\shortmid } \widehat{\mathbf {Z}}ic_{\grave{\imath }\grave{j}}=\, _{\shortmid }R- \, _{\shortmid }\widehat{\mathbf {Z}}, \nonumber \\ \, _{\shortmid }\widehat{\mathbf {Z}}= & {} \mathbf {g}^{\beta \gamma } \, _{\shortmid }\widehat{\mathbf {Z}}ic_{\beta \gamma }=\mathbf {q}^{ \grave{\imath }\grave{j}}\, _{\shortmid }\widehat{\mathbf {Z}}ic_{ \grave{\imath }\grave{j}}\nonumber \\= & {} \, _{h}\widehat{Z}+\, _{v}\widehat{Z},\, _{h}\widehat{Z }=g^{ij}\ \widehat{\mathbf {Z}}ic_{ij},\, _{v}\widehat{Z}=h^{ab}\ \widehat{ \mathbf {Z}}ic_{ab}; \nonumber \\ R= & {} \, _{h}R+\, _{v}R,\ \, _{h}R:=g^{ij}\ R_{ij},\, _{v}R=h^{ab}\ R_{ab}. \end{aligned}$$
(27)

Such values can be computed in explicit form for any class of exact solutions of the nonholonomic Einstein equations (13), when a double \(2+2\) and \(3+1\) splitting is prescribed and the LC-conditions (14) can be imposed additionally.

3.2 Nonholonomc Perelman’s functionals on 3-d hypersurfaces

For standard Ricci flows on a normalized 3-d space-like closed hypersurface \( \ ^{c}\widehat{\Xi }\subset \mathbf {V},\) the normalized Hamilton equations written in a coordinate basis are

$$\begin{aligned}&\partial _{\tau }q_{\grave{\imath }\grave{j}} =-2\, _{\shortmid }R_{\grave{ \imath }\grave{j}}+\frac{2\grave{r}}{5}q_{\grave{\imath }\grave{j}}, \\&q_{\grave{\imath }\grave{j}\mid \tau =0} =q_{\grave{\imath }\grave{j} }^{[0]}[x^{\grave{\imath }}]. \nonumber \end{aligned}$$
(28)

We use the left label “c” for the conditions “compact and closed” and do not emphasize the dependence on space coordinates (writing in brief \(q_{\grave{\imath }\grave{j}}(x^{\grave{\imath }},\tau )=q_{\grave{\imath }\grave{j}}(\tau ))\) if this does not result in ambiguities. In the above formulas, \(\, _{\shortmid }R_{ \grave{\imath }\grave{j}}\) is computed for the Levi-Civita connection \(\, _{\shortmid }\nabla \) of \(q_{\grave{\imath }\grave{j}}(\tau )\) parameterized by a real variable \(\tau ,\) \(0\le \tau <\tau _{0},\) for a differentiable function \(\tau (t).\) The boundary conditions are stated for \(\tau =0\) and the normalizing factor

is introduced in a form so as to preserve the volume of \(\ \ ^{c}\widehat{\Xi },\) i.e. \(\int _{\ ^{c}\Xi }\sqrt{|q_{\grave{\imath }\grave{j}}|}\mathrm{d}\grave{x}^{3}.\) For simplicity, we can find solutions of (1) with \(\grave{r}=0.\)

In order to find explicit solutions of (28) for \(q_{\grave{\imath } \grave{j}}\subset g_{\alpha \beta }\) with \(g_{\alpha \beta }\) defined also as a solution of a 4-d Einstein equations (12), we have to consider a nontrivial \(\grave{r}.\) We can re-write (1) in any nonholonomic basis using the geometric evolution of frame fields, \(\ \partial _{\chi }e_{\grave{\imath }}^{\ \underline{\grave{\imath }}}=q^{ \underline{\grave{\imath }}\underline{\grave{j}}}\, _{\shortmid }R_{\underline{ \grave{j}}\underline{\grave{k}}}e_{\grave{\imath }}^{\ \underline{\grave{k}} }, \) when \(q_{\grave{\imath }\grave{j}}(\tau )=q_{\underline{\grave{\imath }} \underline{\grave{j}}}(\tau )e_{\grave{\imath }}^{\ \underline{\grave{\imath }} }(\tau )e_{\grave{j}}^{\ \underline{\grave{j}}}(\tau )\) for \(e_{\grave{\imath }}(\tau )=e_{\grave{\imath }}^{\ \underline{\grave{\imath }}}(\tau )\partial _{ \underline{\grave{\imath }}}\) and \(e^{j}(\tau )=e_{\ \underline{\grave{j}} }^{j}(\tau )\mathrm{d}x^{\underline{\grave{j}}}.\) There is a unique solution for such systems of linear ODEs for any \(\tau \in [0,\tau _{0}).\)

In nonholonomic variables and for the linear d-connection \(\, _{\shortmid } \widehat{\mathbf {D}},\) the Perelman functionals parameterized in N-adapted form are written

$$\begin{aligned} \, _{\shortmid }\widehat{\mathcal {F}}= & {} \int _{\widehat{\Xi }_{t}}e^{-f}\sqrt{ |q_{\grave{\imath }\grave{j}}|}\mathrm{d}\grave{x}^{3}(\, _{\shortmid }\widehat{R}+|\, _{\shortmid }\widehat{\mathbf {D}}f|^{2}), \end{aligned}$$
(29)
$$\begin{aligned}&\text{ and } \nonumber \\ \, _{\shortmid }\widehat{\mathcal {W}}= & {} \int _{\widehat{\Xi }_{t}}M\sqrt{|q_{ \grave{\imath }\grave{j}}|}\mathrm{d}\grave{x}^{3}[\tau (\, _{\shortmid }\widehat{R}+|\ \, _{\shortmid }^{h}\widehat{\mathbf {D}}f|\nonumber \\&+|\ \, _{\shortmid }^{v}\widehat{ \mathbf {D}}f|)^{2}+f-6], \end{aligned}$$
(30)

where the scaling function f satisfies \(\int _{\widehat{\Xi }_{t}}M\sqrt{ |q_{\grave{\imath }\grave{j}}|}\mathrm{d}\grave{x}^{3}=1\) for \(M=( 4\pi \tau ) ^{-3}e^{-f}\) . The functionals \(\, _{\shortmid }\widehat{\mathcal {F} }\) and \(\, _{\shortmid }\widehat{\mathcal {W}}\) transform into standard Perelman functionals on \(\widehat{\Xi }_{t}\) if \(\, _{\shortmid }\widehat{ \mathbf {D}}\rightarrow \, _{\shortmid }\nabla .\) The W-entropy \(\, _{\shortmid }\widehat{\mathcal {W}}\) is a Lyapunov type non-decreasing functional.

3.3 Nonholonomic Ricci flow evolution equations for 3-d hypersurface metrics

Considering the dependencies of Eqs. (29) and (30) on a smooth parameter \(\chi (\tau )\) for which \(\partial \chi /\partial \tau =-1\) (when, for simplicity, the normalization terms are not included) one can prove the geometric evolution equations for any induced 3-d metric \(\mathbf { q}\) and canonical d-connection \(\, _{\shortmid }\widehat{\mathbf {D}}.\) A subclass of geometric flow models with general relativistic extension from 3-d to 4-d can be elaborated if one of the parameters \(\chi \) or \(\tau \) is taken to be proportional to the time-like fourth coordinate, t.

Applying the variational procedure to \(\, _{\shortmid }\widehat{\mathcal {F}}\) (29) in N-adapted form (see for details [20, 21] and the references therein and for a double nonholonomic splitting), we obtain

$$\begin{aligned}&\partial _{\chi }\mathbf {q}_{\grave{\imath }\grave{j}} =-2( _{\shortmid }\widehat{\mathbf {R}}_{\grave{\imath }\grave{j}}+\, _{\shortmid }\widehat{\mathbf {Z}}ic_{\grave{\imath }\grave{j}}), \end{aligned}$$
(31)
$$\begin{aligned}&\, _{\shortmid }\widehat{\mathbf {R}}_{i\grave{a}} =-\, _{\shortmid }\widehat{\mathbf {Z}}ic_{i\grave{a}}, \nonumber \\&\partial _{\chi }f =-(\, _{\shortmid }\widehat{\Delta }-\, _{\shortmid }^{Z}\widehat{\Delta })f+|(\, _{\shortmid }\widehat{\mathbf {D}}- \, _{\shortmid }\widehat{\mathbf {Z}})f|^{2}-\, _{\shortmid } \widehat{R}-\, _{\shortmid }\widehat{\mathbf {Z}}.\nonumber \\ \end{aligned}$$
(32)

The distortions in such nonholonomic evolution equations are completely determined by \(\mathbf {q}_{\grave{\imath }\grave{j}}\) following Eq. (27). There is also another important property:

$$\begin{aligned}&\partial _{\chi }\, _{\shortmid }\widehat{\mathcal {F}}(q,\, _{\shortmid } \widehat{\mathbf {D}},f)=2\int _{\ \ ^{c}\widehat{\Xi }_{t}}e^{-f}\sqrt{|q_{ \grave{\imath }\grave{j}}|}\mathrm{d}\grave{x}^{3}\nonumber \\&[|\, _{\shortmid }\widehat{\mathbf {R}} _{\grave{\imath }\grave{j}}+\, _{\shortmid }\widehat{\mathbf {Z}}ic_{\grave{ \imath }\grave{j}}+(\, _{\shortmid }\widehat{\mathbf {D}}_{\grave{\imath }}-\, _{\shortmid }\widehat{\mathbf {Z}}_{\grave{\imath }})(\, _{\shortmid }\widehat{ \mathbf {D}}_{\grave{j}}-\, _{\shortmid }\widehat{\mathbf {Z}}_{\grave{j} })f|^{2}], \end{aligned}$$

when \(\int _{\ \ ^{c}\widehat{\Xi }_{t}}e^{-f}\sqrt{|q_{\grave{\imath }\grave{j }}|}\mathrm{d}\grave{x}^{3}\) is constant for a fixed \(\tau \) and \(f(\chi (\tau ))=f(\tau ).\)

The system of equations (31) and (32) is equivalent to (28) up to a certain redefinition of nonholonomic frames and variables. We have to consider (32) as additional constraints because in nonholonomic variables the Ricci d-tensor is (in general) nonsymmetric. If we do not impose such constraints, the geometric evolution goes with nonsymmetric metrics.

3.4 Geometric evolution to 4-d Lorentz configurations as exact solutions in GR

The geometric evolution of metrics of type (19), when \(\mathbf {q} (\tau )\mathbf {\rightarrow } \mathbf {g}(\tau ):=( \mathbf {q}(\tau ) \mathbf {,}\breve{N}(\tau ))\), is described by respective generalizations of the functionals (29) and (30). Considering N-connection adapted foliations \(\widehat{\Xi }_{t}\) parameterized by a spacetime coordinate t,  we introduce such 4-d functionals,

$$\begin{aligned}&\widehat{\mathcal {F}}(\mathbf {q},\, _{\shortmid }\widehat{\mathbf {D}},f) =\int _{t_{1}}^{t_{2}}\mathrm{d}t\ \breve{N}(\tau )\, _{\shortmid }\widehat{\mathcal {F}} (\mathbf {q},\, _{\shortmid }\widehat{\mathbf {D}},f), \end{aligned}$$
(33)
$$\begin{aligned}&\text{ and } \nonumber \\&\widehat{\mathcal {W}}(\mathbf {q},\, _{\shortmid }\widehat{\mathbf {D}},f) =\int _{t_{1}}^{t_{2}}\mathrm{d}t\ \breve{N}(\tau )\ \, _{\shortmid }\widehat{\mathcal {W }}(\mathbf {q},\, _{\shortmid }\widehat{\mathbf {D}},f(t)).\ \end{aligned}$$
(34)

For arbitrary frame transforms on 4-d nonholonomic Lorentz manifolds, such values can be redefined, respectively, in terms of the data \((\mathbf {g}(\tau ), \widehat{\mathbf {D}}(\tau )).\) We obtain

$$\begin{aligned} \widehat{\mathcal {F}}= & {} \int _{t_{1}}^{t_{2}}\int _{\widehat{\Xi }_{t}}e^{- \widehat{f}}\sqrt{|\mathbf {g}_{\alpha \beta }|}\mathrm{d}^{4}u(\widehat{R}+|\widehat{ \mathbf {D}}\widehat{f}|^{2}), \end{aligned}$$
(35)

and

$$\begin{aligned} \widehat{\mathcal {W}}= & {} \int _{t_{1}}^{t_{2}}\int _{\widehat{\Xi }_{t}} \widehat{M}\sqrt{|\mathbf {g}_{\alpha \beta }|}\mathrm{d}^{4}u[\tau (\widehat{R}+|\ ^{h}\widehat{D}\widehat{f}|\nonumber \\&+|\ ^{v}\widehat{D}\widehat{f}|)^{2}+\widehat{f} -8], \end{aligned}$$
(36)

where the scaling function \(\widehat{f}\) satisfies \(\int _{t_{1}}^{t_{2}} \int _{\widehat{\Xi }_{t}}\widehat{M}\sqrt{|\mathbf {g}_{\alpha \beta }|} \mathrm{d}^{4}u=1\) for \(\widehat{M}=( 4\pi \tau ) ^{-3}e^{-\widehat{f}} . \)

We emphasize that the functionals \(\widehat{\mathcal {F}}\) and \(\widehat{ \mathcal {W}}\) for 4-d pseudo-Riemannian metrics are not of entropy type like in the 3-d Riemannian case. Nevertheless, they describe nonlinear general relativistic diffusion-type processes if \(\mathbf {q}\subset \mathbf {g}\ \) and \(_{\shortmid }\widehat{\mathbf {D}}\) \(\subset \widehat{\mathbf {D}}\) are determined by certain lapse and shift functions as certain solutions of the 4-d gravitational equations. This is motivation to construct such functionals using nonholonomic variables.

Equations (31) and (32) can be considered for 4-d configurations with the coefficients determined by the Ricci d-tensors and the distortions will keep the left label “\(\shortmid \)”. Using the formulas \( \mathbf {q}_{\alpha \beta }=\mathbf {g}_{\alpha \beta }+\mathbf {n}_{\alpha } \mathbf {n}_{\beta }\) (see (23)), we get

$$\begin{aligned}&\partial _{\tau }\mathbf {g}_{\alpha \beta } =-2(\, _{\shortmid } \widehat{\mathbf {R}}_{\alpha \beta }+\, _{\shortmid }\widehat{ \mathbf {Z}}ic_{\alpha \beta })-\partial _{\tau }(\mathbf {n}_{\alpha }\mathbf { n}_{\beta }), \nonumber \\&\, _{\shortmid }\widehat{\mathbf {R}}_{ia} =-\, _{\shortmid }\widehat{\mathbf {Z}}ic_{ia},\quad \text{ for } \; \, _{\shortmid }\widehat{ \mathbf {R}}_{\alpha \beta }\quad \text{ with } \; \alpha \ne \beta . \end{aligned}$$
(37)

The term \(\partial _{\tau }(\mathbf {n}_{\alpha }\mathbf {n}_{\beta })\) can be computed in explicit form using the formulas for the geometric evolution of N-adapted frames; see below Eq. (39). The formulas for \( \partial _{\tau }f,\partial _{\tau }(_{\shortmid }\widehat{\mathcal {F}})\) and \(\partial _{\tau }\widehat{\mathcal {W}}\) can be re-written for values with 4-d indices (we omit such formulas in this work).

For 4-d configurations with a corresponding redefinition of the scaling function, \(f\rightarrow \widehat{f},\) and for necessary type N-adapted distributions, we can construct models of the geometric evolution with h- and v-splitting for \(\mathbf {\widehat{\mathbf {D}},}\)

$$\begin{aligned}&\partial _{\tau }g_{ij} =-2\widehat{R}_{ij},\quad \partial _{\tau }g_{ab}=-2 \widehat{R}_{ab}, \nonumber \\&\partial _{\tau }\widehat{f} =-\widehat{\square }\widehat{f}+\left| \widehat{\mathbf {D}}\widehat{f}\right| ^{2}-\ ^{h}\widehat{R}-\ ^{v} \widehat{R}. \end{aligned}$$
(38)

These formulas can be derived from the functional (33), following a similar calculation to that presented in the proof of Proposition 1.5.3 of [12] but in N-adapted form as in Refs. [20, 21]. We have to impose the conditions \(\widehat{R}_{ia}=0\) and \(\widehat{R} _{ai}=0 \) if we want to keep the total metric symmetric under Ricci flow evolution. The general relativistic character of the 4-d geometric flow evolution is encoded in operators like \(\widehat{\square }=\widehat{\mathbf {D }}^{\alpha }\widehat{\mathbf {D}}_{\alpha },\) d-tensor components \(\widehat{R} _{ij}\) and \(\widehat{R}_{ab},\) their scalars \(\ ^{h}\widehat{R}=g^{ij} \widehat{R}_{ij}\) and \(\ ^{v}\widehat{R}=g^{ab}\widehat{R}_{ab}\) with the data \( (g_{ij},g_{ab},\widehat{\mathbf {D}}_{\alpha })\) constrained by the condition to define solutions of certain 4-d Einstein equations.

The evolution with the parameter \(\chi \in [0,\chi _{0})\) of N-adapted frames in a 4-d nonholonomic Lorentz manifold can be computed as

$$\begin{aligned} \ \mathbf {e}_{\alpha }(\chi )=\ \mathbf {e}_{\alpha }^{\ \underline{\alpha } }(\chi ,u)\partial _{\underline{\alpha }}. \end{aligned}$$

Up to frame/coordinate transforms the frame coefficients are

$$\begin{aligned}&\mathbf {e}_{\alpha }^{\ \underline{\alpha }}(\chi ,u)=\left[ \begin{array}{cc} \ e_{i}^{\ \underline{i}}(\chi ,u) &{} -~N_{i}^{b}(\chi ,u)\ e_{b}^{\ \underline{a}}(\tau ,u) \\ 0 &{} \ e_{a}^{\ \underline{a}}(\chi ,u) \end{array} \right] ,\nonumber \\&\mathbf {e}_{\ \underline{\alpha }}^{\alpha }(\chi ,u)\ =\left[ \begin{array}{cc} e_{\ \underline{i}}^{i}=\delta _{\underline{i}}^{i} &{} e_{\ \underline{i} }^{b}=N_{k}^{b}(\chi ,u)\ \ \delta _{\underline{i}}^{k} \\ e_{\ \underline{a}}^{i}=0 &{} e_{\ \underline{a}}^{a}=\delta _{\underline{a} }^{a} \end{array} \right] , \end{aligned}$$

with \(\tilde{g}_{ij}(\chi )=\ e_{i}^{\ \underline{i}}(\chi ,u)\ e_{j}^{\ \underline{j}}(\chi ,u)\eta _{\underline{i}\underline{j}}\) and \(\tilde{g} _{ab}(\chi )=\ e_{a}^{\ \underline{a}}(\chi ,u)\ e_{b}^{\ \underline{b} }(\chi ,u)\eta _{\underline{a}\underline{b}}.\) For \(\eta _{\underline{i} \underline{j}}=\mathrm{diag}[+,+]\) and \(\eta _{\underline{a}\underline{b}}=\mathrm{diag}[1,-1]\) corresponding to the chosen signature of \({\tilde{\mathbf {g}}}_{\alpha \beta }^{[0]}(u),\) we have the evolution equations

$$\begin{aligned} \frac{\partial }{\partial \chi }\mathbf {e}_{\ \underline{\alpha }}^{\alpha }\ =\ \mathbf {g}^{\alpha \beta }~\widehat{\mathbf {R}}_{\beta \gamma }~\ \mathbf {e}_{\ \underline{\alpha }}^{\gamma }. \end{aligned}$$
(39)

Such equations are prescribed for models with geometric evolution determined by the canonical d-connection \(\widehat{\mathbf {D}}.\) Equation (39) can be written in terms of the Levi-Civita connection if distortions of type (27) are considered for the 4-d values determined by exact solutions in GR.

4 Generation of off-diagonal solutions

Let us summarize the anholonomic frame deformation method, AFDM, of constructing generic off-diagonal exact solutions in GR with possible dependencies on all spacetime coordinates (see details and various examples in [42, 47, 50, 51] and the references therein). Using N-adapted \(2+2\) frame and coordinate transforms,

$$\begin{aligned} \mathbf {g}_{\alpha \beta }(x^{i},t)=e_{\ \alpha }^{\alpha ^{\prime }}(x^{i},y^{a})e_{\ \beta }^{\beta ^{\prime }}(x^{i},y^{a})\widehat{\mathbf {g }}_{\alpha ^{\prime }\beta ^{\prime }}(x^{i},y^{a})\quad \text{ and } \\ \Upsilon _{\alpha \beta }(x^{i},t)=e_{\ \alpha }^{\alpha ^{\prime }}(x^{i},y^{a})e_{\ \beta }^{\beta ^{\prime }}(x^{i},y^{a})\widehat{\Upsilon }_{\alpha ^{\prime }\beta ^{\prime }}(x^{i},y^{a}), \end{aligned}$$

for a time-like coordinate \(y^{4}=t\) (\(i^{\prime },i,k,k^{\prime },\ldots =1,2\), and \(a,a^{\prime },b,b^{\prime },\ldots =3,4),\) we can parameterize the metric and effective source in certain adapted forms. We consider

$$\begin{aligned} \mathbf {g}= & {} \mathbf {g}_{\alpha ^{\prime }\beta ^{\prime }}\mathbf {e} ^{\alpha ^{\prime }}\otimes \mathbf {e}^{\beta ^{\prime }}=g_{i}(x^{k})\mathrm{d}x^{i}\otimes \mathrm{d}x^{j}\nonumber \\&+\,\omega ^{2}(x^{k},y^{3},t)h_{a}(x^{k},t) \mathbf {e}^{a}\otimes \mathbf {e}^{a} \nonumber \\= & {} q_{i}(x)\mathrm{d}x^{i}\otimes \mathrm{d}x^{i}+q_{3}(x,y)\mathbf {e}^{3}\otimes \mathbf {e} ^{3}\nonumber \\&-\breve{N}^{2}(x^{k},y^{3},t)\mathbf {e}^{4}\otimes \mathbf {e}^{4}, \\ \mathbf {e}^{3}= & {} \mathrm{d}y^{3}+n_{i}(x^{k},t)\mathrm{d}x^{i},\mathbf {e} ^{4}=\mathrm{d}t+w_{i}(x^{k},t)\mathrm{d}x^{i}. \nonumber \end{aligned}$$
(40)

This ansatz is a general one for the 4-d metric which can be written in the form (19) with

$$\begin{aligned} \breve{N}^{2}(u)=-\omega ^{2}h_{4}=-g_{4}. \end{aligned}$$
(41)

It allows for a straightforward extension of 3-d ansatz to 4-d configurations by introducing the values \(\breve{N}^{2}(x^{k},t)\) and \(w_{i}(x^{k},t)\) in order to generate exact solutions of the Einstein equations. The nontrivial respective N-connection, d-metric, and matter source coefficients are denoted

$$\begin{aligned} \ N_{i}^{3}= & {} \, _{\shortmid }N_{i}^{3}=n_{i}(x^{k},t);N_{i}^{4}=w_{i}(x^{k},t); \nonumber \\ \{g_{i^{\prime }j^{\prime }}\}= & {} \mathrm{diag}[g_{i}],g_{1}=g_{2}=q_{1}=q_{2}=e^{\psi (x^{k})}; \nonumber \\ \{g_{a^{\prime }b^{\prime }}\}= & {} \mathrm{diag}[\omega ^{2}h_{a}],h_{a}=h_{a}(x^{k},t),q_{3}=\omega ^{2}h_{3}; \nonumber \\&\quad \text{ and } \quad \Upsilon _{\alpha \beta }=\mathrm{diag}[\Upsilon _{i};\Upsilon _{a}],\nonumber \\&\quad \text{ for } \quad \Upsilon _{1}=\Upsilon _{2}=\widetilde{\Upsilon } (x^{k}),\Upsilon _{3}=\Upsilon _{4}=\Upsilon (x^{k},t). \nonumber \\ \end{aligned}$$
(42)

The ansatz (40) determines d-metrics of type (7) and (19) as 4-d generalizations of the 3-d hypersurface metric (17) for a nontrivial lapse function (41). The N-adapted coefficients (42) can be very general ones but for the assumption that there are N-adapted frames with respect to which the exact solutions for \(\omega =1\) have the Killing symmetry on \(\partial /\partial y^{3}\). For such configurations, there are N-adapted bases when the geometric and physical values do not depend on the coordinate \(y^{3}.\) For simplicity, we shall consider solutions with one Killing symmetry.Footnote 7

4.1 Decoupling of Einstein equations for nonholonomic \(2+2\) splitting

Let us introduce the values \(\alpha _{i}=h_{3}^{*}\partial _{i}\varpi ,\beta =h_{3}^{*}\ \varpi ^{*},\gamma =( \ln |h_{3}|^{3/2}/|h_{4}|) ^{*}\) determined by the generating function

$$\begin{aligned} \varpi&:=\ln |h_{3}^{*}/\sqrt{|h_{3}h_{4}|}|,\quad \nonumber \\&\text{ we } \text{ shall } \text{ also } \text{ use } \text{ the } \text{ value } \; \Psi :=e^{\varpi }. \end{aligned}$$
(43)

We shall use brief notations for the partial derivatives: \(a^{\bullet }=\partial _{1}a,b^{\prime }=\partial _{2}b,h^{*}=\partial _{4}h=\partial _{t}h.\)

For ansatz (40) written in terms of the data for the generating function (43) and the above coefficients, the nonholonomic Einstein equations (13) transform into a system of nonlinear PDEs with decoupling property, Footnote 8

$$\begin{aligned}&\psi ^{\bullet \bullet }+\psi ^{\prime \prime }=2~\widetilde{\Upsilon },\quad \varpi ^{*}h_{3}^{*}=2h_{3}h_{4}\Upsilon ,\nonumber \\&\quad n_{i}^{**}+\gamma n_{i}^{*}=0,\quad \beta w_{i}-\alpha _{i}=0. \end{aligned}$$
(44)

The unknown functions for this system are \(\psi (x^{i}),h_{a}(x^{k},t),w_{i}(x^{k},t)\) and \(n_{i}(x^{k},t).\)

We can simplify the system (43) and (44) using an important property which allows us to re-define the generating function, \(\Psi =\exp \varpi \longleftrightarrow \widetilde{\Psi }=\exp \widetilde{\varpi },\) and the effective source, \(\Upsilon \longleftrightarrow \Lambda =\mathrm{const},\Lambda \ne 0.\) Such nonlinear transforms are given by the formulas

$$\begin{aligned}&\Lambda (\Psi ^{2})^{*}=|\Upsilon |(\widetilde{\Psi }^{2})^{*}\quad \text{ and } \nonumber \\&\quad \Lambda \Psi ^{2}=\widetilde{\Psi }^{2}|\Upsilon |+\int \mathrm{d}t\widetilde{ \Psi }^{2}|\Upsilon |^{*}. \end{aligned}$$
(45)

For generating off-diagonal inhomogeneous and locally anisotropic cosmological solutions depending on t,  we have to consider generating functions for which \(\Psi ^{*}\ne 0.\) We obtain a system of nonlinear PDEs with effective cosmological constant \(\Lambda ,\)

$$\begin{aligned}&\psi ^{\bullet \bullet }+\psi ^{\prime } =2~\widetilde{\Upsilon },\ \end{aligned}$$
(46)
$$\begin{aligned}&\widetilde{\varpi }^{*}h_{3}^{*}=2h_{3}h_{4}\Lambda , \end{aligned}$$
(47)
$$\begin{aligned}&\ n_{i}^{**}+\gamma n_{i}^{*} =0,\ \ \end{aligned}$$
(48)
$$\begin{aligned}&\varpi ^{*}\ w_{i}-\partial _{i}\varpi =0 \end{aligned}$$
(49)
$$\begin{aligned}&\partial _{k}\omega +n_{k}\partial _{3}\omega +w_{k}\omega ^{*} = \mathbf {e}_{k}\omega =0. \end{aligned}$$
(50)

This system can be solved in very general forms by prescribing \(~\widetilde{ \Upsilon },\) \(\Lambda \) and \(\Psi ,\) or \(\widetilde{\Psi },\) by integrating the equations “step by step”. We obtain

$$\begin{aligned} g_{i}= & {} g_{i}[\psi ,~\widetilde{\Upsilon }]\simeq e^{\psi (x^{k})}\quad \nonumber \\&\text{ as } \text{ a } \text{ solution } \text{ of } \text{2-d } \text{ Poisson/Laplace } \text{ equations } \text{(46) }; \nonumber \\ h_{a}= & {} h_{a}[\widetilde{\Psi },\Lambda ]=h_{a}[\Psi ,\Upsilon ],\quad \nonumber \\&\text{ where } \; h_{3}=\frac{\widetilde{\Psi }^{2}}{4\Lambda }\quad \text{ and } \; h_{4}=\frac{(\widetilde{\Psi }^{*})^{2}}{\Theta }; \nonumber \\ n_{k}= & {} \, _{1}n_{k}+\, _{2}n_{k}\int \mathrm{d}t\ h_{4}/\left( \sqrt{|h_{3}|}\right) ^{3}\nonumber \\= & {} \, _{1}n_{k}+\, _{2}\widetilde{n}_{k}\int \mathrm{d}t\ (\widetilde{\Psi }^{*})^{2}/ \widetilde{\Psi }^{3}\Theta ; \nonumber \\ w_{i}= & {} \partial _{i}\varpi /\varpi ^{*}=\partial _{i}\Psi /\Psi ^{*}=\partial _{i}(\Psi )^{2}/(\Psi ^{2})^{*}\nonumber \\= & {} \partial _{i}\Theta /\Theta ^{*}; \nonumber \\ \omega= & {} \omega [\widetilde{\Psi },\Lambda ]=\omega [\Psi ,\Upsilon ]\;\nonumber \\&\text{ is } \text{ any } \text{ solution } \text{ of } \text{ first } \text{ order } \text{ system } \text{(50) }. \end{aligned}$$
(51)

By \(h_{a}[\Psi ,\Upsilon ],\) we denote the fact that the coefficients \(h_{a}\) depend functionally on two functions, \([\Psi ,\Upsilon ].\)

The solutions (51) contain also integration functions \(\ ^{0}h_{3}(x^{k}),\) \(\, _{1}n_{k}(x^{i})\) and \(\, _{2}n_{k}(x^{i}),\) or \(_{2} \widetilde{n}_{k}(x^{i})=8\, _{2}n_{k}(x^{i})|\Lambda |^{3/2},\) and the generating source

$$\begin{aligned} \Theta =\int \mathrm{d}t\ \Upsilon (\widetilde{\Psi }^{2})^{*}. \end{aligned}$$
(52)

We can satisfy the conditions for \(\omega \) in the second line in (49) if we keep, for simplicity, the Killing symmetry on \(\partial _{i}\) and take, for example, \(\omega ^{2}=|h_{4}|^{-1}.\) Such solutions are constructed in explicit form by solving Eqs. (44) and/or (46)–(50) for certain prescribed values of \(\Psi \) and \( \Upsilon \) and following certain assumptions on initial/boundary/asymptotic conditions, physical arguments on symmetries of solutions, compatibility with observational data etc. It should be emphasized that redefinitions of the generating functions of type (45) allow one to construct exact solutions for general sources \(\Upsilon \) using certain classes of solutions with nontrivial cosmological constant \(\Lambda \).

4.2 The Levi-Civita conditions

The solutions (51) are defined for the canonical d-connection \( \widehat{\mathbf {D}}\) and with respect to N-adapted frames. There are nontrivial coefficients of the nonholonomically induced torsion. We have to subject the d-metric and N-connection coefficients to additional nonholonomic constraints (14) in order to satisfy the torsionless conditions and extract Levi-Civita configurations. For the ansatz (40), such conditions can be written (for details see [42, 47, 50, 51])

$$\begin{aligned}&w_{i}^{*}=(\partial _{i}-w_{i}\partial _{4})\ln \sqrt{|h_{4}|},(\partial _{i}-w_{i}\partial _{4})\ln \sqrt{|h_{3}|}=0,\nonumber \\&\partial _{i}w_{j}=\partial _{j}w_{i},n_{i}^{*}=0,\partial _{i}n_{j}=\partial _{j}n_{i}. \end{aligned}$$
(53)

We must consider additional constraints on the data \((\Psi ,\Upsilon ),\) or \(( \tilde{\Psi },\ \Lambda ),\) and non-zero integration functions \(\, _{1}n_{j}(x^{k})\) but \(\, _{2}n_{k}(x^{i})=0.\)

To generate explicit solutions, we can consider any functional dependence \(H= \widetilde{\Psi }[\Psi ],\) for which

$$\begin{aligned} \mathbf {e}_{i}H=(\partial _{i}-w_{i}\partial _{4})H=\frac{\partial H}{ \partial \Psi }(\partial _{i}-w_{i}\partial _{4})\Psi \equiv 0; \end{aligned}$$

see (49). For instance, \(H=\widetilde{\Psi }=\ln \sqrt{|\ h_{3}|}\) results in \(\mathbf {e}_{i}\ln \sqrt{|\ h_{3}|}=0.\) If we work with classes of generating functions \(\Psi =\check{\Psi }\) for which

$$\begin{aligned} (\partial _{i}\check{\Psi })^{*}=\partial _{i}(\check{\Psi }^{*}), \end{aligned}$$
(54)

we obtain \(w_{i}^{*}=\mathbf {e}_{i}\ln |\check{\Psi }^{*}|.\) For a given functional dependence \(h_{4}[\Psi ,\Upsilon ],\) we can express

$$\begin{aligned} \mathbf {e}_{i}\ln \sqrt{|\ h_{4}|}=\mathbf {e}_{i}[\ln |\check{\Psi }^{*}|-\ln \sqrt{|\ \Upsilon |}] \end{aligned}$$

(we used the property \(\mathbf {e}_{i}\check{\Psi }=0)\). As a result,

$$\begin{aligned} w_{i}^{*}=\mathbf {e}_{i}\ln \sqrt{|\ h_{4}|} \end{aligned}$$

if \(\mathbf {e}_{i}\ln \sqrt{|\ \Upsilon |}=0.\) This is possible for any \( \Upsilon =\mathrm{const}, \) or any effective source expressed as a functional \(\ \Upsilon (x^{i},t)=\ \Upsilon [\check{\Psi }].\)

The conditions \(\partial _{i}w_{j}=\partial _{j}w_{i}\) can be solved by any function \(\check{A}=\check{A}(x^{k},t)\) for which

$$\begin{aligned} w_{i}=\check{w}_{i}=\partial _{i}\check{\Psi }/\check{\Psi }^{*}=\partial _{i}\check{A}. \end{aligned}$$
(55)

This is a system of first order PDEs which allows one to find a function \(\check{ A}[\check{\Psi }]\) if a functional \(\check{\Psi }\) is prescribed. For the second set of N-coefficients, we choose \(\, _{1}n_{j}(x^{k})=\partial _{j}n(x^{k})\) for a function \(n(x^{k}).\)

Finally, we conclude that we can generate off-diagonal torsionless solutions of the Einstein equations (12) by choosing certain subclasses of generating functions and effective sources in (51), when

$$\begin{aligned} \check{\Upsilon }=\Upsilon (x^{i},t)=\ \Upsilon [\check{\Psi } ],\quad w_{i}=\partial _{i}\check{A}[\check{\Psi }],\quad n_{i}=\partial _{i}n, \end{aligned}$$
(56)

and the generating function \(\Psi =\check{\Psi }\) and “associated” \(\check{A}\) are subjected to the conditions (54) and (55).

4.3 General solutions for (non) holonomic Einstein manifolds

Summarizing the results obtained in previous subsections, we can construct the quadratic linear elements for generic off-diagonal metrics defining 4-d Einstein spaces with an effective cosmological constant \(\Lambda \) and nonholonomic deformations to general sources of type (42).

4.3.1 Solutions with nonholonomically induced torsion

The quadratic line elements determined by coefficients (51), with label torsRs,  are written

$$\begin{aligned}&\mathrm{d}s_{torsRs}^{2} =g_{\alpha \beta }(x^{k},t)\mathrm{d}u^{\alpha }\mathrm{d}u^{\beta }=e^{\psi }[(\mathrm{d}x^{1})^{2}+(\mathrm{d}x^{2})^{2}] \nonumber \\&\quad +\,\omega ^{2}\frac{\tilde{\Psi }^{2}}{4\Lambda }\left[ dy^{3}+\left( \, _{1}n_{k}+_{2}\widetilde{n}_{k}\int \mathrm{d}t\frac{(\tilde{\Psi }^{*})^{2}}{ \tilde{\Psi }^{3}\ \Theta }\right) \mathrm{d}x^{k}\right] ^{2}\nonumber \\&\quad +\,\omega ^{2}\frac{( \tilde{\Psi }^{*})^{2}}{\Theta }\ \left[ \mathrm{d}t+\frac{\partial _{i}\Theta }{ \Theta ^{*}}\mathrm{d}x^{i}\right] ^{2}. \end{aligned}$$
(57)

Fixing in (57) an effective cosmological constant \(\Lambda \ne 0,\) an “associated” generating function \(\tilde{\Psi }(x^{k},t),\) for which \( \tilde{\Psi }^{*}\ne 0,\) and a generating source \(\Theta (x^{k},t),\) for which \(\Theta ^{*}\ne 0,\) we generate exact solutions of the nonholonomic Einstein equations (13). Using Eq. (52), we compute the effective source \(\Upsilon =\) \(\Theta ^{*}/( \widetilde{\Psi }^{2})^{*}.\) Such effective values and Eq. (45) determine a generating function \(\Psi :=e^{\varpi }\) (43) computed from \(\Psi ^{2}=\Lambda ^{-1}(\widetilde{\Psi }^{2}|\Upsilon |+\int \mathrm{d}t\widetilde{\Psi }^{2}|\Upsilon |^{*}).\) We can express equivalently the N-adapted coefficients \(h_{a},n_{i}\) and \(w_{i}\) as functionals of \( (\Psi ,\Upsilon )\) following Eq. (51). The h-metric \( e^{\psi }\) is any solution of the 2-d Poisson equation (46) with source \(\widetilde{\Upsilon }.\) The vertical conformal factor \(\omega (x^{i},y^{3},t)\) can be determined as a solution of (50). The subclass of solutions (57) with Killing symmetry on \(\partial _{3}\) can be extracted for \(\omega ^{2}=1.\)

It should be noted that the effective source \(\Upsilon \) determines an effective source (52) up to a class of frame transforms \(e_{\ \alpha }^{\alpha ^{\prime }}(x^{i},y^{a}).\) Such coefficients must be defined from a system of quadratic algebraic equations for any prescribed \(\Upsilon _{\beta }^{\alpha }=[ \ \tilde{\Upsilon }\delta _{j}^{i},\ \Upsilon \delta _{b}^{a}] .\) The coefficients \(e_{\ \alpha }^{\alpha ^{\prime }} \) and the integration functions \(\, _{1}n_{k}(x^{i})\) and \(_{2}\widetilde{n }_{k}(x^{i})\) can be defined in explicit form if we choose the respective boundary/asymptotic conditions, or solve (for another type of solutions) the Cauchy problem.

The solutions (57) can be parameterized in certain forms generating nonholonomically deformed black ellipsoid, black hole, wormhole, solitonic, and inhomogeneous solutions in various MGTs and in GR; see the discussion and examples in [41, 42, 47, 48].

4.3.2 LC-varieties for effective Einstein manifolds

If the generating/effective functions and sources are subjected to the LC-conditions (54)–(56), we obtain the quadratic linear elements

$$\begin{aligned} ds_{LCRs}^{2}= & {} g_{\alpha \beta }(x^{k},t)\mathrm{d}u^{\alpha }\mathrm{d}u^{\beta }=e^{\psi }[(\mathrm{d}x^{1})^{2}+(\mathrm{d}x^{2})^{2}] \nonumber \\&+\,\omega ^{2}[\check{\Psi }]\frac{(\tilde{\Psi }[\check{\Psi }])^{2}}{4\Lambda } [ \mathrm{d}y^{3}+(\partial _{k}n)\ \mathrm{d}x^{k}] ^{2}\nonumber \\&+\,\omega ^{2}[\check{\Psi }] \frac{(\tilde{\Psi }^{*}[\check{\Psi }])^{2}}{\Theta [\check{\Psi }]} \ [ \mathrm{d}t+(\partial _{i}\ \check{A}[\check{\Psi }])\mathrm{d}x^{i}] ^{2}. \end{aligned}$$
(58)

The coefficients of these generic off-diagonal metrics also generate exact solutions of (12) with effective source (52) but with zero torsion. Such solutions can be modeled equivalently in GR using the LC-connection \(\nabla .\)

There are such generating functions and sources when, for instance, black ellipsoid/hole de Sitter configurations can be extracted from (58) in the limit of certain small off-diagonal deformations; for details see [41]. We note that the metrics are generic off-diagonal if the anholonomic coefficients \(\ W_{\alpha \beta }^{\gamma }\) (5) are not zero.

4.4 Off-diagonal deformations of physically important solutions in GR

In this section, we present four classes of exact solutions generated by the AFDM. The are constructed explicit examples of generic off-diagonal metrics (58) with parameterizations of the generating and integration functions and constants, and corresponding sources, when certain “prime” metrics in GR are transformed into “target” metrics with modified physical constants (running type, or effective polarizations) and/or deformed horizons and/or self-consistent interactions, for instance, with gravitational solitonic waves.

4.4.1 Ellipsoid Kerr–de Sitter configurations

Let us consider a parameterization of 4-coordinated like in the Kerr geometry (we cite [39], for a review of physically important solutions, and [42], for nonholonomic deformations of such solutions). The coordinates \(x^{k^{\prime }}=x^{k^{\prime }}(r,\vartheta ),\) \( y^{3}=t,y^{4}=\varphi ,\) when \(u^{\alpha }=(x^{i^{\prime }},t,\varphi ).\) The prime metric is taken to be the Kerr solution written in the so-called Boyer–Linquist coordinates \((r,\vartheta ,\varphi ,t),\) for \(r=m_{0}(1+p \widehat{x}_{1}),\widehat{x}_{2}=\cos \vartheta .\) The parameters pq are related to the total black hole mass, \(m_{0}\) and the total angular momentum, \(am_{0},\) for the asymptotically flat, stationary, and axisymmetric Kerr spacetime. The formulas \(m_{0}=Mp^{-1}\) and \(a=Mqp^{-1}\) when \( p^{2}+q^{2}=1\) imply \(m_{0}^{2}-a^{2}=M^{2}.\) In such variables, the vacuum Kerr solution can be written

$$\begin{aligned} \mathrm{d}s_{[0]}^{2}= & {} (\mathrm{d}x^{1^{\prime }})^{2}+(\mathrm{d}x^{2^{\prime }})^{2}+\overline{A}( \mathbf {e}^{3^{\prime }})^{2}+(\overline{C}-\overline{B}^{2}/\overline{A})( \mathbf {e}^{4^{\prime }})^{2}, \nonumber \\ \mathbf {e}^{3^{\prime }}= & {} \mathrm{d}t+\mathrm{d}\varphi \overline{B}/\overline{A} =\mathrm{d}y^{3^{\prime }}-\partial _{i^{\prime }}(\widehat{y}^{3^{\prime }}+\varphi \overline{B}/\overline{A})\mathrm{d}x^{i^{\prime }},\nonumber \\ \mathbf {e}^{4^{\prime }}= & {} \mathrm{d}y^{4^{\prime }}=\mathrm{d}\varphi . \end{aligned}$$
(59)

We can consider any coordinate functions

$$\begin{aligned}&x^{1^{\prime }}(r,\vartheta ),\ x^{2^{\prime }}(r,\vartheta ),\ y^{3^{\prime }}=t+\widehat{y}^{3^{\prime }}(r,\vartheta ,\varphi )+\varphi \overline{B}/ \overline{A},y^{4^{\prime }}\\&\quad =\varphi ,\ \partial _{\varphi }\widehat{y} ^{3^{\prime }}=-\overline{B}/\overline{A}, \end{aligned}$$

for which \((\mathrm{d}x^{1^{\prime }})^{2}+(\mathrm{d}x^{2^{\prime }})^{2}=\Xi ( \Delta ^{-1}\mathrm{d}r^{2}+\mathrm{d}\vartheta ^{2}) \), and the coefficients are

$$\begin{aligned} \overline{A}= & {} -\overline{\Xi }^{-1}(\overline{\Delta }-a^{2}\sin ^{2}\vartheta ),\quad \\&\overline{B}=\overline{\Xi }^{-1}a\sin ^{2}\vartheta [ \overline{\Delta }-(r^{2}+a^{2})], \\ \overline{C}= & {} \overline{\Xi }^{-1}\sin ^{2}\vartheta [ (r^{2}+a^{2})^{2}-\overline{\Delta }a^{2}\sin ^{2}\vartheta ] ,\quad \text{ and } \\ \overline{\Delta }= & {} r^{2}-2m_{0}+a^{2},\quad \overline{\Xi }=r^{2}+a^{2}\cos ^{2}\vartheta . \end{aligned}$$

The quadratic linear element (59) with prime data

$$\begin{aligned}&\mathring{g}_{1} =1,\mathring{g}_{2}=1,\mathring{h}_{3}=-\rho ^{2}Y^{-1}, \mathring{h}_{4}=Y,\mathring{N}_{i}^{a}=\partial _{i}\widehat{y}^{a}, \\&\quad (\text{ or } \mathring{g}_{1^{\prime }} =1,\mathring{g}_{2^{\prime }}=1, \mathring{h}_{3^{\prime }}=\overline{A},\mathring{h}_{4^{\prime }}=\overline{ C}-\overline{B}^{2}/\overline{A}, \\&\mathring{N}_{i^{\prime }}^{3} =\mathring{n}_{i^{\prime }}=-\partial _{i^{\prime }}(\widehat{y}^{3^{\prime }}+\varphi \overline{B}/\overline{A}), \mathring{N}_{i^{\prime }}^{4}=\mathring{w}_{i^{\prime }}=0) \end{aligned}$$

define solutions of the vacuum Einstein equations parameterized in the form (13) and (14) with zero sources.

We construct a subclass of solutions with rotoid configurations for generating functions \(\ \tilde{\Psi }=e^{\varpi }[\,\widetilde{ \Lambda }/2\ ^{\mu }\tilde{\Lambda }+\underline{\zeta }\sin (\omega _{0}\varphi +\varphi _{0})],\) for \(\tilde{\Psi }=e^{\varpi },\)

$$\begin{aligned} \mathrm{d}s^{2}= & {} e^{\psi (x^{k^{\prime }})}(1+\varepsilon \chi (x^{k^{\prime }}))[(\mathrm{d}x^{1^{\prime }})^{2}+(\mathrm{d}x^{2^{\prime }})^{2}] \nonumber \\&-\frac{e^{2\varpi }}{4|\Lambda |}\overline{A}[1+2\varepsilon \underline{ \zeta }\sin (\omega _{0}\varphi +\varphi _{0})]\nonumber \\&\times [\mathrm{d}y^{3^{\prime }}+\left( \partial _{k^{\prime }}\ ^{\eta }n(x^{i^{\prime }})-\partial _{k^{\prime }}( \widehat{y}^{3^{\prime }}+\varphi \frac{\overline{B}}{\overline{A}})\right) \mathrm{d}x^{k^{\prime }}]^{2} \nonumber \\&+\frac{(\varpi ^{*})^{2}}{\ \Lambda }\left( \overline{C}-\frac{\overline{B} ^{2}}{\overline{A}}\right) \nonumber \\&\times \bigg [1+\varepsilon \bigg (\frac{\Lambda }{\widetilde{\lambda }} \partial _{\varphi }\varpi \,\nonumber \\&+2\partial _{4}\varpi \underline{\zeta }\sin (\omega _{0}\varphi +\varphi _{0})+2\omega _{0}\ \underline{\zeta } \cos (\omega _{0}\varphi +\varphi _{0})\bigg )\bigg ] \nonumber \\&\times [\mathrm{d}\varphi +(\partial _{i^{\prime }}\ \widetilde{A}+\varepsilon \partial _{i^{\prime }}\ \ ^{1}\check{A})\mathrm{d}x^{i^{\prime }}]^{2}. \end{aligned}$$
(60)

Such metrics have a Killing symmetry on \(\partial /\partial t\) and are completely defined by a generating function \(\varpi (x^{k^{\prime }},\varphi )\) and the source and \(\ \Lambda .\) They describe \(\varepsilon \)-deformations of Kerr–de Sitter black holes into ellipsoid configurations with effective (polarized) cosmological constants determined by possible geometric flows and generic off-diagonal interactions. If the zero torsion conditions are satisfied (like in (58)), such metrics can be modeled in GR.

4.4.2 Nonholonomically deformed wormhole configurations

Let us consider a different class of stationary configurations which define a wormhole solutions and their off-diagonal deformations. We begin with a diagonal prime wormhole metric,

$$\begin{aligned} \mathbf {\mathring{g}}= & {} \mathring{g}_{i}(x^{k})\mathrm{d}x^{i}\otimes \mathrm{d}x^{i}+ \mathring{h}_{a}(x^{k})\mathrm{d}y^{a}\otimes \mathrm{d}y^{a} \\= & {} [1-b(r)/r]^{-1}\mathrm{d}r\otimes \mathrm{d}r+r^{2}\mathrm{d}\theta \otimes \mathrm{d}\theta \\&-\,e^{2B(r)}\mathrm{d}t\otimes \mathrm{d}t+r^{2}\sin ^{2}\theta \mathrm{d}\varphi \otimes \mathrm{d}\varphi , \end{aligned}$$

where B(r) and b(r) are called, respectively, the red-shift and form functions. On holonomic wormhole solutions, we refer the reader to [52,53,54,55]. There the local coordinates \(u^{\alpha }=(r,\theta ,t,\varphi )\) are used. The radial coordinate has the range \(r_{0}\le r<a.\) The minimum value \(r_{0}\) is for the wormhole throat. The constant a is the distance at which the interior spacetime joins to an exterior vacuum solution (\(a\rightarrow \infty \) for specific asymptotically flat wormhole geometries). The coefficients of the diagonal stress-energy tensor

$$\begin{aligned} \mathring{T}_{~\nu }^{\mu }=\mathrm{diag}[~^{r}p=\tau (r),~^{\theta }p\!=\!p(r),~^{\varphi }p=p(r),~^{t}p=\rho (r)] \end{aligned}$$

are subjected to certain conditions in order to generate wormhole solutions of the Einstein equations in GR.

Wormhole metrics are constructed to possess conformal symmetry determined by a vector \(\mathbf {X}=\{X^{\alpha }(u)\},\) when considering the Lie derivative,

$$\begin{aligned} X^{\alpha }\partial _{\alpha }\mathring{g}_{\mu \nu }+\mathring{g}_{\alpha \nu }\partial _{\mu }X^{\alpha }+\mathring{g}_{\alpha \mu }\partial _{\nu }X^{\alpha }=\sigma \mathring{g}_{\mu \nu }, \end{aligned}$$

where \(\sigma =\sigma (u)\) is the conformal factor. A class of such solutions is parameterized by

$$\begin{aligned} B(r)= & {} \frac{1}{2}\ln (C^{2}r^{2})-\kappa \\&\times \int r^{-1}( 1-b(r)/r) ^{-1/2}\mathrm{d}r,~b(r)=r[1-\sigma ^{2}(r)], \\ \tau (r)= & {} \frac{1}{\kappa ^{2}r^{2}}(3\sigma ^{2}-2\kappa \sigma -1),~p(r)\\= & {} \frac{1}{\kappa ^{2}r^{2}}(\sigma ^{2}-2\kappa \sigma +\kappa ^{2}+2r\sigma \sigma ^{\bullet }), \\ \rho (r)= & {} \frac{1}{\kappa ^{2}r^{2}}(1-\sigma ^{2}-2r\sigma \sigma ^{\bullet }). \end{aligned}$$

These data generate “diagonal” wormhole configurations determined by “exotic” matter because the null energy condition (NEC), \(\mathring{T}_{\mu \nu }k^{\mu }k^{\nu }\ge 0\) (\(k^{\nu }\) is any null vector), is violated.

In this section, we analyze wormhole configurations which match the interior geometries to an exterior de Sitter one which (in general) can be also determined by an off-diagonal metric. The exotic matter and effective matter configurations are considered to be restricted to spatial distributions in the throat neighborhood which limit the dimension of locally isotropic and/or anisotropic wormhole to be not arbitrarily large. The Schwarzschild–de Sitter (SdS) metric,

$$\begin{aligned} \mathrm{d}s^{2}=q^{-1}(r)(\mathrm{d}r^{2}+r^{2}\ \mathrm{d}\theta ^{2})+r^{2}\sin ^{2}\theta \ \mathrm{d}\varphi ^{2}-q(r)\ \mathrm{d}t^{2}, \end{aligned}$$

can be re-parameterized for any \((x^{1}(r,\theta ),x^{2}(r,\theta ),y^{3}=\varphi ,y^{4}=t)\) in order to express \(q^{-1}(r)(\mathrm{d}r^{2}+r^{2}\ \mathrm{d}\theta ^{2})=e^{\mathring{\psi }{(x^{k})}}[(\mathrm{d}x^{1})^{2}+(\mathrm{d}x^{2})^{2}].\) Such a metric defines two real static solutions of the Einstein equations with cosmological constant \(\Lambda \) if \(M<1/3\sqrt{|\Lambda |},\) for \(q(r)=1-2 \overline{M}(r)/r,\overline{M}(r)=M+\Lambda r^{3}/6,\) where M is a constant mass parameter. For diagonal configurations, we can identify \( \Lambda \) with the effective cosmological constant.

The next step is to consider conformal, ellipsoid, and/or solitonic/toroidal deformations related in certain limits to the SdS metric written in the form

$$\begin{aligned} ~_{\Lambda }\mathbf {g}&=\mathrm{d}\xi \otimes \mathrm{d}\xi +r^{2}(\xi )\ \mathrm{d}\theta \otimes \mathrm{d}\theta +r^{2}(\xi )\sin ^{2}\theta \ \mathrm{d}\varphi \\&\otimes \mathrm{d}\varphi -q(\xi )\ \mathrm{d}t\otimes \ \mathrm{d}t. \end{aligned}$$

Local coordinates \(x^{1}=\xi =\int \mathrm{d}r/\sqrt{\left| q(r)\right| },x^{2}=\vartheta ,y^{3}=\varphi ,y^{4}=t,\) are used for a system of h-coordinates when \((r,\theta )\rightarrow (\xi ,\vartheta )\) with \(\xi \) and \(\vartheta \) of length dimension. The data for this primary metric are written

$$\begin{aligned}&\mathring{g}_{i}=\mathring{g}_{i}(x^{k})=e^{\mathring{\psi }{(x^{k})}},\quad \mathring{h}_{3}=r^{2}(x^{k})\sin ^{2}\theta (x^{k}),\\&\quad \mathring{h}_{4}=-q(r{ (x^{k})}),\quad \mathring{w}_{i}=0,\quad \mathring{n}_{i}=0. \end{aligned}$$

Off-diagonal deformations on a small parameter \(\varepsilon \) of the wormhole metrics are described by quadratic elements of type (58), when the generating and integration functions are chosen

$$\begin{aligned} \mathbf {ds}^{2}= & {} e^{\tilde{\psi }(\widetilde{\xi },\theta )}(\mathrm{d}\widetilde{ \xi }^{2}+\ \mathrm{d}\vartheta ^{2})-\frac{e^{2\widetilde{\varpi }}}{4\Lambda } [1+\varepsilon \overline{\chi }_{3}(\widetilde{\xi },\varphi )]e^{2B( \widetilde{\xi })}\nonumber \\&\times \,[\mathrm{d}t+\partial _{\widetilde{\xi }}(\ ^{\eta }n+\varepsilon \overline{n})~\mathrm{d}\widetilde{\xi }+\partial _{\vartheta }(\ ^{\eta }n+\varepsilon \overline{n})~\mathrm{d}\vartheta ]^{2} \nonumber \\&+\,\frac{[\partial _{\varphi }\widetilde{\varpi }]^{2}}{\Lambda } \left( 1+\varepsilon \frac{\partial _{\varphi }[\overline{\chi }_{3}\widetilde{ \varpi }]}{\partial _{\varphi }\widetilde{\varpi }}\right) \ ^{0}\overline{h} _{4}\nonumber \\&\times \,[\mathrm{d}\varphi +\partial _{\widetilde{\xi }}(\ ^{\eta }\widetilde{A} +\varepsilon \overline{A})\mathrm{d}\widetilde{\xi }+\partial _{\vartheta }(\ ^{\eta } \widetilde{A}+\varepsilon \overline{A})\mathrm{d}\vartheta ]^{2}. \nonumber \\ \end{aligned}$$
(61)

We prescribe such generating functions \(\widetilde{\varpi }\) and effective source \(\Lambda \) so that the effective polarization functions can be approximated \(\widetilde{\eta }_{a}\simeq 1\) and \(^{\eta }\widetilde{A}\) and \(\ ^{\eta }n\) are “almost constant”, with respect to certain systems of radial coordinates. For such conditions, the metric (61) mimics small rotoid wormhole like configurations with off-diagonal terms and possible geometric flow modifications of the diagonal coefficients; see [56] for further details. It is possible to choose such integration functions and constants so that this class of stationary solutions define wormhole-like metrics depending generically on three space coordinates with self-consistent “embedding” in an off-diagonal GR background, for \(\ ^{0} \overline{h}_{3}=r^{2}(\widetilde{\xi })\sin ^{2}\theta (\widetilde{\xi } ,\vartheta ),\) where \(\ \widetilde{\xi }=\int \mathrm{d}r/\sqrt{|1-b(r)/r|}\) and \(B( \widetilde{\xi })\) are determined by a wormhole metric.

4.4.3 Solitonic waves for inhomogeneous cosmological solutions

For simplicity, we can consider solutions of type (58) generated by a nonlinear radial (solitonic, with left s-label) generating function \( \tilde{\Psi }=\ ^{s}\tilde{\Psi }(r,t)=4\arctan e^{q\sigma (r-vt)+q_{0}}\) and construct a metric

$$\begin{aligned} \mathbf {ds}^{2}= & {} e^{\psi (r,\theta )}(\mathrm{d}r^{2}+\ \mathrm{d}\theta ^{2})+~\frac{\ ^{s} \tilde{\Psi }^{2}}{4\ \Upsilon }\mathring{h}_{3}(r,\theta )\mathrm{d}\varphi ^{2}\nonumber \\&- \frac{(\partial _{t}\ ^{s}\tilde{\Psi })^{2}}{\Upsilon \ ^{s}\tilde{\Psi }^{2}} \mathring{h}_{4}(r,t)[\mathrm{d}t+(\partial _{r}\ \widetilde{A})\mathrm{d}r]^{2}. \end{aligned}$$
(62)

In this metric, for simplicity, we fixed \(n(r,\theta )=0\) and consider \(\widetilde{A}(r,t)\) to be defined as a solution of \(\ ^{s}\tilde{\Psi } ^{\bullet }/\ ^{s}\tilde{\Psi }^{*}=\partial _{r}\ \widetilde{A}\) and \( \mathring{h}_{a}\) are given by homogeneous cosmology model data. The generating function is just a 1-soliton solution of the sine-Gordon equation \(\ ^{s}\tilde{\Psi }^{**}-\ ^{s}\tilde{\Psi }^{\bullet \bullet }+\sin \ ^{s}\tilde{\Psi }=0\). For any class of small polarizations with \(\eta _{a}\sim 1),\) we can consider the source \((\ ^{m}\Upsilon +\ ^{\alpha }\Upsilon )\) to be polarized by \(\ ^{s}\tilde{\Psi }^{-2}\) when \( h_{3}\sim \mathring{h}_{3}\) and \(h_{4}\sim \mathring{h}_{4}(\ ^{s}\tilde{\Psi }^{*})^{2}/\ ^{s}\tilde{\Psi }^{4}\) with an off-diagonal term \(\partial _{r}\ \widetilde{A}\) resulting in a stationary solitonic Universe. If we consider that \((\partial _{\widehat{R}}\widehat{f})^{-1}=\ ^{s}\tilde{\Psi } ^{-2}\), we can model \(\widehat{f}\)-interactions via off-diagonal interactions and “gravitational polarizations”.

In the absence of matter, the off-diagonal cosmology is completely determined by \(\ \Upsilon \) related to an effective cosmological constant induced by geometric flows. Such configurations can be determined alternatively using distribution of matter fields when contributions from massive gravity are with small anisotropic polarization; for details see [57,58,59,60].

4.4.4 Off-diagonal deformations of FLRW metrics and gravitational solitonic waves

Using a redefined time-like coordinate \(\widehat{t},\) when \(t=t(x^{i}, \widehat{t})\), \(\sqrt{|h_{4}|}\partial t/\partial \widehat{t}\), for a scale factor \({\widehat{a}}(x^{i},\widehat{t}),\) the d-metric (62) can be represented in the form

$$\begin{aligned} \mathrm{d}s^{2}=\widehat{a}^{2}(x^{i},\widehat{t})[\eta _{i}(x^{k},\widehat{t} )(\mathrm{d}x^{i})^{2}+\widehat{h}_{3}(x^{k},\widehat{t})(\mathbf {e}^{3})^{2}-( \widehat{\mathbf {e}}^{4})^{2}], \end{aligned}$$
(63)

where \(\eta _{i}=\widehat{a}^{-2}e^{\psi },\widehat{a}^{2}\widehat{h} _{3}=h_{3},\mathbf {e}^{3}=\mathrm{d}y^{3}+\partial _{k}n~\mathrm{d}x^{k},\widehat{\mathbf {e}} ^{4}=\mathrm{d}\widehat{t}+\sqrt{|h_{4}|}(\partial _{i}t+w_{i}).\) This is a non-stationary solution defining inhomogeneous cosmological spacetimes.

We can model small off-diagonal deformations of the Friedmann–Lemaître–Roberstson–Walker (FLRW) metric parameterized by a small parameter \( \varepsilon ,\) with \(0\le \varepsilon <1,\) when

$$\begin{aligned}&\eta _{i}\simeq 1+\varepsilon \chi _{i}(x^{k},\widehat{t}),\quad \partial _{k}n\simeq \varepsilon \widehat{n}_{i}(x^{k}),\\&\quad \quad \sqrt{|h_{4}|}(\partial _{i}t+w_{i})\simeq \varepsilon \widehat{w}_{i}(x^{k},\widehat{t}). \end{aligned}$$

It is possible to choose such generating functions and sources when \( \widehat{a}(x^{i},\widehat{t})\rightarrow \) \(\widehat{a}(t),\widehat{h} _{3}(x^{i},\widehat{t})\rightarrow \widehat{h}_{3}(\widehat{t})\) etc. This results in new classes of solutions even in the diagonal limits because of the generic nonlinear and the nonholonomic character of off-diagonal systems in GR and geometric flow evolution theories. For \(\varepsilon \rightarrow 0\) and \(\widehat{a}(x^{i},\widehat{t})\rightarrow \) \(\widehat{a}(t),\) we obtain scaling factors which are very different from those in FLRW cosmology. They mimic such cosmological models with redefined interaction parameters and possible small off-diagonal deformations of the cosmological evolution both in GR, MGTs, and Ricci flow evolution; see [59, 60].

A metric (62) defines solitonic waves along the coordinate \(x^{1}\) if we take

$$\begin{aligned} \check{\Psi }=\ ^{s}\Psi (x^{1},t)=4\arctan e^{q_{1}q_{2}(x^{1}-\xi t)+q_{0}}, \end{aligned}$$

for certain constants \(q_{0},q_{1}\) and \(q_{2}=1/\sqrt{|1-\xi ^{2}|},\) and put, for simplicity, \(n=0.\) The function \(\widetilde{A}\) is a solution of \( \ ^{s}\Psi ^{*}\widetilde{A}^{\bullet }=\ ^{s}\Psi ^{\bullet }.\) For such conditions, the generating function induces 1-soliton waves as a solution of the sine-Gordon equation,

$$\begin{aligned} \ ^{s}\Psi ^{**}-\ ^{s}\Psi ^{\bullet \bullet }+\sin (\ ^{s}\Psi )=0, \end{aligned}$$

in the v-subspace together with self-consistent off-diagonal propagation of such waves via \(w_{i}=\widetilde{A}^{\bullet }.\) More general 3-d solitonic gravitational waves are possible which can propagate self-consistently in Minkowski spacetime or in a FLRW background. To generate such solutions we need to take a generating function \(\ ^{s}\Psi (x^{1},x^{2},t)\) which is a solution of the Kadomtsev–Petviashvili (KdP) equation [61,62,63].

If we consider that in the ansatz of type (58) \(y^{3}=t\) (time-like coordinate) and \(y^{4}=\varphi \) (an angular space coordinate), we construct stationary solutions with Killing symmetry on \(\partial _{3};\) see the details in [41, 42]. Such solutions in GR and geometric flow theory include as particular cases black ellipsoid configurations, Taub NUT configurations, and other types of Coulomb-like gravitational fields. In particular, for small \(\varepsilon \)-deformations they include black hole and black ellipsoid solutions in GR and certain classes of modified theories. Such gravitational configurations are characterized by different thermodynamic values in a W-thermodynamic model of geometric flows and their stationary Ricci soliton configurations.

5 Nonholonomic thermodynamics of gravitational fields

There is a standard theory of the thermodynamics of the black hole (BH) solutions originally elaborated by Bekenstein–Hawking constructions for stationary solutions in gravity theories (for a review, see [39] and the references therein). Following the formalism of double fibrations, the BH thermodynamics can be derived for a very special example of Lyapunov-type functionals on 3-d hypersurfaces with further \(2+1\) splitting (or \(2+1+1\) holonomic fibrations with horizon configurations and corresponding asymptotic conditions). In general, we cannot apply the BH thermodynamics to study general gravitational configurations, for instance, models of inhomogeneous cosmology and/or other spacetime configurations in GR and MGTs. Certain ideas were proposed to formulate more general thermodynamic descriptions of gravitational fields; see [40] (for the so-called CET model) and the references therein.

In this section, we develop a W-functional thermodynamic model on a 3-d closed space-like hypersurface in GR, following the standard theory of Ricci flows. In general, such an approach is very different from that considered in standard BH thermodynamics. As a matter of principle, it is possible to establish certain equivalence conditions (with very special holonomic vacuum gravity configurations or diagonal de Sitter BH solutions) when the W-functional approach will give the same result as in Bekenstein–Hawking BH thermodynamics. The goal of this section is to elaborate on models on nonholonomic thermodynamics of gravitational fields using 4-d generalizations of Perelman’s W-functionals. We shall also speculate how such constructions can be related to the CET approach if certain additional conditions are imposed. In this work, we restrict our approach to general relativistic flows of the 3-d to the 4-d configuration to a special class of models when the effective temperature \(\beta ^{-1}(t)\) depends, in general, on the time-like coordinate t.

5.1 General relativistic models of W-entropy and geometric flow evolution

We emphasize that Perelman’s functional \(\, _{\shortmid }\mathcal {W}\) (30) is in a sense analogous to minus the entropy. Similar values can be considered for various metric compatible nonholonomic Ricci flows and (non) commutative, fractional derivative and other type modifications [8, 9, 20, 21]. This allows us to associate certain analogous thermodynamical values characterizing a (non) holonomic modified Ricci flow evolution of the metrics with a local Euclidean structure and generalized connections. Using \(\, _{\shortmid }\mathcal {W}\) for \(( \mathbf {q},\, _{\shortmid }\widehat{\mathbf {D}})\) and its 4-d generalization \(\widehat{\mathcal {W}}\) in the form (34), or (36) for \((\ \mathbf {g},\widehat{\mathbf {D}}),\) we can construct relativistic geometric evolution models with generalized Hamilton equations of type (38).

For 4-d general relativistic configurations, we cannot provide a standard statistical thermodynamic interpretation. We shall have to elaborate on relativistic hydrodynamical type generalizations (see Sect. 5.3). Here we emphasize that we can always characterize the geometric flows by analogous thermodynamic systems on corresponding families of 3-d closed hypersurfaces \(\ \widehat{\Xi }_{t}\) using a nonholonomically deformed entropy \(\, _{\shortmid }\widehat{\mathcal {W}}.\) Let us consider the standard partition function,

$$\begin{aligned} \breve{Z}=\exp \left\{ \int _{\widehat{\Xi }_{t}}M\sqrt{|q_{\grave{\imath } \grave{j}}|}\mathrm{d}\grave{x}^{3}[-\breve{f}+n]~\right\} , \end{aligned}$$

for the conditions stated for definition of (29) and (30). This allows us to compute the main thermodynamical values for the Levi-Civita connection \(\, _{\shortmid }\nabla \) [1, 12] and \(n=3.\) A statistical model can be elaborated for any prescribed partition function \(Z=\int \exp (-\beta E)\mathrm{d}\omega (E)\) for a corresponding canonical ensemble at temperature \(\beta ^{-1}\), being defined by the measure taken to be the density of states \(\omega (E).\) The standard thermodynamical values are computed for the average energy, \(\mathcal {E}=\ \langle E\rangle :=-\partial \log Z/\partial \beta ,\) the entropy \(S:=\beta \langle E\rangle +\log Z\), and the fluctuation \(\sigma :=\langle ( E-\langle E\rangle ) ^{2}\rangle =\partial ^{2}\log Z/\partial \beta ^{2}.\)

To elaborate the constructions in N-adapted form we consider a family of \( q_{\grave{\imath }\grave{j}}(\tau (\chi ))\), with \(\partial \tau /\partial \chi =-1\) being a real re-parametrization of \(\chi .\) For 4-d, it can be considered as a time-like parameter, \(\chi \sim t).\) We change \(\, _{\shortmid }\nabla \rightarrow \, _{\shortmid }\widehat{\mathbf {D}}\) as it is determined by the distortions (27); see similar constructions in [20]. By a corresponding re-scaling \(\breve{f}\rightarrow \tilde{f}\) and \(\tau \rightarrow \tilde{\tau }(t)\) (such a re-scaling is useful if we want to compare thermodynamical values for different linear connections in GR), we compute

$$\begin{aligned}&\, _{\shortmid }\widehat{\mathcal {E}} =-\tilde{\tau }^{2}\int _{\widehat{\Xi }_{t}}M\sqrt{|q_{\grave{\imath }\grave{j}}|}\mathrm{d}\grave{x}^{3}\left( \, _{\shortmid }\widehat{R}+|\, _{\shortmid }\widehat{\mathbf {D}}\tilde{f}|^{2}- \frac{3}{\tilde{\tau }}\right) , \nonumber \\&\, _{\shortmid }\widehat{S} =-\int _{\widehat{\Xi }_{t}}M\sqrt{|q_{\grave{ \imath }\grave{j}}|}\mathrm{d}\grave{x}^{3}[ \tilde{\tau }( \ \, _{\shortmid }\widehat{R}+|\, _{\shortmid }\widehat{\mathbf {D}}\tilde{f} |^{2}) +\tilde{f}-6], \nonumber \\&\, _{\shortmid }\widehat{\sigma } = 2\ \tilde{\tau }^{4}\int _{\widehat{\Xi } _{t}}M\sqrt{|q_{\grave{\imath }\grave{j}}|}\mathrm{d}\grave{x}^{3}\left[ |\, _{\shortmid } \widehat{\mathbf {R}}_{\grave{\imath }\grave{j}}+\, _{\shortmid }\widehat{ \mathbf {D}}_{\grave{\imath }}\, _{\shortmid }\widehat{\mathbf {D}}_{\grave{j}} \tilde{f}-\frac{1}{2\tilde{\tau }}q_{\grave{\imath }\grave{j}}|^{2}\right] .\nonumber \\ \end{aligned}$$
(64)

Such formulas can be considered for 4-d configurations considering the lapse function \(\breve{N}=1\) for N-adapted Gaussian coordinates but in such cases it will be more difficult to compute in explicit form using standard forms of solutions of 4-d physically important equations.

Using the thermodynamical equation (64), we can compute the corresponding average energy, entropy, and fluctuations for the evolution both on redefined parameter \(\tau \) and on a time-like parameter t of any family of closed hypersurfaces all determined by \(\, _{\shortmid }\widehat{\mathbf {D}}\),

$$\begin{aligned} \widehat{\mathcal {E}}(\tau )= & {} \int _{t_{1}}^{t_{2}}\mathrm{d}t\breve{N}(\tau )\, _{\shortmid }\widehat{\mathcal {E}}(\tau ),\ \widehat{\mathcal {S}}(\tau )\nonumber \\= & {} \int _{t_{1}}^{t_{2}}\mathrm{d}t\breve{N}(\tau )\, _{\shortmid }\widehat{S}(\tau ),\ \widehat{\Sigma }(\tau )=\int _{t_{1}}^{t_{2}}\mathrm{d}t\breve{N}(\tau )\, _{\shortmid }\widehat{\sigma }(\tau ). \nonumber \\ \end{aligned}$$
(65)

These formulas are related by distortion equations (27) with corresponding values determined by \(\, _{\shortmid }\nabla \),

$$\begin{aligned} \ ^{\nabla }\mathcal {E}(\tau )= & {} \int _{t_{1}}^{t_{2}}\mathrm{d}t\breve{N}(\tau )\, _{\shortmid }^{\nabla }\mathcal {E}(\tau ),\ ^{\nabla }\mathcal {S}(\tau )\\= & {} \int _{t_{1}}^{t_{2}}\mathrm{d}t\breve{N}(\tau )\, _{\shortmid }^{\nabla }S(\tau ),\ ^{\nabla }{\Sigma }(\tau )\\= & {} \int _{t_{1}}^{t_{2}}\mathrm{d}t\breve{N}(\tau )\, _{\shortmid }^{\nabla }\sigma (\tau ). \end{aligned}$$

Such values with a “hat” or left label \(\nabla \) are different; we have \(\, _{\shortmid }\widehat{\mathbf {D}}\rightarrow \nabla \) for certain topologically nontrivial configurations.

5.2 W-thermodynamic values for exact solutions in GR

For nonholonomic 4-d Lorentz–Ricci solitonic equations defined by systems of type (57), or (58), the generating function \(\Psi (x^{k},t)\) (or \(\widetilde{\Psi }(x^{k},t),\) or \(\check{\Psi }(x^{k},t)\) for torsionless configurations) and the effective source \(\Upsilon (x^{k},t)\) can be prescribed in some forms not depending one on another. The vertical conformal factor \(\omega \) can be an arbitrary function on \((x^{k},y^{3},t)\) subjected to the conditions (50). Such values should be subjected to an additional constraint if we consider that the 4-d solutions for a d-metric \(\mathbf {g}\) encode also 3-d d-metrics \(\mathbf {q}\). The 3-d metric is parameterized in the form (40) with a Ricci flow evolution determined by the Hamilton equations (28) on a relativistic parameter. To construct explicit solutions, we consider the parameter \(\chi =t\) for 3-d evolution equations written with respect to N-adapted frames.

5.2.1 Modified 3-d Ricci flows with induced nonholonomic torsion

On a 4-d Lorentz–Ricci solitonic space determined by a quadratic line element (57), the 3-d Ricci flow evolution equation (28) is written

$$\begin{aligned} \mathbf {q}_{\grave{\imath }\grave{j}}^{*}=-2\, _{\shortmid }\widehat{ \mathbf {R}}_{\grave{\imath }\grave{j}}+\frac{2\grave{r}}{5}\mathbf {q}_{\grave{ \imath }\grave{j}}, \nonumber \end{aligned}$$

for \(\mathbf {q}_{\grave{\imath }\grave{j}}^{*}=\partial _{t}\mathbf {q}_{ \grave{\imath }\grave{j}},\) where \(\, _{\shortmid }\widehat{\mathbf {R}}_{ \grave{\imath }\grave{j}}=\{\widehat{\mathbf {R}}_{1}^{1}\mathbf {q}_{ij}, \widehat{\mathbf {R}}_{3}^{3}\mathbf {q}_{3}\}.\) In non-explicit form, \( \partial _{t}\) is related to the partial derivation on temperature \(\beta ^{-1}(t)\), which in this work can be related to the time-like variable. For components with \(i,j,k\ldots =1,2,\) in a \(3+1\) nonholonomic distribution, these evolution equations, with \(q_{ij}^{*}=0\) and \(q_{3}^{*}\ne 0,\) decouple in the form

$$\begin{aligned}&_{\shortmid }\widehat{\mathbf {R}}_{j}^{i} =\delta _{j}^{i}\frac{\grave{r}}{ 5}=\delta _{j}^{i}\widetilde{\Upsilon }(x^{k}), \nonumber \\&\partial _{t}\ln |\omega ^{2}h_{3}| =-2\widehat{\mathbf {R}}_{3}^{3}+\frac{2 \grave{r}}{5}, \end{aligned}$$
(66)

for \(\widehat{\mathbf {R}}_{3}^{3}=\Upsilon (x^{k},t).\) The h-source \( \widetilde{\Upsilon }(x^{k})\) can be considered as a normalizing factor for the 3-d Ricci flows when \(\grave{r}=\widetilde{\Upsilon }.\) Taking \(h_{3}= \frac{\tilde{\Psi }^{2}}{4\Lambda }\) and introducing the first equation in (66) into the second one, we obtain

$$\begin{aligned} \partial _{t}\ln |\widetilde{\Psi }|=\widetilde{\Upsilon }-\Upsilon -\partial _{t}\ln |\omega |. \end{aligned}$$
(67)

We conclude that 4-d solitons of the Einstein equations describe also 3-d solutions of the normalized N-adapted Hamilton equations if the associated generating function \(\tilde{\Psi },\) the effective sources \(\widetilde{ \Upsilon }-\Upsilon \) and v-conformal factor \(\omega \) are subject to the conditions (67). We can use this for definition of a self-consistent parametrization of the effective source, \(\Upsilon =\) \( \widetilde{\Upsilon }-(\ln |\omega \widetilde{\Psi }|)^{*}.\) Introducing this value in (45), we can find this nonlinear symmetry for re-defining generating functions:

$$\begin{aligned} \Psi ^{2}= & {} \Lambda ^{-1}\left[ \widetilde{\Psi }^{2}|\widetilde{\Upsilon } -(\ln |\omega \widetilde{\Psi }|)^{*}|\right. \\&\left. +\int \mathrm{d}t\widetilde{\Psi }^{2}| \widetilde{\Upsilon }-(\ln |\omega \widetilde{\Psi }|)^{*}|^{*} \right] . \end{aligned}$$

As particular cases, we can consider \(\widetilde{\Upsilon }=\mathrm{const}\) and \( \omega =1\) in order to generate a self-consistent 3-d Ricci flow evolution on a 4-d effective Einstein spaces. In such cases, all geometric and physical objects have the Killing symmetry on \(\partial _{3}.\)

By definition of the generating effective source (52), \(\Theta ^{*}=\ \Upsilon (\widetilde{\Psi }^{2})^{*}.\) We can write Eq. (67) using values \((\omega ,\widetilde{\Psi },\widetilde{ \Upsilon },\Theta ^{*}),\)

$$\begin{aligned} \Theta ^{*}=(\widetilde{\Psi }^{2})^{*}\ \left[ \widetilde{\Upsilon }- \frac{\omega ^{*}}{\omega }\right] -(\widetilde{\Psi }^{*})^{2}. \end{aligned}$$

We conclude that we can model Ricci flows of 3-d metrics \(q_{\grave{\imath } \grave{j}}\) in N-adapted form on effective 4-d Einstein configurations by imposing additional constraints on the generating functions, effective sources, or effective generating sources. Such geometric evolutions are characterized by nontrivial nonholonomic torsion completely defined by \( \mathbf {g.}\)

5.2.2 Modified 3-d Ricci flows of LC-configurations

Torsionless 4-d Einstein configurations (58) are determined by generating functions \(\tilde{\Psi }[\check{\Psi }]\) and sources \(\check{ \Upsilon }=\Upsilon (x^{i},t)=\ \Upsilon [\check{\Psi }]\) subject to the conditions (54)–(56). The 3-d Ricci flow evolution in N-adapted form is described by (67) redefined for generating functions and sources subjected to the condition

$$\begin{aligned} \partial _{t}\ln |\widetilde{\Psi }[\check{\Psi }]|=\widetilde{\Upsilon }- \check{\Upsilon }-\partial _{t}\ln |\omega |. \end{aligned}$$

The LC-geometric flow configurations are characterized by the effective source \( \check{\Upsilon }=\) \(\widetilde{\Upsilon }-(\ln |\omega \widetilde{\Psi }[ \check{\Psi }]|)^{*}\) and an effective generating source \(\check{\Theta } ^{*}=\ \check{\Upsilon }(\widetilde{\Psi }^{2}[\check{\Psi }])^{*}.\)

5.2.3 N-adapted 3-d Ricci flows on exact solutions in GR

Exact solutions of 3-d Hamilton-like equations (66) are considered for 4-d solutions (58). The evolution equation (67) is modified,

$$\begin{aligned} \partial _{t}\ln |\widetilde{\Psi }|=\widetilde{\Upsilon }-\Upsilon -\partial _{t}\ln |\omega |. \end{aligned}$$

The generalization for additional geometric flow effective source and effective generating source is given, respectively, by

$$\begin{aligned} \ \Upsilon =\ \widetilde{\Upsilon }-(\ln |\omega \widetilde{\Psi }|)^{*}\quad \text{ and } \quad \Theta ^{*}=\ \Upsilon (\widetilde{\Psi }^{2})^{*}. \end{aligned}$$

The generic off-diagonal contributions to non-Riemannian evolution models have nontrivial nonholonomic torsion. LC-configurations can be extracted by respective functional dependencies \(\widetilde{\Psi }^{2}[ \check{\Psi }]\) and effective sources and effective generating sources, respectively, \(\ \check{\Upsilon }=\ \) \(\widetilde{\Upsilon }-(\ln |\omega \widetilde{\Psi }[\check{\Psi }]|)^{*}\) and \(\check{\Theta }^{*}=\ \check{\Upsilon }(\widetilde{\Psi }^{2}[\check{\Psi }])^{*}.\) One of the fundamental consequences of such nonholonomic evolution theories is that various massive, modified, and GR effects can be modeled by nonholonomic constraints even any value of effective sources \(\Upsilon \) can be zero for certain configurations on a 3-d hypersurface \(\Xi _{0}.\)

5.2.4 Relativistic thermodynamic values for N-adapted 3-d modified Ricci flows

One of motivations to find 3-d Ricci flow solutions of Hamilton equations (66) embedded in 4-d spacetimes in GR is that such solutions are characterized by analogous thermodynamical models which can be generalized for associated 4-d modified Lorentz-Ricci solitons, i.e. effective Einstein spaces. We note that the lapse function \(N(u)=-\omega ^{2}h_{4}=-g_{4}\) (41) is contained in explicit form in the integration measure for Eq. (64) if we do not work in N-adapted Gaussian coordinates. For 3-d thermodynamical values, we obtain

$$\begin{aligned}&\, _{\shortmid }\widehat{\mathcal {E}} =\tilde{\tau }^{2}\int _{\widehat{ \Xi }_{t}}\tilde{M}\omega ^{2}h_{4}\sqrt{|q_{1}q_{2}q_{3}|}\mathrm{d}\grave{x} ^{3}\left( \, _{\shortmid }\widehat{R}+|\, _{\shortmid }\widehat{ \mathbf {D}}\tilde{f}|^{2}-\frac{3}{\tilde{\tau }}\right) , \nonumber \\&\, _{\shortmid }\widehat{S} =\int _{\widehat{\Xi }_{t}}\tilde{M}\omega ^{2}h_{4}\sqrt{|q_{1}q_{2}q_{3}|}\mathrm{d}\grave{x}^{3}[ \tilde{\tau }( _{\shortmid }\widehat{R}+|\, _{\shortmid }\widehat{\mathbf {D}} \tilde{f}|^{2}) \!+\!\tilde{f}\!-\!6] , \nonumber \\&\, _{\shortmid }\widehat{\sigma } =-2\ \tilde{\tau }^{4}\int _{\widehat{\Xi } _{t}}\tilde{M}\omega ^{2}h_{4}\sqrt{|q_{1}q_{2}q_{3}|}\mathrm{d}\grave{x}^{3}\nonumber \\&\quad \quad \quad \times \left[ |\, _{\shortmid }\widehat{\mathbf {R}}_{\grave{\imath }\grave{j}}+\, _{\shortmid } \widehat{\mathbf {D}}_{\grave{\imath }}\, _{\shortmid }\widehat{\mathbf {D}}_{ \grave{j}}\tilde{f}-\frac{1}{2\tilde{\tau }}q_{\grave{\imath }\grave{j}}|^{2}\right] , \end{aligned}$$
(68)

up to any parametric function \(\tilde{\tau }(t)\) in \(\tilde{M}=( 4\pi \tilde{\tau }) ^{-3}e^{-\tilde{f}}\) with any \(\tilde{\tau }(t)\) for \( \partial \tilde{\tau }/\partial t=-1\) and \(\chi >0\) Taking the respective 3-d coefficients of a solution (57), or (58) [or any solution of type (60), (61), (62), (63)], and prescribing a closed \(\widehat{\Xi }_{0}\) we can compute such values for any closed \(\widehat{\Xi }_{t}.\) We have to fix an explicit N-adapted system of reference and scaling function \(\tilde{f}\) in order to find certain explicit values for corresponding average energy, entropy, and fluctuations for evolution on a time-like parameter t of any family of closed hypersurfaces. We can decide if certain solutions with an effective Lorentz–Ricci soliton source and/or with contributions from additional MGT sources may be more convenient thermodynamically than other configurations.

Let us compute the values (68) for systems (46)–(50) with solutions (57) with \(q_{1}=q_{2}=e^{\psi },q_{3}= \tilde{\Psi }^{2}/4\Lambda ,h_{4}=(\tilde{\Psi }^{*})^{2}/\ ^{F}\Theta .\) Using the nonlinear symmetry (45) with effective \(\Lambda ,\) for \( \, _{\shortmid }\widehat{\mathbf {R}}_{\grave{\imath }\grave{j}}=\Lambda \mathbf {q}_{\grave{\imath }\grave{j}}\) and \(\, _{\shortmid } \widehat{R}=3\Lambda ;\) taking, for simplicity, \(\ \tilde{f}=0,\) and re-defining \(\frac{3}{2(4\pi )^{3}}\omega ^{2}\rightarrow \omega ^{2},\) we compute

$$\begin{aligned}&\, _{\shortmid }\widehat{\mathcal {E}} =\tilde{\tau }^{2}(t)\int _{\widehat{ \Xi }_{t}}\frac{\tilde{\Psi }(\tilde{\Psi }^{*})^{2}}{\ \Theta \tilde{\tau } ^{3}(t)}e^{\psi }\omega ^{2}\mathrm{d}\grave{x}^{3}\left[ \sqrt{|\Lambda |}-\frac{1}{ \tilde{\tau }(t)\sqrt{|\Lambda |}}\right] , \nonumber \\&\, _{\shortmid }\widehat{S} =\frac{1}{2}\int _{\widehat{\Xi }_{t}}\frac{ \tilde{\Psi }(\tilde{\Psi }^{*})^{2}}{\Theta \tilde{\tau }^{3}(t)}e^{\psi }\omega ^{2}\mathrm{d}\grave{x}^{3}\left[ \tilde{\tau }(t)\sqrt{|\Lambda |}-\frac{2}{ \sqrt{|\Lambda |}}\right] , \nonumber \\&\, _{\shortmid }\widehat{\sigma } =-\ \tilde{\tau }^{4}(t)\int _{\widehat{\Xi }_{t}}\frac{\tilde{\Psi }(\tilde{\Psi }^{*})^{2}}{\Theta \tilde{\tau } ^{3}(t)}e^{\psi }\omega ^{2}\mathrm{d}\grave{x}^{3}\nonumber \\&\quad \quad \times \left[ \sqrt{|\Lambda |}-\frac{1}{2 \tilde{\tau }(t)\sqrt{|\Lambda |}}\right] ^{2}, \end{aligned}$$
(69)

where \(\ \Theta =\int \mathrm{d}t\ \ \Upsilon (\widetilde{\Psi }^{2})^{*}\).

It is possible to define and compute thermodynamic like values (69) for generic off-diagonal solutions in GR as we explained for solutions (58) and corresponding functionals \(\tilde{\Psi }[\check{ \Psi }],\Theta [\check{\Psi }]\) and \(\omega [\check{\Psi }].\) We emphasize that in both cases with zero, or non-zero nonholonomic torsion, the above formulas for thermodynamical values are defined for a non-zero effective \(\widetilde{\Lambda }\) and non-zero source \(\Upsilon .\) For (effective) vacuum configurations, such formulas have to be computed using corresponding classes of off-diagonal solutions.

The entities (68) and/or (69) can be computed for 4-d configurations determined by modified Lorentz–Riemann solitons, \(\widehat{ \mathcal {E}}=\int _{t_{1}}^{t_{2}}\mathrm{d}t\, _{\shortmid }\widehat{\mathcal {E}},\) \( \widehat{\mathcal {S}}=\int _{t_{1}}^{t_{2}}\mathrm{d}t\, _{\shortmid }\widehat{S},\) \( \widehat{{\Sigma }}=\int _{t_{1}}^{t_{2}}\mathrm{d}t\, _{\shortmid }\widehat{\sigma }\) determined by \(_{\shortmid }\widehat{\mathbf {D}}.\) In a similar form, we can use the distortion formulas (27) and compute \(\ ^{\nabla }\mathcal {E} =\int _{t_{1}}^{t_{2}}\mathrm{d}t\, _{\shortmid }^{\nabla }\mathcal {E},\ ^{\nabla } \mathcal {S}=\int _{t_{1}}^{t_{2}}\mathrm{d}t\, _{\shortmid }^{\nabla }S,\ ^{\nabla }{ \Sigma }=\int _{t_{1}}^{t_{2}}\mathrm{d}t\, _{\shortmid }^{\nabla }\sigma \) for \(\, _{\shortmid }\nabla .\) Such values are positively different for \(\, _{\shortmid }\widehat{\mathbf {D}}\rightarrow \, _{\shortmid }\nabla \) for certain topologically nontrivial configurations. As a result, we can analyze if a nonholonomic configuration with N-adapted \(_{\shortmid }\widehat{ \mathbf {D}}\) may be more, or less, convenient thermodynamically than a similar holonomic one as determined by \(\, _{\shortmid }\nabla .\)

Finally, we note that geometric thermodynamics values (68) are defined both for (modified) black hole solutions and inhomogeneous cosmological solutions. The physical meaning of such a thermodynamical approach is very different from that of standard black hole thermodynamics in GR. Nevertheless, certain criteria for equivalent modeling can be analyzed in various MGTs; see [42, 56,57,58,59,60, 62].

5.3 General relativistic hydrodynamics and thermodynamics for geometric flows of gravitational fields

The thermodynamic quantities (64) and (65) are formulated in terms of variables with coefficients computed with respect to N-adapted frame of reference which allows extensions of 3-d thermodynamic values to 4-d ones. This is a hidden aspect of models of non-relativistic hydrodynamics and an apparent property of relativistic theories; see modern approaches in [32,33,34,35,36,37]. In this section, we speculate on a model of geometric flows with local thermodynamic equilibrium for which there exists a single distinguished velocity field \(\mathbf {v}_{\alpha }\) (if it is convenient, we can take \(\mathbf {v} _{\alpha }=\mathbf {n}_{\alpha }),\) and a corresponding N-adapted rest frame, where all flow evolutions of the related physical quantities are described by simple formulas.

For a fluid model with particle production, we characterize general relativistic Ricci flows by the particle number density d-vector \(\widehat{ \mathcal {N}}^{\alpha },\) the entropy density d-vector \(\ \widehat{\mathcal {S }}^{\alpha }\), and the effective energy-momentum density tensor \(\widehat{ \mathcal {T}}^{\alpha \beta }.\) Using the LC-connection \(\nabla ,\) we postulate covariant forms for conservation of the particle number and the effective energy momentum,

$$\begin{aligned} \nabla _{\alpha }\widehat{\mathcal {N}}^{\alpha }(u^{\gamma })=0\quad \text{ and } \quad \nabla _{\alpha }\ \widehat{\mathcal {T}}^{\alpha \beta }(u^{\gamma })=0. \end{aligned}$$
(70)

Using the distortion relations \(\widehat{\mathbf {D}}=\nabla +\widehat{\mathbf {Z}} \) (10), we can compute nonholonomic deformations of such values when \(\widehat{\mathbf {Z}}\) determines the respective sources induced by nonholonomic torsion. This is typical for models of nonholonomic continuous mechanics and generalized hydrodynamic/fluid models. We note that the effective entropy is not conserved. It is supposed that the (general) relativistic effective entropy is zero only in thermodynamic equilibrium:

$$\begin{aligned} \nabla _{\alpha }\widehat{\mathcal {S}}^{\alpha }(u^{\gamma })\ge 0. \end{aligned}$$

With the help of the N-adapted velocity field \(\mathbf {v}^{\alpha }(u^{\gamma }) \) and velocity orthogonal projection operator \(\pi ^{\alpha \beta }=\delta ^{\alpha \beta }-\mathbf {v}^{\alpha }\mathbf {v}^{\beta },\) we introduce necessary values for an effective hydrodynamical model:

$$\begin{aligned} \begin{array}{ll} \widehat{n}:=\widehat{\mathcal {N}}^{\alpha }\mathbf {v}_{\alpha } &{} \text{ is } \text{ the } \text{ particle } \text{ number } \text{ density }; \\ \widehat{j}^{\alpha }:=\pi ^{\alpha \beta }\widehat{\mathcal {N}}_{\beta } &{} \text{ is } \text{ the } \text{ diffusion } \text{ current }; \\ \ \widehat{S}:=\widehat{\mathcal {S}}^{\alpha }\mathbf {v}_{\alpha } &{} \text{ is } \text{ the } \text{ entropy } \text{ density }; \\ \widehat{\mathcal {J}}^{\alpha }:=\pi ^{\alpha \beta }\widehat{\mathcal {S}} _{\beta } &{} \text{ is } \text{ the } \text{ entropy } \text{ current }; \\ \widehat{\mathcal {E}}:=\mathbf {v}_{\alpha }\ \widehat{\mathcal {T}}^{\alpha \beta }\mathbf {v}_{\beta } &{} \text{ is } \text{ the } \text{ energy } \text{ density }; \\ \widehat{p}^{\alpha }:=\pi _{\ \gamma }^{\alpha }\ \widehat{\mathcal {T}} ^{\gamma \beta }\mathbf {v}_{\beta } &{} \text{ is } \text{ the } \text{ momentum } \text{ density }\\ &{}\text{(i.e. } \text{ the } \text{ energy } \text{ current) }; \\ \widehat{\mathcal {P}}^{\alpha \beta }:=\pi _{\ \alpha ^{\prime }}^{\alpha }\pi _{\ \beta ^{\prime }}^{\beta }\ \widehat{\mathcal {T}}^{\alpha ^{\prime }\beta ^{\prime }} &{} \text{ is } \text{ the } \text{ pressure } \text{ tensor }. \end{array} \end{aligned}$$

Hats are used in order to emphasize that necessary coefficients are defined in N-adapted form. As a matter of principle, all constructions can be redefined in arbitrary frames of reference (when symbols “lose” their hats). In addition to these local N-adapted rest frame quantities, we consider such values:

$$\begin{aligned} \begin{array}{ll} \text{ the } \text{ effective } \text{ energy-momentum } \text{ vector } &{} \widehat{\mathcal {E}} ^{\alpha }:=\mathbf {v}_{\beta }\ \widehat{\mathcal {T}}^{\alpha \beta }; \\ \text{ the } \text{ energy-momentum } \text{ current } \text{ density } &{} \widehat{\mathcal {B}}^{\alpha \beta }:=\pi _{\ \gamma }^{\alpha }\ \widehat{\mathcal {T}}^{\gamma \beta }; \\ \text{ momentum } \text{ density } &{} \widehat{p}^{\alpha }:=\pi _{\ \gamma }^{\alpha }\ \widehat{\mathcal {E}}^{\gamma }; \\ \text{ the } \text{ energy } \text{ current } &{} \widehat{p}^{\alpha }:=\mathbf {v}_{\beta } \widehat{\mathcal {B}}^{\alpha \beta }. \end{array} \end{aligned}$$

Using the above formulas, we compute the local rest N-adapted frame densities determined by general relativistic Ricci flows as

$$\begin{aligned}&\widehat{\mathcal {N}}^{\alpha } =\widehat{n}\mathbf {v}^{\alpha }+\widehat{j }^{\alpha }, \nonumber \\&\widehat{\mathcal {S}}^{\alpha } =\widehat{S}\mathbf {v}^{\alpha }+\widehat{ \mathcal {J}}^{\alpha }, \nonumber \\&\ \widehat{\mathcal {T}}^{\alpha \beta } =\widehat{\mathcal {E}}^{\alpha } \mathbf {v}^{\beta }+\widehat{\mathcal {B}}^{\alpha \beta }=\widehat{\mathcal {E }}\mathbf {v}^{\alpha }\mathbf {v}^{\beta }+\widehat{p}^{\alpha }\mathbf {v} ^{\beta }+\mathbf {v}^{\alpha }\widehat{p}^{\beta }+\widehat{\mathcal {P}} ^{\alpha \beta }.\nonumber \\ \end{aligned}$$
(71)

These formulas conveniently express the particle number density for the 4-d d-vector and the energy-momentum density d-tensor with the help of local rest frame quantities relative to a velocity field. Such formulas can be related to a 3-d hypersurface Perelman thermodynamic model (64) if \(\widehat{S}=\, _{\shortmid }\widehat{S}\) and \(\widehat{\mathcal {E}}=\, _{\shortmid }\widehat{\mathcal {E}}.\) For \(\mathbf {v}^{\beta }=\mathbf {n} ^{\beta }\) determined by an exact solution in GR, the densities (71) can be normalized to respective values in (65).

Using the canonical d-connection \(\widehat{\mathbf {D}}\) and the above formulas, the balance equations (70) can be written in the form

$$\begin{aligned} (\widehat{\mathbf {D}}_{\alpha }-\widehat{\mathbf {Z}}_{\alpha })\widehat{ \mathcal {N}}^{\alpha }&=\mathbf {v}^{\alpha }(\widehat{\mathbf {D}}_{\alpha }- \widehat{\mathbf {Z}}_{\alpha })\widehat{n}\\&\quad +\widehat{n}(\widehat{\mathbf {D}} _{\alpha }-\widehat{\mathbf {Z}}_{\alpha })\mathbf {v}^{\alpha }+(\widehat{ \mathbf {D}}_{\alpha }-\widehat{\mathbf {Z}}_{\alpha })\widehat{j}^{\alpha } \end{aligned}$$

and

$$\begin{aligned}&(\widehat{\mathbf {D}}_{\beta }-\widehat{\mathbf {Z}}_{\beta })\widehat{ \mathcal {T}}^{\alpha \beta } =\mathbf {v}^{\gamma }(\widehat{\mathbf {D}} _{\gamma }-\widehat{\mathbf {Z}}_{\gamma })\widehat{\mathcal {E}}^{\alpha }\\&\qquad + \widehat{\mathcal {E}}^{\alpha }(\widehat{\mathbf {D}}_{\gamma }-\widehat{ \mathbf {Z}}_{\gamma })\mathbf {v}^{\gamma }+(\widehat{\mathbf {D}}_{\beta }- \widehat{\mathbf {Z}}_{\beta })\widehat{\mathcal {B}}^{\alpha \beta } \\&\quad =\mathbf {v}^{\alpha }\mathbf {v}^{\gamma }(\widehat{\mathbf {D}}_{\gamma }- \widehat{\mathbf {Z}}_{\gamma })\widehat{\mathcal {E}}+\widehat{\mathcal {E}} \mathbf {v}^{\alpha }(\widehat{\mathbf {D}}_{\gamma }-\widehat{\mathbf {Z}} _{\gamma })\mathbf {v}^{\gamma }\\&\qquad +\,\mathbf {v}^{\gamma }(\widehat{\mathbf {D}} _{\gamma }-\widehat{\mathbf {Z}}_{\gamma })\widehat{p}^{\alpha }+\widehat{p} ^{\alpha }(\widehat{\mathbf {D}}_{\gamma }-\widehat{\mathbf {Z}}_{\gamma }) \mathbf {v}^{\gamma } \\&\qquad +\,\widehat{\mathcal {E}}\mathbf {v}^{\gamma }(\widehat{\mathbf {D}}_{\gamma }- \widehat{\mathbf {Z}}_{\gamma })\mathbf {v}^{\alpha }+\mathbf {v}^{\alpha }( \widehat{\mathbf {D}}_{\gamma }-\widehat{\mathbf {Z}}_{\gamma })\widehat{p} ^{\gamma }\\&\qquad +\,\widehat{p}^{\gamma }(\widehat{\mathbf {D}}_{\gamma }-\widehat{ \mathbf {Z}}_{\gamma })\mathbf {v}^{\alpha }+(\widehat{\mathbf {D}}_{\beta }- \widehat{\mathbf {Z}}_{\beta })\widehat{\mathcal {P}}^{\alpha \beta }. \end{aligned}$$

These formulas are written in N-adapted form. In general, they contain non-dissipative and dissipative components of general relativistic Ricci flows. In terms of the LC-connection \(\nabla =\widehat{\mathbf {D}}-\widehat{ \mathbf {Z}},\) we can eliminate certain sources of nonholonomic torsion dissipation but in LC-variables the formulas cannot be decoupled in general form.

Let us discuss the fact and the conditions when the particle number density d-vector and the energy-momentum d-tensor split into a non-dissipative and a dissipative part. Such parts are easily distinguished in a local rest N-adapted frame when the non-dissipative particle number density d-vector \(\widehat{\mathcal {N}}_{0}^{\alpha }\) is parallel to the d-velocity and the non-dissipative energy-momentum d-tensor \( \widehat{\mathcal {T}}^{\gamma \beta }\) is diagonal. We have a particular case of balance equations (71) when

$$\begin{aligned} \widehat{\mathcal {N}}_{0}^{\alpha }=\widehat{n}_{0}\mathbf {v}^{\alpha }\quad \text{ and } \quad \widehat{\mathcal {T}}_{0}^{\alpha \beta }=\widehat{\mathcal {E}} _{0}\mathbf {v}^{\alpha }\mathbf {v}^{\beta }-p_{0}\mathcal {\pi }^{\alpha \beta }, \nonumber \end{aligned}$$

with an effective static (scalar) pressure determined by the state equation of state of the effective fluid. Such approximations are possible for gravitational configurations without singularities. For instance, black holes may have a non-zero entropy (even determined in a different form from an approach with W-entropy) together with a non-zero particle production. In such a case, we must introduce a term like \(\widehat{\mathcal {S}} _{0}^{\alpha }=\widehat{S}_{0}\mathbf {v}^{\alpha }+\widehat{\mathcal {J}} _{0}^{\alpha }\) with an entropy current \(\widehat{\mathcal {J}}_{0}^{\alpha }\) determined by geometric flows and \(\widehat{S}_{0}\) identified with the standard black hole entropy. Therefore the diffusion current density \( \widehat{j}^{\alpha }\) consists of the dissipative part of the particle number (and the momentum density/energy) currents. This defines the difference of the total and the equilibrium pressure \(\Pi ^{\alpha \beta }=\widehat{ \mathcal {P}}^{\alpha \beta }+p_{0}\mathcal {\pi }^{\alpha \beta }.\) The effective viscous pressure determines the dissipative parts of the energy momentum.

It should be mentioned that a splitting into some dissipative and non-dissipative parts of certain physical quantities is related to and depends on the local rest frame (how it is chosen and N-adapted). Additionally to the velocity field of the geometric flow continuum, this requires a particular thermostatics to determine the static pressure. However, the dissipation of thermodynamical systems (in our approach, by general relativistic Ricci flows) is principally defined by the entropy production. This is implicitly related to the background thermostatics, which is a concept for systems in local equilibrium. On the other hand, the rest frame is not determined in the case of dissipation, neither by any special form of the physical quantities. For geometric flows related to 4-d exact solutions in GR, one may consider to extend the Perelman thermodynamic approach by additional construction with the velocity field. We have to address nonholonomic kinetic theory and possible Finsler-like modifications of GR, relativistic thermodynamics, and diffusion as in Refs. [7,8,9]. Such models are planned to be elaborated in our further work.

The general relativistic Ricci flow evolution is related to problems of stability and causality, similar to those in relativistic hydrodynamics. This involves a detailed analysis of the Second Law as in Refs. [35,36,37]. We can assume the possibility of acceleration independent entropy production when the local rest frame entropy density is a function of type \(\widehat{S}(\widehat{\mathcal {E}},\widehat{n})=\widehat{S }(\sqrt{\widehat{\mathcal {E}}^{\alpha }\widehat{\mathcal {E}}_{\alpha }}, \widehat{n})\) for the quantities defined above. Then it is possible to write a Gibbs-like relation,

$$\begin{aligned} \mathrm{d}\widehat{\mathcal {E}}+\widehat{\mathcal {E}}^{-1}\widehat{p}_{\gamma }\mathrm{d} \widehat{p}^{\gamma }=\beta ^{-1}\mathrm{d}\widehat{S}+\widehat{\mu }\mathrm{d}\widehat{n}, \end{aligned}$$

for an effective temperature \(\beta \) and effective particle production potential \(\widehat{\mu }.\) Such an approach to relativistic thermodynamics of Ricci flows eliminates the generic instability which exists in the original Eckart theory and the stability conditions are independent of any flow frames; see [28, 29]. Such conditions are the non-negativity of the geometric flow transport coefficients and the effective thermodynamic stability and the concavity of the entropy density. This is in contrast to the former variants of Eckart theory (which is unstable), or to Israel–Stewart theory. In the latter case, there are several complicated conditions; for details see [33, 35, 36].

Finally, the main conclusion is that we can characterize the gravitational field equations by Perelman’s W-entropy on any closed 3-d hypersurface. In such a case, a natural statistical thermodynamic model can be associated. To extend the approach in a 4-d general relativistic form it is necessary to elaborate on modified relativistic thermodynamic and hydrodynamic theories, which is the main purpose of this work. The constructions can be performed in explicit form and related to general classes of exact solutions in GR if we formulate all gravitational field, geometric evolution, and effective relativistic thermodynamic theories in certain N-adapted nonholonomic variables. In the next section, we show how this can be connected to exact solutions.

5.4 Parameterizations for the CET model

A thermodynamic model is relativistic if it is derived for the energy momentum conservation equations

$$\begin{aligned} \widehat{\mathbf {D}}_{\beta }(\mathbf {v}_{\alpha }\mathbf {T}^{\alpha \beta })=\widehat{\mathbf {D}}_{\beta }(\mathbf {v}_{\alpha })\mathbf {T}^{\alpha \beta }-\mathbf {v}_{\alpha }\widehat{\mathbf {J}}^{\alpha }, \end{aligned}$$
(72)

considering the heat flow \(\mathbf {v}_{\alpha }\widehat{\mathbf {J}}^{\alpha } \) into an effective fluid with \(\widehat{\mathbf {J}}^{\alpha }=-\widehat{ \mathbf {D}}_{\beta }\mathbf {T}^{\alpha \beta }.\) Footnote 9 For perfect (pressureless) matter, \(\mathbf {T}_{\alpha \beta }=p\mathbf {g}_{\alpha \beta }+(\rho +p) \mathbf {v}_{\alpha }\mathbf {v}_{\beta },\) where \(\mathbf {v}_{\alpha }\mathbf { v}^{\alpha }=-1\) , \(\rho \), and p are, respectively, the density and pressure. The 4-velocity of the fluid can be taken as \(v^{\alpha }=(0,0,0,1) \) in a certain co-moving N-adapted frames when \(\mathbf {T} _{~\beta }^{\alpha }=\mathrm{diag}[0,0,0,-\rho ].\) We can consider also fluids with nontrivial momentum density \(\mathbf {q}^{\alpha }\) and anisotropic (tracefree) pressure \(\pi _{\alpha \beta },\) when

$$\begin{aligned} \mathbf {T}_{\alpha \beta }=p\mathbf {h}_{\alpha \beta }+\rho \mathbf {v} _{\alpha }\mathbf {v}_{\beta }+\mathbf {q}_{\alpha }\mathbf {v}_{\beta }+ \mathbf {v}_{\alpha }\mathbf {q}_{\beta }+\pi _{\alpha \beta }. \end{aligned}$$
(73)

Defining \(\widehat{\mathbf {\Theta }}:=(\mathbf {v}^{\gamma }\widehat{\mathbf {D }}_{\gamma }v)/v\) for a spatial domain \(\upsilon \) with volume \( V=\int \nolimits _{\upsilon }v,\) the entropy is given by \(S=\int \nolimits _{ \upsilon }\vartheta ,\) when (72) and (73) lead to

$$\begin{aligned} \widehat{\mathcal {T}}\ \mathbf {v}^{\gamma }\ \widehat{\mathbf {D}}_{\gamma }\vartheta:= & {} \mathbf {v}^{\beta }\widehat{\mathbf {D}}_{\beta }(\rho v)+p \mathbf {v}^{\gamma }\widehat{\mathbf {D}}_{\gamma }v\\= & {} v(\mathbf {v}_{\alpha } \widehat{\mathbf {J}}^{\alpha }-\widehat{\mathbf {D}}_{\beta }\mathbf {q} ^{\beta }-\mathbf {q}^{\alpha }\mathbf {v}^{\gamma }\widehat{\mathbf {D}} _{\gamma }\mathbf {v}_{\alpha }-\sigma ^{\alpha \beta }\pi _{\alpha \beta }). \end{aligned}$$

The variables \(\vartheta \) and \(\widehat{\mathcal {T}},\) represent the point-wise entropy and temperature. If we define \(\widehat{\mathcal {T}}\) independently, we get \(S=\int \nolimits _{\upsilon }\vartheta .\) Footnote 10

We are seeking to construct a Ricci flow and the gravitational analog of the fundamental laws of the thermodynamics in the form

$$\begin{aligned} ~^{g}\mathcal {T}\mathrm{d}~^{g}S=\mathrm{d}~^{g}U+~^{g}p\mathrm{d}V, \end{aligned}$$
(74)

where \(~^{g}\mathcal {T},~^{g}S,~^{g}U\) and \(~^{g}p\) are, respectively, the effective temperature, entropy, internal energy, isotropic pressure; V is the spatial volume. For the vacuum gravitational fields, we consider the quantities

$$\begin{aligned} \widehat{\mathbf {W}}= & {} \frac{1}{4}(\widehat{\mathbf {E}}_{\alpha }^{\quad \beta }\widehat{\mathbf {E}}_{\quad \beta }^{\alpha }+\widehat{\mathbf {H}}_{\alpha }^{\quad \beta }\widehat{\mathbf {H}}_{\quad \beta }^{\alpha }),\nonumber \\ \widehat{\mathbf {J}}_{\alpha }= & {} \frac{1}{2}[\widehat{\mathbf {E}},\widehat{ \mathbf {H}}]_{\alpha },\quad \mathbf {v}^{\gamma }\widehat{\mathbf {D}} _{\gamma }\widehat{\mathbf {W}}+\widehat{\mathbf {D}}^{\alpha }\widehat{ \mathbf {J}}_{\alpha }\simeq 0, \end{aligned}$$

where \(\widehat{\mathbf {W}}\) is the ‘super-energy density’ and the ’super-Poynting vector’ \(\widehat{\mathbf {J}}_{\alpha }=-\mathbf {h}_{\alpha }^{~\delta }\widehat{\mathbf {T}}_{\delta \beta \gamma \tau }\mathbf {v} ^{\beta }\mathbf {v}^{\gamma }\mathbf {v}^{\tau }\) is determined by the Bel–Robinson tensor (for holonomic configurations; for details see [40, 64,65,66]),

$$\begin{aligned} \widehat{\mathbf {T}}_{\alpha \beta \gamma \tau }:=\frac{1}{4}(\widehat{ \mathbf {C}}_{\varepsilon \alpha \beta \varphi }\widehat{\mathbf {C}}_{\ \gamma \tau }^{\varepsilon \quad \varphi }+\ ^{*}\widehat{\mathbf {C}} _{\varepsilon \alpha \beta \varphi }\ ^{*}\widehat{\mathbf {C}}_{\ \gamma \tau }^{\varepsilon \quad \varphi }), \end{aligned}$$

when the dual Weyl tensor is \(\ ^{*}\widehat{\mathbf {C}}_{\alpha \beta \varphi \tau }:=\frac{1}{2}\eta _{\alpha \beta \varepsilon \delta }\widehat{ \mathbf {C}}_{\ \varphi \tau }^{\varepsilon \delta }.\) There is a two-index ‘square-root’ [67], \(\mathfrak {t}_{\alpha \beta },\) constructed as a solution of

$$\begin{aligned} \widehat{\mathbf {T}}_{\alpha \beta \gamma \tau }&=\mathfrak {t}_{(\alpha \beta }\mathfrak {t}_{\gamma \tau )}-\frac{1}{2}\mathfrak {t}_{\varepsilon (\alpha } \mathfrak {t}_{\beta }^{\ \varepsilon }\mathbf {g}_{\gamma \tau )}\\&\quad \,\, +\frac{1}{24} \left[ \mathfrak {t}_{\varepsilon \mu }\mathfrak {t}^{\varepsilon \mu }+\frac{1}{2}( \mathfrak {t}_{\beta }^{\ \varepsilon })^{2}\right] \mathbf {g}_{(\alpha \beta } \mathbf {g}_{\gamma \tau )}. \end{aligned}$$

The solutions of these equations depend on the algebraic type (following Petrov’s classification) of spacetime and related classes of solutions of the Einstein equations; for details see [40, 68]. There are variants such as

$$\begin{aligned} \mathfrak {t}_{\alpha \beta }=\left\{ \begin{array}{ll} 3\epsilon |\widehat{\Psi }_{2}|(\mathbf {m}_{(\alpha }{\bar{\mathbf {m}}}_{\beta )}+\mathbf {l}_{\alpha }\mathbf {k}_{\beta }), &{} \text{ Petrov } \text{ type } \text{ N, } \text{ similar } \text{ to } \text{ pure } \text{ radiation; } \\ \epsilon |\widehat{\Psi }_{4}|\mathbf {k}_{\alpha }\mathbf {k}_{\beta }, &{} \text{ Petrov } \text{ type } \text{ D, } \text{ Coulomb-like } \text{ gravitational } \text{ configurations; } \\ \text{ various } \text{ forms }, &{} \text{ other } \text{ Petrov } \text{ types, } \text{ factorized } \text{ as } \text{ D } \text{ or } \text{ N, } \text{ or } \text{ more } \text{ complicated }, \end{array} \right. \end{aligned}$$

where \(\epsilon =\pm 1,\) \(\widehat{\Psi }_{2}=\widehat{\mathbf {C}} _{\varepsilon \alpha \beta \varphi }\mathbf {k}^{\varepsilon }\mathbf {m} ^{\alpha }{\bar{\mathbf {m}}}^{\beta }\mathbf {l}^{\varphi }\) and \(\widehat{\Psi }_{4}=\widehat{\mathbf {C}}_{\varepsilon \alpha \beta \varphi }{\bar{\mathbf {m}} }^{\varepsilon }\mathbf {l}^{\alpha }{\bar{\mathbf {m}}}^{\beta }\mathbf {l} ^{\varphi }.\) The first two cases above contain very interesting examples corresponding to stationary black hole solutions and the cases of scalar perturbations of FLRW geometries and their off-diagonal deformations [41].

Using the respective \(\mathfrak {t}_{\alpha \beta }\) and \(|\widehat{\Psi }_{2}|= \sqrt{2\widehat{\mathbf {W}}/3},\) we can construct a thermodynamic model (74) for Coulomb-like gravitational fields with effective energy-momentum tensor of type (73),

$$\begin{aligned} \mathfrak {t}_{\alpha \beta }\simeq \varkappa \ ^{g}\mathbf {T}_{\alpha \beta }=\hat{\alpha }\sqrt{2\widehat{\mathbf {W}}/3}[\mathbf {x}_{\alpha }\mathbf {x} _{\beta }+\mathbf {y}_{\alpha }\mathbf {y}_{\beta }+2(\mathbf {v}_{\alpha } \mathbf {v}_{\beta }-\mathbf {z}_{\alpha }\mathbf {z}_{\beta })], \end{aligned}$$

for a constant \(\hat{\alpha }\) to be determined from certain experimental data and/or other theoretic considerations. The corresponding effective energy density and pressure are \(\varkappa \ ^{g}\rho =2\hat{\alpha }\sqrt{2 \widehat{\mathbf {W}}/3}\ge 0\) and \(\ ^{g}p=0.\ \)We can take \(\ ^{g}q_{\alpha }=0\) and write the effective fundamental thermodynamic equation for \(\ ^{g}S=\int \nolimits _{\upsilon }\ ^{g}\vartheta \) in the presence of perfect matter field in the form

$$\begin{aligned}&~^{g}\mathcal {T\ }\mathbf {v}^{\gamma }\widehat{\mathbf {D}}_{\gamma }(~^{g}\vartheta )=\mathbf {v}^{\gamma }\widehat{\mathbf {D}}_{\gamma }(\ ^{g}\rho v)\\&\quad =-v\sigma _{\alpha \beta }[\pi ^{\alpha \beta }+\varkappa (\rho +p)\widehat{\mathbf {E}}^{\alpha \beta }/2\hat{\alpha }\sqrt{6\widehat{\mathbf { W}}}]. \end{aligned}$$

For wave-like gravitational fields \(|\Psi _{4}|=2\sqrt{\widehat{\mathbf {W}}} .\) Working with plane wave geometries in Kundt’s class of solutions (for holonomic configurations; for details see [40]), the effective thermodynamic quantities

$$\begin{aligned}&\varkappa \ ^{g}\rho =2\hat{\beta }\sqrt{\widehat{\mathbf {W}}},\ ^{g}p=\ ^{g}\rho /3,\quad \varkappa \ ^{g}q_{\alpha }=2\hat{\beta }\sqrt{\widehat{\mathbf { W}}}\mathbf {z}_{\alpha }, \\&\quad 6\varkappa \ ^{g}\pi _{\alpha \beta } =-2\hat{\beta }\sqrt{\widehat{\mathbf {W}}}(\mathbf {x}_{\alpha }\mathbf {x}_{\beta }+\mathbf { y}_{\alpha }\mathbf {y}_{\beta }-2\mathbf {z}_{\alpha }\mathbf {z}_{\beta }) \end{aligned}$$

can be taken as

$$\begin{aligned} \mathfrak {t}_{\alpha \beta }\simeq \varkappa \ ^{g}\mathbf {T}_{\alpha \beta }=2\hat{\beta }\sqrt{\widehat{\mathbf {W}}}\mathbf {k}_{\alpha }\mathbf {k} _{\beta } \end{aligned}$$

with a constant \(\hat{\beta }\) (in general, \(\hat{\beta }\ne \hat{\alpha }).\) In the presence of a perfect fluid, the resulting fundamental thermodynamic equation for \(\ ^{g}S=\int \nolimits _{\upsilon }\ ^{g}\vartheta \) is

$$\begin{aligned}&~^{g}\mathcal {T\ }\mathbf {v}^{\gamma }\widehat{\mathbf {D}}_{\gamma }(~^{g}\vartheta ) =\mathbf {v}^{\gamma }\widehat{\mathbf {D}}_{\gamma }(\ ^{g}\rho v)+\ ^{g}p\mathbf {v}^{\gamma }\widehat{\mathbf {D}}_{\gamma }v \\&\quad =-v[\sigma _{\alpha \beta }\ ^{g}\pi ^{\alpha \beta }+\mathbf {g}^{\alpha \beta }\widehat{\mathbf {D}}_{\beta }(\ ^{g}\mathbf {q}_{\alpha })+\ ^{g} \mathbf {q}_{\alpha }\mathbf {v}^{\gamma }\widehat{\mathbf {D}}_{\gamma } \mathbf {v}^{\alpha }]\\&\qquad -\,\varkappa (\rho +p)v\sigma _{\alpha \beta }\widehat{ \mathbf {E}}^{\alpha \beta }/4\hat{\beta }\sqrt{\widehat{\mathbf {W}}}. \end{aligned}$$

There are alternative definitions of the gravitational temperature which depend on the type of solutions considered, for instance, in a black hole or cosmological model. Usually, one postulates the expression

$$\begin{aligned} ~^{g}\mathcal {T}=8|\mathbf {k}^{\gamma }\mathbf {l}^{\beta }\widehat{\mathbf {D} }_{\gamma }\mathbf {v}_{\beta }|/\varkappa \ =4|\mathbf {z}^{\beta }\mathbf {v} ^{\gamma }\widehat{\mathbf {D}}_{\gamma }\mathbf {v}_{\beta }+\widehat{H} +\sigma _{\alpha \beta }\mathbf {z}^{\alpha }\mathbf {z}^{\beta }|/\varkappa , \end{aligned}$$
(75)

where \(\widehat{H}=\widehat{\Theta }/3\) is the isotropic Hubble rate. Such a formula reproduces in appropriate limits the formulas from quantum field theory in curved spacetimes, black hole thermodynamics, de Sitter spaces, or Unruh temperature etc. It can be defined in MGTs with effective modeling by off-diagonal deformations of Einstein spaces.

The smooth function f on a 3-dimensional closed hypersurface can be considered as a function determining the natural log of partition function \( Z=\int e^{-\beta E}\mathrm{d}\omega (E),\) with density of states measure \(\omega (E),\)

$$\begin{aligned} \log Z=\int _{\Sigma _{t}}[-f+3/2](4\pi \tau )^{-3/2}e^{-f}\mathrm{d}\check{v}. \end{aligned}$$
(76)

To find thermodynamic values in explicit form we consider the log of the partition function \(\log Z\) (76) on a closed 3-d hypersurface \( \Sigma _{t}\) of volume \(V=\int \nolimits _{\upsilon }v,\) when the entropy is given by \(S=\) \(\int \nolimits _{\upsilon }\vartheta \) and the constant \(\tau =\beta ^{-1}\) is treated as an effective temperature. In order to reproduce relativistic Ricci flows effective thermodynamics models for gravitational fields with temperature (75), the corresponding Perelman type entropy

$$\begin{aligned} S=\left( 1-\beta \frac{\partial }{\partial \beta }\right) \log Z=-\int \nolimits _{\upsilon }[\tau ({\breve{\mathbf {R}}}+|{\breve{\mathbf {D}}}f|^{2})+f-3]\mathrm{d}\check{v} \end{aligned}$$
(77)

can be considered as a functional \(~^{g}S(S)\) when

$$\begin{aligned} ~^{g}\mathcal {T}\mathrm{d}~^{g}S=\tau \mathrm{d}S. \end{aligned}$$
(78)

In this way, we transform the thermodynamical variables for \( ~^{g}U(~^{g}S,V)=U(S,V)\) and the fundamental law of thermodynamics (74) is re-written in the form

$$\begin{aligned} \mathrm{d}U=\tau \mathrm{d}S-~^{g}p\mathrm{d}V, \end{aligned}$$

where \(p=~^{g}p\) and the relation \(\partial S/\partial \tau =\sigma ^{2}/\tau ^{3}\) follows from the definition of the W-entropy.

In explicit form, the effective entropy of the gravitational fields and thermodynamic transforms should be defined differently, for instance, for the Coulomb and/or wave-like gravitational fields, when the respective constants \(\hat{\alpha }\) and \(\hat{\beta }\) are chosen so as to obtain equivalent models for \(~^{g}S\) and/or S. Prescribing any values for independent data \( ~^{g}S\) and \(~^{g}\mathcal {T}\), we can construct S and \(\log Z\) using, respectively, Eqs. (77) and (78). We denote the solution of the last equation in the form \(\log ~^{g}Z.\) Inverting Eq. (76), we can find a corresponding function \(f=~^{g}f\) describing the geometric evolution of a configuration of gravitational fields characterized by certain relativistic thermodynamics data \( (~^{g}S,~^{g}\mathcal {T}).\)

In [1], the W-entropy was defined by analogy to statistical thermodynamics but a microscopic description was not considered. In statistical mechanics, the partition function \(Z=\sum \nolimits _{i}e^{-\beta E_{i}}\) for a thermodynamic system is \(\Gamma [\varphi _{i},E_{i},\beta ]\) in a canonical ensemble with the \(E_{i}\) being the energies associated with corresponding “microstates” \(\varphi _{i}\) at inverse temperature \(\beta \) (we consider summation on all available “microstates”). For non-discrete microstates, one considers integrals over the spaces of microstates \(\Omega ,\) \(Z=\int \nolimits _{\upsilon }e^{-\beta E(\omega )}\mathrm{d}\omega ,\) with associated energy functions \(E:\Omega \rightarrow \mathbb {R }\) and \(\mathrm{d}\omega \) being the density of states measure on \(\Omega .\) The entropy S (77) of \(\Gamma \) was computed for equilibrium states at temperature \(\tau .\) In the relativistic case, we can consider such constructions for an effective theory for the evolution of certain 3-d metrics embedded into 4-d spacetime for which a local gravitational equilibrium exists for a time-like variable t. As a thermodynamic system, the analogous \(\Gamma \) is characterized by a full set of “macrostates” describing its large-scale properties and equivalently modeled by the data \(~^{g}S\) and \(~^{g} \mathcal {T}\).

It should be noted that Perelman called f the dilaton field and used also an F-functional, \(F=-\tau \log Z,\) with a “first variation” formula which “can be found in the literature on string theory, where it describes the low energy effective action; \(\ldots \)” see the end of Sect. 1 in [1]. Recently, the analogy with string theory in the Polyakov formulation [69] was exploited in Ref. [22], where a general scheme of defining partition functions associated to relevant geometric flows was proposed. The Lin approach is based on the functional determinant of differential operators and has well-defined microstates as members of a functional space. For each microstate, it associates the “Dirichlet energy” which in turn is associated to the underlying operator. Using double, \(3+1\), and \(2+2\) fibrations, the corresponding assumptions and data (on analogous macroscopic thermodynamical values associated to solutions of the Einstein equations, additional smooth functions etc.) we can determine the operator for energy determination or macrostates.

6 Final remarks and conclusions

Perelman’s proof of the Poincaré conjecture provided a fundamental result on the topological structure of the Universe. He elaborated a general schematic procedure in geometric analysis and speculated on a number of perspectives for applications in modern physics. A statistical and thermodynamic analogy was developed for geometric evolution scenarios using Lyapunov-type functionals. The approach was based on the concept of W-entropy. It was supposed from the very beginning [1] that geometric flows may have certain implications for black hole physics and string theory but the original theory of Ricci evolution flows was formulated in a non-relativistic form. So, Perelman’s ideas could not be developed in a framework of theories with geometric flow evolution of 3-d Riemannian metrics.

This article has a very exact goal: to understand the W-entropy in a wide range of general relativistic geometric flows and analyze possible connections to other formulations of gravitational thermodynamics describing general nonhomogeneous cosmological solutions and black hole configurations. To this end, we need to provide a satisfactory relativistic account of analogous statistical thermodynamical models associated to geometric flows of exact solutions in GR resulting in effective polarizaton (running) of the fundamental constant and generic off-diagonal effects. This will remain a challenge with many interesting directions in MGTs and modern accelerating cosmology emerging from the work presented here. An interesting problem is that, for instance, certain scenarios in modern cosmology can be modeled equivalently by generic off-diagonal interactions in GR and/or by analogous solutions in a MGT. Nevertheless, such gravitational field configurations have very different properties in a framework of geometric flow theory, considering its self-similar configurations as (relativistic) Ricci solitons. The corresponding analogous thermodynamics is determined, in general, by different types of physical values. This motivates our approach to study general relativistic models of Ricci flow evolution of exact solutions in gravity theories.

In this paper, we have shown in explicit form how the W-entropy can be generalized for relativistic geometric evolution of solutions of the Einstein equations. This allows for a statistical thermodynamic description of gravitational interactions, their (fractional) diffusion and kinetic processes [7,8,9]. The marriage of the analogous Perelman’s thermodynamics and GR is an issue of very active debates on physical meaning of such constructions and new mathematics. A temperature-like evolution parameter can be considered also for relativistic generalizations but we cannot address the problem as a specific type heat propagation in the context of relativistic dynamics of a non-perfect fluid like in various approaches to relativistic thermodynamics, nonlinear diffusion and relativistic kinetic theory [32,33,34,35,36,37]. In theories with geometric flows, we consider the evolution of metrics and fundamental geometric objects (in general, these can be certain generalized connections, curvatures, and torsions). Similar to heat fluxes, we can consider non-relativistic gradient flows and consider an unbounded speed of effective thermal disturbances (fluctuations of tensor fields). This is due to the parabolic nature of PDEs describing the heat transport. In the relativistic cases, we get not only an inconvenience but indeed a fundamental problem. Various issues were studied related to this problem in relativistic hydrodynamical and dissipative theories [25,26,27,28,29, 31] (for a review, see [38]).

Such constructions, in a macroscopic relativistic thermodynamic formulation [40], seem to be important for the structure formation of Universes, black hole thermodynamics, and the characterization of nonlinear waves, locally anisotropic and/inhomogeneous gravitational configurations. It is possible to formulate an underlying operator formalism and microscopic approach in the spirit of Polyakov’s approach to string theory. This scheme requires a method of constructing generic off-diagonal solutions depending on all spacetime variables [41, 42] which drive in nonlinear parametric form the relativistic Ricci flows evolution of the 3-d hypersurface metrics in a 4-d spacetime.

We conclude that the anholonomic frame deformation method involved in our constructions can be applied to a wide range in GR and MGTs [59, 60, 70,71,72,73,74] and generalized geometric evolution models [19,20,21]. Such constructions are more general and different from those elaborated in [75, 76] where Ricci flows were studied in connection with the standard black hole thermodynamics and possible extensions to quantum diffusion and informational entropy.