Black hole solutions in the quadratic Weyl conformal geometric theory of gravity

Abstract

We consider numerical black hole solutions in the Weyl conformal geometry and its associated conformally invariant Weyl quadratic gravity. In this model, Einstein gravity (with a positive cosmological constant) is recovered in the spontaneously broken phase of Weyl gravity after the Weyl gauge field ($\omega _{\mu }$) becomes massive through a Stueckelberg mechanism and it decouples. As a first step in our investigations, we write down the conformally invariant gravitational action, containing a scalar degree of freedom and the Weyl vector. The field equations are derived from the variational principle in the absence of matter. By adopting a static spherically symmetric geometry, the vacuum field equations for the gravitational, scalar, and Weyl fields are obtained. After reformulating the field equations in a dimensionless form, and by introducing a suitable independent radial coordinate, we obtain their solutions numerically. We detect the formation of a black hole from the presence of a Killing horizon for the timelike Killing vector in the metric tensor components, indicating the existence of the singularity in the metric. Several models corresponding to different functional forms of the Weyl vector are considered. An exact black hole model corresponding to a Weyl vector having only a radial spacelike component is also obtained. The thermodynamic properties of the Weyl geometric type black holes (horizon temperature, specific heat, entropy, and evaporation time due to Hawking luminosity) are also analyzed in detail.

1 Introduction

One of the fundamental results in theoretical relativistic astrophysics indicates that compact objects having mass functions greater than 3–4 solar masses must be black holes [1]. Black hole solutions were among the first to have been considered in the framework of general relativity, with the first vacuum solution of the Einstein gravitational field equations obtained more than 100 years ago [2]. However, the first firm observational evidence for the existence of black holes was found relatively recently, beginning with the astronomical discoveries of the 1970s. One of the first black hole candidates was observed in a binary system consisting of the supergiant star HD 226868 associated with the high-mass X-ray binary object Cygnus (Cyg) X-1. The determination of the mass function $f\left( M_x\right)$, giving the mass $M_x$ of the compact object Cyg X-1 in terms of the companion star $M_c$, and of the inclination angle i already showed that $M_x>4M_{\odot }$ [3], indicating that Cyg X-1 may be a stellar mass black hole. In the Milky Way alone, the total number of stellar mass black holes (in binaries or isolated) is estimated to be of the order of 100 million [4]. On the other hand, compact objects consisting of exotic matter (axion stars, Bose-Einstein condensate stars, quark stars, etc.) have many properties very similar to those of the stellar type black holes [5, 6]. Hence, finding effective methods for distinguishing between different types of compact objects is a fundamental problem in present-day astronomy and astrophysics.

An important and intriguing class of black holes is represented by the supermassive black holes, which are located at the center of every massive galaxy [7]. A supermassive black hole accretes during its lifetime a significant amount of matter, which leads to the formation of an accretion disk around the compact object. For a review of the supermassive black hole properties, see [8]. If the supermassive black hole is in the accreting state, it is called an active galactic nucleus. The most important observational evidence for the existence of supermassive black holes is provided by the very-long baseline interferometry (VLBI) imaging of molecular $\mathrm{H_2O}$ masers. The VLBI imaging, using Doppler shift measurements that assume a Keplerian type motion of the masering source, led to the very precise estimation of the mass of the black hole at the center of the active galaxy NGC4258, which was obtained as $3.6 \times 10^7M_{\odot }$ [9]. The supermassive black hole closest to Earth is the variable infrared, radio and X-ray source Sagittarius (Sgr) A*, located at the center of the Milky Way [10,11,12]. The investigation of the orbital motion of the stars around Sgr A* lead to further confirmation of the predictions of the general theory of relativity. An important recent advance in black hole physics was the observation by the Event Horizon Telescope collaboration of the first image of the supermassive black hole M87* [13,14,15]. These observations are consistent with the Kerr-like nature of the Sgr A* black hole, but still deviations from the predictions of general relativity cannot be excluded a priori.

Hence, finding black hole solutions is essential for the theoretical understanding and the observational testing of gravitational theories. The first vacuum solution of the general relativistic field equations [2], found by Karl Schwarzschild in 1916, was obtained by considering a static, spherically symmetric central compact object. There are a large number of exact vacuum and black hole solutions in general relativity and in its extensions. For a review of the exact solutions of the Einstein gravitational field equations, see [16]. One effective method to test black hole properties, as well as the nature of the gravitational force, is by using the electromagnetic emissivity properties of thin accretion disks that form around compact astrophysical type objects [17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]. For a review of the observational possibilities of testing black hole geometries by using electromagnetic radiation, see [37].

Black hole solutions have also been investigated in modified theories of gravity. Some investigations in this field can be found in the works [38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75]. Static spherically symmetric solutions of the gravitational field equations obtained by including four-derivative terms of the form $\int {R_{\mu \nu }R^{\mu \nu }\sqrt{-g}d^4x}$ and $\int {R^2\sqrt{-g}d^4x}$, respectively, were considered in [38]. The static, linearized solutions of the field equations are combinations of the Newtonian and Yukawa potentials, while the acceptable static metric solutions in the full nonlinear theory are regular at the origin. Static black hole solutions in Einstein’s gravity with additional quadratic curvature terms were obtained in [43] with the use of a Lichnerowicz type theorem that simplifies the analysis through the condition of a vanishing Ricci scalar curvature. The existence of black hole solutions was proven by using numerical methods, and their thermodynamic properties were also discussed.

A family of exact black hole solutions on a static spherically symmetric background in second-order generalized Proca theories with derivative vector-field interactions coupled to gravity was obtained in [50]. The intrinsic vector-field derivative interactions give rise to a secondary hair, induced by nontrivial field profiles. The deviation from general relativity is most significant around the horizon. The properties of black holes on a static and spherically symmetric background in U(1) gauge-invariant scalar-vector-tensor theories with second-order equations of motion were studied in [62]. It was shown that in shift-symmetric theories invariant under the shift of scalar $\phi \rightarrow \phi +c$, new hairy black hole solutions exist. Vacuum static spherically symmetric solutions in the hybrid metric-Palatini gravity theory, which is a combination of the metric and Palatini f(R) formalisms unifying local constraints at the solar system level and the late-time cosmic acceleration, were considered in [67] by adopting the scalar-tensor representation of the theory, in which the scalar-tensor definition of the potential can be represented as the solution of a Clairaut differential equation. To study the behavior of the metric functions and of the scalar field, a numerical approach was adopted. The static and spherically symmetric solutions in a gravity theory that extends the standard Hilbert–Einstein action with a Lagrangian constructed from a three-form field $A_{\alpha \beta \gamma }$, which is related to the field strength and a potential term, were investigated in [68]. The field equations were obtained explicitly for a static and spherically symmetric geometry in vacuum. Several models corresponding to different functional forms of the three-field potential, including the Higgs and exponential types, were considered.

The effects of the quantum fluctuations on the spherically symmetric static black hole solutions were investigated in [69] by considering a set of modified field equations and by assuming that the quantum correction tensor is given by the coupling of a scalar field to the metric tensor. The formation of a black hole was detected from the presence of a singularity in the metric tensor components. Several models corresponding to different functional forms of the scalar field potential were considered. The thermodynamic properties of the black hole solutions were also investigated. Several classes of exact and perturbative spherically symmetric solutions in f(T, B) gravity theories were considered in [72]. General methods and strategies that allow one to find spherically symmetric static vacuum solutions in modified teleparallel theories of gravity, such as generalized Bianchi identities, were also presented. Charged black hole and wormhole solutions in the Einstein–Maxwell system in the presence of a massless, real scalar field were classified and studied in [73]. Using a conformal transformation, the static, spherically symmetric possible structures in the minimally coupled system were analyzed. Besides wormholes and naked singularities, only a restricted class of black holes does exist, exhibiting a horizon with an infinite surface and a timelike central singularity.

The birth of general relativity did have a deep effect not only on theoretical and observational physics, but also on mathematics. Soon after Einstein proposed his gravitational field equations, Weyl [76, 77] proposed a geometric extension of the Riemannian geometry (on which Einstein’s theory is based), in which the covariant divergence of the metric tensor does not vanish identically, but satisfies a relation of the form $\nabla _{\lambda }g_{\alpha \beta }=Q_{\lambda \alpha \beta }$, where $Q_{\lambda \alpha \beta }$ is the nonmetricity tensor. In the initial formulation of its geometry, Weyl assumed that $Q_{\lambda \alpha \beta }=\omega _{\lambda }g_{\alpha \beta }$, where $\omega _{\lambda }$ was assumed to be the electromagnetic potential. By using this geometric interpretation, Weyl tried to construct a unified theory of gravitation and of the electromagnetic field, which was severely criticized by Einstein. Due to this criticism, this unified theory was later completely abandoned. However, Weyl’s geometry has many attractive features, and it did open some new perspectives in physics, mainly due to the concept of conformal invariance, on which Weyl geometry is based. For a review of Weyl geometry and its applications, see [78].

The present-day physical applications of Weyl geometry are obtained by adopting the essential assumption that all physical laws must satisfy the principle of conformal invariance. A local conformal transformation is defined as $\mathrm{{d}}{\tilde{s}}^2=\Sigma ^n(x)g_{\alpha \beta }\mathrm{{d}}x^{\alpha }\mathrm{{d}}x^{\beta }={\tilde{g}}_{\alpha \beta }\mathrm{{d}}x^{\alpha }\mathrm{{d}}x^{\beta }$, where $\Sigma (x)$ is the conformal factor and n is the Weyl charge. A generalization of Weyl’s gravity was proposed by Dirac [79, 80], and this approach was further investigated in [81, 82]. In Weyl–Dirac type theories, ordinary baryonic matter is generated at the very beginning of the cosmological evolution as a direct consequence of the presence of the Dirac gauge function. On the other hand, in the late universe, Dirac’s gauge function creates the dark energy that triggers the transition to a de Sitter type accelerating phase. In [83] it was suggested that the conformal invariance of the string theory has its origin in the worldvolume diffeomorphism invariance of the membrane theory, a result obtained by postulating an ansatz for the dimensional reduction that leaves the conformal factor as an “integration constant.” For the role of conformal symmetry in field theory and quantum gravity, see [84, 85], respectively.

Conformally invariant gravitational theories can be formulated using actions constructed with the help of the Weyl tensor $C_{\alpha \beta \gamma \delta }$, and which can be obtained from the action $S=-(1/4)\int {C_{\alpha \beta \gamma \delta }C^{\alpha \beta \gamma \delta }\sqrt{-g}\mathrm{{d}}^4x}$ [86,87,88,89,90,91]. Gravitational theories constructed with the help of such actions are called conformally invariant, or Weyl gravity type theories. Another interesting application of Weyl geometry is related to the f(Q) type gravity theories, in which it is assumed that the basic quantity describing the gravitational field is the nonmetricity [92,93,94,95,96,97]. f(Q) or f(Q, T) type gravity theories may represent an attractive alternative to standard general relativity and to the standard $\Lambda$CDM cosmological model.

The idea of conformal invariance also plays a central role in the conformal cyclic cosmology (CCC) theoretical model [98], in which the universe exists as a collection of eons, which are geometric structures corresponding to time-oriented spacetimes. The conformal compactification of eons has spacelike null infinities. Different aspects of the CCC model were investigated in [99,100,101]. In [102] it was proposed that the conformal symmetry is an exact symmetry of a nature that is spontaneously broken, and that it could be as important as the Lorentz invariance of physical laws. The breaking of the conformal invariance could provide a mechanism that would allow the physical understanding of the small-scale structure of the gravitational interaction and also of the Planck-scale physics. A gravitational theory assuming that local conformal symmetry is an exact but spontaneously broken symmetry was introduced in [103].

A conformally invariant Weyl type gravity developed by fully using the geometrical structures of Weyl geometry was investigated in both metric and Palatini formulations in [104,105,106,107,108,109,110,111,112,113]. The physical implications of the theory for the understanding of the elementary particle physics, as well as its cosmological aspects related to the very early evolution of the universe, were studied in detail. The Weyl geometry-inspired quadratic action has a spontaneous symmetry breaking through the Stueckelberg mechanism, which leads to the important result that the Weyl gauge field becomes massive. Hence, in this approach, one re-obtains the standard Hilbert–Einstein action of general relativity in the presence of a positive cosmological constant. Moreover, a Proca type term for the massive Weyl gauge field also appears in the action [104]. A Weyl-geometric invariant Lagrangian without ghosts, given by

$$\begin{aligned} L= & {} \sqrt{-g}\left\{ - \frac{\xi _j}{2} \left[ \frac{1}{6}\, \, \phi ^2_j\,R+ g^{\mu \nu } \, \partial _\mu \phi _j\, \partial _\nu \phi _j\right] \right. \nonumber \\{} & {} \left. +(1+\xi _j)\, \frac{1}{2} g^{\mu \nu } {{\tilde{D}}}_\mu \phi _j\,{{\tilde{D}}}_\nu \phi _j -V(\phi _j)\right\} , \end{aligned}$$

(1)

was investigated in [105], where a potential term $V(\phi _j)$ for the scalars $\phi _j$ was also considered, with V given by a homogeneous function of the form $V(\phi _j)=\phi _k^4\, V(\phi _j/\phi _k), k=\textrm{fixed}$.

It turns out that in the Weyl geometric approach, a successful description of the early universe inflation can be achieved if one of the scalar fields is considered as the inflaton field [104,105,106]. Inflation in the Weyl geometry-inspired gravity in the presence of a scalar field gives similar results as in the Starobinsky model [114], which is recovered in the limit of the vanishing of the non-minimal coupling term [106]. Hence, Weyl geometry could play a fundamental role in the early universe, once one assumes that the effective theory at short distances is conformally invariant. Moreover, in [107] it was shown that after a particular Weyl gauge transformation (or gauge fixing), Weyl conformal geometry acquires a geometric Stueckelberg mechanism, which is naturally built in the theory. This mechanism is broken spontaneously, and leads to the appearance of the Riemannian geometry. On the other hand, the Stueckelberg mechanism conserves the total number of the degrees of freedom, only rearranging them.

In the Palatini formalism, quadratic Weyl gravity with Lagrangian $R^2+R_{\mu \nu }^2$ was studied in [108]. In this approach, the connection and the metric are considered as independent variables, with the action having a gauged scale symmetry. For non-minimally coupled Higgs-like fields, the theory can describe inflation well. In [109], a comparative study of inflation in the original Weyl quadratic gravity and in a theory obtained by replacing the Weyl connection by its Palatini counterpart was considered. After the Weyl vector becomes massive via the Stueckelberg mechanism, the Planck scale, the metricity condition, and the Einstein–Proca action of the Weyl field arise in the broken phase. For large Higgs fields, inflation is also possible. In [111], the cosmological dynamics of the Weyl geometry-inspired quadratic gravitational model was investigated in detail.

The coupling of matter to geometry in Weyl geometry-inspired conformal quadratic gravity, by assuming a coupling term of the form $L_m{\tilde{R}}^ 2$, where $L_m$ is the ordinary matter Lagrangian and ${\tilde{R}}$ is the Weyl scalar, was investigated in [113]. The coupling explicitly satisfies the conformal invariance of the theory. By expressing ${\tilde{R}}^2$ with the help of an auxiliary scalar field and of the Weyl scalar, the gravitational action can be linearized, leading in the Riemannian space to a conformally invariant $f\left( R,L_m\right)$ type theory, with the matter Lagrangian non-minimally coupled to the Ricci scalar. The cosmological implications of the theory were also considered for the case of a flat, homogeneous and isotropic geometry. The model can give a good description of the observational data for the Hubble function up to a redshift of the order of $z \approx 3$.

It is the goal of the present work to investigate black hole solutions in the Weyl type geometric gravity theories by following the approach introduced and developed in [104,105,106,107,108,109,110,111,112,113]. We assume, as a starting point of our investigations, that the background geometry of the spacetime is of Weyl type. Moreover, we also require the conformal invariance of any gravitational theory. To implement these requirements, we adopt the simplest possible model of a conformally invariant Weyl geometric gravity model, by assuming that the Lagrangian density is constructed from the square of the Weyl scalar, and of the electromagnetic type tensor $F_{\mu \nu }$ only. This Lagrangian density can be linearized with the help of an auxiliary scalar field, and finally it can be written as an Einstein–Proca type Lagrangian in the ordinary Riemannian space, in the presence of a non-minimal coupling between the scalar field and the Ricci scalar, and the Weyl vector, respectively.

As the next step in our study, after obtaining from the considered variational principle the gravitational field equations, we consider a spherically symmetric static metric and write down the corresponding vacuum field equations. Due to their mathematical complexity and their strong nonlinear nature, the differential equations describing the vacuum solutions of the theory can be generally investigated only numerically. We consider three classes of models, corresponding to the three possible choices of the Weyl vector field, assumed to have only one timelike component, only one radial spacelike component, and both components. The field equations are rewritten in a dimensionless form, and, after introducing a new independent radial coordinate, defined as the inverse of r, we investigate their solutions numerically. Thus, three distinct classes of black hole models are obtained. We detect the presence of an event horizon from the appearance of the singularities in the metric tensor components. An exact solution of the gravitational field equations, obtained in the presence of the radial component of the Weyl vector, is also obtained. The thermodynamic properties of the Weyl geometric type black holes, including their horizon temperature, specific heat, and entropy, as well as the evaporation time due to Hawking luminosity, are also considered in detail by using numerical methods for all considered cases.

The present paper is organized as follows. We review the foundations of the Weyl geometry, we introduce the gravitational action, and we obtain the general field equations of the Weyl geometric gravity in Sect. 2. The spherically symmetric static field equations in the Weyl geometric model of gravity are obtained in Sect. 3 for three distinct choices of the functional form of the Weyl vector $\omega ^{\mu }$. Numerical black hole solutions in Weyl geometric gravity, corresponding to the different choices of the Weyl vector, are considered in Sect. 4. An exact solution of the static spherically symmetric field equations in the presence of a Weyl vector having only a radial spacelike component is also obtained. The thermodynamic properties of the obtained black hole solutions are considered in Sect. 5. We discuss our results and conclusions in Sect. 6. In this paper we use a system of units with $c = 1$.

2 Gravitational field equations and equations of motion in Weyl conformal gravity

We begin the present section with a brief discussion of the basic elements of Weyl geometry. Then we proceed to introduce the simplest Weyl geometry-inspired conformally invariant gravitational action. After linearizing the quadratic Weyl action, we obtain the gravitational field equations describing the dynamical evolution of the metric, the auxiliary scalar field, and the Weyl vector in their general form.

2.1 A quick introduction to Weyl geometry

A basic property of a manifold is the existence of an intrinsic metric that gives the distance between two infinitesimally close points. In Riemannian geometry, the distance is preserved under the parallel transport with respect to the metric-compatible connection, called the Levi-Civita connection, or the Christoffel symbols, if we refer to its components. Shortly after Einstein and Hilbert obtained the correct form of the general relativistic field equations, with Hilbert introducing the variational principle for their derivations, Weyl [76, 77] proposed the existence of a conformal transformation in every point of the spacetime manifold. In the case of the metric, such a conformal transformation takes the form ${\tilde{g}}_{\mu \nu } =\Sigma ^{n}(x)g_{\mu \nu }$, where n is the Weyl charge. In the following we will consider only the case $n=1$.

In Weyl geometry, under parallel transport, the length $\ell$ of a vector varies when it is parallelly transported by an infinitesimal displacement between the points $x^\mu$ and $x^\mu + \delta x^\mu$. The change in the length of a vector is given by

$$\begin{aligned} \delta \ell = \ell \omega _\mu \delta x^\mu , \end{aligned}$$

(2)

where $\omega _\mu$ is the Weyl vector field. In the general case, one can introduce the nonmetricity $Q_{\lambda \mu \nu }$ of the Weyl geometry, via the covariant derivative of the metric tensor, according to the definition

$$\begin{aligned}&{\tilde{\nabla }}_\lambda g_{\mu \nu } = - \alpha \omega _\lambda g_{\mu \nu }\equiv Q_{\lambda \mu \nu }, \end{aligned}$$

(3)

where $\alpha$ is the Weyl gauge coupling constant. For the connection of the Weyl geometry, from the nonmetricity condition (3) we obtain the expression

$$\begin{aligned}&{\tilde{\Gamma }}^\lambda _{\mu \nu } = \Gamma ^\lambda _{\mu \nu } + \frac{1}{2} \alpha \Big [ \delta ^\lambda _\mu \omega _\nu + \delta ^\lambda _\nu \omega _\mu - g_{\mu \nu } \omega ^\lambda \Big ], \end{aligned}$$

(4)

where $\Gamma ^\lambda _{\mu \nu }$ is the standard Levi-Civita connection associated with the metric g. In the following, all geometrical and physical quantities defined in the Weyl geometry will be denoted by a tilde.

By using the Weyl connection, one can easily construct the curvature tensor ${\tilde{R}}_{\mu \nu \sigma }^{\lambda }$, defined as

$$\begin{aligned} {\tilde{R}}_{\mu \nu \sigma }^{\lambda }=\partial _{\nu }{\tilde{\Gamma }}_{\mu \sigma }^{\lambda }-\partial _{\sigma }{\tilde{\Gamma }}_{\mu \nu }^{\lambda }+ {\tilde{\Gamma }}_{\rho \nu }^{\lambda }{\tilde{\Gamma }}_{\mu \sigma }^{\rho }- {\tilde{\Gamma }}_{\rho \sigma }^{\lambda }{\tilde{\Gamma }}_{\mu \nu }^{\rho }, \end{aligned}$$

(5)

its first contraction,

$$\begin{aligned} {\tilde{R}}_{\mu \nu }={\tilde{R}}_{\mu \lambda \nu }^{\lambda },{\tilde{R}} =g^{\mu \sigma }{\tilde{R}}_{\mu \sigma }, \end{aligned}$$

(6)

and its second contraction (the Weyl scalar), respectively, given by

$$\begin{aligned} {{\tilde{R}}} = g^{\mu \nu } \Big ( \partial _\rho {\tilde{\Gamma }}^\rho _{\mu \nu } - \partial _\nu {\tilde{\Gamma }}^\rho _{\mu \rho } + {\tilde{\Gamma }}^\rho _{\mu \nu } {\tilde{\Gamma }}^\sigma _{\rho \sigma } - {\tilde{\Gamma }}^\sigma _{\mu \rho }\tilde{ \Gamma }^\rho _{\nu \sigma } \Big ), \end{aligned}$$

(7)

or, equivalently, by

$$\begin{aligned} {\tilde{R}}=R-3n\alpha \nabla _{\mu }\omega ^{\mu }-\frac{3}{2}\left( n\alpha \right) ^{2}\omega _{\mu }\omega ^{\mu }, \end{aligned}$$

(8)

where R is the Ricci scalar defined in the Riemannian geometry.

Another important geometrical and physical quantity, the Weyl tensor, is defined according to

$$\begin{aligned} {\tilde{C}}_{\mu \nu \rho \sigma }^{2}=C_{\mu \nu \rho \sigma }^{2}+\frac{3}{2} \left( \alpha n\right) ^{2}{\tilde{F}}_{\mu \nu }^{2}, \end{aligned}$$

(9)

where $C_{\mu \nu \rho \sigma }$ is the Weyl tensor as introduced in the standard Riemannian geometry [115]. $C_{\mu \nu \rho \sigma }^{2}$ can be obtained in terms of the Riemann and Ricci tensors and of the Ricci scalar as

$$\begin{aligned} C_{\mu \nu \rho \sigma }^{2}=R_{\mu \nu \rho \sigma }R^{\mu \nu \rho \sigma }-2R_{\mu \nu }R^{\mu \nu }+\frac{1}{3}R^{2}. \end{aligned}$$

(10)

If we apply a conformal transformation with a conformal factor $\Sigma$ in one point, the variances of the metric tensor, the Weyl field, and a scalar field $\phi$ are given by

$$\begin{aligned} {\hat{g}}_{\mu \nu } = \Sigma g_{\mu \nu }, {\hat{\omega }} _\mu = \omega _\mu - \frac{1}{\alpha } \partial _\mu \ln \Sigma , {\hat{\phi }} = \Sigma ^{-\frac{1}{2 }} \phi . \end{aligned}$$

(11)

With the help of the Weyl vector, we can construct the Weyl field strength $F_{\mu \nu }$ of $\omega _\mu$, defined as

$$\begin{aligned} {\tilde{F}}_{\mu \nu } = {\tilde{\nabla }}_{[\mu } \omega _{\nu ]} = \nabla _{[\mu } \omega _{\nu ]} = \partial _{[\mu } \omega _{\nu ]}=\partial _{\mu }\omega _\nu -\partial _\nu \omega _\mu . \end{aligned}$$

(12)

2.2 Action and field equations

In the following, we consider the simplest conformally invariant Lagrangian density that can be constructed in Weyl geometry, and which is given by [106,107,108,109,110]

$$\begin{aligned} L_0=\left[ \, \frac{1}{4!}\,\frac{1}{\xi ^2}\,{{\tilde{R}}}^2 - \frac{1}{4}\, F_{\mu \nu }^{\,2} \right] \sqrt{-g}, \end{aligned}$$

(13)

with perturbative coupling $\xi < 1$. In order to extract the scalar degree of freedom of the above Lagrangian, and to linearize it, in $L_0$ we replace ${\tilde{R}}^2$ with ${\tilde{R}}^2\rightarrow 2 \phi _0^2\,{{\tilde{R}}}-\phi _0^4$, where $\phi _0$ is an auxiliary scalar field. The new Lagrangian density is equivalent to the initial one, since by using the solution $\phi _0^2={{\tilde{R}}}$ of the equation of motion of $\phi _0$ in the new $L_0$, we recover Eq. (13). Hence, we obtain a mathematically equivalent Weyl geometric Lagrangian,

$$\begin{aligned} L_0=\sqrt{-g} \left[ \frac{1}{12}\frac{1}{\xi ^2}\,\phi _0^2\,{{\tilde{R}}} -\frac{1}{4} \,F_{\mu \nu }^2-\frac{\phi _0^4}{4!\,\xi ^2} \right] . \end{aligned}$$

(14)

This is the simplest possible gravitational Lagrangian density with Weyl gauge symmetry and conformal invariance fully implemented. $L_0$ has a spontaneous breaking to an Einstein–Proca Lagrangian of the Weyl gauge field. On the other hand, we would like to point out that the simplest gravitational action with conformal symmetry, realized locally, is the conformal gravity model based on the $C_{\alpha \beta \gamma \delta }^2$ action [86,87,88,89,90,91], with the Weyl tensor $C_{\alpha \beta \gamma \delta }^2$ defined in the standard Riemannian geometry and given by Eq.  (10). Note that the tensor defined by Eq. (10) is different from the Weyl geometric tensor as introduced in Eq. (9). In the simple model of conformal Weyl gravity, no Weyl gauge field $\omega _\mu$ or scalar $\phi$ is present. However, in four dimensions, the theory is still invariant under local Weyl gauge transformations. Generally, a physical system is called conformally invariant if it satisfies the condition that the variation of the action $S\left[ g_{\mu \nu },\phi \right]$ with respect to the group of the conformal transformations is zero, $\delta _cS\left[ g_{\mu \nu },\phi \right] =\int {d^n x\left( \delta S/\delta \phi \right) \delta _c\phi }=0$ [116]. Another type of transformation is given by the Weyl rescaling, implying the simultaneous pointwise transformations of the metric and of the fields ${\tilde{g}}_{\mu \nu }(x)=e^{2\sigma (x)}g_{\mu \nu }$ and ${\tilde{\phi }}=e^{-\Delta \sigma (x)}\phi (x)$, under which the action transforms as $\delta _\sigma S\left[ g_{\mu \nu },\phi \right] =\int {d^n x\sigma \left[ 2\left( \delta S/\delta g_{\mu \nu }\right) g_{\mu \nu }-\Delta _n\left( \delta S/\delta \phi \right) \phi \right] }$ [116]. For an in-depth discussion of the conformal and Weyl invariance, see [116].

By replacing in Eq. (14) ${\tilde{R}}$ as given by Eq. (8), and by performing a gauge transformation that allows the redefinition of the variables, we obtain an action, defined in the Riemannian space, invariant under conformal transformation, given by [106,107,108]

$$\begin{aligned} {\mathcal {S}}= & {} \int d^4x \sqrt{-g} \Bigg [ \frac{1}{12} \frac{\phi ^2}{\xi ^2} \Big ( R - 3\alpha \nabla _\mu \omega ^\mu - \frac{3}{2} \alpha ^2 \omega _\mu \omega ^\mu \Big ) \nonumber \\{} & {} - \frac{1}{4!}\frac{\phi ^4}{\xi ^2} - \frac{1}{4} F_{\mu \nu } F^{\mu \nu } \Bigg ]. \end{aligned}$$

(15)

The variation of the action Eq. (15) with respect to the metric tensor gives the field equation

$$\begin{aligned}{} & {} \frac{\phi ^{2}}{\xi ^{2}}\Big (R_{\mu \nu }-\frac{1}{2}Rg_{\mu \nu } \Big )\nonumber \\{} & {} \quad -\frac{3\alpha }{2\xi ^{2}}\Big (\omega ^{\rho }\nabla _{\rho }\phi ^{2}g_{\mu \nu }-\omega _{\nu }\nabla _{\mu }\phi ^{2}-\omega _{\mu }\nabla _{\nu }\phi ^{2}\Big )\nonumber \\{} & {} \quad +\frac{3\alpha ^{2}}{4\xi ^{2}}\phi ^{2}\Big (\omega _{\rho }\omega ^{\rho }g_{\mu \nu }-2\omega _{\mu }\omega _{\nu }\Big ) -6F_{\rho \mu }F_{\sigma \nu }g^{\rho \sigma }\nonumber \\{} & {} \quad +\frac{3}{2}F_{\rho \sigma }^{2}g_{\mu \nu }+\frac{1}{4\xi ^{2}}\phi ^{4}g_{\mu \nu }+\frac{1}{\xi ^{2}} \Big (g_{\mu \nu }\Box -\nabla _{\mu }\nabla _{\nu }\Big )\phi ^{2}=0.\nonumber \\ \end{aligned}$$

(16)

The trace of Eq. (16) gives

$$\begin{aligned} \Phi R+3\alpha \omega ^{\rho }\nabla _{\rho }\Phi -\Phi ^{2}-\frac{3}{2} \alpha ^{2}\Phi \omega _{\rho }\omega ^{\rho }-3\Box \Phi =0, \end{aligned}$$

(17)

where we have denoted $\Phi \equiv \phi ^{2}$. Variation of the action given by Eq. (15) with respect to the scalar field $\phi$ gives the equation of motion of $\phi$,

$$\begin{aligned} R-3\alpha \nabla _{\rho }\omega ^{\rho }-\frac{3}{2}\alpha ^{2}\omega _{\rho }\omega ^{\rho }-\Phi =0. \end{aligned}$$

(18)

From Eqs. (17) and (18) it immediately follows that

$$\begin{aligned} \Box \Phi -\alpha \nabla _{\rho }(\Phi \omega ^{\rho })=0. \end{aligned}$$

(19)

The variation of Eq. (15) with respect to $\omega _{\mu }$ gives the equation of motion of the Weyl vector as

$$\begin{aligned} 4\xi ^{2}\nabla _{\nu }F^{\mu \nu }+\alpha ^{2}\Phi \omega ^{\mu }-\alpha \nabla ^{\mu }\Phi =0. \end{aligned}$$

(20)

Applying to both sides of the above equation the operator $\nabla _{\mu }$, we obtain Eq. (19), a result that can be seen as a consistency check of the theory.

3 The spherically symmetric vacuum field equations

In the following, we adopt a static and spherically symmetric geometry in which all quantities depend only on the radial coordinate, r, in a system of coordinates defined according to $\left( t,r,\theta , \varphi \right)$. For mathematical convenience we will use geometric units in all equations. In the following, we denote by a prime the derivative with respect to the radial coordinate r. Hence, in the adopted geometry, the spacetime interval becomes

$$\begin{aligned} \mathrm{{d}}s^{2}=e^{\nu (r)}\mathrm{{d}}t^{2}-e^{\lambda (r)}\mathrm{{d}}r^{2}-r^{2}\mathrm{{d}}\Omega ^{2}, \end{aligned}$$

(21)

where we have denoted $\mathrm{{d}}\Omega ^{2}=\mathrm{{d}}\theta ^{2}+\sin ^{2}\theta \mathrm{{d}}\varphi ^{2}$. Thus, the metric tensor components are given by

$$\begin{aligned} g_{\mu \nu }=\textrm{diag}(e^{\nu (r)},-e^{\lambda (r)},-r^{2},-r^{2}\sin ^{2}\theta ). \end{aligned}$$

(22)

If we consider a spherically symmetric configuration, the third and the fourth components of the Weyl vector vanish identically. Hence, the Weyl vector can be represented in a general form as $\omega _{\mu }=\left( \omega _{0},\omega _{1},0,0\right)$.

In the following, we will write down the field equations resulting from the three possibilities allowed by the various choices of the Weyl vector, and we also present some of their consequences.

3.1 The case $\omega _{\mu }=\left( 0,\omega _{1},0,0\right)$

In the first case we are considering, the temporal component, $\omega _{0}$, of the Weyl vector can be taken as zero, and thus $\omega _{\mu }$ is given by $\omega _{\mu }=(0,\omega _{1},0,0)$. For this choice of $\omega _{\mu }$, we have $F_{\mu \nu }\equiv 0$. Hence, Eq. (20) immediately gives

$$\begin{aligned} \Phi '=\alpha \Phi \omega _1. \end{aligned}$$

(23)

By taking into account that

$$\begin{aligned} \Box \Phi = \frac{1}{\sqrt{-g}}\frac{\partial }{\partial x^{\alpha }}\left( \sqrt{-g} g^{\alpha \beta }\frac{\partial \Phi }{\partial x^{\beta }}\right) \end{aligned}$$

(24)

and

$$\begin{aligned} \nabla _{\alpha }\omega ^{\alpha }=\frac{1}{\sqrt{-g}}\frac{\partial }{\partial x^{\alpha }}\left( \sqrt{-g}\omega ^{\alpha }\right) , \end{aligned}$$

(25)

respectively, Eq. (19) becomes

$$\begin{aligned}{} & {} \frac{1}{\sqrt{-g}}\frac{\mathrm{{d}}}{\mathrm{{d}}r}\left( \sqrt{-g}g^{11}\frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r}\right) -\alpha \omega ^{1}\frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r} -\frac{\alpha \Phi }{\sqrt{-g}}\frac{\mathrm{{d}}}{\mathrm{{d}}r} \left( \sqrt{-g}\omega ^{1}\right) =0.\nonumber \\ \end{aligned}$$

(26)

In order to make Eqs. (23) and (26) consistent, we need to impose the gauge condition on the Weyl vector, which generally can be formulated as

$$\begin{aligned} \nabla _{\alpha }\omega ^{\alpha }=\frac{1}{\sqrt{-g}}\frac{\partial }{ \partial x^{\alpha }}\left( \sqrt{-g}\omega ^{\alpha }\right) =0, \end{aligned}$$

(27)

giving for the present choice of $\omega _{\mu }$ the relation

$$\begin{aligned} \frac{1}{\sqrt{-g}}\frac{\mathrm{{d}}}{\mathrm{{d}}r}\left( \sqrt{-g}\omega ^{1} \right) =0 \end{aligned}$$

(28)

or

$$\begin{aligned} \omega ^{1}=\frac{C_{1}}{\sqrt{-g_\mathrm{{rad}}}}=C_{1}\frac{e^{-{(\nu +\lambda )}/2} }{r^{2}} \end{aligned}$$

(29)

and

$$\begin{aligned} \omega _{1}=g_{11}\omega ^{1}=-C_{1}\frac{e^{-\left( \nu -\lambda \right) /2}}{r^{2}}, \end{aligned}$$

(30)

respectively, where $C_{1}$ is an arbitrary integration constant, and we have denoted $\sqrt{-g_\mathrm{{rad}}}=e^{\left( \nu +\lambda \right) /2}r^{2}$. By taking into account the gauge condition and the explicit form of $\omega ^{1}$, Eq. ( 26) becomes

$$\begin{aligned} \frac{1}{\sqrt{-g}}\frac{\mathrm{{d}}}{\mathrm{{d}}r}\left( \sqrt{-g}g^{11}\frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r}\right) -\alpha \omega ^{1}\frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r}=0 \end{aligned}$$

(31)

or, equivalently,

$$\begin{aligned} \frac{1}{\sqrt{-g_\mathrm{{rad}}}}\frac{\mathrm{{d}}}{\mathrm{{d}}r}\left( \sqrt{-g_\mathrm{{rad}}}g^{11}\frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r }\right) -\alpha \frac{C_{1}}{\sqrt{-g_\mathrm{{rad}}}}\frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r}=0, \end{aligned}$$

(32)

giving

$$\begin{aligned} \sqrt{-g_\mathrm{{rad}}}g^{11}\frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r}-\alpha C_{1}\Phi =C, \end{aligned}$$

(33)

where C is an arbitrary integration constant which for consistency with Eq. (23) must be taken as zero. Hence, we re-obtain Eq. (23), $\Phi ^{\prime }=\alpha \omega _{1}\Phi$.

By taking into account the explicit expression of $\omega ^{1}$ as given by Eq. (29), we obtain for $\Phi$ the differential equation

$$\begin{aligned} \frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r}=-\alpha C_{1}\frac{e^{-\left( \nu -\lambda \right) /2}}{r^{2}}\Phi . \end{aligned}$$

(34)

Alternatively, from Eqs. (31) and (23) we obtain a nontrivial dynamical equation for $\Phi (r)$,

$$\begin{aligned} \frac{\mathrm{{d}}}{\mathrm{{d}}r}\left( \sqrt{-g}g^{11}\frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r} \right) - \sqrt{-g} g^{11} \frac{1}{\Phi } \left( \frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r}\right) ^2 = 0. \end{aligned}$$

(35)

Next, the gravitational field equations give the relations

$$\begin{aligned}{} & {} -1+e^{\lambda }-\frac{1}{4}e^{\lambda }r^{2}\Phi -\frac{2r\Phi ^{\prime }}{\Phi }+ \frac{3r^{2}}{4}\frac{\Phi ^{\prime 2}}{\Phi ^{2}}+r\lambda ^{\prime }\nonumber \\{} & {} \qquad +\frac{ r^{2}\lambda ^{\prime }}{2}\frac{\Phi ^{\prime }}{\Phi }-\frac{r^{2}\Phi ^{\prime \prime }}{\Phi }=0, \end{aligned}$$

(36)

$$\begin{aligned}{} & {} 1-e^{\lambda }+\frac{1}{4}e^{\lambda }r^{2}\Phi +\frac{2r\Phi ^{\prime }}{\Phi }+\frac{ 3r^{2}}{4}\frac{\Phi ^{\prime 2}}{\Phi ^{2}}+r\nu ^{\prime }\left( 1 +\frac{r}{2}\frac{\Phi ^{\prime }}{\Phi }\right) \nonumber \\{} & {} \quad =0 \end{aligned}$$

(37)

and

$$\begin{aligned}{} & {} 2(\nu ^{\prime }-\lambda ^{\prime })+(4-2r\lambda ^{\prime }+2r\nu ^{\prime })\frac{\Phi ^{\prime }}{\Phi }\nonumber \\{} & {} \quad +r\left( e^{\lambda }\Phi +4\frac{\Phi ^{\prime \prime }}{\Phi }-3\frac{\Phi ^{\prime 2}}{\Phi ^{2}}-\lambda ^{\prime }\nu ^{\prime }+\nu ^{\prime 2}+2\nu ^{\prime \prime }\right) =0,\nonumber \\ \end{aligned}$$

(38)

respectively.

It can be easily proved that Eq. (38) is a linear combination of the other two field equations. By adding Eqs. (36) and (37), we obtain

$$\begin{aligned} \frac{\Phi ''}{\Phi }-\frac{3}{2}\frac{\Phi ^{\prime 2}}{\Phi ^2}-\frac{\nu '+\lambda '}{2}\frac{\Phi '}{\Phi }-\frac{\nu '+\lambda '}{r}=0. \end{aligned}$$

(39)

3.2 The case $\omega _{\mu }=\left( \omega _0,0,0,0\right)$

We now consider the case in which it is possible to neglect all the effects of the spatial components $\omega _1$ in the Weyl vector. Thus, only the $\omega _0$ component of the Weyl vector is nonzero. Therefore, Eq. (20) immediately gives

$$\begin{aligned}&-\frac{\alpha ^2}{4\xi ^2} e^\lambda \Phi \omega _0 + \frac{-4+r \lambda ^\prime + r \nu ^\prime }{2r}\omega _0^\prime - \omega ^{\prime \prime }_0 = 0, \end{aligned}$$

(40)

and

$$\begin{aligned} \Phi ^\prime = 0, \end{aligned}$$

(41)

respectively. Hence, the gauge condition imposed on the Weyl vector field in the former section is now fulfilled trivially. In this case, the gravitational field equations give the following relations

$$\begin{aligned}{} & {} -1 + e^\lambda + \frac{1}{4} e^\lambda r^2 \Phi - \frac{3}{4} \alpha ^2 e^{\lambda -\nu } r^2 \omega _0^2 + r \lambda ^\prime + \frac{3\xi ^2 e^{-\nu } r^2 \omega _0^{\prime 2}}{\Phi }\nonumber \\{} & {} \quad = 0, \end{aligned}$$

(42)

$$\begin{aligned}{} & {} 1 - e^\lambda - \frac{1}{4} e^\lambda r^2 \Phi - \frac{3}{4} \alpha ^2 e^{\lambda -\nu } r^2 \omega _0^2 + r \nu ^\prime - \frac{3\xi ^2 e^{-\nu } r^2 \omega _0^{\prime 2}}{\Phi }\nonumber \\{} & {} \quad =0. \end{aligned}$$

(43)

It is also easy to prove that Eq. (18) gives a trivial constraint for $\Phi$, if Eq. (41) is used to express $\Phi$.

3.3 The case $\omega _{\mu }=\left( \omega _0,\omega _1,0,0\right)$

In this case, Eq. (20) gives

$$\begin{aligned} \Phi '=\alpha \Phi \omega _1 \end{aligned}$$

(44)

and

$$\begin{aligned}&- \frac{\alpha ^2}{4\xi ^2} e^\lambda \Phi \omega _0 + \frac{-4+r \lambda ^\prime + r \nu ^\prime }{2r}\omega _0^\prime - \omega ^{\prime \prime }_0 = 0, \end{aligned}$$

(45)

respectively.

By substituting $\omega _1$ from Eq. (44), one can simplify the field equations Eq. (16) as

$$\begin{aligned}&-1+e^{\lambda }+\frac{1}{4}e^{\lambda }r^{2}\Phi -\frac{2r\Phi ^{\prime }}{\Phi }+ \frac{3r^{2}}{4}\frac{\Phi ^{\prime 2}}{\Phi ^{2}}+r\lambda ^{\prime }\nonumber \\&\quad +\frac{ r^{2}\lambda ^{\prime }}{2}\frac{\Phi ^{\prime }}{\Phi }-\frac{r^{2}\Phi ^{\prime \prime }}{\Phi }-\frac{3}{4}\alpha ^2 e^{\lambda -\nu }r^2\omega _0^2+\frac{3\xi ^2}{\Phi }e^{-\nu }r^2\omega ^{\prime 2}_0=0 \end{aligned}$$

(46)

and

$$\begin{aligned}{} & {} 1-e^{\lambda }-\frac{1}{4}e^{\lambda }r^{2}\Phi +\frac{2r\Phi ^{\prime }}{\Phi }+\frac{ 3r^{2}}{4}\frac{\Phi ^{\prime 2}}{\Phi ^{2}}+r\nu ^{\prime }\nonumber \\{} & {} \quad +\frac{r^{2}\nu ^{\prime }}{2}\frac{\Phi ^{\prime }}{\Phi }-\frac{3}{4}\alpha ^2 e^{\lambda -\nu }r^2\omega _0^2-\frac{3\xi ^2}{\Phi }e^{-\nu }r^2\omega ^{\prime 2}_0=0, \end{aligned}$$

(47)

respectively. By adding Eqs. (46) and (47), one obtains

$$\begin{aligned} \frac{\Phi ''}{\Phi }-\frac{3}{2}\frac{\Phi ^{\prime 2}}{\Phi ^2}-\frac{\nu '+\lambda '}{2}\frac{\Phi '}{\Phi }-\frac{\nu '+\lambda '}{r}+\frac{3}{2}\alpha ^2 e^{\lambda -\nu }\omega _0^2=0. \end{aligned}$$

(48)

3.4 Asymptotic and near-horizon behavior

We now consider the asymptotic and near-horizon behavior of the field equations of the geometric Weyl gravity. In particular, the asymptotic values of the metric tensor components and of the scalar and Weyl vector fields allow us to fix the conditions at infinity for the numerical integration and analysis of the field equations. We will consider independently the three cases corresponding to the three different choices of the Weyl vector considered in the previous subsection. We assume that at infinity, the metric can either be asymptotically flat, satisfying the conditions

$$\begin{aligned} \lim _{r\rightarrow \infty }\nu (r)=0, \lim _{r\rightarrow \infty }\lambda (r)=0 \end{aligned}$$

(49)

and

$$\begin{aligned} \lim _{r\rightarrow \infty }\nu '(r)=0, \lim _{r\rightarrow \infty }\lambda '(r)=0, \end{aligned}$$

(50)

respectively, or of the de Sitter type, so that

$$\begin{aligned} \lim _{r\rightarrow \infty }e^{\nu (r)}=\lim _{r\rightarrow \infty }e^{-\lambda }=1-\frac{r^2}{\mu ^2} \end{aligned}$$

(51)

and

$$\begin{aligned} \lim _{r\rightarrow \infty }\nu '(r)=-\lim _{r\rightarrow \infty }\lambda '(r)=-\frac{2r/\mu ^2}{1-r^2/\mu ^2}, \end{aligned}$$

(52)

respectively, where $\mu$ is a constant. In both cases, at infinity, the metric tensor components satisfy the conditions $\lim _{r\rightarrow \infty }\left( \nu (r)+\lambda (r)\right) =0$, and $\lim _{r\rightarrow \infty }\left( \nu '(r)+\lambda '(r)\right) =0$. For a correct and complete definition of the asymptotic limits of the components of the metric tensor, one should also consider some bounds of these functions, such as that for large r, they fall faster than 1/r.

In the general case, the static spherically symmetric field equations of the Weyl geometric gravity are too complicated to be solved analytically. Therefore, to obtain their solutions, we must resort to numerical methods. In our investigations, we assume that the field equations have a black hole solution, whose presence is indicated by the presence of a horizon at a radius $r = r_0 > 0$, where the metric functions $e^{\nu }$ and $e^{\lambda }$ become singular. Then, near the horizon, the metric functions can be approximated by their Taylor expansions [38, 43],

$$\begin{aligned} e^{\nu }=A_1\left( r-r_0\right) +A_2\left( r-r_0\right) ^2+A_3\left( r-r_0\right) ^3+.... \end{aligned}$$

(53)

and

$$\begin{aligned} e^{-\lambda }=B_1\left( r-r_0\right) +B_2\left( r-r_0\right) ^2+B_3\left( r-r_0\right) ^3+...,\nonumber \\ \end{aligned}$$

(54)

respectively, where $A_i$ and $B_i$, $i=1,2,3,...$ are constants that can be determined recursively after substitution into the field equations.

3.4.1 Models with radial component of the Weyl vector

a. Asymptotic limits. In the presence of a Weyl vector having only a radial component, in the limit of large r, for both the asymptotically flat and de Sitter cases, Eq. (39) takes the form

$$\begin{aligned} \Phi \Phi ''=\frac{3}{2} \Phi ^{\prime 2}, \end{aligned}$$

(55)

with the solution

$$\begin{aligned} \Phi (r)= \frac{K_1}{\left( r+K_2\right) ^2}, \end{aligned}$$

(56)

where $K_1$ and $K_2$ are arbitrary integration constants. Hence, for the scalar and the vector fields, we obtain the general asymptotic conditions

$$\begin{aligned} \lim _{r\rightarrow \infty }\Phi (r)= & {} 0, \lim _{r\rightarrow \infty }\Phi ' (r)=0, \nonumber \\ \lim _{r\rightarrow \infty }\omega _1 (r)= & {} -\frac{2}{\alpha }\lim _{r\rightarrow \infty }\frac{1}{r+K_2}=0. \end{aligned}$$

(57)

In the asymptotically flat case, Eq. (36) gives $-16K_2r+\left( K_1-4\right) r^2\approx 0$, which is identically satisfied for $K_2=0$ and $K_1=4$, respectively. For the de Sitter type asymptotic behavior, we find $-16K_2r+r^2\left[ -12K_2^2+\left( K_1-4\right) \alpha ^2\right] /\alpha ^2\approx 0$, a condition that is again satisfied for $K_2=0$ and $K_1=4$.

However, in the asymptotic limit, Eq. (39) also has the solution $\Phi =\textrm{constant}$, implying that at infinity, the scalar field can take arbitrary values. By taking into account the relation between the scalar field and the radial component of the Weyl vector, as given by Eq. (23), it follows that $\lim _{r\rightarrow \infty }\omega _1(r)=0$; that is, the Weyl vector vanishes at infinity, independently of the asymptotic values of $\Phi$.

b. Near-horizon behavior. In order to consider the near-horizon behavior of the model, we also assume that near $r_0$, the metric tensor components and the scalar field admit Taylor expansions of the form [38, 43]

$$\begin{aligned} f_i(r)=K_{i1}\left( r-r_0\right) +K_{i2}\left( r-r_0\right) ^2+K_{i3}\left( r-r_0\right) ^3+...,\nonumber \\ \end{aligned}$$

(58)

where $f_i=\left\{ e^{\nu }, e^{\lambda }, \Phi \right\}$, while $K_{ij}$, $i,j=1,2,3,...$ are constants.

As a simple example of the near-horizon behavior of Weyl geometric black holes in the presence of the radial component of the Weyl vector, we assume that the metric tensor components can be represented near the singularity as

$$\begin{aligned} e^{\nu (r)}=\left( r-r_0\right) ^2\left[ K_{11}+K_{12}\left( r-r_0\right) \right] ^2 \end{aligned}$$

(59)

and

$$\begin{aligned} e^{\lambda (r)}=\left( r-r_0\right) ^2\left[ K_{21}+K_{22}\left( r-r_0\right) \right] ^2, \end{aligned}$$

(60)

respectively, where $K_{ij}$, $i,j=1,2$ are constants. For the near-horizon behavior of the scalar field, we adopt the functional form

$$\begin{aligned} \Phi (r)=\Phi _0e^{-A/r}r^B\left[ C+D\left( r-r_0\right) \right] ^{-B}, \end{aligned}$$

(61)

where $\Phi _0,A,B,C,D$ are constants. By substituting the above expressions of metric tensor components and of the scalar field into Eq. (34), it follows that the equation is identically satisfied if the unknown coefficients satisfy the conditions

$$\begin{aligned}{} & {} C=K_{11},D=K_{12}, \end{aligned}$$

(62)

$$\begin{aligned}{} & {} A\left( C-r_0D\right) +\alpha C_1\left( K_{21}-r_0K_{22}\right) =0 \end{aligned}$$

(63)

and

$$\begin{aligned} AD+B\left( C-r_0D\right) +C_1K_{22}\alpha =0, \end{aligned}$$

(64)

respectively, giving the relations

$$\begin{aligned} K_{21}= & {} \frac{B r_0 \left( K_{12} r_0-K_{11}\right) -A K_{11}}{\alpha C_1}, \end{aligned}$$

(65)

$$\begin{aligned} K_{22}= & {} -\frac{A K_{12}+B\left( K_{11}- K_{12} r_0\right) }{\alpha C_1}. \end{aligned}$$

(66)

By substituting the above representations of the metric tensor and the scalar field into the field equations Eq. (36) and (37), one obtains a system of ordinary nonlinear algebraic equations for the determination of the parameters of the solution. The resulting equations can be generally solved only approximately, due to their extremely complicated character. The series solutions can be extended to any order of approximation. Hence, at least in principle, one can obtain a complete power series solution of the Weyl geometric field equations near the black hole horizon.

An alternative approximate solution can be obtained by neglecting the nonlinear term $\Phi ^{\prime 2}/\Phi ^2$ in the scalar field equation (48) and in the gravitational field equations (36) and (37). In this case, we obtain a system of algebraic nonlinear equations that can be solved to recursively obtain the values of the coefficients in the Taylor series expansions of the physical and geometrical quantities. Thus, a series solution of the static spherically symmetric field equations of Weyl geometric gravity in the presence of the radial component of the Weyl vector can be constructed in both its linearized and nonlinear versions.

3.4.2 Models with a temporal component of the Weyl vector

a. Asymptotic limits. In the case in which the Weyl vector has only a nonzero temporal component $\omega _0$, the scalar field is a constant, $\Phi =\Phi _0=\textrm{constant}$. In the asymptotic limit of the flat spacetimes, with $\lambda =0$ and $\nu '+\lambda '=0$, Eq. (40) becomes

$$\begin{aligned} \omega _0^{\prime \prime }+\frac{2}{r}\omega _0^{\prime }+\frac{\alpha ^2\Phi _0}{4\xi ^2}\omega _0=0, \end{aligned}$$

(67)

with the general solution given by

$$\begin{aligned} \omega _0=\frac{c_1\cos ( \sqrt{A} r)}{r}+\frac{c_2 \sin (\sqrt{A} r)}{r}, \end{aligned}$$

(68)

where we have denoted $A=\alpha ^2\Phi _0/4\xi ^2$, and $c_1$ and $c_2$ are constants of integration. Thus, we obtain

$$\begin{aligned} \lim _{r\rightarrow \infty }\omega _0=0, \lim _{r\rightarrow \infty }\omega _0^{\prime }=0. \end{aligned}$$

(69)

By assuming that, asymptotically, the metric is de Sitter, with $e^{\lambda }=\left( 1-r^2/\mu ^2\right) ^{-1}$, Eq. (40) takes the form

$$\begin{aligned} \omega _0^{\prime \prime }+\frac{2}{r}\omega _0^{\prime }+\frac{A}{1-r^2/\mu ^2}\;\omega _0=0, \end{aligned}$$

(70)

with the general solution given by

$$\begin{aligned} \omega _{0}(r)= & {} c_{2}\,_{2}F_{1}\left( \delta _{-},\delta _{+};\frac{3}{2}; \frac{r^{2}}{\mu ^{2}}\right) \nonumber \\{} & {} -\frac{ic_{1}\mu \,_{2}F_{1}\left( -\delta _{+},-\delta _{-};\frac{1}{2}; \frac{r^{2}}{\mu ^{2}}\right) }{r}, \end{aligned}$$

(71)

where $c_{1}$ and $c_{2}$ are arbitrary integration constants, $_{2}F_{1}\left( a,b;c;z\right)$ is the hypergeometric function, and we have denoted

$$\begin{aligned} \delta _{-}=\frac{1}{4}\left( 1-\sqrt{1+4A\mu ^{2}}\right) ,\delta _{+}= \frac{1}{4}\left( 1+\sqrt{1+4A\mu ^{2}}\right) .\nonumber \\ \end{aligned}$$

(72)

Since $\omega _{0}$ is a real quantity, we need to take $c_{1}=0$ in the solution, thus obtaining

$$\begin{aligned} \omega _{0}(r)=c_{2}\,_{2}F_{1}\left( \delta _{-},\delta _{+};\frac{3}{2}; \frac{r^{2}}{\mu ^{2}}\right) . \end{aligned}$$

(73)

Here, the hypergeometric function $_{2}F_{1}(a,b;c;z)$, representing the regular solution of the hypergeometric differential equation, is defined for $|z|<1$ by a power series of the form $_{2}F_{1}(a,b;c;z)=\sum _{k=0}^{\infty }{\left[ (a)_k(b)_k/(c)_k\right] z^k/k!}$, where $(q)_k$ is the (rising) Pochhammer symbol [117]. For c not a negative integer, the function $_{2}F_{1}(a,b;c;z)$ converges for all of $|z|<1$ and, if $\textrm{Re}(c-a-b)>0$, also on the unit circle $|z|=1$. Different solutions of the hypergeometric differential equations can also be derived for other values of z, not restricted to the range $|z|<1$, and these solutions are valid in different regions of the complex plane.

Thus, by taking into account that the radius of convergence of the hypergeometric function in the expression (73) of $\omega _0(r)$ is $r^{2}<\mu ^{2}$, we find

$$\begin{aligned} \lim _{r\rightarrow \mu }\omega _{0}(r)=-\frac{2\sqrt{\pi }c_{2}}{A\mu ^{2}\Gamma \left( \delta _{-}\right) \Gamma \left( \delta _{+}\right) }, \end{aligned}$$

(74)

where $\Gamma (z)=\int _{0}^{\infty }{t^{z-1}e^{-t}\mathrm{{d}}t}$ is the Euler gamma function. The value $r=\mu$ of the radial coordinate defines a cosmological horizon for the present model.

For the limit at infinity of the derivative of $\omega _0$, we obtain

$$\begin{aligned} \lim _{r\rightarrow \mu }\omega _{0}^{\prime }(r)=-\frac{1}{3}Ac_{2}\mu \,_{2}F_{1}\left( 1+\delta _{-},1+\delta _{+};\frac{5}{2};1\right) . \end{aligned}$$

(75)

In the case of an asymptotically flat geometry, with the use of the conditions (69), the field equations (42) as evaluated at infinity give $\Phi =0$. Hence, it follows that if one assumes the presence of a nonzero scalar field, the asymptotic limit of the metric cannot be flat.

3.4.3 Models with both temporal and radial components of the Weyl vector

a. Asymptotic limits. Finally, we consider the behavior at infinity of the static spherically symmetric Weyl geometric models in the presence of both temporal and radial components of the Weyl field. In the case of the asymptotically flat geometry, the coupled system of equations describing the behavior of the scalar field and the temporal component of the Weyl vector are given by

$$\begin{aligned} \omega _0^{\prime \prime }+\frac{2}{r}\omega _0^{\prime }+\frac{\alpha ^2}{4\xi ^2}\Phi \omega _0=0 \end{aligned}$$

(76)

and

$$\begin{aligned} \frac{\Phi ^{\prime \prime }}{\Phi }-\frac{3}{2}\frac{\Phi ^{\prime 2}}{\Phi ^2}+\frac{3}{2}\alpha ^2\omega _0^2=0, \end{aligned}$$

(77)

respectively. By neglecting the nonlinear term $\Phi \omega _0$ in Eq. (76), we find

$$\begin{aligned} \omega _0^{\prime }(r)=\frac{c_4}{r^2},\omega _0(r)=c_3-\frac{c_4}{r}. \end{aligned}$$

(78)

Hence,

$$\begin{aligned} \lim _{r\rightarrow \infty }\omega _0^{\prime }=0, \lim _{r\rightarrow \infty }\omega _0=c_3. \end{aligned}$$

(79)

Substitution on Eq. (77) gives for $\Phi (r)$ the equation

$$\begin{aligned} \frac{\Phi ^{\prime \prime }}{\Phi }-\frac{3}{2}\frac{\Phi ^{\prime 2}}{\Phi ^2}+\frac{3}{2}\alpha ^2\left( c_3-\frac{c_4}{r}\right) ^2=0. \end{aligned}$$

(80)

By neglecting the nonlinear term $\Phi ^{\prime 2}/\Phi ^2$, and in the limit of large r, Eq. (80) becomes

$$\begin{aligned} \Phi ^{\prime \prime }+\frac{3}{2}\alpha ^2c_3^2\Phi =0, \end{aligned}$$

(81)

with the general solution given by

$$\begin{aligned} \Phi (r)=C_5\cos \left( \sqrt{\frac{3}{2}}\alpha c_3r+C_6\right) , \end{aligned}$$

(82)

where $C_5$ and $C_6$ are arbitrary integration constants. Hence, at least in the approximation considered, the scalar field has an oscillatory behavior at infinity. For the Weyl vector component, we obtain

$$\begin{aligned} \omega _1(r)\approx -\frac{3}{2}\alpha c_3\tan \left( \sqrt{\frac{3}{2}}\alpha c_3r+C_6\right) . \end{aligned}$$

(83)

Thus, the component of the Weyl vector has an oscillatory behavior at infinity.

In the case of an asymptotic de Sitter geometry, the scalar field equation (48) can be approximated as

$$\begin{aligned} \Phi ^{\prime \prime }+\frac{3}{2}\frac{\alpha ^2c_3^2}{\left( 1-r^2/\mu ^2\right) ^2}\Phi =0, \end{aligned}$$

(84)

and it has the general solution

$$\begin{aligned}{} & {} \Phi (r)= c_5 (r-\mu )^{\frac{1}{4} \left( \sqrt{4-6 \alpha ^2 c_3^2 \mu ^2}+2\right) }\nonumber \\{} & {} \quad (r+\mu )^{\frac{1}{4} \left( 2-\sqrt{4-6 \alpha ^2 c_3^2 \mu ^2}\right) }\nonumber \\{} & {} \quad -\frac{c_6 (r-\mu )^{\frac{1}{4} \left( 2-\sqrt{4-6 \alpha ^2 c_3^2 \mu ^2}\right) } (r+\mu )^{\frac{1}{4} \left( \sqrt{4-6 \alpha ^2 c_3^2 \mu ^2}+2\right) }}{\mu \sqrt{4-6 \alpha ^2 c_3^2 \mu ^2}},\nonumber \\ \end{aligned}$$

(85)

where $c_5$ and $c_6$ are integration constants. In the limit $r\rightarrow \zeta$, the scalar field tends to zero. Similarly, in the limit of large r, the derivative of the scalar field tends to zero. However, we would like to point out that the above results are approximate, but even so they indicate that both arbitrary nonzero and very small (zero) numerical values are possible for the scalar and Weyl vector fields at infinity.

4 Black hole solutions in Weyl geometric conformal gravity

In order to simplify the mathematical formalism, we now introduce the following representation for $e^{-\lambda }$

$$\begin{aligned}&e^{-\lambda } = 1-\frac{2GM(r)}{c^{2}r}=1-\frac{2nGM_{\odot }m(r)}{c^{2}r}, \end{aligned}$$

(86)

where

$$\begin{aligned} m(r) = \frac{M(r)}{n M_{\odot }}, \end{aligned}$$

(87)

and $M_{\odot }$ denotes the mass of the sun. M(r) and m(r) are effective masses incorporating the effects of both the scalar field and the Weyl vector, while n is a positive integer.

We now define a group of dimensionless quantities $(\eta ,m,\nu ,\psi ,\zeta ,\Omega _0,\Theta _0)$, given by

$$\begin{aligned}&r = \frac{2r_g}{\eta },\;\; \psi = r_g^2\Phi , \;\; \zeta = r_g^3 \frac{\mathrm{{d}}\Phi }{\mathrm{{d}}r},\;\; \Omega _0 = \alpha r_g \omega _0, \nonumber \\&\Theta _0 = \xi r_g^2 \frac{\mathrm{{d}}\omega _0}{\mathrm{{d}}r},\;\; e^{-\lambda } = 1- m(\eta )\eta , \end{aligned}$$

(88)

where we have denoted $r_g=nGM_{\odot }/c^{2}$. Hence, for the derivative with respect to $\eta$, we obtain

$$\begin{aligned} \frac{\mathrm{{d}}}{\mathrm{{d}}r} = -\frac{\eta ^2}{2r_g} \frac{\mathrm{{d}}}{\mathrm{{d}}\eta } \end{aligned}$$

(89)

and

$$\begin{aligned}&\frac{\mathrm{{d}}\lambda }{\mathrm{{d}}r} = -\frac{\eta ^2}{2r_g} \frac{\mathrm{{d}}\lambda }{\mathrm{{d}}\eta } = - \frac{\eta ^2}{2r_g}\frac{\frac{\mathrm{{d}}m(\eta )}{\mathrm{{d}}\eta }\eta + m(\eta )}{1-m(\eta )\eta }, \end{aligned}$$

(90)

respectively. We also define a new constant $\gamma$ as $\gamma = \alpha /\xi$.

4.1 Black hole solutions with radial component of the Weyl vector field

We will first consider black hole solutions in which the Weyl vector has only a radial component $\omega _1$. In this case it is possible to find an exact solution of the field equations. Numerical black hole solutions are also obtained for a specific set of initial conditions of the Weyl vector and the scalar field at infinity.

4.1.1 Exact black hole solutions

In the following, we will first look for exact black hole solutions satisfying the condition

$$\begin{aligned} \nu (r)+\lambda (r) =0, \forall r>0. \end{aligned}$$

(91)

Then, Eq. (39) immediately gives for $\Phi$ the equation

$$\begin{aligned} \Phi ''=\frac{3}{2}\frac{\Phi ^{\prime 2}}{\Phi }, \end{aligned}$$

(92)

with the general solution given by

$$\begin{aligned} \Phi (r)=\frac{C_1}{\left( r+C_2\right) ^2}, \end{aligned}$$

(93)

where $C_1$ and $C_2$ are arbitrary constants of integration. Equation (36) can be reformulated as

$$\begin{aligned}{} & {} -1+\frac{\mathrm{{d}}}{\mathrm{{d}}r}\left( re^{-\lambda }\right) -\frac{1}{4}r^2\Phi +2re^{-\lambda }\frac{\Phi '}{\Phi }-\frac{3}{4}r^2e^{-\lambda }\frac{\Phi ^{\prime 2}}{\Phi ^2}\nonumber \\{} & {} -\frac{r^2\lambda 'e^{-\lambda }}{2}\frac{\Phi '}{\Phi }+r^2e^{-\lambda }\frac{\Phi ''}{\Phi }=0. \end{aligned}$$

(94)

By taking into account the identity

$$\begin{aligned} -\frac{r^2\lambda 'e^{-\lambda }}{2}\frac{\Phi '}{\Phi }=-\frac{re^{-\lambda }}{2}\frac{\Phi '}{\Phi }+\frac{r}{2}\frac{d}{dr}\left( re^{-\lambda }\right) \frac{\Phi '}{\Phi }, \end{aligned}$$

(95)

Eq. (94) becomes

$$\begin{aligned}{} & {} \left( 1+\frac{r}{2}\frac{\Phi '}{\Phi }\right) \frac{d}{dr}\left( re^{-\lambda }\right) + \left( \frac{3}{2}\frac{\Phi '}{\Phi }-\frac{3}{4}r\frac{\Phi ^{\prime 2}}{\Phi ^2}+r\frac{\Phi ''}{\Phi }\right) \left( re^{-\lambda }\right) \nonumber \\{} & {} \quad -\frac{1}{4}r^2\Phi -1=0. \end{aligned}$$

(96)

With the use of the expression of $\Phi (r)$ as given by Eq. (93), we obtain

$$\begin{aligned}{} & {} \left( 1-\frac{r}{C_2+r}\right) u'(r)+3\left( \frac{ r}{(C_2+r)^2}-\frac{1}{C_2+r}\right) u(r)\nonumber \\{} & {} \quad -\frac{C_1 r^2}{4 (C_2+r)^2}-1=0, \end{aligned}$$

(97)

where we have denoted $u=re^{-\lambda }$. Equation (97) has the exact general solution

$$\begin{aligned} re^{-\lambda }= & {} \frac{r^2 \left( 12 C_3 C_2^2-C_1-4\right) }{4 C_2}+r \left( 3 C_3 C_2^2-\frac{C_1}{4}-2\right) \nonumber \\{} & {} +C_3 C_2^3-\frac{1}{12} (C_1+12) C_2+C_3 r^3, \end{aligned}$$

(98)

where $C_3$ is an integration constant. Depending on the values of the integration constants $C_1$, $C_2$ and $C_3$, the metric (58) can take two distinct forms. If we assume the condition

$$\begin{aligned} 3 C_3C_2^2-\frac{C_1}{4}-2=1 \end{aligned}$$

(99)

or, equivalently, $C_3C_2^3-C_1C_2/12=C_2$, the metric (98) becomes

$$\begin{aligned} e^{-\lambda }=e^{\nu }=1+\frac{2}{C_2}r+C_3r^2, \end{aligned}$$

(100)

and it represents a generalization of the static cosmological de Sitter solution. The solution is not asymptotically flat. If one imposes the condition

$$\begin{aligned} 3 C_3C_2^2-\frac{C_1}{4}-2=2\beta \end{aligned}$$

(101)

where $\beta$ is a constant, the metric tensor components (98) become

$$\begin{aligned} e^{-\lambda }=e^{\nu }=2\beta +\frac{1+2\beta }{C_2}r-\frac{C_2\left( 1-2\beta \right) }{3}\frac{1}{r}+C_3r^2. \end{aligned}$$

(102)

By assuming, in analogy with the Schwarzschild metric, that $C_2\left( 1-2\beta \right) /3=r_g$, the metric tensor components (98) take the form

$$\begin{aligned} e^{-\lambda }=e^\nu =2\beta +\frac{1-4\beta ^2}{3}\frac{r}{r_g}-\frac{r_g}{r}+C_3r^2. \end{aligned}$$

(103)

If we denote $1-2\beta =\delta$, the metric (103) can be written as

$$\begin{aligned} e^{-\lambda }=e^\nu =1-\delta +\frac{\delta (2-\delta )}{3}\frac{r}{r_g}-\frac{r_g}{r}+C_3r^2. \end{aligned}$$

(104)

In Fig.  1 we have plotted the behavior of the metric tensor $e^\nu$ as a function of $r/r_g$ for different values of the constants $c_3\equiv C_3r_g^2$ and $\delta$, respectively. In each of these cases a singularity in the metric tensor component does appear, indicating the formation of an event horizon and the presence of a black hole.

Fig. 1

Variation of the metric tensor component $e^{\nu }$ for the exact solution of the Weyl geometric gravity in the presence of the radial component of the Weyl vector only, for $c_3=0.01$, and different values of $\delta$ (left panel), and for $\delta =1$, and different values of $c_3$ (right panel)

It is interesting to note that very similar solutions of the vacuum field equations have been obtained in the framework of other modified gravity theories. In particular, in conformal Weyl gravity [86,87,88,89,90,91], for a static spherically symmetric metric having the standard form $\mathrm{{d}}s^2=-B(r)\mathrm{{d}}t^2+B^{-1}(r)\mathrm{{d}}r^2+r^2\mathrm{{d}}\Omega$, the theory admits vacuum solutions of the form [86]

$$\begin{aligned} B(r)=1-3\beta \gamma -\frac{\beta \left( 2-3\beta \gamma \right) }{r}+\gamma r+kr^2, \end{aligned}$$

(105)

where $\beta$, $\gamma$ and k are arbitrary integration constants [86]. With the use of the metric (105) in [86], it was suggested that Weyl conformal gravity can explain the flat rotation curves of galaxies without invoking the presence of dark matter. Another interesting modified gravity theory, the dRGT massive gravity theory [118,119,120], with action

$$\begin{aligned} S = \int \mathrm{{d}}^4 x \sqrt{-g} \frac{M_{Pl}^2}{2} \left[ R + m_g^2 {\mathcal {U}}(g,f) \right] + S_m \,, \end{aligned}$$

(106)

where $M_{Pl}$ is the reduced Planck mass, R is the Ricci scalar, $m_g$ is the graviton mass, and ${\mathcal {U}}$ is the self-interacting potential of the gravitons, also admits a static spherically symmetric solution of the form (see [120], and references therein)

$$\begin{aligned} B(r) = 1 - \frac{2 G m(r)}{r} - \frac{\Lambda r^2}{3} + \gamma r + \zeta , \end{aligned}$$

(107)

where m(r) is the mass within radius r, and $\Lambda \equiv - 3 m^2_g (1 + \alpha + \beta )$, $\gamma \equiv - m^2_g C (1 + 2 \alpha + 3 \beta )$, and $\zeta \equiv m_g^2 C^2 (\alpha + 3 \beta )$.

Here, C, $\gamma$, and $\zeta$ are constants depending on the graviton mass $m_g$, and $\Lambda$ corresponds to the cosmological constant, while $\alpha =3\alpha _3+1$ and $\beta =4\alpha _4-(1-\alpha )/3$ are coefficients related to the decomposition of the self-interacting potential of the graviton as ${\mathcal {U}}={\mathcal {U}}_2+\alpha _3{\mathcal {U}}_3+\alpha _4{\mathcal {U}}_4$, with ${\mathcal {U}}_2 \equiv [{\mathcal {K}}]^2 - [{\mathcal {K}}^2]$, ${\mathcal {U}}_3 \equiv [{\mathcal {K}}]^3 - 3 [{\mathcal {K}}][{\mathcal {K}}^2] + 2 [{\mathcal {K}}^3]$, and ${\mathcal {U}}_4 \equiv [{\mathcal {K}}]^4 - 6 [{\mathcal {K}}]^2 [{\mathcal {K}}^2] + 3[{\mathcal {K}}^2]^2 + 8 [{\mathcal {K}}][{\mathcal {K}}^3] - 6 [{\mathcal {K}}^4]$, where ${\mathcal {K}}^{\mu }_{\nu } \equiv \delta ^{\mu }_{\nu } - \sqrt{g^{\mu \lambda } \partial _{\lambda } \varphi ^a \partial _{\nu } \varphi ^b f_{ab}}$, and we have denoted $[{\mathcal {K}}]={\mathcal {K}}_{\mu }^{\mu }$. Here, $g_{\mu \nu }$ is the physical metric, $f_{\mu \nu }$ is a reference (fiducial) metric, and $\varphi ^a$ are the Stueckelberg fields [120]. For $m_g = 0$, we recover the standard Schwarzschild solution of general relativity.

The existence of mathematically identical vacuum solutions in several distinct gravity theories raises the question of the possible universality of metrics of the form (104), representing the most general extension of the Schwarzschild metric.

4.1.2 Numerical black hole solutions

We will now proceed to obtaining numerical black hole solutions in Weyl geometric gravity in the presence of a Weyl vector field having only a radial component. In a static and spherically symmetric configuration, Eqs. (35)–(37) become

$$\begin{aligned}{} & {} \frac{\mathrm{{d}}m}{\mathrm{{d}}\eta } = -\frac{7(1-\eta m)\zeta ^3+\eta (5-4\eta m)\zeta ^2\psi +\zeta (3\eta ^3 m-4\eta ^2-3\psi )\psi ^2-\eta \psi ^4}{\eta ^3\psi (\zeta +\eta \psi )^2},\end{aligned}$$

(108)

$$\begin{aligned}{} & {} \frac{\mathrm{{d}}\nu }{\mathrm{{d}}\eta } = \frac{\zeta (1-\eta m)(3\zeta +4\eta \psi )-(\eta ^3 m+\psi )\psi ^2}{\eta ^2\psi (1-\eta m)(\zeta +\eta \psi )},\end{aligned}$$

(109)

$$\begin{aligned}{} & {} \frac{\mathrm{{d}}\psi }{\mathrm{{d}}\eta } = -\frac{2\zeta }{\eta ^2},\end{aligned}$$

(110)

$$\begin{aligned}{} & {} \frac{\mathrm{{d}}\zeta }{\mathrm{{d}}\eta } = -\frac{\zeta ^2(1-\eta m)(5\zeta +2\eta \psi )+\zeta (3\eta ^3 m-4\eta ^2-\psi )\psi ^2}{\eta ^2\psi (1-\eta m)(\zeta +\eta \psi )}, \end{aligned}$$

(111)

Note that the above system of ordinary, strongly nonlinear differential equations is independent of the coupling constants $\alpha$ and $\xi$. We integrate the system from infinity, corresponding to $\eta =\lim _{r\rightarrow \infty }1/r=0$, up to the event horizon of the black hole. We find a numerical black hole solution by detecting a singularity in metric tensor components $e^{\nu }$ and $e^{\lambda }$. The singularities in the $e^{\nu }$ and $e^{\lambda }$ metric tensor components are indicated by the fulfillment of the conditions

$$\begin{aligned} \left. e^{\nu (\eta )}\right| _{\eta =\eta _{\text {hor}}}=0,\;\left. e^{-\lambda (\eta )}\right| _{\eta =\eta _{\text {hor}}}=1-\left. \left[ m(\eta )\eta \right] \right| _{\eta =\eta _{\text {hor}}}=0,\nonumber \\ \end{aligned}$$

(112)

where $\eta _{\text {hor}}$ is the horizon of the black hole. For a Schwarzschild black hole, the position of the event horizon corresponds to $\eta _{\text {hor}}=1$.

a. The initial conditions. In order to numerically integrate the gravitational field equations in the variable $\eta$, we need to fix the initial conditions at $\eta =0$, corresponding to the asymptotic values of the scalar and vector fields, and of the metric. As we have seen in the discussion of the asymptotic limits of this model, at infinity, the values of the scalar field and its derivative can be taken either as zero, so that $\psi (0)=0$ and $\zeta (0)=0$, or as having some arbitrary numerical values. As for the metric tensor components, we assume that the metric either can be either asymptotically flat, corresponding to $\nu (0)=\lambda (0)=0$, or can have a de Sitter form, $\left. e^{-\lambda }\right| _{\eta =\eta _{0}}=\left. \left[ 1-\left( 4r_g^2/\mu ^2\right) \left( 1/\eta ^2\right) \right] \right| _{\eta =\eta _{0}}$, giving $\left. m(\eta )\right| _{\eta =\eta _{0}}=\left( 4r_g^2/\mu ^2\right) \left( 1/\eta ^3\right) _{\eta =\eta _0}=m_0$.

Hence, in the de Sitter case, we fix the initial values of $(m,\nu )$ at the physical infinity $\eta =\eta _0$ as $m\left( \eta _0\right) =m_0$ and $\nu =\left. \ln \left( 1-m(\eta )\eta \right) \right| _{\eta =\eta _0}=\ln \left( 1-m_0\eta _0\right)$.

Fig. 2

Variation of the metric tensor components $e^\nu$ (upper left panel) and $e^{-\lambda }$ (upper right panel), of $\psi$ (lower left panel), and of $\zeta$ (lower right panel) as a function of $\eta$, for a Weyl geometric black hole in the presence of the radial component of the Weyl vector only, for $\zeta (0)=1\times 10^{-35}$, and for different values of $\psi (0)$, presented in the legends of the figures

To obtain the black hole solutions, we vary the initial values of the dimensionless scalar field $\psi (0)$ and its derivative $\zeta (0)$, and, through the numerical integration of the field equations, we obtain the variations of the metric tensor components, the scalar field, and the Weyl vector field as functions of the dimensionless (inverse) radial coordinate $\eta$.

b. Numerical results. In the following, we consider numerical black hole solutions obtained by numerically integrating the gravitational field equations by assuming the presence of the spatial component of the Weyl vector only. We consider two types of models, obtained by fixing the initial value of $\zeta$ at infinity and varying the values of the scalar field, and by fixing the scalar field at infinity and varying its derivative $\zeta (0)$.

Models with fixed $\zeta (0)$ and varying $\psi (0)$ . In Fig. 2, we present the results of the numerical integration of the field equations in the presence of a radial component of the Weyl vector only, obtained by fixing the initial value of the derivative of the scalar field $\zeta$ at infinity as $\zeta (0)=1\times 10^{-35}$, and by varying the initial values of the scalar field $\psi (0)$. We restrict our analysis to the case of the flat asymptotic conditions. As one can see from the upper panels of Fig. 2, the metric tensor components are decreasing functions of $\eta$, and they become singular for a finite value of the radial coordinate, indicating the formation of a black hole. The position of the event horizon is very sensitive to the small variations of the initial conditions. The scalar field is a decreasing function of $\eta$, taking only positive values, while $\zeta$ increases with increasing $\eta$.

Selected values of $\eta _{\textit{hor}}$ corresponding to different values of $\psi (0)$, and for a fixed $\zeta (0)$, are presented in Table 1. As one can see from Table 1, the position of the event horizon of the black hole is strongly dependent on the values of the scalar field at infinity, and large variations of its location are possible, ranging from half to twice of the Schwarzschild values. Thus, for example, for $\psi (0)=10^{-15}$, the physical radius of the event horizon, $r_\mathrm{{hor}}$, is located at $r_\mathrm{{hor}}\approx 0.54$, while $r_\mathrm{{hor}}\approx 2.10$ for $\psi (0)=6\times 10^{-15}$.

Table 1 Variation of the position of the event horizon $\eta _{\text {hor}}$ of the Weyl geometric black hole with radial component of the Weyl vector for $\zeta (0)=1\times 10^{-35}$, and different initial values of $\psi (0)$

Models with fixed $\psi (0)$ and varying $\zeta (0)$ . In Fig. 3 we present the results of the numerical integration of the field equations in the presence of a radial component of the Weyl vector only, obtained by fixing the initial value of the scalar field $\psi (0)$ at infinity as $\psi (0)=1\times 10^{-15}$, and by varying the initial value of its derivative $\zeta (0)$.

Fig. 3

Variation of the metric tensor components $e^\nu$ (upper left panel) and $e^{-\lambda }$ (upper right panel), of $\psi$ (lower left panel), and $\zeta$ (lower right panel) as a function of $\eta$, for a Weyl geometric black hole in the presence of the radial component of the Weyl vector only, for $\psi (0)=1\times 10^{-15}$, and for different values of $\zeta (0)$, presented in the legends of the figures

As one can see from the upper panels of Fig. 3, a singularity does appear in the metric tensor components, corresponding to the formation of a black hole. The position of the event horizon is strongly dependent on the derivatives of the scalar field at the (physical) infinity. The scalar field is a decreasing function of $\eta$, while its derivative monotonically increases towards the event horizon.

The variation of $\eta _{\textit{hor}}$ with respect to the changes in the values of the derivative of the scalar field at infinity $\zeta (0)$, for a fixed $\psi$, are presented in Table 2.

Table 2 Variation of the position of the event horizon $\eta _{\text {hor}}$ of the Weyl geometric black hole with radial component of the Weyl vector for $\psi (0)=1\times 10^{-15}$, and different initial values $\zeta (0)$

As one can see from Table 2, the position of the event horizon of the black hole decreases as a function of $\eta$ (increases as a function of r) with increasing initial values of the derivative of the scalar field $\zeta (0)$. Thus, for $\zeta (0)=0.8\times 10^{-35}$, the physical radius of the event horizon has the value $r_\mathrm{{hor}}\approx 0.45$, while $r_\mathrm{{hor}}\approx 0.63$ for $\zeta (0)=1.5\times 10^{-35}$. Hence, higher values of $\zeta (0)$ generate black holes having larger radii. However, we would like to point out that in these considered examples, the event horizons of the black holes are located at much smaller radii than the event horizons of their standard general relativistic counterparts.

Fig. 4

Variations of the metric tensor components $e^\nu$ (upper left panel), of the effective mass $\eta m$ (upper right panel), of $\Omega _0$ (lower left panel), and of $\Theta _0$ (lower right panel) as a function of $\eta$ in the presence of the temporal component of the Weyl vector only for $\psi (0)=1\times 10^{-15}$, $\Omega _0(0)=1\times 10^{-11}$, $\Theta _0(0)=1\times 10^{-18}$, and for different values of $\gamma$, presented in the legends of the figures

4.2 Black hole solutions with temporal component of the Weyl vector field

We now consider a second class of Weyl type geometric black holes, in which the Weyl vector has only a temporal component, with $\omega _{\mu }$ given by $\omega _{\mu }=\left( \omega _0, 0,0,0\right)$. In a static and spherically symmetric gravitational field configuration, and in the presence of a Weyl vector with a temporal component only, the gravitational field equations Eqs. (40)–(43) become

$$\begin{aligned} \frac{\mathrm{{d}}m}{\mathrm{{d}}\eta }= & {} \frac{\psi }{\eta ^4} +\frac{12(1-\eta m)\Theta _0^2-3\psi \Omega _0^2}{e^\nu \eta ^4 \psi }, \end{aligned}$$

(113)

$$\begin{aligned} \frac{\mathrm{{d}}\nu }{\mathrm{{d}}\eta }= & {} - \frac{e^\nu \psi (\psi +\eta ^3m)+12(1-\eta m)\Theta _0^2-3\psi \Omega _0^2}{e^\nu \eta ^3 \psi (1-\eta m)},\nonumber \\ \end{aligned}$$

(114)

$$\begin{aligned} \frac{\mathrm{{d}}\psi }{\mathrm{{d}}\eta }= & {} 0,\end{aligned}$$

(115)

$$\begin{aligned} \frac{\mathrm{{d}}\Omega _0}{\mathrm{{d}}\eta }= & {} - \frac{2\gamma \Theta _0}{\eta ^2},\end{aligned}$$

(116)

$$\begin{aligned} \frac{\mathrm{{d}}\Theta _0}{\mathrm{{d}}\eta }= & {} \frac{-\gamma \eta \psi \Omega _0+\Theta _0(4\eta ^2(1-\eta m)-6e^{-\nu }\Omega _0^2)}{2\eta ^3(1-\eta m)}, \end{aligned}$$

(117)

4.2.1 The initial conditions

In this case the model depends on six parameters, the initial conditions at infinity of m and $\nu$, describing the metric tensor components, as well as of the scalar field $\psi$, of the Weyl vector $\Omega _0$, and of its derivative $\Theta _0$. Moreover, a coupling constant $\gamma$ is also present in the model. In the following, we investigate numerically only the asymptotically flat case, thus assuming $\nu (0)=\lambda (0)=0$. As one can see from Eqs. (69), asymptotically, the Weyl vector temporal component and its derivative tend to zero, and hence we have $\Omega _0(0)=0$ and $\Theta _0(0)=0$, respectively. For the initial value of m we adopt the condition $m(0)=0$. The scalar field $\psi$ is a constant, and its value can be fixed arbitrarily. From the numerical point of view, we investigate the effects of the variation of the coupling constant and of the numerical value of the scalar field on the position of the event horizon.

4.2.2 Numerical results

We now present the numerical results of the numerical integration of the field equations in the presence of a temporal component of the Weyl vector only.

Fig. 5

Variations of the metric tensor components $e^\nu$ (upper left panel), and of the effective mass $\eta m$ (upper right panel), of $\Omega _0$ (lower left panel), and of $\Theta _0$ (lower right panel) as a function of $\eta$ of a Weyl black hole in the presence of the temporal component of the Weyl vector only, for $\gamma =100$, $\Omega _0 (0)=1\times 10^{-11}$, $\Theta _0(0)=1\times 10^{-20}$, and for different values of $\psi (0)$, presented in the legends of the figures

a. Varying the value of the coupling constant $\gamma$. As a first example of black hole solutions with a temporal Weyl vector, we consider models in which only the coupling parameter $\gamma$ varies, with all other parameters fixed. The results of the numerical integration of the field equations (113)–(117) are represented in Fig. 4.

As one can see from Fig. 4, the metric tensor components present a singular behavior, indicating the formation of the event horizon. There is a strong dependence on the numerical values of $\gamma$ of $\eta _{\text {hor}}$. The behavior and the numerical values of the effective mass are also strongly dependent on the coupling constant $\gamma$, indicating an increase in the effective mass when approaching the event horizon of the black holes. The variation of the Weyl vector temporal component $\Omega _0$ is represented in the lower right panel of Fig. 4. As a function of $\eta$, the Weyl vector field is a linearly decreasing function, while its derivative is an increasing function of the inverse of the radial variable. The effects of the variation of the coupling constant $\gamma$ on the behavior of the fields are significant.

The positions of the event horizon of the black holes are presented in Table 3 for different values of $\gamma$, with all other initial conditions fixed as in Fig. 4.

Table 3 The position of the event horizon of the Weyl geometric type black holes in the presence of the temporal component of the Weyl vector only for different values of $\gamma$, and with the other parameters fixed as in the plots in Fig. 4

There is a significant dependence of the position of the event horizon of the black hole on the numerical value of the coupling $\gamma$, with the position of the physical radius $r_\mathrm{{hor}}$ of the event horizon of the black hole decreasing with increasing $\gamma$. Thus, for $\gamma =100$, $r_\mathrm{{hor}}\approx 1.27$, while for $\gamma =500$, $r_{hor}\approx 0.95$. Hence, Weyl geometric black holes with radii larger than their Schwarzschild counterparts can be obtained for small values of the coupling constant $\gamma$.

b. Varying the value of the scalar field. We now consider black hole models in Weyl geometric type gravity in the presence of the temporal component of the Weyl vector obtained by varying the numerical values of the constant scalar field $\psi$, while keeping all the other numerical values at infinity of the physical and geometrical quantities fixed. The variations with respect to $\eta$ of the metric tensor components and of the effective mass are represented for this case in Fig. 5. As one can see from Fig. 5, for this choice of parameters, the position of the event horizon is strongly dependent on the values of the scalar field, which as a constant acts as a cosmological background. The Weyl vector $\Omega _0$ is linearly decreasing with respect to $\eta$ (linearly increasing as a function of r), while its derivative is a monotonically increasing function of $\eta$ (monotonically decreasing with respect to r).

The positions of the event horizon are presented for different initial values of the scalar field $\psi (0)$ in Table 4.

Table 4 Variation of the position of the event horizon of the Weyl geometric black holes in the presence of a temporal component of the Weyl vector for different initial values of the scalar field $\psi _0(0)$. The numerical values of the other quantities are fixed as in the plots in Fig. 5

Fig. 6

Variations of the metric tensor component $e^\nu$ (upper left panel), of the effective mass $\eta m$ (upper right panel), of $\Omega _0$ (lower left panel), and of $\Theta _0$ (lower right panel) as a function of $\eta$ for a Weyl black hole in the presence of a temporal component of the Weyl vector only, for $\gamma =100$, $\psi (0)=1\times 10^{-15}$, $\Theta _0(0)=1\times 10^{-20}$, and for different values of $\Omega _0(0)$

There is a significant dependence of the position of the event horizon on the values of the scalar field, with $\eta _{\text {hor}}$ decreasing with the increase in the numerical values of $\psi (0)$. The physical radius $r_\mathrm{{hor}}$ increases with the increase in $\psi _0(0)$, from $r_\mathrm{{hor}}\approx 0.33$ ($psi (0)=10^{-15}$) to $r_\mathrm{{hor}}\approx 1.33$ for $\psi (0)=4\times 10^{-15}$.

c. Varying the initial value $\Omega _0(0)$ of the temporal component of the Weyl vector. We now consider Weyl geometric type black holes obtained by varying only the initial value of $\Omega _0$, while fixing the values at infinity of the other physical and geometrical quantities. The variations of the metric tensor components and of the effective mass are presented in Fig. 6. As one can see from the figure, the presence of singularities in the metric tensor components indicates the formation of an event horizon and, therefore, the presence of a black hole. The position of the singularities is strongly dependent on the values at infinity of the Weyl vector. This dependence is also apparent in the behavior of the effective mass m, which increases with increasing $\eta$. The variations of the temporal component of the Weyl vector is represented in the lower right panel in Fig. 6. The Weyl field decreases linearly towards the event horizon, and its values significantly depend on the asymptotic value of the field. The derivative of the Weyl vector is an increasing function of $\eta$ (a decreasing function of r).

The explicit values of the positions of the event horizon are presented in Table 5. As one can see from Table 5, there is a very strong dependence of $\eta _\mathrm{{hor}}$ on $\Omega _0(0)$. The event horizon location $\eta _{\text {hor}}$ decreases with increasing $\Omega _0(0)$, thus indicating an important effect of the Weyl vector on the global properties of the black holes. On the other hand, the physical radius of the event horizon $r_\mathrm{{hor}}$ increases from $r_\mathrm{{hor}}\approx 0.33$ ($\Omega _0(0)=10^{-11}$), to $r_\mathrm{{hor}}\approx 1.27$, for $\Omega _0(0)=6\times 10^{-10}$.

Table 5 Variation of the position of the event horizon for a Weyl geometric black hole in the presence of the temporal component of the Weyl vector $\Omega _0$ with respect to $\Omega _0(0)$. The numerical values of the other quantities are fixed as in the plots in Fig. 6

4.3 Black hole solutions in the presence of both temporal and radial components of the Weyl vector

We now consider the general case in which both components of the Weyl vector are nonzero, and therefore $\omega _{\mu }=\left( \omega _0,\omega _1,0,0\right)$. Then, the field equations (44) and (45) can be reformulated in a dimensionless form as a first-order dynamical system given by

$$\begin{aligned} \frac{\mathrm{{d}}m}{\mathrm{{d}}\eta }= & {} \frac{1}{\eta ^4 \psi ^2 (\zeta + \eta \psi )} \Big (12e^{-\nu }\Theta _0^2\psi (1-\eta m)(2\zeta +\eta \psi )\nonumber \\{} & {} -\zeta ^2(1-\eta m)(4\zeta +5\eta \psi )\nonumber \\{} & {} +\zeta \psi ^2(\eta ^3m+2\psi )+\eta \psi ^4 \Big ), \end{aligned}$$

(118)

$$\begin{aligned} \frac{\mathrm{{d}}\nu }{\mathrm{{d}}\eta }= & {} \frac{1}{\eta ^2 \psi (1-\eta m) (\zeta + \eta \psi )} \nonumber \\{} & {} \times \Big ( \zeta (1-\eta m)(3\zeta +4\eta \psi )-(\eta ^3m+\psi )\psi ^2\nonumber \\{} & {} -3e^{-\nu }\psi \big (4(1-\eta m)\Theta _0^2+\psi \Omega _0^2\big )\Big ), \end{aligned}$$

(119)

$$\begin{aligned} \frac{\mathrm{{d}}\psi }{\mathrm{{d}}\eta }= & {} - \frac{ 2\zeta }{\eta ^2}, \end{aligned}$$

(120)

$$\begin{aligned} \frac{\mathrm{{d}}\zeta }{\mathrm{{d}}\eta }= & {} -\frac{\zeta }{\eta ^3\psi ^2} \bigg [ 2\zeta (\zeta +2\eta \psi )-12e^{-\nu }\psi \Theta _0^2\nonumber \\{} & {} -\frac{\psi ^2}{1-\eta m}\big (\eta ^2(2-\eta m)-1\big ) \bigg ], \end{aligned}$$

(121)

$$\begin{aligned} \frac{\mathrm{{d}}\Omega _0}{\mathrm{{d}}\eta }= & {} - \frac{2\gamma \Theta _0}{\eta ^2}, \end{aligned}$$

(122)

$$\begin{aligned} \frac{\mathrm{{d}}\Theta _0}{\mathrm{{d}}\eta }= & {} \frac{1}{2\eta ^3\psi ^2(1-\eta m)(\zeta + \eta \psi )} \nonumber \\{} & {} \times \Big ( 2\Theta _0(1-\eta m)(12e^{-\nu }\zeta \Theta _0^2\psi -2\zeta ^3-\eta \zeta ^2\psi \nonumber \\{} & {} +2\eta ^3\psi ^3)-6e^{-\nu }\eta \Theta _0\psi ^3\Omega _0^2\nonumber \\{} & {} +2\eta ^2\Theta _0\zeta \psi ^2(4-3\eta m)\nonumber \\{} & {} +\gamma \psi ^3\Omega _0(\zeta +\eta \psi )+2\Theta _0\zeta \psi ^3 \Big ), \end{aligned}$$

(123)

4.3.1 The initial conditions

For this model, the general solution of the system depends on six parameters $\left( \gamma , m(0),\psi (0),\zeta (0),\Omega _{0}(0),\Theta _{0}(0)\right)$, representing the numerical values of the coupling constant $\gamma$, and of the initial conditions at infinity. As we have already discussed, the scalar field and the $\omega _1$ component have an oscillatory behavior for large r, and hence at infinity they do not converge to a single value. On the other hand, the $\omega _0$ component of the field becomes an arbitrary integration constant at infinity. However, we will select only the set of initial conditions that is consistent with the previously analyzed particular cases, and which lead to astrophysical black holes similar to the general relativistic Schwarzschild ones, with similar values for the positions of the event horizon. These conditions again imply very small initial values of the scalar and Weyl vector fields and their derivatives.

4.3.2 Numerical results

We now consider the results of the numerical integration of the static, spherically symmetric field equations of Weyl geometric gravity, obtained by varying the numerical values of the coupling constant and the initial conditions at infinity.

a. Varying the coupling constant $\gamma$. We consider first classes of numerical black hole solutions obtained by varying the coupling constant $\gamma$ only, while keeping the initial conditions at infinity of the physical and geometrical quantities fixed. The variations of the metric tensor components and of the effective mass are represented for different values of $\gamma$ in Fig. 7. The formation of an event horizon is indicated by the presence of the singularities in the metric tensor components. For the adopted values of the coupling constant, there is a strong dependence of the position of $\eta _\mathrm{{hor}}$ on $\gamma$. On the other hand, the effective mass m of the Weyl black hole is a monotonically increasing function of $\eta$, also showing a stronger dependence on the Weyl couplings.

Fig. 7

Variation of the metric tensor component $e^\nu$ (upper left panel), of $\eta m$ (upper right panel), of $\Theta _0$ (middle left panel), of $\Omega _0$ (middle right panel), of $\psi$ (lower left panel), and of $\zeta$ (lower right panel) as a function of $\eta$ for Weyl geometric black holes in the presence of both temporal and radial components of the Weyl vector for $\psi (0)=1\times 10^{-15}$, $\zeta (0)= 1\times 10^{-28}$, $\Omega _0(0)=1\times 10^{-11}$, $\Theta _0(0)=1\times 10^{-18}$, and for different values of the coupling constant $\gamma$, presented in the legends of the figures

The positions of the event horizon of the Weyl black holes for different values of $\gamma$ and fixed initial conditions are presented in Table 6. The modification of the coupling constant on a large range of values has only a relatively weak effect on the position of the event horizon as compared to the Schwarzschild case. However, the increase in $\gamma$ leads to a decrease in the value of the physical radius of the event horizon, from $r_{hor}\approx 1.28$ ($\gamma =100$) to $r_{hor}\approx 0.95$, corresponding to $\gamma =500$.

b. Varying the initial values of the scalar field $\psi (0)$ and of the temporal component of the Weyl vector $\Omega _0(0))$. We now consider numerical black hole solutions obtained by varying the initial values of the scalar field and the temporal component of the Weyl vector, with all other quantities fixed. The variation of the position of the event horizon is presented in Fig. 8.

The positions of the event horizon of the Weyl black holes for different values of $\psi (0)$ and fixed initial conditions and values of the coupling constant $\gamma$ are presented in Table 7. As one can see from the table, very small variations of the initial values of the scalar field can induce very significant changes in the position of the event horizon. Interestingly enough, increasing the value of the scalar field at infinity leads to a decrease in the inverse of the radius of the horizon of the black hole, leading to smaller values of $\eta _\mathrm{{hor}}$ as compared to the Schwarzschild case. The physical radius of the event horizon increases with increasing $\psi (0)$, from $r_\mathrm{{hor}}\approx 0.33$, corresponding to $\psi (0)=10^{_15}$, to $r_\mathrm{{hor}}\approx 1.33$ for $\psi (0)=4\times 10^{-15}$.

Table 6 Location of the event horizon of the Weyl geometric black holes in the presence of both temporal and radial components of the Weyl vector, for $\psi (0)=1\times 10^{-15}$, $\zeta (0)= 1\times 10^{-28}$, $\Omega _0(0)=1\times 10^{-11}$, and $\Theta _0(0)=1\times 10^{-18}$, respectively, and for different values of the coupling constant $\gamma$

Fig. 8

Variation of the position of the event horizon of the Weyl geometric black hole in the presence of both temporal and radial components of the Weyl vector as a function of the initial values of the scalar field $\psi _0(0)$ and of the temporal component of the Weyl vector $\Omega _0(0)$. To obtain the figure, we have assumed $\gamma =10^2$, $\Theta _{0}(0)=10^{-20}$, and $\zeta (0)=10^{-28}$, respectively

Table 7 Location of the event horizon of the Weyl geometric black holes for $\gamma =100$, $\zeta (0)= 1\times 10^{-28}$, $\Omega _0(0)=1\times 10^{-11}$, and $\Theta _0(0)=1\times 10^{-20}$, respectively, and for different values of $\psi (0)$

The positions of the event horizon of the Weyl black holes for different values of $\Omega _0(0)$ and fixed initial conditions are presented for $\gamma =100$ in Table 8. Increasing the values of $\Omega _0(0)$ leads again to a decrease in the inverse of the radius of the Weyl geometric black hole, and to an increase in the physical radius of the event horizon $r_\mathrm{{hor}}=1/\eta _\mathrm{{hor}}$. While for $\Omega _0(0)=10^{-11}$ the event horizon of the black hole is located at $\eta _\mathrm{{hor}}\approx 3$, $r_\mathrm{{hor}}=1/\eta _\mathrm{{hor}}\approx 0.33$, for $\Omega _0(0)=6\times 10^{-10}$ the position of the event horizon of the black hole has a value of $\eta _\mathrm{{hor}}\approx 0.8$, leading to a physical radius of the event horizon of the Weyl black hole larger than the value of the horizon of the Schwarzschild black hole, $r_\mathrm{{hor}}=1/\eta _\mathrm{{hor}}\approx 1.25$.

Table 8 Location of the event horizon of the Weyl geometric black holes in the presence of both temporal and radial components of the Weyl vector for $\gamma =100$, $\zeta (0)= 1\times 10^{-28}$, $\psi (0)=1\times 10^{-15}$, and $\Theta _0(0)=1\times 10^{-20}$, respectively, and for different values of $\Omega _0(0)$

c. Varying the initial values of the derivatives of the scalar field $\zeta (0)$ , and of the temporal component of the Weyl vector, $\Theta _0 (0)$. We now consider models obtained by varying the initial values of the derivatives of the scalar field, $\zeta (0)$, and of the derivative of the temporal component of the Weyl vector $\Theta _0(0)$ only, while keeping the initial conditions at infinity of the other physical and geometrical quantities fixed. The variations of the positions of the position of the event horizon are represented in Fig. 9. Our results indicate the formation of an event horizon due to the presence of the singularities in the metric tensor components.

Fig. 9

Variation of the position of the event horizon of the Weyl geometric black hole in the presence of both temporal and radial components of the Weyl vector as a function of the initial values of the derivatives of the scalar field $\psi _0(0)$ and of the temporal component of the Weyl vector $\Theta _0(0)$. To obtain the figure we have assumed $\gamma =10^2$, $\psi (0)=10^{-15}$, and $\Omega _0(0)=10^{-11}$, respectively

The positions of the event horizon of the Weyl black holes for different $\zeta (0)$ values and fixed initial conditions are presented in Table 9.

Table 9 Location of the event horizon of the Weyl geometric black holes in the presence of both temporal and radial components of the Weyl vector, for $\gamma =100$, $\Omega _0(0)= 1\times 10^{-11}$, $\psi (0)=1\times 10^{-15}$, and $\Theta _0(0)=1\times 10^{-20}$, and for different values of $\zeta (0)$

With the increase in $\zeta (0))$, the position of the event horizon increases in the physical radial distance variable r from $r_\mathrm{{hor}}\approx 0.34$ to $r_\mathrm{{hor}}\approx 1.375$, indicating the possibility for creation of both smaller and larger black holes as compared to their general relativistic counterparts.

5 Thermodynamics of the Weyl geometric black holes

In our analysis of the vacuum field equations of the linear representation of the quadratic Weyl geometric gravity, we have assumed that the metric tensor components $e^{\nu }$ and $e^{\lambda }$, and consequently the effective mass function m, all depend only on the radial coordinate r. Therefore, the geometry of the spacetime is static, and moreover, a timelike Killing vector $t^{\mu }$ always does exist [121, 122].

a. Surface gravity of the Weyl black holes. For a static black hole that has a Killing horizon, the surface gravity ${\tilde{\kappa }}$ is generally defined according to [121, 122]

$$\begin{aligned} t^{\mu }\nabla _{\mu }t^{\nu }=t^{\nu }{\tilde{\kappa }}. \end{aligned}$$

(124)

By adopting a static, spherically symmetric black hole geometry given by

$$\begin{aligned} \mathrm{{d}}s^2=-{\tilde{\sigma }} ^2 (r)f(r)c^2\mathrm{{d}}t^2+\frac{\mathrm{{d}}r^2}{f(r)}+r^2\mathrm{{d}}\Omega ^2, \end{aligned}$$

(125)

where ${\tilde{\sigma }}$ and f are functions of the radial coordinate only, and after suitably normalizing the Killing vector $t^{\mu }$ as $t^{\mu }=\left( 1/{\tilde{\sigma }}_{\infty },0,0,0\right)$, the surface gravity of the black hole is given by [122]

$$\begin{aligned} {\tilde{\kappa }}=\left( \frac{{\tilde{\sigma }} _\mathrm{{hor}}}{{\tilde{\sigma }} _{\infty }}\right) \frac{c^4}{4GM_{hor}}\left. \left[ 1-\frac{2GM'(r)}{c^2}\right] \right| _\mathrm{{hor}}, \end{aligned}$$

(126)

where the subscript hor specifies that the calculation of all physical quantities must be done on the outer apparent horizon of the black hole. If ${\tilde{\sigma }} \equiv 1$, and $M=\textrm{constant}$, then the expression of the surface gravity of a Schwarzschild black hole, ${\tilde{\kappa }}=c^4/4GM_{\text {hor}}$ [121], is re-obtained.

5.1 Hawking temperature of the quadratic Weyl geometric black holes.

The Hawking temperature $T_\mathrm{{BH}}$ of the black hole is defined according to [121, 122]

$$\begin{aligned} T_\mathrm{{BH}}=\frac{\hbar }{2\pi ck_B} {\tilde{\kappa }}, \end{aligned}$$

(127)

where by $k_B$ we have denoted Boltzmann’s constant. In the system of dimensionless variables defined in Eq. (88), the temperature of the black hole is obtained as

$$\begin{aligned} T_\mathrm{{BH}}=T_\mathrm{{H}}\frac{1}{m\left( \eta _{\text {hor}}\right) }\left. \left( 1+\eta ^2\frac{\mathrm{{d}}m}{\mathrm{{d}}\eta }\right) \right| _{\eta =\eta _\mathrm{{hor}}}, \end{aligned}$$

(128)

where we have introduced the notation

$$\begin{aligned} T_\mathrm{{H}}=\frac{\hbar c^3}{8\pi Gk_BnM_{\odot }}, \end{aligned}$$

(129)

corresponding to the Hawking temperature of the standard general relativistic Schwarzschild black hole. The variation of the dimensionless horizon temperature

$$\begin{aligned} \theta =\frac{T_\mathrm{{BH}}}{T_\mathrm{{H}}} \end{aligned}$$

(130)

of the black holes in quadratic Weyl geometric gravity is presented for selected values of the model parameters and for the general case, with the Weyl vector possessing both temporal and radial components, in Fig. 10.

Fig. 10

Variation of the dimensionless black hole horizon temperature for a Weyl geometric gravity black hole in the presence of both temporal and radial components of the Weyl vector for $\theta$ for $\Theta _{0}(0)=10^{-18}$, $\Omega _{0}(0)=10^{-11}$, and $\zeta (0)=10^{-28}$, respectively, and different values of $\psi (0)$ (left panel), and for $\psi (0)=10^{-15}$, and different values of $\zeta (0)$ (right panel). In both panels, the values of $\gamma$ are varied continuously

As one can see from Fig. 10, the temperature of the horizon of the Weyl black holes is generally higher than that of their general relativistic counterparts. This is also related to the wider range of event horizon positions as compared to the Schwarzschild black hole case. The Hawking temperature is a monotonically increasing function of the position of the event horizon, and it is strongly dependent on the initial values of the scalar field and its derivatives.

5.2 Specific heat of the Weyl black holes

Another important physical quantity characterizing the thermodynamic properties of the black holes is their specific heat $C_\mathrm{{BH}}$, which can be obtained from the definition [121, 122]

$$\begin{aligned} C_\mathrm{{BH}}= & {} \frac{\mathrm{{d}}M}{\mathrm{{d}}T_\mathrm{{BH}}}=\left. \frac{\mathrm{{d}}M}{\mathrm{{d}}r}\frac{\mathrm{{d}}r}{\mathrm{{d}}T_\mathrm{{BH}}}\right| _{r=r_\mathrm{{hor}}} \nonumber \\= & {} \frac{nM_{\odot }}{T_\mathrm{{H}}} \left. \frac{\mathrm{{d}}m\left( \eta \right) }{\mathrm{{d}}\eta }\frac{\mathrm{{d}}\eta }{\mathrm{{d}}\theta }\right| _{\eta =\eta _\mathrm{{hor}}}, \end{aligned}$$

(131)

where we have denoted $C_\mathrm{{H}}=nM_{\odot }/T_\mathrm{{H}}$. The variations of the dimensionless specific heat $C_\mathrm{{eff}}=C_\mathrm{{BH}}/C_\mathrm{{H}}$ of the Weyl black holes as a function of the dimensionless horizon radius in quadratic Weyl geometric gravity are represented for some selected values of the model parameters in Fig. 11.

Fig. 11

Variation of the dimensionless specific heat $C_\mathrm{{eff}}=C_\mathrm{{BH}}/C_\mathrm{{H}}$ of the Weyl geometric black holes in the presence of both temporal and radial components of the Weyl vector for $\Theta _{0}(0)=10^{-18}$, $\Omega _{0}(0)=10^{-11}$, and $\zeta (0)=10^{-28}$, respectively, and for different values of $\psi (0)$ (left panel) and for $\psi (0)=10^{-15}$, and different values of $\zeta (0)$ (right panel). In both panels, the values of $\gamma$ are varied continuously

Similarly to the Hawking temperature, the numerical values of $C_\mathrm{{eff}}$ depend on the initial conditions at infinity of the scalar field and its derivative. The specific heat of the Weyl geometric black holes is a rapidly decreasing function of $\eta _\mathrm{{hor}}$.

5.3 The entropy of the Weyl geometric black holes

The entropy $S_\mathrm{{BH}}$ of the Weyl geometric black hole is generally obtained as [121, 122]

$$\begin{aligned} S_\mathrm{{BH}}= & {} \int _{\infty }^{r_{\text {hor}}}{\frac{\mathrm{{d}}M}{T_\mathrm{{BH}}}}=\int _{\infty }^{r_{\text {hor}}}{\frac{1}{T_\mathrm{{BH}}}\frac{\mathrm{{d}}M}{\mathrm{{d}}r}\mathrm{{d}}r}. \end{aligned}$$

(132)

In the set of the dimensionless variables considered in the present study, we have

$$\begin{aligned} S_\mathrm{{BH}}\left( \eta _\mathrm{{hor}}\right) =C_\mathrm{{H}}\int _0^{\eta _\mathrm{{hor}}}{\frac{1}{\theta \left( \eta \right) }\frac{\mathrm{{d}}m\left( \eta \right) }{\mathrm{{d}}\eta }\mathrm{{d}}\eta }. \end{aligned}$$

(133)

In the following we denote $S_\mathrm{{eff}}=S_\mathrm{{BH}}\left( \eta _\mathrm{{hor}}\right) /C_\mathrm{{H}}$. The variation as a function of the dimensionless horizon radius $\eta _\mathrm{{hor}}$ of the entropy $S_\mathrm{{eff}}$ of the black holes in the Weyl geometric gravity theory with only the radial component in the Weyl vector is represented as a function of the dimensionless horizon radius for different values of $\psi (0)$ and $\zeta (0)$, in Fig. 12.

Fig. 12

Variation of the dimensionless entropy $S_\mathrm{{eff}}=S_\mathrm{{BH}}/C_H$ of the Weyl geometric black holes in the presence of both temporal and radial components of the Weyl vector as a function of the position of the event horizon for $\Theta _{0}(0)=10^{-18}$, $\Omega _{0}(0)=10^{-11}$, and $\zeta (0)=10^{-28}$, respectively, and for different values of $\psi (0)$ (left panel), and for $\psi (0)=10^{-15}$, and different values of $\zeta (0)$ (right panel). In both panels, the values of $\gamma$ are varied continuously

The black hole entropy decreases with the increase in the radius of the event horizon. While there is a significant dependence on the values at infinity of the scalar field, the dependence of the entropy on the derivatives of the field is weak, and no significant change in its values occurs due to the modifications of $\zeta (0)$.

Fig. 13

Variation of the evaporation time $\tau _\mathrm{{eff}}=\tau _\mathrm{{BH}}/\tau _\mathrm{{H}}$ of the Weyl geometric black holes in the presence of both temporal and radial components of the Weyl vector as a function of the position of the event horizon for $\Theta _{0}(0)=10^{-18}$, $\Omega _{0}(0)=10^{-11}$, and $\zeta (0)=10^{-28}$, and for different values of $\psi (0)$ (left panel), and for $\psi (0)=10^{-15}$, and different values of $\zeta (0)$ (right panel). In both panels, the values of $\gamma$ are varied continuously

5.4 The luminosity of the Weyl black holes

The formation and the evaporation of a spherically symmetric black hole in conformal gravity was investigated in [123] by considering the collapse of a spherically symmetric thin shell of radiation, leading to the formation of a singularity-free non-rotating black hole. The black hole has the same Hawking temperature as a Schwarzschild black hole with the same mass, and it completely evaporates in either a finite or an infinite time, depending on the ensemble. The evaporation process of a spherical neutral anti-de Sitter black hole in four-dimensional conformal gravity, where the equations of states are branched, was investigated in [124].

The luminosity of the Weyl geometric black holes due to the Hawking evaporation processes can be obtained according to the relation [121, 122]

$$\begin{aligned} L_\mathrm{{BH}}=-\frac{\mathrm{{d}}M}{\mathrm{{d}}}=-\sigma A_\mathrm{{BH}}T_\mathrm{{BH}}^4, \end{aligned}$$

(134)

where $\sigma$ is a parameter that depends on the considered model, while $A_\mathrm{{BH}}=4\pi r_\mathrm{{hor}}^2$ is the area of the event horizon. For the black hole evaporation time $\tau$, we thus have

$$\begin{aligned} \tau= & {} \int _{t_\mathrm{{in}}}^{t_\mathrm{{fin}}}{\mathrm{{d}}t}=-\frac{1}{4\pi \sigma }\int _{t_\mathrm{{in}}}^{t_\mathrm{{fin}}}{\frac{\mathrm{{d}}M}{r_{\text {hor}}^2T_\mathrm{{BH}}^4}}, \end{aligned}$$

(135)

where $t_{\textrm{in}}$ and $t_{\textrm{fin}}$ represent the initial and the final times considered for the evaporation process. Equivalently, we can obtain the evaporation times of the Weyl geometric black hole in the form

$$\begin{aligned} \tau _\mathrm{{BH}}\left( \eta _\mathrm{{hor}}\right) = \tau _\mathrm{{H}}\int _0^{\eta _\mathrm{{hor}}}{\frac{1}{\eta ^2\theta ^4\left( \eta \right) }\frac{\mathrm{{d}}m\left( \eta \right) }{\mathrm{{d}}\eta }\mathrm{{d}}\eta }, \end{aligned}$$

(136)

where we have denoted

$$\begin{aligned} \tau _\mathrm{{H}}=\frac{c^4}{8\pi G^2\sigma nM_{\odot }T_\mathrm{{BH}}^4}. \end{aligned}$$

(137)

The variations of the dimensionless Hawking evaporation time $\tau _\mathrm{{eff}}=\tau _\mathrm{{BH}}/\tau _\mathrm{{H}}$ of black holes in the Weyl geometric gravity as a function of the position of the event horizon are presented in Fig. 13. In this figure we have varied the values of $\gamma$ continuously, and kept the values of the other parameters fixed.

The Hawking evaporation times of the Weyl geometric black holes show a strong dependence on the initial conditions at infinity of both the scalar field and its derivative. They are rapidly decreasing functions of the radius of the event horizon, and they become negligibly small for a sufficient quantity of massive black holes.

6 Discussion and final remarks

The idea that conformal symmetry is an exact but spontaneously broken symmetry [102, 103] has many attractive features. If Einstein’s field equations are reformulated by strictly imposing conformal symmetry, then the conformal component of the metric field can be treated as a dilaton field with only renormalizable interactions. Hence, this condition imposes some strong constraints on the gravitational theories, and they are equivalent to demanding regularity of the action as the dilaton field variable tends to zero. Moreover, such a procedure can turn a black hole into a regular, topologically trivial soliton with no singularities, horizons, or firewalls [103].

In the present work, we have considered the possible existence of black hole type structures in the framework of the simplest conformally invariant gravitational theory, constructed ab initio in a Weyl geometry. The simplest such gravitational action contains the quadratic Weyl scalar and the electromagnetic type scalar, defined in a way similar to the Faraday tensor in standard electromagnetic theory, with the Weyl vector playing the role of the electromagnetic four-potential. This theory intrinsically contains a scalar degree of freedom, and, once reformulated in the standard Riemannian geometry, it can be linearized in the Ricci scalar. Hence, geometric conformally invariant Weyl gravity is equivalent in the Riemannian space to a scalar-vector-tensor theory, in which, besides the metric tensor, the gravitational properties are determined by a scalar and a tensor field, with both having purely geometric origins.

To investigate the physical properties of the scalar-vector-tensor Weyl geometric theory, in the present work we have considered one of the simplest possible physical and geometrical situations, namely, the case of the vacuum static and spherically symmetric systems. Even by adopting this simple theoretical model, the gravitational field equations are extremely complicated. Hence, in order to obtain solutions of the field equations, one must make extensive use of numerical methods. In order to do so in a computationally efficient way, one must first rewrite the static spherically symmetric Weyl geometric gravitational field equations in vacuum in a dimensionless form suitable for numerical integration. This can be achieved by introducing as an independent variable the inverse of the radial coordinate. Moreover, we have reformulated the gravitational field equations as a first-order dynamical system. In this representation, the numerical integration procedure is significantly simplified. In analogy with the Schwarzschild black hole solution of standard general relativity, we have also represented the metric tensor component $e^{\lambda }$ in terms of an effective mass. Hence, in order to obtain a numerical solution of the field equations of the Weyl geometric gravitational theory, we need to give the initial conditions at infinity of the metric tensor components, the scalar, and the vector fields, and their derivatives, respectively. As for the metric tensor components, we have assumed the condition of the asymptotically flat geometry, while the numerical values at infinity of the components of the scalar and vector fields and their derivatives have been chosen in an arbitrary way. We consider that the presence of a singularity in the field equations—or, more precisely, of a singular point in the behavior of the metric tensor components—indicates the formation of an event horizon and, consequently, indicates the presence of a black hole type object. The total mass of the black hole can be determined from the effective mass appearing in the radial metric tensor component, and it can be interpreted as containing the standard (baryonic) mass of the black hole to which the contributions from the scalar and Weyl type vector fields are added.

Generally, in spherical symmetry, only two components of the Weyl vector survive, namely, the temporal and the radial components, respectively. Based on this result, we have considered the solutions of the gravitational field equations of the quadratic Weyl geometric gravity in three cases, corresponding to the Weyl vector having a radial component only, a temporal component only, and in the presence of both components. As a first general conclusion of our study, we have detected the formation of an event horizon in all three cases. Consequently, these results indicate the formation of black holes within the quadratic geometric Weyl gravity theory. The position of the event horizon of the black holes depends significantly on the numerical values of the scalar and vector fields and their derivatives at infinity (representing the initial conditions for our dynamical systems). These results show the existence of an interesting relation between the asymptotic values of the scalar and Weyl fields, and the black hole properties. For example, in the case of the presence of the temporal component of the Weyl vector only, for particular values of the vector field and its derivative at infinity, the physical position $r_\mathrm{{hor}}$ of the event horizon of the black hole can be located at distances of the order of 0.33–1.27 of the standard Schwarzschild radius, a result that indicates the possibility of the formation of a large variety of Weyl type black holes, including the presence of more compact black holes having the same mass but a smaller radius than the standard Schwarzschild black hole of general relativity. On the other hand, black holes with the same mass but with radii larger than the Schwarzschild radius can also exist.

As a general rule, the position of the event horizon of the Weyl geometric black holes is also dependent on the coupling parameter $\gamma$ of the model. Thus, there is a multi-parametric dependence of the Weyl geometric black hole properties on the geometric couplings and on the asymptotic conditions at infinity of the metric and scalar and vector fields. We would also like to point out that in our numerical investigations, we could not detect any case of the formation of a naked singularity with the singularity located at the center $r=0$. In all the considered numerical cases, the black holes are hidden beyond an event horizon. But, of course, the possible existence in the present conformally invariant quadratic Weyl geometric model of naked singularities or topologically trivial soliton type structures cannot be excluded a priori. The numerical detection of these types of objects requires a significant extension of the range of values of the parameter space and more detailed numerical investigations of the dynamical systems describing static spherically symmetric quadratic Weyl geometric gravity structures in vacuum.

In the present work, we have presented only the numerical results obtained by considering the asymptotically flat conditions for the metric tensor at infinity. Numerical solutions corresponding to the de Sitter solutions can also be easily obtained. However, if the de Sitter condition is assumed to be of cosmological type, the deviations from the asymptotically flat case are negligible. But if one assumes a significant difference from flatness at infinity, several distinct types of black holes can be obtained. In particular, solutions corresponding to very large “super-massive black holes” can also be constructed by imposing appropriate initial conditions on the mass function m. One such example of super-massive black hole solutions of the Weyl geometric field equations in the presence of the radial component of the Weyl vector is presented in Table 10 for $m(0)=0.4$. In the coordinate $\eta$, these black holes have an extremely small value of the position of the event horizon, of the order of $10^{-5}$, much smaller than the position of the inverse of the event horizon of the Schwarzschild black holes. However, the physical radius of these particular types of black holes is very large, ranging from $r_{hor}\approx 10^4$ ($\eta _\mathrm{{hor}}=9.96\times 10^{-5}$), to $r_\mathrm{{hor}}\approx 3.69\times 10^4$ ($\eta _\mathrm{{hor}}=2.71\times 10^{-5}$). Such Weyl type structures may be considered as possible alternatives to the super-massive black holes of standard general relativity.

Table 10 Variation of the position of the event horizon $\eta _{\text {hor}}$ of the Weyl geometric black hole with radial component of the Weyl vector only for $\zeta (0)=1\times 10^{-30}$, $m(0)=0.4$ and different initial values of $\psi$

Even though most of our investigations have been performed by using a numerical approach, an exact solution of the field equations has also been obtained using analytical methods. In the case of a Weyl vector having only a radial component, if the metric functions satisfy the condition $\nu +\lambda =0$, an exact solution of the Weyl geometric gravity field equations can be obtained, with the metric given by a function of the form $e^{\nu }=e^{-\lambda }=A-B/r-Cr+Dr^2$, with A, B, C, and D constants that do not depend on the initial conditions at infinity. Thus, apart from the Schwarzschild term B/r, the metric tensor components contain two new terms proportional to r and $r^2$, respectively. We could not find similar simple analytical solutions for the cases of the presence of the temporal component of the Weyl only, or for the general case corresponding to the two-component Weyl vector. The analytical solutions are extremely useful in the study of the physical properties of black holes, and in particular for the study of the dynamics and motion of matter particles around them. They can also be applied for the investigation of the electromagnetic properties of thin accretion disks that may be present around black holes. The exact solution obtained could also help in discriminating quadratic Weyl geometric black holes from standard general relativistic black holes, and for obtaining observational constraints on the coupling constants and on the Weyl vector components.

We have also investigated in some detail the thermodynamic properties of the numerical black hole solutions of the quadratic Weyl geometric gravity. One extremely interesting physical property of black holes is their Hawking temperature, an essential parameter that has important theoretical implications. As compared to the Hawking temperature of the standard general relativistic Schwarzschild black holes, the horizon temperature of the quadratic Weyl geometric gravity black holes has a strong dependence on the initial conditions at infinity of the scalar and the Weyl vector fields. As one can observe from Fig. 10, the decrease in the physical horizon radius $r_{hor}$ (corresponding to an increase in the coordinate $\eta$) leads to an increase in the black hole temperature, a property that is specific to Weyl geometric black holes. In the case of the specific initial conditions for the scalar and Weyl vector fields considered in Fig. 10, the increase in the temperature is several times higher as compared to the Schwarzschild general relativistic case. Similar behavior appears for the specific heat, entropy, and evaporation time, respectively, of quadratic Weyl geometric black holes, with all these physical quantities strongly dependent on the initial conditions of the scalar and the Weyl vector field at infinity. The specific heat and the entropy of the Weyl black holes decrease with decreasing $r_\mathrm{{hor}}$ (increasing $\eta _\mathrm{{hor}}$). The black hole evaporation times may be very different for Weyl type black holes than Schwarzschild black holes, decreasing with increasing $r_\mathrm{{hor}}$. However, we would like to point out that our results on the thermodynamics of Weyl gravity black holes are obtained for a limited range of initial conditions at infinity for the scalar field and the Weyl vector, and therefore they may be considered as having a qualitative interpretation only. But, even within this limited level of numerical investigation, they provide an indicator of the complex physical behavior of the Weyl black holes, and the interesting and novel properties related to them.

The no-hair theorem is an important result in black hole physics [125,126,127,128]. In short, this theorem asserts that asymptotically flat black holes cannot have external nontrivial scalar fields possessing a non-negative field potential $V (\phi )$. It would be interesting to consider the no-hair theorem for Weyl geometric black holes. The preliminary, and mostly numerical, results obtained in the present work seem to point towards the conclusion that the no-hair theorem in its standard formulation may not be valid for quadratic Weyl black holes. All the black hole solutions we have considered have an asymptotically flat geometry, and scalar and vector fields exist around them. However, the question of whether these results follow from the particular choice of model parameters (coupling constant and initial conditions at infinity) or whether they are essential properties of the theory deserves further studies and investigations.

Quadratic Weyl gravity black holes have greater variability with respect to their basic properties, leading to more complicated external dynamics, as compared with the Schwarzschild black holes of general relativity. These richer properties follow from the presence of the scalar and vector degrees of freedom, resulting in very complicated and strongly nonlinear field equations. The effects associated with the scalar field and Weyl vector degrees of freedom could also lead to some specific astrophysical imprints and signatures, whose observational detection could open new perspectives on the possibilities of testing Weyl geometry and its associated gravitational theory on cosmic scales. The possible observational/astrophysical implications of the existence of quadratic Weyl geometric black holes will be considered in a future investigation.