Surrounded Vaidya solution by cosmological fields

Abstract

In the present work, we study the general surrounded Vaidya solution by the various cosmological fields and its nature describing the possibility of the formation of naked singularities or black holes. Motivated by the fact that real astrophysical black holes as non-stationary and non-isolated objects are living in non-empty backgrounds, we focus on the black hole subclasses of this general solution describing a dynamical evaporating-accreting black holes in the dynamical cosmological backgrounds of dust, radiation, quintessence, cosmological constant-like and phantom fields, the so called “surrounded Vaidya black hole”. Then, we analyze the timelike geodesics associated with the obtained surrounded black holes and we find that some new correction terms arise relative to the case of Schwarzschild black hole. Also, we address some of the subclasses of the obtained surrounded black hole solution for both dynamical and stationary limits. Moreover, we classify the obtained solutions according to their behaviors under imposing the positive energy condition and discuss how this condition imposes some severe and important restrictions on the black hole and its background field dynamics.

1 Introduction

Nowadays, we know that black holes are not just a mathematically possible solution to Einstein’s field equations, rather they seem to be some realistic astrophysical objects. It is more than a decade that we have obtained good evidences indicating that most of the galaxies, as our Milky Way, host many stellar active black holes, as well as a super-massive active black hole, in their centers. On the other hand, due to black hole evaporation [1, 2] and accretion–absorbtion processes [3, 4], it is accepted that the mass and other parameters of black holes are not fixed, and should change with time. Therefore, generally speaking, real black holes are non-stationary, and the stationary black holes such as Schwarzschild and Reissner–Nordström are only ideal models. Thus, the study of non-stationary black holes is meaningful and so motivating in the exploration of real black holes. There are a lot of research works on general dynamical black holes and their properties, see the works of Ashtekar and Krishnan [5, 6] and Hayward [7] as instances. Actually, studying black holes from astrophysical point of view and by astrophysicists has been originated in recent decades due to the dramatic increase in the number of black hole candidates from the sole candidate Cygnus X-1. This study needs a deeper understanding of the black hole physics and especially the black hole radiation by astrophysicists and relativists. Hawking used a quantum field theoretical approach to explore the black hole radiation, for the first time [1, 2]. Afterwards, some models to describe the classical essence of this radiation in a language which is free from the usual quantum field theoretic tools and is more familiar to the astrophysicists and relativists, have been introduced. For instance, the Vaidya solution [8,9,10] has provided a simple classical model for the black hole radiation and has been vastly investigated in this regard [11,12,13,14,15,16,17], see also [18,19,20,21,22,23,24,25,26,27,28,29] for more studies. In fact, the Vaidya solution is one of the non-static solutions of the Einstein field equations and can be regarded as a generalization of the static Schwarzschild black hole solution. This solution is characterized by a dynamical mass function depending on the retarded time coordinate u, i.e $m=m(u)$ and an ingoing/outgoing flow $\sigma (u,r)$. Thus, it can be implemented as a classical model for a dynamical black hole which is effectively evaporating or accreting, regarding its effective flow direction. On the other hand, the Vaidya solution has been used for studying the process of spherical symmetric gravitational collapse and as a testing ground for the cosmic censorship conjecture [30,31,32,33], see also [34] where a possible astrophysical application of the model for describing the energy source of gamma-ray bursts is discussed. These studies also are motivated by the time dependant mass parameter of this solution along with its outgoing radiation flow during the collapse ending by a naked singularity or a black hole. It is shown that if the outgoing flux diverges, the back-reaction will prevent the formation of naked singularity [35]. The observable sign of the formation of a naked singularity, by the collapse process, appears to be the burst of a radiation possessing a non-thermal spectrum, as the Cauchy horizon is approached [36]. Indeed, this is in contrast to the slow evaporation of a black hole via black-body spectrum of the Hawking radiation [1, 2]. Then, it would be important to carefully investigate Vaidya solution to better understanding of real dynamical black holes or the typical signs of naked singularities and to explore if there are any astrophysical objects whose properties resemble those of a naked singularity [36]. The Vaidya solution was generalized to the charged case known as the Bonnor–Vaidya solution [37], see also its application for example in [38,39,40,41,42]. Also, a generalisation of the Vaidya solution is introduced in [43]. This generalisation is based on the fact that the total supporting energy–momentum tensor of spacetime, constructed from type I and type II energy–momentum tensors [44], is linear in term of the mass function. Consequently, any linear superposition of particular solutions to the Einstein field equations will also a solution. Then, using this approach, we can construct more general solutions such as the Bonnor–Vaidya [37], Vaidya–de Sitter [45], radiating dyon solution [46], Bonnor–Vaidya–de Sitter [40, 47,48,49] and the Husain solution [50].

On the other hand, a new exact static solution to Einstein field equations has been recently introduced by Kiselev [51]. Actually, the Kiselev solution is nothing but the static generalization of the Schwarzschild solution to include a non-empty cosmological background, especially well known for the quintessence background. This generalization is well motivated by the fact that black holes in real world are not isolated and are not embedded in empty backgrounds. The black hole solutions coupled to matter fields, such as Kiselev solution, are of interest in studying astrophysical distorted black holes [52,53,54,55], as well as in exploring the no hair theorems [56,57,58,59]. Indeed, a crucial assumption for the no-hair theorem is that the black hole is isolated, i.e., the spacetime is asymptotically flat and contains no other sources. However, in real world astrophysical situations this requirement is not fulfilled, for examples, for black holes in binary systems, for black holes surrounded by plasma, or black holes having an accretion disk or jets in their vicinity. All these situations indicate that a black hole may put on different types of wigs. For these cases, the standard no-hair theorem for the isolated black holes can be questioned, see for examples [57, 60]. In a recent research, the authors of [61] discussed on distinguishing rotating Kiselev black hole from naked singularity using spin precession of test gyroscope. In general, since black holes possess strong gravitational attraction such that their nearby matter, even light, cannot escape from their gravitational field, they cannot be observed directly and there are some different ways to detect them in binary systems as well as at the centers of their host galaxies. The most promising way is the accretion process. In the language of astrophysics, the accretion is defined as the inward flow of matter fields surrounding a compact object, such as black holes and neutron stars, due to the gravitational attraction. Then, the process of accretion into black holes is one of the most interesting research fields in relativistic astrophysics [62,63,64,65,66]. This process may be described by a perfect fluid coupled to general relativity representing a plasma which obeys the equations of ideal or resistive magnetohydrodynamics or a fluid coupled to radiation. Such accretion processes along with their detailed physical descriptions, can be found in [67] and references therein, see also [68,69,70,71,72,73,74,75]. On the other hand, there are also other kind of accretion processes related to the black holes surrounded by exotic matter fields as potential models of dark energy, whose existence and features are motivated by the problems in the standard model of cosmology. A number of theoretical and observational studies confirmed that our universe in its early stages experienced an inflation process while it is undergoing an accelerated expansion in the late time. In order to explain these events, an energy component, known as the dark energy, is required to be introduced to the framework of general theory of relativity. The cosmological constant is a leading candidate for dark energy while there are other proposals including the dynamical scalar fields such as quintessence and phantom fields. In the Bousso’s work [76], one finds that “Q-space exhibits thermodynamic properties similar to those of the de Sitter horizon. The horizon radius in Q-space grows linearly with time, and consequently the temperature slowly decreases. We find that this behavior is consistent with the first law of thermodynamics: the temperature and entropy respond appropriately to the flux of quintessence stress-energy across the horizon”, for a cosmological setup, where “Q-space” stands for quintessence-space. This important result along with the observational data confirming a dark energy fluid responsible for the accelerating expansion of Universe with the equation of state parameter $\omega <- \ 1/3$, has motivated the community to study in detail the black hole solutions in the quintessence background. For some recent studies of Kiselev black holes, see [77] for its generalization to rotating case, [78,79,80] for quasinormal modes and Hawking radiation, [81,82,83,84] for thermodynamical studies, [85,86,87,88,89,90] for trajectories and particle dynamics around this black hole, [91] for accretion process and [92] for gravitational lensing among the others. One should also note that the Kiselev solution can be implemented for more generic backgrounds of dust, radiation, quintessence, cosmological constant and phantom fields as well as for any realistic combination of these cosmological fields. Then, by the presence of such fields around the black holes, one may have interest to explore some interesting facts such as whether black holes have hair or scalar wigs [93], how black holes affect these cosmological surrounding fields and what are the consequences or what are the influences of these surrounding fields on the features, behaviors and abundance of black holes. In this regard, one may find the reference [94] as a good review including various scenarios of accretion process into black holes, see also [69, 95] for charged black hole accretion. Among the all of the accretion processes, the most interesting one are related to those that the accretion of the surrounding fields enforcing a black hole to shrink. These surrounding field include the scalar fields or fluid violating the weak energy condition, i.e $\rho >0$ and $\rho +p>0$ [94]. Specific scenarios involving the accretion of phantom energy have shown that the black-hole area decreases with the accretion [96,97,98,99]. For example, in [96], it is shown that black holes will gradually vanish as the universe approaches a cosmological big rip state. The big rip scenario for a cosmos occurs when its filling dark energy is the phantom energy with $p< -\rho$. In this scenario, the cosmological phantom field disrupts finally all bounded objects of the universe up to sub-nuclear scales. For the test-field approximation, one may find the accretion process of a scalar field violating the energy conditions leading the decrease in the black holes area in [98, 100]. Moreover, the shrink of the black hole area through the accretion of a phantom scalar field has been confirmed in full nonlinear general relativity [101, 102]. In this regard, the shrink of the black hole area by the accretion of a potentially surrounding field is an interesting phenomena in the sense that it can be an alternative process for black hole evaporation through the Hawking radiation or even be an auxiliary for speeding up it. One physical explanation for a black hole mass diminishing may be is that accreting particles of a phantom scalar field have a total negative energy [103]. Similar particles possessing negative energies are created through the Hawking radiation process and also in the energy extraction process from a black hole by the Penrose mechanism. The effect of phantom-like dark energy onto a charged Reissner–Nordström black hole is studied in [104] and it is found that accretion is possible only through the outer horizon. On the other hand, for scalar fields regarding the energy conditions, there is a possibility indicating that the accretion of a scalar field can be partial such that the amount of accreted scalar field depends on features of the incident wave packet, i.e. the wave number and the width of the packet. This has been studied both in the test-field approximation [105] and in full general relativity [101, 102]. In this line, some studies in the test-field limit indicate that a scalar field can also be sustained by a black hole without being accreted [106].

In the present work, following the approach of [51, 107] introduced for the static black holes, and motivated by the facts that real astrophysical black holes are neither stationary nor isolated and are not embedded in empty backgrounds, we wish to find a more realistic dynamical solution for the classical description of the evaporating-accreting black holes in generic dynamical backgrounds. The organization of the paper is as follows. In Sect. 2, we introduce the general surrounded Vaidya solution, its nature describing the possibility of the formation of naked singularities or black holes, interaction of its possible black holes with their backgrounds as well as its timelike geodesic analysis in the general form. Then, in Sects. 3–7, we investigate in detail the special classes of this solution as the surrounded Vaidya black hole by the dust, radiation, quintessence, cosmological constant-like and phantom fields, respectively. The paper ends with a conclusion, in Sect. 8.

2 The general surrounded Vaidya solutions

In this section, we are looking for the general surrounded Vaidya solutions by the approach of [51, 107]. Then, we consider the general spherical symmetric spacetime metric in the form of

$$\begin{aligned} ds^{2}=-f(u,r)du^2+2\epsilon dudr+r^2d\Omega ^2,\quad \epsilon =\pm 1, \end{aligned}$$

(1)

where $d\Omega ^2=d\theta ^2+sin^2\theta d\phi ^2$ is the two dimensional unit sphere and f(u, r) is a generic metric function depending on both of the advanced/retarded time coordinate u and the radial coordinate r. The cases, $\epsilon =- \ 1$ and $\epsilon =+ \ 1$ represent the outgoing and ingoing flows corresponding to the effectively evaporating and accreting Vaidya black hole solutions, respectively. Using the metric (1), we obtain nonvanishing components of the Einstein tensor as

$$\begin{aligned}&{G^{0}}_{0}={G^{1}}_{1}=\epsilon G_{01}=\epsilon G_{10}=\frac{1}{r^2} (f^{\prime }r-1+f ),\nonumber \\&{G^{1}}_{0}=\epsilon G_{00}+fG_{01}=-\frac{\dot{f}}{r},\nonumber \\&{G^{2}}_{2}=\frac{1}{r^2}G_{22}=\frac{1}{r^2}\left( rf^{\prime }+\frac{1}{2}r^2 f^{\prime \prime }\right) ,\nonumber \\&{G^{3}}_{3}=\frac{1}{r^2 sin^2 \theta }G_{33}=\frac{1}{r^2}\left( rf^{\prime }+\frac{1}{2}r^2 f^{\prime \prime }\right) , \end{aligned}$$

(2)

where dot and prime signs represent the derivatives with respect to the time coordinate u and the radial coordinate r, respectively. Then, the total energy–momentum supporting this spacetime should have the following non-diagonal form

$$\begin{aligned} {T^{\mu }}_{\nu }= \begin{pmatrix}{T^{0}}_{0} &{} 0 &{} 0 &{} 0 \\ {T^{1}}_{0} &{} {T^{1}}_{1} &{} 0 &{} 0 \\ 0 &{} 0 &{} {T^{2}}_{2} &{} 0 \\ 0 &{} 0 &{} 0 &{} {T^{3}}_{3}\\ \end{pmatrix}, \end{aligned}$$

(3)

where also must obey the symmetries in Einstein tensor ${G^{\mu }}_{\nu }$. With respect to the field equations in (2), the equalities ${G^{0}}_{0}={G^{1}}_{1}$ and ${G^{2}}_{2}={G^{3}}_{3}$ require ${T^{0}}_{0}={T^{1}}_{1}$ and ${T^{2}}_{2}={T^{3}}_{3}$, respectively. Then, for the nature of the Vaidya solution in the presence of a dynamical background, one can consider a total energy–momentum tensor supporting the Einstein field equations in the following form

$$\begin{aligned} {T^{\mu }}_{\nu }={\tau ^{\mu }}_{\nu }+{\mathcal {T}^{\mu }}_{\nu }, \end{aligned}$$

(4)

where ${\tau ^{\mu }}_\nu$ is the energy–momentum tensor associated to the Vaidya null radiation-accretion as

$$\begin{aligned} {\tau ^{\mu }}_{\nu }=\sigma k^{\mu }k_{\nu }, \end{aligned}$$

(5)

such that $\sigma =\sigma (u,r)$ is the measure of the energy flux or the energy density of the outgoing radiation-ingoing accretion flow [108] and $k_\mu ={\delta ^{0}}_{\mu }$ is a null vector field while ${\mathcal {T}^{\mu }}_{\nu }$ is the energy–momentum tensor of the surrounding fluid defined as in [51]

$$\begin{aligned} {\mathcal {T}^{0}}_{0}= & {} -\rho _s(u,r),\nonumber \\ {\mathcal {T}^{i}}_{j}= & {} -\rho _{s}(u,r)\alpha \left[ -(1+3\beta )\frac{r_i r^j}{r_n r^n}+\beta {\delta ^{i}}_{j}\right] , \end{aligned}$$

(6)

where subscript “s” stands for the surrounding field which can be a dust, radiation, quintessence, cosmological constant, phantom field or even any complex field constructed by the combination of these fields.¹ This form of energy–momentum for the surrounding fluid is implying that the spatial profile of the Vaidya solution surrounding energy–momentum tensor is proportional to the time component, describing the dynamical energy density $\rho _s(u,r)$, with the arbitrary constant parameters $\alpha$ and $\beta$ depending the internal structure of the surrounding fields. The isotropic averaging over the angles results in [51]

$$\begin{aligned} \langle {\mathcal {T}^{i}}_{j}\rangle =\frac{\alpha }{3}\rho _{s}(u,r){\delta ^{i}}_{j}=p_{s}(u,r){\delta ^{i}}_{j}, \end{aligned}$$

(7)

since we considered $\langle r^{i}r_{j}\rangle =\frac{1}{3}{\delta ^{i}}_{j}r_n r^n$. Then, we have the barotropic equation of state for the surrounding field as

$$\begin{aligned} p_s(u,r)=\omega _s \rho _s(u,r), \quad \omega _s=\frac{1}{3}\alpha , \end{aligned}$$

(8)

where $p_s(u,r)$ and $\omega _s$ are the dynamical pressure and the constant equation of state parameter of the surrounding field, respectively.² Thus, regarding the Einstein tensor components in (2) and the total energy–momentum tensor given by the Eqs. (3)–(6), we have ${\mathcal {T}^{0}}_{0}={\mathcal {T}^{1}}_{1}$ and ${\mathcal {T}^{2}}_{2}={\mathcal {T}^{3}}_{3}$. These exactly provide the so called principle of additivity and linearity considered in [51] in order to determine the free parameter $\beta$ of the energy momentum-tensor ${\mathcal {T}^{\mu }}_{\nu }$ of the surrounding field as

$$\begin{aligned} \beta =-\frac{1+3\omega _s}{6\omega _s}. \end{aligned}$$

(9)

Then, by substituting $\alpha$ and $\beta$ parameters in (8) and (9) into (6), the non-vanishing components of the surrounding energy–momentum tensor ${\mathcal {T}^{\mu }}_{\nu }$ will be

$$\begin{aligned}&{\mathcal {T}^{0}}_{0}={\mathcal {T}^{1}}_{1}=-\rho _s(u,r),\nonumber \\&{\mathcal {T}^{2}}_{2}={\mathcal {T}^{3}}_{3}=\frac{1}{2}\left( 1+3\omega _s\right) \rho _s(u,r). \end{aligned}$$

(10)

Now, by having the Einstein tensor components and the corresponding total energy–momentum tensor ${T^{\mu }}_{\nu }$, one can obtain the associated field equations. Then, the ${G^{0}}_{0}={T^{0}}_{0}$ and ${G^{1}}_{1}={T^{1}}_{1}$ components of the Einstein field equations give the following differential equation

$$\begin{aligned} \frac{1}{r^2} (f^{\prime }r-1+f )=-\rho _s. \end{aligned}$$

(11)

Similarly, the ${G^{1}}_{0}={T^{1}}_{0}$ component leads to

$$\begin{aligned} -\frac{\dot{f}}{r}=\epsilon \sigma , \end{aligned}$$

(12)

and ${G^{2}}_{2}={T^{2}}_{2}$ and ${G^{3}}_{3}={T^{3}}_{3}$ components read as

$$\begin{aligned} \frac{1}{r^2} (rf^{\prime }+\frac{1}{2}r^2 f^{\prime \prime } )=\frac{1}{2}(1+3\omega _{s} )\rho _{s}. \end{aligned}$$

(13)

Thus, we see that there are three unknown dynamical functions f(u, r), $\sigma (u,r)$ and $\rho _s(u,r)$ which can be determined analytically by the above three differential equations. Simultaneous solving the differential equations (11) and (13), one obtains the following solution for the metric function

$$\begin{aligned} f(u,r)=1-\frac{2M(u)}{r}-\frac{N_s(u)}{{r}^{{3\omega _s +1}} }, \end{aligned}$$

(14)

with the energy density $\rho _s (u,r)$ of the surrounding field in the form of

$$\begin{aligned} \rho _s (u,r)=- \frac{3\omega _s N_{s}(u)}{{r}^{{3(\omega _s +1)}}}, \end{aligned}$$

(15)

where M(u) and $N_s(u)$ are integration coefficients representing the Vaidya dynamical mass and the surrounding dynamical field structure parameter, respectively.

On the other hand, respecting to the weak energy condition imposing the positivity of any kind of energy density of the surrounding field, i.e $\rho _s\ge 0$, demands

$$\begin{aligned} \omega _s N_{s}(u)\le 0. \end{aligned}$$

(16)

This implies that for the surrounding field with a positive equation of state parameter $\omega _s$, it is needed to have $N_{s}(u)\le 0$ and conversely for a negative $\omega _s$, it is required to have $N_{s}(u)\ge 0$. Then, this condition determines the gravitational nature of the term associated to surrounding field in the metric function f(u, r).

Regarding the metric function (14), the spacetime metric (1) reads as

$$\begin{aligned} ds^{2}=-\left( 1-\frac{2M(u)}{r}-\frac{N_{s}(u)}{{r}^{{3\omega _s +1}} }\right) du^2+2\epsilon dudr+r^2d\Omega ^2, \end{aligned}$$

(17)

representing an effectively evaporating-accreting Vaidya spacetime in a dynamical background. One may realize the following distinct subclasses of this general solution as

• The solution by setting $f=f(u,r)$ and $\rho _s=\rho _s(r)$ in the field equations (11)–(13).

  These considerations lead to $M=M(u)$ and $N_s=constant$ in the metric function f(u, r) and $\sigma \ne 0$ for the black hole’s radiation density. In this case, there is no dynamics in the surrounding field and consequently there is no accretion to the black hole. Indeed, this case represents an evaporating black hole solution with $\epsilon =-1$ in a static background. Then, the evaporating black hole in an empty background, i.e $\rho _s=0$ [8, 9], and (anti)-de Sitter space, i.e $\rho _s = \rho _{\Lambda }=constant$ [45, 109, 110], are special subclasses of our general solution. Some interesting physical features of these solutions can be found in the references [11, 12, 31, 42, 111,112,113,114,115,116,117,118,119].

• The solution by setting $f=f(r)$ and $\rho _s=\rho _s(r)$ in the field equations (11)–(13).

  These considerations lead to $M=constant$, $N_s=constant$ in the metric function and $\sigma =0$ for the radiation-accretion density. This case represents a non-dynamical back hole in a static background and consequently, there are no accretion and evaporation. The Schwarzschild black hole as well as its generalization to (anti)-de Sitter background are two special subclasses of our general solution. For a general background, not just the (anti)-de Sitter background, it is interesting that using the following coordinate transformation

  $$\begin{aligned} du=dt+\frac{\epsilon dr}{1-\frac{2M}{r}-\frac{N_{s}}{{r}^{{3\omega _s +1}} }}, \end{aligned}$$

  (18)

  one can obtain the general static solution of the Schwarzschild black hole surrounded by a surrounding field as

  $$\begin{aligned} ds^2=-\left( 1-\frac{2M}{r}-\frac{N_s}{r^{3\omega _s +1}} \right) dt^2 +\frac{dr^2}{ 1-\frac{2M}{r}-\frac{N_s}{r^{3\omega _s +1}} } +r^2 d\Omega ^2, \end{aligned}$$

  (19)

  which was found by Kiselev [51]. Then, the Kiselev solution also can be obtained as a subclass of our general dynamical solution (17) in the stationary limit.

• The solution for $\epsilon =+1$ with changing the background field parameters as $\omega _s\rightarrow \frac{1}{3}(2k-1)$ and $N_s(u)\rightarrow -\frac{2g(u)}{2k-1}$.

  By this considerations, we recover the Husain solution describing a null fluid collapse [50] as

  $$\begin{aligned} ds^{2}=-\left( 1-\frac{2M(u)}{r}+\frac{2g(u)}{(2k-1){r}^{{2k}} }\right) du^2+2\epsilon dudr+r^2d\Omega ^2, \end{aligned}$$

  (20)

  with the energy density

  $$\begin{aligned} \rho _s (u,r)=\frac{2g(u)}{{r}^{{2k+2}}}. \end{aligned}$$

  (21)

This solution and it various applications are widely studied in the literature, see for instances [120, 121] and [122, 123] where a barotropic equation of state is considered for the collapse study. There is a difference in the method obtaining the solutions in the present work and in [50] as well as in the other mentioned works. The solution (20), as in [50], is obtained by the “pre-imposed” equation of state $p=k\rho ^a$, whereas in our approach, the effective equation of state is resulting from the isotropic averaging over the angles for the surrounding field distribution. Our approach is motivated by the present anisotropy in the Einstein tensor components (2) and the corresponding total energy–momentum tensor (3), such that the surrounding fluid behaves effectively as a perfect fluid with the effective (averaged) equation of state $p_s(u,r)=\omega _s \rho (u,r)$, see (7). As the advantage of this averaging method, one can substitute for $\omega _s$ the same known cosmological field equation of state parameters $\frac{1}{3},0,-1,-\frac{2}{3}$ and $-\frac{4}{3}$ for the radiation, dust, cosmological, quintessence and phantom fields, respectively, when the black hole is embedded in these cosmological backgrounds. Substituting the same values of cosmological parameters for k, through $p=k\rho ^a$ even for $a=1$, in (20) gives different solutions with respect to (17) for the general dynamical case as well as for the known static solution in [51] in the stationary limit, by doing a similar transformation to (18). For example, throwing a bunch of dust with the mass of $M_{dust}(=g)$ to the black hole with mass M, one expects a resulting metric for the final black hole as $f(r)=1-\frac{2M_{eff}}{r}$ where $M_{eff}=M+M_{dust}$, whereas substituting $k=0$ in the metric (20) gives $f(r)=1-\frac{2M}{r}-2M_{dust}$ which seems to be incorrect due to the gravitational potential form of the final black hole and also the dimensional consideration. One also realizes that, as we will see in the next sections of the paper, there is a possibility of the formation of both the naked singularities and black holes in different backgrounds for the solution with $\epsilon =-1$.

2.1 The analysis of naked singularity or black hole formations

In order to investigate the formation of naked singularity or black hole associated to the obtained solution (17), we follow the approach of [123]. The equation for the radial null geodesics using the metric (1), or (17), can be obtained by setting $ds^{2}=0$ and $d\Omega _{2}^{2}=0$ as

$$\begin{aligned} \frac{du}{dr}=\frac{2\epsilon }{f(u,r)}. \end{aligned}$$

(22)

This system has a singularity at $r=0,~u=0$. Defining the function X as $X=\frac{u}{r}$ gives us the possibility of studying the limiting behavior of X as we approach the singularity located at $r=0,~u=0$, along the radial null geodesics. Denoting this limiting value of X by $X_{0}$, we have

(23)

Using the metric function (14) in (23), we obtain

(24)

Now, following the method of [123] for our case, we consider $M(u)=m u$ and $N_s(u)=n u^{3\omega _s+1}$, where m and n are constants. Thus, using (24), we obtain the following algebraic equation in terms of $X_{0}$

$$\begin{aligned} n X_{0}^{3\omega _s+2}+2m X_{0}^{2}-X_{0}+2\epsilon =0. \end{aligned}$$

(25)

A black hole will be formed if one obtains only non-positive solutions of this equation. However, if we find a positive real root for (25), then this system describes a naked singularity and consequently provides counterexamples for the cosmic censorship conjecture by Penrose [30]. It is difficult to find exact solutions for $X_{0}$ in (25) for the generic values of $n, m, \epsilon$ and $\omega _s$ parameters. However, as a result, one can find that there are possibilities of the formation of both the naked singularities and black holes for the various backgrounds of dust, radiation, quintessence, cosmological constant-like and phantom backgrounds for some particular ranges of m and n parameters. We postpone the detailed study of this equation for the mentioned backgrounds, to the Sects. 3–7.

2.2 The analysis of the black hole-background field interactions

Because in this work, we are mainly interested in the possible interactions between a dynamical black hole and its surrounding background, hence regarding the possibility of formation of black holes as mentioned in the previous subsection and as we will see in the Sects. 3–7, here we consider only the case that black holes are formed and we analyze in detail the general radiation-accretion profile for the corresponding systems and classify the possible situations under the positive energy condition.

Substituting the metric function (14) in the Eq. (12) gives the radiation-accretion density of the effectively evaporating-accreting black hole as

$$\begin{aligned} \sigma (u,r)=\epsilon \left( \frac{2\dot{M}(u)}{r^{2}} +\frac{\dot{N}_{s}(u)}{r^{3\omega _s +2}} \right) , \end{aligned}$$

(26)

where the first and second terms in RHS are the radiation-accretion density corresponding to the mass change of the black hole and the dynamics of the surrounding field, respectively. This shows that for construction of a realistic effectively evaporating-accreting black hole model, one needs to implement such a solution including a dynamical black hole in a dynamical background described by the energy–momentum (10). Considering (26), the following points can be realized.

• By turning off the background field dynamics, i.e. $\dot{N}_s(u)=0$, we recover the energy flux associated to the mass change of the central black hole corresponding to the original Vaidya solution [8, 9]. See [108] for more discussion on the properties of the original Vaidya solution.

• For the background field possessing $\omega _s>0$, if $\dot{M}(u)$ and $\dot{N}_s(u)$ have a same order of magnitude, the surrounding background field contribution to the total density $\sigma (u,r)$ is dominant near the black hole while at far distances from the black hole it decreases faster than the contribution of the black hole mass changing term. In contrast, for the background field possessing $\omega _s<0$, the surrounding background field contribution is dominant at large distances while the black hole contribution is dominant near the black hole itself. Then, from the astrophysical point of view, the detected amount of the radiation-accretion density by the observer not only depends on the distance from the black hole but also depends on the nature of background field.

Considering the positive energy density condition (by the weak energy condition) on the total radiation-accretion density $\sigma (u,r)$ in (26) requires

$$\begin{aligned} \epsilon \left( \frac{2\dot{M}(u)}{r^{2}}+\frac{\dot{N}_{s}(u)}{r^{3\omega _s +2}}\right) \ge 0. \end{aligned}$$

(27)

This inequality confines the dynamical behaviours of the black hole and its background field at any time and distance (u, r). In the case of a static background, as in the Vaidya’s original solution [8, 9], it is required that $\epsilon$ and $\dot{M}(u)$ have the same signs to have positive energy density. This shows that for a radiating black hole with $\dot{M}(u)<0$ we have $\epsilon =-1$ which represents the outgoing null flow, while for an accreting black hole it is required to have $\epsilon =+1$, representing the ingoing null flow. In the presence of the background dynamics, it is not mandatory that $\epsilon$ and $\dot{M}(u)$ take the same signs and the satisfaction of the positive energy density condition can be achieved even by their opposite signs depending on the background field parameters $\dot{N}_s(u)$ and $\omega _s$. Based on the relation (27), the dynamical behaviour of the background field is governed by

$$\begin{aligned} {\left\{ \begin{array}{ll}\dot{N}_{s}(u)\le -2\,{r}^{3\omega _s }\,\dot{M}(u), &{} \epsilon =-1, \\ \\ \dot{N}_{s}(u)\ge -2\,{r}^{3\omega _s }\,\dot{M}(u), &{}\epsilon =+1. \end{array}\right. } \end{aligned}$$

(28)

Then, at any distance r from the black hole, the background field must obey the above conditions. One astrophysical importance of such a physical constraint is that the observer knows the dynamical range of the background field at any distance and that, prior to any observation, he knows how to include or remove the background field contribution if he is only interested in black hole’s contribution, or vice versa. Interestingly, for the special case of $\dot{N}_{s}(u)= -2\,{r}^{3\omega _s }\,\dot{M}(u)$, there is no pure radiation-accretion density, i.e $\sigma (u,r)=0$. This case corresponds to two possible physical situations. The first one is related to the situation where for any particular distance $r_0$, the background $\dot{N}(u)$ and black hole $\dot{M}(u)$ behave such that their contributions cancel out each others leading to $\sigma (u,r_0)=0$. The second situation is related to the case that for the given dynamical behaviors of the black hole and its background, one can always find the particular time dependent distance

$$\begin{aligned} r_{*}(u)=\left( -\frac{\dot{N}_{s}(u)}{2\dot{M}(u)}\right) ^{\frac{1}{3\omega _{s}}}, \end{aligned}$$

(29)

possessing zero energy density $\sigma (u,r_{*}(u))$. For the case of constant rates of $\dot{N}_{s}(u)$ and $\dot{M}(u)$, the distance $r_{*}$ is fixed to a particular value. To have a particular distance at which the density $\sigma (u,r_{*})$ is zero, the positivity of $r_*$ also requires that $\dot{M}(u)$ and $\dot{N}_{s}(u)$ have opposite signs. For the cases in which $r_*$ is not positive, the lack of a positive real value radial coordinate is interpreted as follows: the radiation-accretion density $\sigma (u,r)$ never and nowhere vanishes.

In the case of being the positive radial coordinate $r_*$, for the given radiation-accretion behaviors of the black hole and its surrounding field, i.e $\dot{M}(u)$ and $\dot{N}_{s}(u)$, it is possible to find a distance at which we have no any radiation-accretion energy density contribution. In other words, it turns out that the rate of outgoing radiation energy density of the black hole is exactly balanced by the rate of ingoing absorption rate of surrounding field at the distance $r_{*}$ and vice versa. Beyond or within this particular distance, the various general situations can be realized in the Tables 1 and 2 for the black hole (BH) and its surrounding field (SF). One practical importance of (29) for an astrophysicist is that a particle detector at this distance will detect vanishing radiation-accretion density.

Table 1 General BH and SF parameters for $\epsilon =-1$

Table 2 General BH and SF parameters for $\epsilon =+1$

Then, regarding these tables and Eq. (27), we find the following results.

• The cases possessing negative values of $r_*$ (the cases I, IV, V and VIII) mean that the radiation-accretion density does not vanish somewhere and forever. Among these cases, the ones which have positive $\sigma (u,r)$ are only physical, i.e the cases IV and VIII for $\epsilon =-1$, and I and V for $\epsilon =+1$. Then, one realizes that how the weak energy condition causes in practice the nonphysical events to be hidden to an astrophysicist aiming to investigate a black hole and his surrounding field.

• The remaining positive values of $r_*$, corresponding to a zero radiation-accretion density, are physically viable and their corresponding physical processes are listed in the last column. These properties are determined according to the behaviours of the parameters $\epsilon$, $\omega _s$, and quantities $\dot{M}(u)$, $\dot{N}_{s}(u)$, and $\sigma (u,r)$. Those values of $r_*$ corresponding to the negative energy density $\sigma (u,r)$ represent no physical situation about the evaporation–absorption or accretion. The real features of those regions are hidden by the weak energy condition. Then, it is physically reasonable to do any astrophysical experiment in the regions respecting the energy condition.

• For $\omega _s>-\frac{2}{3}$, the particular distance $r_*$, where $\sigma (u,r)$ vanishes, corresponds to two possible cases as $r_{*}(u)=\left( -\frac{\dot{N}_{s}(u)}{2\dot{M}(u)}\right) ^{\frac{1}{3\omega _{s}}}$ and $r_*=\infty$. In the first case, for $-\frac{2}{3}<\omega _s<0$ with $|\dot{N}_s(u)|\ll |\dot{M}(u)|$ and for $\omega _s\ge 0$ with $|\dot{M}(u)|\ll |\dot{N}_{s}(u)|$, we have $r_*\rightarrow \infty$. This means that the first situation indicates that black hole evolves very faster than its background while the second indicates that black hole evolves very slow relative to its background. By satisfaction of these dynamical conditions to hold $r_*\rightarrow \infty$, the positive energy density is respected everywhere in the spacetime. Then, in practice, an astrophysical observer can detect a radiation-accretion density resulting from the interaction of the black hole with its surrounding field even at far distances, in which for $-\frac{2}{3}<\omega _s<0$ and $\omega _s\ge 0$ the main contribution in the detected radiation-accretion density belongs to the black hole and surrounding field, respectively. In other cases, the positive energy density will be respected in some regions while violated beyond those regions.

• For $\omega _s\le -\frac{2}{3}$, the particular distance $r_*$ is given as $r_{*}(u)=\left( -\frac{\dot{N}_{s}(u)}{2\dot{M}(u)}\right) ^{\frac{1}{3\omega _{s}}}$. Then, for a rapidly evolving black hole relative to its background, i.e $|\dot{N}_s(u)|\ll |\dot{M}(u)|$, we have $r_*\rightarrow \infty$. This case implies an evolving black hole in an almost static background in which the positive energy condition is respected everywhere in this spacetime. Then, for $\omega _s\le -\frac{2}{3}$ representing a dark energy fluid, an astrophysicist finds that it is the black hole which has the main contribution in the radiation-accretion density.

2.3 Timelike geodesics for the surrounded black holes

The geodesics for our metric (1), or (17), will all lie on a plane due to the spherical symmetry in which for the sake of simplicity, one can choose $\theta =\pi /2$. The geodesic equations for the above spacetime metric can be derived by varying the following action

$$\begin{aligned} I=\int \mathcal {L}d\tau =\frac{1}{2}\int \left( -f(u,r)\overset{*}{u}^2+2\epsilon \overset{*}{u} \overset{*}{r}+r^2 \overset{*}{\varphi }^2 \right) d\tau , \end{aligned}$$

(30)

where the star sign denotes the derivative with respect to the proper time $\tau$. Then, we have the following three equations

$$\begin{aligned} \overset{*}{\varphi }=\frac{L}{r^2}, \end{aligned}$$

(31)

and

$$\begin{aligned} -\frac{1}{2}f^{\prime }\overset{*}{u}^2 +r\overset{*}{\varphi }^2 -\epsilon \overset{**}{u}=0, \end{aligned}$$

(32)

and

$$\begin{aligned} \epsilon \overset{**}{r}=\frac{1}{2}\dot{f} \overset{*}{u}^2+f\overset{**}{u}+f^{\prime }\overset{*}{u}\overset{*}{r}, \end{aligned}$$

(33)

for $\varphi , r$ and u variables respectively, where L is the conserved angular momentum per unit mass and dot and prime signs denote the derivative with respect to u and r, respectively.³ Using (31) in (32), one finds

$$\begin{aligned} f\overset{**}{u}=\epsilon f\frac{L^2}{r^3}-\frac{1}{2}\epsilon f f^\prime \overset{*}{u}^2. \end{aligned}$$

(34)

On the other, using the timelike geodesics condition as $g_{\mu \nu }\dot{x}^{\mu }\dot{x}^{\nu }=-1$, one finds

$$\begin{aligned} f^\prime \overset{*}{r} \overset{*}{u}=-\frac{1}{2}\epsilon f^\prime +\frac{1}{2}\epsilon ff^\prime -\frac{1}{2}\epsilon f^\prime \frac{L^2}{r^2} \overset{*}{u}^2, \end{aligned}$$

(35)

where the Eq. (31) has been used. Then, by substituting (34) and (35) in (33), we arrive at the following general equation of motion in term of the metric function for the radial coordinate

$$\begin{aligned} \overset{**}{r}=\frac{1}{2}\epsilon \dot{f}\overset{*}{u}^2-\frac{1}{2}f^\prime -\frac{1}{2}f^\prime \frac{L^2}{r^2}+f\frac{L^2}{r^3}. \end{aligned}$$

(36)

Then, using the metric function $f(u,r)=1-\frac{2M(u)}{r}-\frac{N_s(u)}{{r}^{{3\omega _s +1}}}$, this equation takes the following form

$$\begin{aligned} \overset{**}{r}= & {} -\frac{M(u)}{r^2}+\frac{L^2}{r^3}-\frac{3M(u)L^{2}}{r^4}\nonumber \\&-\frac{(3\omega _s+1)N(u)}{2r^{3\omega _s+2}}-\frac{3(\omega _s+1)N(u)L^{2}}{2r^{3\omega _s+4}}\nonumber \\&+\frac{1}{2}\epsilon \dot{f}\overset{*}{u}^2. \end{aligned}$$

(37)

Then, one realizes the following three interesting points.

• 1. The terms in the first line are exactly the same as that of the standard Schwarzschild black hole in which the first term represents the Newtonian gravitational force, the second term represents a repulsive centrifugal force and the third term is the relativistic correction of the Einstein GR which accounts for the perihelion precession.

• 2. The terms in the second line are new correction terms due to the presence of the background field which surrounds the Vaidya black hole, in which its first term is similar to the term of gravitational potential in the first line, while its second term is similar to the relativistic correction of GR. Then, regarding (37) one realizes that for the more realistic non-empty backgrounds, the geodesic equation of any object depends strictly not only on the mass of the central object of the system and the conserved angular momentum of the orbiting body, but also on the background field nature. The new correction terms may be small in general in comparison to their Schwarzschild counterparts (the first and third terms in the first line). However, one can show that there are possibilities that these terms are comparable to them. Then, in order to find a situation where these forces are comparable to the Newtonian gravitational force and the GR correction term in (37), we define the distances $D_{s_{1}}$ and $D_{s_{2}}$ which satisfy $|\frac{a_{s_1}}{a_{N}}|\simeq 1$ and $|\frac{a_{s_2}}{a_{L}}|\simeq 1$, respectively, where $a_N$, $a_L$ are the Newtonian and the relativistic correction accelerations, respectively, and $a_{s_1}$ and $a_{s_2}$ are defined as

  $$\begin{aligned} a_{s_1}=\frac{(3\omega _s+1)N(u)}{2r^{3\omega _s+2}},\quad a_{s_2}=\frac{3(\omega _s+1)N(u)L^{2}}{2r^{3\omega _s+4}}. \end{aligned}$$

  (38)

  Then, the distances $D_{s_{1}}$ and $D_{s_{2}}$ will be given by

  $$\begin{aligned} D_{s_{1}}^{3\omega _s}\!=\left( \frac{|(3\omega _s \!+\!1)N_s(u)|}{2M(u)}\right) ,\,\,\, D_{s_{2}}^{3\omega _s}\!=\left( \frac{|(\omega _s \!+\!1)N_s(u)|}{2M(u)}\right) \!. \end{aligned}$$

  (39)

  We give the detailed study of these particular distances for the various cosmological backgrounds, in the Sects. 3–7.

• 3. The new correction term in the third line is also a non-Newtonian gravitational force originated from the dynamics of black hole and its surrounding field. It is associated with the radiation-accretion power of the black hole and its surrounding field.⁴ Calling this acceleration as the induced acceleration $a_{i}$, where the subscript i stands for “induced”, we have

  $$\begin{aligned} a_{i}=\frac{1}{2}\epsilon \dot{f}\overset{*}{u}^{2} =-\epsilon \left( \frac{\dot{M}(u)}{r} +\frac{\dot{N}(u)}{2r^{3\omega _s+1}}\right) \overset{*}{u}^{2}, \end{aligned}$$

  (40)

  in which, following Lindquist et al. [124], one can define the generalized “total apparent flux” as $\mathcal {A}_{F}=\epsilon \left( \dot{M}(u) +\frac{\dot{N}(u)}{2r^{3\omega _s}}\right) {\overset{*}{u}}^{2}=\mathfrak {L}+\frac{\mathfrak {N}}{2r^{3\omega _{s}}}$ where $\mathfrak {L}$ and $\mathfrak {N}$ are the apparent fluxes associated to the black hole and its surrounding field radiation-accretion rates, respectively. Using these definitions (40), takes the following form

  $$\begin{aligned} a_{i}=-\frac{\mathfrak {L}}{r}-\frac{\mathfrak {N}}{2r^{3\omega _s+1}}. \end{aligned}$$

  (41)

  As mentioned in [124], this new correction term may be small in general in comparison to the Newtonian term. However, one can show that there are possibilities that these two terms are comparable. Then, in order to find a situation where this induced force is comparable to the Newtonian gravitational force in (37), we define the distance R which satisfies $a_i\simeq a_N$, where $a_N$ is the Newtonian gravitational acceleration. Then, this distance will be given by the solutions of the following equation for different values of M, $\omega _s$ and apparent fluxes $\mathfrak {L}$ and $\mathfrak {N}$ as

  $$\begin{aligned} \mathfrak {L}R^{3\omega _s} +\frac{1}{2}\mathfrak {N}\simeq MR^{3\omega _s-1}. \end{aligned}$$

  (42)

  Finding the general solutions to this equation in terms of the generic $\mathfrak {L}, \mathfrak {N}, M$ and $\omega _s$ parameters is not simple. However, one can find that there are possible solutions for the various backgrounds of dust, radiation, quintessence, cosmological constant-like and phantom fields for some particular ranges of the parameters. We give the detailed study of this equation for the mentioned backgrounds, in the Sects. 3–7.

3 Evaporating-accreting Vaidya black hole surrounded by the dust field

3.1 Naked singularity or black hole formation analysis

For this case, the Eq. (25) takes the following form

$$\begin{aligned} (n+2m) X_{0}^{2}-X_{0}+2\epsilon =0. \end{aligned}$$

(43)

Then, one can obtain the following set of solutions to (43)

$$\begin{aligned} X_{01}= & {} \frac{1-\sqrt{1+16m+8n}}{2(2m+n)},\nonumber \\ X_{02}= & {} \frac{1+\sqrt{1+16m+8n}}{2(2m+n)},\quad \epsilon =-1,\nonumber \\ X_{01}= & {} \frac{1-\sqrt{1-16m-8n}}{2(2m+n)},\nonumber \\ X_{02}= & {} \frac{1+\sqrt{1-16m-8n}}{2(2m+n)},\quad \epsilon =+1. \end{aligned}$$

(44)

Thus, one finds that some particular conditions on the parameters m and n are required for having positive or negative solutions. In Fig. 1, we have plotted the solutions of (43) for some typical ranges of m and n parameters. Then, regarding this figure, one realizes the possibility of the formation of both naked singularities and black holes in the dust background depending on the value of parameters.

Fig. 1

The variation $X_0$ versus typical values of the m and n parameters in (43) for the dust background

3.2 Black hole-dust background field interactions

For the dust surrounding field, we set the equation of state parameter of the dust field as $\omega _d=0$ [51, 126]. Then, the metric (17) takes the following form

$$\begin{aligned} ds^{2}=-\left( 1-\frac{2M(u)+N_d(u)}{r}\right) du^2+2\epsilon dudr+r^2d\Omega ^2, \end{aligned}$$

(45)

where $N_{d}(u)$ denotes the normalization parameter for the dust field surrounding the back hole, with the dimension of $\left[ N_d \right] =l$ where l denotes the length. It is seen that the effectively radiating-accreting black hole in the dust background appears as an effectively radiating-accreting black hole with an effective mass $M_{eff}(u)=2M(u)+N_d(u)$. In this case, the presence of new mass term changes the thermodynamics, causal structure and Penrose diagrams just up to a re-scaling in the original Vaidya solution.

The radiation-accretion density in the dust background is given by

$$\begin{aligned} \sigma (u,r)=\epsilon \left( \frac{2\dot{M}(u)+\dot{N}_{d}(u)}{r^{2}}\right) . \end{aligned}$$

(46)

For the Vaidya’s original solution in an empty background, i.e $N_d(u)=0$, or even in a static background, i.e $\dot{N}_d(u)=0$, the positive energy density condition, i.e $\sigma (u,r)\ge 0$, requires that $\epsilon$ and $\dot{M}(u)$ always have the same signs. This means that for $\epsilon =+ \ 1$, M(u) is a monotone increasing mass function while for the case of $\epsilon =- \ 1$, M(u) is a monotone decreasing mass function. In our general solution for the Vaidya black hole in the dust background, the condition $\sigma (u,r)\ge 0$ imposed on (46) is satisfied for more general situations indicated in the Table 3.

Table 3 BH and its surrounding dust field parameters for $\epsilon =\pm 1$. For these cases, the positive energy condition is satisfied everywhere in spacetime

Interestingly, for the special case of $\dot{N}_{d}(u)= -\,2\dot{M}(u)$, there is no pure radiation-accretion density, i.e $\sigma (u,r)=0$, and the energy–momentum tensor (4) will be diagonalized. This means that the black hole and its surrounding background completely cancel out the effects of each others. For $\dot{N}_d(u)\ne -\,2\dot{M}(u)$, regarding (46), we find that for $r_*\rightarrow \infty$, the radiation-accretion density vanishes, i.e $\sigma (u,r)\rightarrow 0$. This means that for the effective emission case, the out going radiation can penetrate through the dust background so far from the black hole and for the effective accretion case by the black hole, the black hole affect its so far surrounding objects. Regrading the conditions in the Table 3 for $\epsilon =-1$ and $\epsilon =+1$, the behaviour of radiation-accretion density $\sigma$ in (46) is plotted for some typical values of $\dot{M}$ and $\dot{N}_{d}$ in the Fig. 2. Using these plots, one can compare the radiation-accretion density values for the various situations.

Fig. 2

Left figure: The radiation-accretion density $\sigma$ versus the distance r for some typical constant values of $\dot{M}$ and $\dot{N}_{d}$ for $\epsilon =-1$ in the dust background. The four upper cases and the four lower cases correspond to the conditions $|\dot{N}_{d}(u)|\ge |2\,\dot{M}(u)|$ and $|\dot{N}_{d}(u)|\le |2\,\dot{M}(u)|$, respectively. By these conditions, it is clear that $\sigma (r)$ is a decreasing function but is positive, and consequently the positive energy condition is satisfied everywhere in spacetime. Right figure: the radiation-accretion density $\sigma$ versus the distance r for some typical constant values of $\dot{M}$ and $\dot{N}_{d}$ for $\epsilon =+ \ 1$ in the dust background. The four upper cases and the four lower cases correspond to the conditions $|\dot{N}_{d}(u)|\le |2\,\dot{M}(u)|$ and $|\dot{N}_{d}(u)|\ge |2\,\dot{M}(u)|$, respectively. By these conditions, it is clear that $\sigma (r)$ is a decreasing function but is positive, and consequently the positive energy condition is satisfied everywhere in spacetime

3.3 Timelike geodesics for the black hole in the dust field background

For this case, we have $D_{s_1}=D_{s_2}$ and both the particular situations associated with $|\frac{a_{s_1}}{a_{N}}|\simeq 1$ and $|\frac{a_{s_2}}{a_{L}}|\simeq 1$ are met for $M(u)=\frac{|N_d(u)|}{2}$ in the whole spacetime. In Fig. 3, we have plotted the possibility of being these particular situations for some typical ranges of M(u) and $N_d(u)$ parameters. Then, one realizes the possibility of equality of the Newtonian force as well as GR correction terms to the corresponding dust background field contributions.

Fig. 3

The variation of $D_{s_1}$ and $D_{s_2}$ versus typical values of the M(u) and $N_d(u)$ parameters for the dust background

Also, for this case, the Eq. (42) associated with $a_i\simeq a_N$ takes the following form

$$\begin{aligned} \mathfrak {L} +\frac{1}{2}\mathfrak {N}\simeq MR^{-1}. \end{aligned}$$

(47)

One can find the following solution to (47)

$$\begin{aligned} R\simeq \frac{2M}{2\mathfrak {L}+\mathfrak {N}}. \end{aligned}$$

(48)

Then, one realizes that how this particular distance depends on the parameters $\mathfrak {L}, \mathfrak {N}$ and M. In the Fig. 4, we have plotted the solutions of (47) for some typical ranges of $\mathfrak {L}$ and $\mathfrak {N}$ parameters. This figure indicates that depending the parameter values, there are locations where the induced force, resulting from the radiation-accretion phenomena in the dust background, is equal to the Newtonian gravitational force.

Fig. 4

The variation of R versus typical values of the $\mathfrak {L}$ and $\mathfrak {N}$ parameters in (47) for the dust background. We have set $M=1$ without loss of generality

4 Evaporating-accreting vaidya black hole surrounded by the radiation field

4.1 Naked singularity or black hole formation analysis

For this case, the Eq. (25) takes the following form

$$\begin{aligned} nX_{0}^3 +2m X_{0}^{2}-X_{0}+2\epsilon =0. \end{aligned}$$

(49)

Then, we obtain the following solutions to (49)

$$\begin{aligned}&X_{01}=-\frac{2m}{3n}-\frac{-4m^2 -3n}{3n\Delta }+\frac{\Delta }{3n}, \nonumber \\&X_{02}=-\frac{2m}{3n}+\frac{(1+i\sqrt{3})(-4m^2 -3n)}{6n\Delta }-\frac{(1-i\sqrt{3})\Delta }{6n}, \nonumber \\&X_{03}=-\frac{2m}{3n}+\frac{(1-i\sqrt{3})(-4m^2 -3n)}{6n\Delta }-\frac{(1+i\sqrt{3})\Delta }{6n},\nonumber \\ \end{aligned}$$

(50)

where $\Delta$ is given by

$$\begin{aligned}&\Delta =\Delta _{-}=\bigg (-8m^3-9mn+27n^2 \nonumber \\&\quad +3\sqrt{3}\sqrt{-m^2n^2-16m^3n^2-n^3-18mn^3+27n^4}\bigg )^{\frac{1}{3}}, \quad \epsilon =-1,\nonumber \\&\Delta =\Delta _{+}=\bigg (-8m^3-9mn-27n^2 \nonumber \\&\quad +3\sqrt{3}\sqrt{-m^2n^2+16m^3n^2-n^3+18mn^3+27n^4}\bigg )^{\frac{1}{3}}, \quad \epsilon =+1.\nonumber \\ \end{aligned}$$

(51)

Then, one finds that some particular conditions are needed on the parameters m and n for having positive or negative solutions. In Fig. 5, we have plotted the solutions of (49) for some typical ranges of m and n parameters. This figure indicates the possibility of the formation of both the naked singularities and black holes in the radiation background depending on the value of parameters.

Fig. 5

The variation of $X_0$ versus typical values of the m and n parameters in (49) for the radiation background

4.2 Black hole-radiation background field interactions

For the radiation surrounding field, we set the equation of state parameter of the radiation field as $\omega _r=\frac{1}{3}$ [51, 126]. Then, the metric (17) takes the following form

$$\begin{aligned} ds^{2}=-\left( 1-\frac{2M(u)}{r}-\frac{N_{r}(u)}{r^2}\right) du^2+2\epsilon dudr+r^2d\Omega ^2, \end{aligned}$$

(52)

where $N_{r}(u)$ is the normalization parameter for the radiation field surrounding the black hole, with the dimension of $[N_r ]=l^2$. Regarding the positive energy condition on the surrounding radiation field, represented by the relation (16), it is required that $N_r(u)\leqslant 0$. Then, by defining the positive parameter $\mathcal {N}_r(u)=-N_r(u)$, we have

$$\begin{aligned} ds^{2}=-\left( 1-\frac{2M(u)}{r}+\frac{\mathcal {N}_{r}(u)}{r^2}\right) du^2+2\epsilon dudr+r^2d\Omega ^2. \end{aligned}$$

(53)

This metric looks like a radiating charged Vaidya black, namely the Bonnor–Vaidya black hole [37], with the dynamical charge $Q(u)=\sqrt{\mathcal {N}_r(u)}$, see also [46] for the radiating dyon solution. This result can be interpreted as the positive contribution of the characteristic feature of the surrounding radiation field to the effective charge term of the Vaidya black hole with the $\frac{1}{r^2}$ gravitational contribution. The appearance of an effective charge in the black hole solution changes the causal structure and Penrose diagrams of this black hole solution in comparison to the neutral Vaidya black holes. A similar effect in the causal structure of spacetime happens when one adds charge to the static Schwarzschild black hole leading to Reissner–Nordström black hole. Then, turning off the background radiation field which surrounds the dynamical Vaidya black hole is equal to turning off the charge in the static Reissner–Nordström case.

In this case, the total radiation-accretion density is given by

$$\begin{aligned} \sigma (u,r)=\epsilon \left( \frac{2\dot{M}(u)}{r^{2}}-\frac{\mathcal {\dot{N}}_{r}(u)}{r^3} \right) . \end{aligned}$$

(54)

Then, we see that there is no positive $r_{*}(u)$ for $\dot{M}(u)$ and $\mathcal {\dot{N}}_r(u)$ having opposite signs, and consequently $\sigma (u,r)$ never vanishes except at infinity. But as $r_*\rightarrow \infty$, the radiation-accretion density again vanishes, i.e $\sigma (u,r_*)\rightarrow 0$. This means that for the emission case, the out going radiation can penetrate through the radiation background so far from the black hole and for the accretion case by the black hole, the black hole affects its so far surrounding radiation filed. The positivity condition of $\sigma (u,r)$ is satisfied everywhere for the situations present in the Table 4.

Table 4 BH and its surrounding radiation field parameters for $\epsilon =\pm 1$. For these cases, the positive energy condition is satisfied everywhere in spacetime. For any other behaviour of the $\dot{M}(u)$ and $\mathcal {\dot{N}}_r(u)$ parameters, the positive energy condition will be violated

Regrading the Table 4, the behaviour of radiation-accretion density $\sigma$ in (54) is plotted for some typical values of $\dot{M}$ and $\mathcal {\dot{N}}_r$ in Fig. 6. Using these plots, one can compare the radiation-accretion densities for the various situations.

Fig. 6

Left figure: the radiation-accretion density $\sigma$ versus the distance r for some typical constant $\dot{M}$ and $\mathcal {\dot{N}}_{r}$ values for $\epsilon =-1$ in the radiation background. Here, $\sigma (r)$ is a decreasing function but is positive, and consequently the positive energy condition is satisfied everywhere in spacetime. Right figure: the radiation-accretion density $\sigma$ versus the distance r for some typical constant $\dot{M}$ and $\mathcal {\dot{N}}_{r}$ values for $\epsilon =+1$ in the radiation background. Here, $\sigma (r)$ is a decreasing function but is positive, and consequently the positive energy condition is satisfied everywhere in spacetime

4.3 Timelike geodesics for the black hole in the radiation field background

For this case, the distances $D_{s_{1}}$ and $D_{s_{2}}$ associated with $|\frac{a_{s_1}}{a_{N}}|\simeq 1$ and $|\frac{a_{s_2}}{a_{L}}|\simeq 1$, respectively, will be given by

$$\begin{aligned} D_{s_{1}}=\frac{|N_r(u)|}{M(u)},\quad D_{s_{2}}=\frac{2|N_r(u)|}{3M(u)}. \end{aligned}$$

(55)

In Fig. 7, we have plotted the location of these particular distances versus some typical ranges of M(u) and $N_r(u)$ parameters. Then, one realizes the possibility of equality of the Newtonian force and GR correction terms to the corresponding radiation background field contributions.

Fig. 7

The variation of $D_{s_1}$ (yellow plot) and $D_{s_2}$ (red plot) versus typical values of the M(u) and $N_r(u)$ parameters for the radiation background

Moreover, for this case, the Eq. (42) associated with $a_i\simeq a_N$ takes the following form

$$\begin{aligned} \mathfrak {L} R+\frac{1}{2}\mathfrak {N}\simeq M. \end{aligned}$$

(56)

Then, one can find the following solution to (56)

$$\begin{aligned} R\simeq \frac{2M-\mathfrak {N}}{2\mathfrak {L}}. \end{aligned}$$

(57)

It is seen that how this particular distance depends on the parameters $\mathfrak {L}, \mathfrak {N}$ and M. In Fig. 8, we have plotted the solutions of (56) for some typical ranges of $\mathfrak {L}$ and $\mathfrak {N}$ parameters. This figure shows that depending the parameter values, there are locations where the induced force, resulting from the radiation-accretion phenomena in the radiation background, is equal to the Newtonian gravitational force.

Fig. 8

The variation of R versus typical values of the $\mathfrak {L}$ and $\mathfrak {N}$ parameters in (47) for the radiation background. We have set $M=1$ without loss of generality

5 Evaporating-accreting vaidya black hole surrounded by the quintessence field

5.1 Naked singularity or black hole formation analysis

For this case, the Eq. (25) takes the following form

$$\begin{aligned} n +2m X_{0}^{2}-X_{0}+2\epsilon =0. \end{aligned}$$

(58)

Then, one can find the solutions as

$$\begin{aligned} X_{01}= & {} \frac{1-\sqrt{1+16m-8mn}}{4m},\nonumber \\ X_{02}= & {} \frac{1+\sqrt{1+16m-8mn}}{4m},\quad \epsilon =-1,\nonumber \\ X_{01}= & {} \frac{1-\sqrt{1-16m-8mn}}{4m},\nonumber \\ X_{02}= & {} \frac{1+\sqrt{1-16m-8mn}}{4m},\quad \epsilon =+1. \end{aligned}$$

(59)

Similarly, some particular conditions are required on the parameters m and n for having positive or negative solutions. In Fig. 9, we have plotted the solutions of (58) for some typical ranges of m and n parameters. Then, regarding this figure, one realizes the possibility of the formation of both naked singularities and black holes in the quintessence background depending on the value of parameters.

Fig. 9

The variation $X_0$ versus typical values of the m and n parameters in (9) for the quintessence background

5.2 Black hole-quintessence background field interactions

In the cosmological context, the quintessence filed is known as the simplest scalar field dark energy model without having theoretical problems such as Laplacian instabilities or ghosts. The energy density and the pressure profile of the quintessence filed are generally considered to vary with time and depend on the scalar field and the potential, which are given by $\rho = \frac{1}{2}\dot{\phi }^2 + V(\phi )$ and $p = \frac{1}{2}\dot{\phi }^2 -V(\phi )$, respectively. Then, the associated equation of state parameter for quintessence field lies in the range $-1<\omega _q<-\frac{1}{3}$. The static Schwarzschild black hole solution surrounded by a quintessence field was found by Kiselev [51]. This solution was generalized to the charged case and studied in [136,137,138,139].

For the quintessence surrounding field, we set the equation of state parameter of quintessence field as $\omega _q=-\frac{2}{3}$ [51, 126]. Then, the metric (17) takes the following form

$$\begin{aligned} ds^{2}=-\left( 1-\frac{2M(u)}{r}-N_{q}(u)r\right) du^2+2\epsilon dudr+r^2d\Omega ^2, \end{aligned}$$

(60)

where $N_q(u)$ is the normalization parameter for the quintessence field surrounding the black hole, with the dimension of $[N_{q}] = l^{-1}$. This result shows a non-trivial contribution of the characteristic feature of the surrounding quintessence field to the metric of the Vaidya black hole. The presence of the background quintessence filed changes the causal structure and Penrose diagrams of this black hole solution in comparison to the black hole in an empty background. A rather similar effect happens when one immerses an static Schwarzschild in a (anti)-de Sitter background with the difference that here the spacetime tends asymptotically to quintessence rather than (anti)-de Sitter asymptotics.

Table 5 BH and its surrounding quintessence field parameters for $\epsilon =\pm 1$. For the quintessence background, the positive energy condition may be completely or partially respected regarding to the above situations

Fig. 10

Left figure: the radiation-accretion density $\sigma$ versus the distance r for some typical constant $\dot{M}$ and $\dot{N}_{q}$ values for $\epsilon =-1$ in the quintessence background. In the four upper cases, the accretion density is an increasing function from negative to the positive values. In the four lower cases, the radiation density is a decreasing function decreases from positive values to negative values. Then, for a dynamical quintessence background, if the condition $|\dot{N}_q(u)|\ll |\dot{M}(u)|$ is not met, the positive energy condition is violated in some regions of spacetime. Right figure: the radiation-accretion density $\sigma$ versus the distance r for some typical constant $\dot{M}$ and $\dot{N}_{q}$ values for $\epsilon =+1$ in the quintessence background. In the upper panel, the accretion density is a decreasing function from positive to the negative values. In the lower panel, the radiation density is an increasing function from negative values to positive values. Then, for a dynamical quintessence background, if the condition $|\dot{N}_q(u)|\ll |\dot{M}(u)|$ is not met, the positive energy condition is violated in some regions of spacetime

Regarding the positive energy condition for the quintessence background, represented by the relation (16), it is required to have $N_q(u)\geqslant 0$. The radiation density is given by

$$\begin{aligned} \sigma (u,r)=\epsilon \left( \frac{2\dot{M}(u)}{r^{2}}+\dot{N}_{q}(u)\right) . \end{aligned}$$

(61)

Then, the dynamical behaviour of the background quintessence field is governed by

$$\begin{aligned} {\left\{ \begin{array}{ll}\dot{N}_{q}(u)\le -\,\frac{2}{r^2}\,\dot{M}(u), &{} \epsilon =-1, \\ \\ \dot{N}_{q}(u)\ge -\frac{2}{r^2}\,\dot{M}(u), &{}\epsilon =+1. \end{array}\right. } \end{aligned}$$

(62)

Consequently, at any distance r from the black hole, the surrounding quintessence field must obey the above conditions. Interestingly, for the special case of $\dot{N}_{q}(u)= -\frac{2\dot{M}(u)}{r^2}$, there is no pure radiation-accretion density, i.e $\sigma (u,r)=0$. This case corresponds to two possible physical situations. The first one is related to the situation where observer can be located at any distance r such that the quintessence background’s and black hole’s contributions cancel out each others leading to $\sigma (u,r)=0$ for a moment or even a period of time. Then, it is required that for an evaporating black hole, we have an equal absorbing quintessence background or for an accreting black hole we have an equal accreted quintessence background. The second situation is related to the case that for the given dynamical behaviors of the black hole and its quintessence background, one can find the particular distance $r_{*}=\sqrt{-\frac{2\dot{M}(u)}{\dot{N}_{q}(u)}}$ possessing zero energy density. For $|\dot{N}_q(u)|\ll |\dot{M}(u)|$, we have $r_*\rightarrow \infty$. This indicates that for an evolving black hole in an almost static quintessence background, the positive energy condition is satisfied everywhere. Also, the positivity of $r_*$ also requires that $\dot{M}(u)$ and $\dot{N}_q(u)$ have opposite signs. Then, if one realize the black and its surrounding quintessence filed behaviors, i.e $\dot{M}(u)$ and $\dot{N}_q(u)$ values, he can find a distance at which we have no any radiation-accretion energy density contribution. Based on these possibilities, the various situations in the Table 5 can be realized.

Then, regarding this table, the positive values of $r_*$ are physically viable and their corresponding physical processes are listed in the last column. These properties are determined according to the behaviours of the parameters $\epsilon$, $\omega _q$, and quantities $\dot{M}(u)$, $\dot{N}_q(u)$, and $\sigma (u,r)$. Those values of $r_*$ corresponding to the negative energy density $\sigma (u,r)$ represent no physical situation about the evaporation–absorption or accretion. The real features of those regions are hidden by the weak energy condition. In the reference [125], the accretion into a static Kiselev black hole with a static exterior spacetime surrounded by a quintessence field without the back-reaction effect is studied. The obtained results in [125] are implying that the accretion rate and the critical points depend on the background quintessence parameter $N_q$. Then, these features deserve to be incorporated in astrophysical studies of the accretion processes.

Regrading the Table 5, the behaviour of radiation-accretion density $\sigma$ in (61) is plotted for some typical values of $\dot{M}$ and $\dot{N}_q$ in Fig. 10. Using these plots, one can compare the radiation-accretion densities for the various situations.

5.3 Timelike geodesics for the black hole in the quintessence field background

For this case, the distances $D_{s_{1}}$ and $D_{s_{2}}$ associated with $|\frac{a_{s_1}}{a_{N}}|\simeq 1$ and $|\frac{a_{s_2}}{a_{L}}|\simeq 1$, respectively, are given as

$$\begin{aligned} D_{s_{1}}^2=\frac{2M(u)}{|-N_q(u)|},\quad D_{s_{2}}^2=\frac{6M(u)}{|N_q(u)|}. \end{aligned}$$

(63)

In Fig. 11, we have plotted the location of these particular distances versus some typical ranges of the black hole mass M(u) and background quintessence field $N_q(u)$ parameters. Then, one finds that there are possibilities for the equality of the Newtonian force and GR correction terms to the corresponding quintessence background field contributions.

Fig. 11

The variation of $D_{s_1}$ (yellow plot) and $D_{s_2}$ (red plot) versus typical values of the M(u) and $N_q(u)$ parameters for the quintessence background

The Eq. (42) associated with $a_{i}\simeq {a_{N}}$ for this case takes the following form

$$\begin{aligned} \mathfrak {L} R^{-2}+\frac{1}{2}\mathfrak {N}\simeq MR^{-3}. \end{aligned}$$

(64)

Then, we obtain the following solutions

$$\begin{aligned} R_1\simeq & {} - \frac{2\mathfrak {L}}{3^{\frac{1}{3}} \left( 9M\mathfrak {N}^{2}+\sqrt{3}\sqrt{8\mathfrak {L}^{3}\mathfrak {N}^{3}} +27M^2 \mathfrak {N}^4 \right) ^{\frac{1}{3}}}\nonumber \\&+ \frac{\left( 9M\mathfrak {N}^{2}+\sqrt{3}\sqrt{8\mathfrak {L}^{3}\mathfrak {N}^{3}} +27M^2 \mathfrak {N}^4 \right) ^{\frac{1}{3}}}{3^{\frac{2}{3}}\mathfrak {N}},\nonumber \\ R_2\simeq & {} \frac{(1+i\sqrt{3})\mathfrak {L}}{3^{\frac{1}{3}} \left( 9M\mathfrak {N}^{2}+\sqrt{3}\sqrt{8\mathfrak {L}^{3}\mathfrak {N}^{3}} +27M^2 \mathfrak {N}^4 \right) ^{\frac{1}{3}}}\nonumber \\&- \frac{(1-i\sqrt{3})\left( 9M\mathfrak {N}^{2}+\sqrt{3}\sqrt{8\mathfrak {L}^{3}\mathfrak {N}^{3}} +27M^2 \mathfrak {N}^4 \right) ^{\frac{1}{3}}}{2\times 3^{\frac{2}{3}}\mathfrak {N}},\nonumber \\ R_3\simeq & {} \frac{(1-i\sqrt{3})\mathfrak {L}}{3^{\frac{1}{3}} \left( 9M\mathfrak {N}^{2}+\sqrt{3}\sqrt{8\mathfrak {L}^{3}\mathfrak {N}^{3}} +27M^2 \mathfrak {N}^4 \right) ^{\frac{1}{3}}}\nonumber \\&- \frac{(1+i\sqrt{3})\left( 9M\mathfrak {N}^{2}+\sqrt{3}\sqrt{8\mathfrak {L}^{3}\mathfrak {N}^{3}} +27M^2 \mathfrak {N}^4 \right) ^{\frac{1}{3}}}{2\times 3^{\frac{2}{3}}\mathfrak {N}}.\nonumber \\ \end{aligned}$$

(65)

Thus, one finds that the location of this particular distance depends on the parameters $\mathfrak {L}, \mathfrak {N}$ and M. In Fig. 12, we have plotted the solutions of (64) for some typical ranges of $\mathfrak {L}$ and $\mathfrak {N}$ parameters. This figure shows that depending the parameter values, there are locations where the induced force, resulting from the radiation-accretion phenomena in the quintessence background, is equal to the Newtonian gravitational force.

Fig. 12

The variation of R versus typical values of the $\mathfrak {L}$ and $\mathfrak {N}$ parameters in (47) for the quintessence background. We have set $M=1$ without loss of generality

6 Evaporating-accreting vaidya black hole surrounded by the cosmological constant

6.1 Naked singularity or black hole formation analysis

For this case, the Eq. (25) takes the following form

$$\begin{aligned} nX_{0}^{-1} +2m X_{0}^{2}-X_{0}+2\epsilon =0. \end{aligned}$$

(66)

Then, one can find the following solutions to (66)

$$\begin{aligned}&X_{01}=\frac{1}{6m}-\frac{-1-12m}{3\times 2^{\frac{2}{3}}m\Delta }+\frac{\Delta }{6\times 2^{\frac{1}{3}}m}, \nonumber \\&X_{02}=\frac{1}{6m}+\frac{(1+i\sqrt{3})(-1-12m)}{6\times 2^{\frac{2}{3}}m\Delta }-\frac{(1-i\sqrt{3})\Delta }{12\times 2^{\frac{1}{3}}m}, \nonumber \\&X_{03}=\frac{1}{6m}+\frac{(1-i\sqrt{3})(-1-12m)}{6\times 2^{\frac{2}{3}}m\Delta }-\frac{(1+i\sqrt{3})\Delta }{12\times 2^{\frac{1}{3}}m},\nonumber \\ \end{aligned}$$

(67)

where $\Delta$ is given by

$$\begin{aligned} \Delta= & {} \Delta _{-}=\left( 2+36m-108m^{2}n\right. \nonumber \\&\left. +\sqrt{4(-1-12m)^3+(2+36m-108m^{2}n)^{2}}\right) ^{\frac{1}{3}}, \quad \epsilon =-1,\nonumber \\ \Delta= & {} \Delta _{+}=\left( 2-36m-108m^{2}n\right. \nonumber \\&\left. +\sqrt{4(-1+12m)^3+(2-36m-108m^{2}n)^{2}}\right) ^{\frac{1}{3}}, \quad \epsilon =+1.\nonumber \\ \end{aligned}$$

(68)

Then, one see that some particular conditions on the parameters m and n are required for having positive or negative solutions. In Fig. 13, we have plotted the solutions of (66) for some typical ranges of m and n parameters. Then, regarding this figure, one realizes the possibility of the formation of the both the naked singularities and black holes in the cosmological constant-like background depending on the value of parameters.

Fig. 13

The variation $X_0$ versus typical values of the m and n parameters in (66) for the cosmological constant-like background

6.2 Black hole-cosmological background field interactions

For the cosmological constant-like surrounding field, we set the equation of state parameter of the cosmological field as $\omega _c=-1$ [51, 126]. Then, the metric (17) takes the following form

$$\begin{aligned} ds^{2}=-\left( 1-\frac{2M(u)}{r}-N_c(u)r^2\right) du^2+2\epsilon dudr+r^2d\Omega ^2, \end{aligned}$$

(69)

where $N_c(u)$ is the normalization parameter for the cosmological field surrounding the black hole, with the dimension of $[N_c(u)] = l^{-2}$. This result indicates the non-trivial contribution of the characteristic feature of the surrounding cosmological constant to the metric of the Vaidya black hole. The presence of the background cosmological field changes the causal structure and Penrose diagrams of this black hole solution in comparison to the black hole in an empty background. This is similar to the case of the static Schwarzschild black hole in a static de Sitter background such that the Penrose diagram changes from Schwarzschild to Schwarzschild-(anti) de Sitter. Then, in our case, the Penrose diagram changes from Vaidya to Vaidya–de Sitter case with dynamical cosmological causal boundaries.

Regarding the positive energy condition for this case, represented by the relation (16), it is required to have $N_c(u)\geqslant 0$. In this case, $N_{c}(u)$ plays the role of a positive dynamical cosmological constant. Then, this case may describes the dynamical black holes in more general cosmological scenarios considering a time varying cosmological term, which have been recently proposed in the literature. The main purpose of these cosmological models is to provide an explanation for the recent accelerating phase of the universe [127,128,129,130,131,132,133,134,135]. These models are well known as the $\Lambda (t)$, where t is the cosmic time. For the case of $N_c=constant=\Lambda$, we recover the solution of the Vaidya black hole embedded in a de Sitter space obtained in [45]. The evolutionary behaviour of such an evaporating black hole including the structures, locations and dynamics of the apparent and event horizons are studied in [109].

In this case, the radiation-accretion density is given by

$$\begin{aligned} \sigma (u,r)=\epsilon \left( \frac{2\dot{M}(u)}{r^{2}}+\dot{N}_{c}(u)r\right) . \end{aligned}$$

(70)

Then, the dynamical behaviour of the background cosmological constant-like field is governed by

$$\begin{aligned} {\left\{ \begin{array}{ll}\dot{N}_{c}(u)\le -\frac{2}{r^3}\,\dot{M}(u), &{} \epsilon =-1, \\ \\ \dot{N}_{c}(u)\ge -\frac{2}{r^3}\,\dot{M}(u), &{}\epsilon =+1. \end{array}\right. } \end{aligned}$$

(71)

Consequently, at any distance r from the black hole, the surrounding cosmological field must obey the above conditions. Similar to the previous solution, for the special case of $\dot{N}_{c}(u)= -\frac{2\dot{M}(u)}{r^3}$, there is no pure radiation-accretion density, i.e $\sigma (u,r)=0$. This case corresponds to two possible physical situations. The first one is related to the situation where observer can be located at any distance r such that the cosmological background’s and black hole’s contributions cancel out each others leading to $\sigma (u,r)=0$ for a moment or even a period of time. Then, it is required that for a radiating black hole, we have an equal absorbing cosmological background or for an accreting black hole we have an equal accreted cosmological background field. The second situation is related to the case that for the given dynamical behaviors of the black hole and its cosmological background, one can find the particular distance $r_{*}=\left( -\frac{2\dot{M}(u)}{\dot{N}_{c}(u)}\right) ^{\frac{1}{3}}$ possessing zero energy density. For $|\dot{N}_c(u)|\ll |\dot{M}(u)|$, we have $r_*\rightarrow \infty$. This indicates that for an evolving black hole in an almost static cosmological background, the positive energy condition is respected everywhere. Here also, the positivity of $r_*$ also guarantees that $\dot{M}(u)$ and $\dot{N}_c(u)$ have opposite signs. Then, if one realize the black and its surrounding cosmological filed behaviors, i.e $\dot{M}(u)$ and $\dot{N}_c(u)$ values, he can find a distance at which we have no any radiation-accretion energy density contribution. Based on these possibilities, the various situations in the Table 6 can be realized.

Table 6 BH and its surrounding cosmological field parameters for $\epsilon =\pm 1$. For the cosmological background, the positive energy condition may be completely or partially respected regarding to the above situations

Regrading the Table 6, the behaviour of radiation-accretion density $\sigma$ in (70) is plotted for some typical values of $\dot{M}$ and $\dot{N}_c$ in Fig. 14. Using these plots, one can compare the radiation-accretion densities for the various situations.

Fig. 14

Left figure: the radiation-accretion density $\sigma$ versus the distance r for some typical constant values $\dot{M}$ and $\dot{N}_{c}$ values for $\epsilon =-1$ in the cosmological background. In the four upper cases, the accretion density is an increasing function from negative to positive values. In the four lower cases, the radiation density is a decreasing function from positive values to negative values. Then, for a dynamical cosmological background, if the condition $|\dot{N}_c(u)|\ll |\dot{M}(u)|$ is not met, the positive energy condition is violated in some regions of spacetime. Right figure: the radiation-accretion density $\sigma$ versus the distance r for some typical constant values $\dot{M}$ and $\dot{N}_{c}$ values for $\epsilon =+ \ 1$ in the cosmological background. In the four upper cases, the accretion density is a decreasing function from positive to negative values. In the four lower cases, the radiation density is an increasing function from negative to positive values. Then, for a dynamical cosmological background, if the condition $|\dot{N}_c(u)|\ll |\dot{M}(u)|$ is not met, the positive energy condition is violated in some regions of spacetime

6.3 Timelike geodesics for the black hole in the cosmological constant-like field background

For this case, the distances $D_{s_{1}}$ and $D_{s_{2}}$ associated with $|\frac{a_{s_1}}{a_{N}}|\simeq 1$ and $|\frac{a_{s_2}}{a_{L}}|\simeq 1$, respectively, read as

$$\begin{aligned} D_{s_{1}}^3=\frac{M(u)}{|-N_q(u)|},\quad D_{s_{2}}\rightarrow \infty . \end{aligned}$$

(72)

The case of $D_{s_2}\rightarrow \infty$ is resulting from the fact that, in contrast to black hole itself, the cosmological constant-like field does not couple to angular momentum L, see $a_{s_2}$ in (38). This shows that there is no similar effect to the GR correction term for the cosmological constant-like field. In Fig. 15, we have plotted the location of the particular distance $D_{s_1}$ for some typical ranges of the black hole mass M(u) and background cosmological constant-like field $N_c(u)$ parameters. Then, one finds that there are possibilities for the equality of the Newtonian force to the corresponding cosmological constant-like background field contributions.

Fig. 15

The variation of $D_{s_1}$ versus typical values of the M(u) and $N_c(u)$ parameters for the cosmological constant-like background

For this case, the Eq. (42) associated with $a_i\simeq a_N$ takes the form of

$$\begin{aligned} \mathfrak {L} R^{-3}+\frac{1}{2}\mathfrak {N}\simeq MR^{-4}. \end{aligned}$$

(73)

Then, we arrive at the following solutions

$$\begin{aligned}&R_1\simeq \frac{1}{2}\Delta ^{\frac{1}{2}}-\frac{1}{2}\sqrt{-\Delta -\frac{4\mathfrak {L}}{\mathfrak {N}\Delta ^{\frac{1}{2}}}},\nonumber \\&R_2\simeq \frac{1}{2}\Delta ^{\frac{1}{2}}+\frac{1}{2}\sqrt{-\Delta -\frac{4\mathfrak {L}}{\mathfrak {N}\Delta ^{\frac{1}{2}}}},\nonumber \\&R_3\simeq - \frac{1}{2}\Delta ^{\frac{1}{2}}-\frac{1}{2}\sqrt{-\Delta +\frac{4\mathfrak {L}}{\mathfrak {N}\Delta ^{\frac{1}{2}}}},\nonumber \\&R_4\simeq - \frac{1}{2}\Delta ^{\frac{1}{2}}+\frac{1}{2}\sqrt{-\Delta +\frac{4\mathfrak {L}}{\mathfrak {N}\Delta ^{\frac{1}{2}}}}, \end{aligned}$$

(74)

where $\Delta$ is

$$\begin{aligned} \Delta= & {} -\frac{4\times 2^{\frac{2}{3}} M}{3^{\frac{1}{3}} \left( 9\mathfrak {L}^{2}\mathfrak {N}+\sqrt{3}\sqrt{27\mathfrak {L}^{4}\mathfrak {N}^{2}} +128M^3 \mathfrak {N}^3 \right) ^{\frac{1}{3}}}\nonumber \\&+\frac{2^{\frac{1}{3}}\left( 9\mathfrak {L}^{2}\mathfrak {N}+\sqrt{3}\sqrt{27\mathfrak {L}^{4}\mathfrak {N}^{2}} +128M^3 \mathfrak {N}^3 \right) ^{\frac{1}{3}}}{3^{\frac{2}{3}}\mathfrak {N}}. \end{aligned}$$

(75)

Again, we see that how the solutions of this particular distance depends on the parameters $\mathfrak {L}, \mathfrak {N}$ and M. In Fig. 16, we have plotted the solutions of (73) for some typical ranges of $\mathfrak {L}$ and $\mathfrak {N}$ parameters. This figure shows that depending the parameter values, there are locations where the induced force, resulting from the radiation-accretion phenomena in the cosmological background, is equal to the Newtonian gravitational force.

Fig. 16

The variation of R versus typical values of the $\mathfrak {L}$ and $\mathfrak {N}$ parameters in (47) for the cosmological constant-like background. We have set $M=1$ without loss of generality

7 Evaporating-accreting Vaidya black hole surrounded by the phantom field

7.1 Naked singularity or black hole formation analysis

For this case, the Eq. (25) takes the following form

$$\begin{aligned} nX_{0}^{-2} +2m X_{0}^{2}-X_{0}+2\epsilon =0. \end{aligned}$$

(76)

Then, one can find the following solutions to this equation

$$\begin{aligned}&X_{01}=\!\frac{1}{8m}-\frac{1}{2}\Delta ^{\frac{1}{2}}-\frac{1}{2}\sqrt{\frac{3}{16m^2} +\frac{6}{3m}-\Delta -\frac{\frac{1}{8m^3}-\frac{2}{m^2}}{4\Delta ^{\frac{1}{2}}}}, \nonumber \\&X_{02}=\!\frac{1}{8m}-\frac{1}{2}\Delta ^{\frac{1}{2}}+\frac{1}{2}\sqrt{\frac{3}{16m^2} +\frac{6}{3m}-\Delta -\frac{\frac{1}{8m^3}-\frac{2}{m^2}}{4\Delta ^{\frac{1}{2}}}}, \nonumber \\&X_{03}=\!\frac{1}{8m}+\frac{1}{2}\Delta ^{\frac{1}{2}}-\frac{1}{2}\sqrt{\frac{3}{16m^2} +\frac{6}{3m}-\Delta +\frac{\frac{1}{8m^3}-\frac{2}{m^2}}{4\Delta ^{\frac{1}{2}}}}, \nonumber \\&X_{04}=\!\frac{1}{8m}+\frac{1}{2}\Delta ^{\frac{1}{2}}+\frac{1}{2}\sqrt{\frac{3}{16m^2} +\frac{6}{3m}-\Delta +\frac{\frac{1}{8m^3}-\frac{2}{m^2}}{4\Delta ^{\frac{1}{2}}}},\nonumber \\ \end{aligned}$$

(77)

where $\Delta$ is given by

$$\begin{aligned} \Delta= & {} \Delta _{-}=\frac{1}{16 m^2}+\frac{2}{3 m}\nonumber \\&+\frac{2\times 2^{\frac{1}{3}} (1+6 m n)}{3 m \left( -16\!+\!27n\!+\!288 m n\!+\!\sqrt{(-16\!+\!27 n\!+\!288 m n)^2-\!4 (4\!+\!24 m n)^3}\right) ^{\frac{1}{3}}} \nonumber \\&+\frac{\left( \!-16\!+\!27n\!+\!288 m n\!+\!\sqrt{(16\!+\!27 n\!+\!288 m n)^2\!-\!4 (4\!+\!24 m n)^3}\right) ^{\frac{1}{3}}}{6\times 2^{\frac{1}{3}}m}, \,\, \epsilon =-1,\nonumber \\ \Delta= & {} \Delta _{+}=\frac{1}{16 m^2}-\frac{2}{3 m}\nonumber \\&+\frac{2\times 2^{\frac{1}{3}} (1+6 m n)}{3 m \left( 16\!+\!27n\!-\!288 m n\!+\!\sqrt{(16\!+\!27n\!-\!288 m n)^2\!-\!4 (4\!+\!24 m n)^3}\right) ^{\frac{1}{3}}} \nonumber \\&+\frac{\left( 16\!+\!27n\!-\!288 m n\!+\!\sqrt{(16\!+\!27n\!-\!288 m n)^2\!-\!4 (4\!+\!24 m n)^3}\right) ^{\frac{1}{3}}}{6\times 2^{\frac{1}{3}}m}, \quad \epsilon =+1.\nonumber \\ \end{aligned}$$

(78)

Then, similar to the previous cases, some particular conditions on the parameters m and n are required for having positive or negative solutions. In Fig. 17, we have plotted the solutions of (76) for some typical ranges of m and n parameters. Regarding this figure, we find the possibility of the formation of the both the naked singularities and black holes in the phantom background depending on the value of parameters.

Fig. 17

The variation $X_0$ versus typical values of the m and n parameters in (76) for the phantom background

7.2 Black hole-phantom background field interactions

For the phantom surrounding field, we set the equation of state parameter of phantom field as $\omega _p=-\frac{4}{3}$ [126]. Then, the metric (17) takes the following form

$$\begin{aligned} ds^{2}=-\left( 1-\frac{2M(u)}{r}-N_p(u)r^3\right) du^2+2\epsilon dudr+r^2d\Omega ^2, \end{aligned}$$

(79)

where $N_p(u)$ is the normalization parameter for the phantom field surrounding the black hole, with the dimension of $[N_{p}(u)] = l^{-3}$. Similarly, this result is interpreted as the non-trivial contribution of the characteristic feature of the surrounding phantom field to the metric of the Vaidya black hole. The presence of the background phantom filed changes the causal structure and Penrose diagrams of this black hole solution in comparison to the Vaidya black hole in an empty background.

Regarding the weak energy condition for this case, represented by the relation (16), it is required to have $N_p(u)\geqslant 0$. In this case, the radiation-accretion density is given by

$$\begin{aligned} \sigma (u,r)=\epsilon \left( \frac{2\dot{M}(u)}{r^{2}}+\dot{N}_{p}(u)r^{2}\right) . \end{aligned}$$

(80)

Then, the dynamical behaviour of the background field is governed by

$$\begin{aligned} {\left\{ \begin{array}{ll}\dot{N}_{p}(u)\le -\frac{2}{r^4}\,\dot{M}(u), &{} \epsilon =-1, \\ \\ \dot{N}_{p}(u)\ge -\frac{2}{r^4}\dot{M}(u), &{}\epsilon =+1. \end{array}\right. } \end{aligned}$$

(81)

Consequently, at any distance r from the black hole, the surrounding phantom field must obey the above conditions. Similar to the previous solutions, for the special case of $\dot{N}_{p}(u)= -\frac{2\dot{M}(u)}{r^4}$, there is no pure radiation-accretion density, i.e $\sigma (u,r)=0$. This case corresponds to two possible physical situations. The first one is related to the situation where observer can be located at any distance r such that the phantom background’s and black hole’s contributions cancel out each others leading to $\sigma (u,r)=0$ for a moment or even a period of time. Then, it is required that for a radiating black hole, we have an equal absorbing phantom background or for an accreting black hole we have an equal accreted phantom background. The second situation is related to the case that for the given dynamical behaviors of the black hole and its phantom background, one can find the particular distance $r_{*}=\left( -\frac{2\dot{M}(u)}{\dot{N}_{p}(u)}\right) ^{\frac{1}{4}}$ possessing zero energy density. Similarly, for $|\dot{N}_p(u)|\ll |\dot{M}(u)|$, we have $r_*\rightarrow \infty$. This indicates that for an evolving black hole in an almost static phantom background, the positive energy condition is satisfied everywhere. Also, the positivity of $r_*$ also requires that $\dot{M}(u)$ and $\dot{N}_p(u)$ must have opposite signs. Then, if one realize the black and its surrounding phantom filed behaviors, i.e $\dot{M}(u)$ and $\dot{N}_p(u)$ values, he can find a distance at which we have no any radiation-accretion energy density contribution. Based on these possibilities, the various situations in the Table 7 can be realized.

Table 7 BH and its surrounding phantom field parameters for $\epsilon =\pm 1$. For the phantom background, the positive energy condition may be completely or partially respected regarding to the above situations

Fig. 18

Left figure: the radiation-accretion density $\sigma$ versus the distance r for some typical constant $\dot{M}$ and $\dot{N}_{p}$ values for $\epsilon =- \ 1$ in the phantom background. In the four upper cases, the accretion density is an increasing function from negative to positive values. In the four lower cases, the radiation density is a decreasing function from positive to negative values. Then, for a dynamical phantom background, if the condition $|\dot{N}_p(u)|\ll |\dot{M}(u)|$ is not met, the positive energy condition is violated in some regions of spacetime. Right figure: the radiation-accretion density $\sigma$ versus the distance r for some typical constant $\dot{M}$ and $\dot{N}_{p}$ values for $\epsilon =+ \ 1$ in the phantom background. In the four upper cases, the accretion density is a decreasing function from positive to negative values. In the four lower cases, the radiation density is an increasing function from negative to positive values. Then, for a dynamical phantom background, if the condition $|\dot{N}_p(u)|\ll |\dot{M}(u)|$ is not met, the positive energy condition is violated in some regions of spacetime

Specific scenarios involving the accretion of phantom energy and resulting in the area decrease of black hole [96,97,98,99] are related to the first case in the above table. For example, in [96], it is shown that the black holes will gradually vanish as the universe approaches a cosmological big rip state with a phantom field. One should note that the astrophysically “observed” infall of quintessence/phantom fields onto black holes are not detected till now. In the present work, we have just introduced a new dynamical solution to the Einstein field equations which can provide a classical model for the “possible” black hole accretion and evaporation (presumably very tiny) in different cosmological backgrounds. Such theoretical studies of the accretion of exotic fields to the black holes are well motivated by cosmology in which these exotic fields can be responsible for the current accelerating expansion of the universe. Regrading the Table 7, the behavior of radiation-accretion density $\sigma$ in (80) is plotted for some typical values of $\dot{M}$ and $\dot{N}_p$ in Fig. 18. Using these plots, one can compare the radiation-accretion densities for the various situations.

7.3 Timelike geodesics for the black hole in the phantom field background

For this case, the distances $D_{s_{1}}$ and $D_{s_{2}}$ associated with $|\frac{a_{s_1}}{a_{N}}|\simeq 1$ and $|\frac{a_{s_2}}{a_{L}}|\simeq 1$, respectively, are given as

$$\begin{aligned} D_{s_{1}}^4=\frac{2M(u)}{3|-N_p(u)|}, \quad D_{s_{2}}^4=\frac{6M(u)}{|-N_p(u)|}. \end{aligned}$$

(82)

In Fig. 19, we have plotted the location of these particular distances for some typical ranges of the black hole mass M(u) and background phantom field $N_p(u)$ parameters. Then, one realizes the possibilities of the equality of the Newtonian force and GR correction terms to the corresponding phantom background field contributions.

Moreover, the Eq. (42) for this case takes the following form

$$\begin{aligned} \mathfrak {L} R^{-4}+\frac{1}{2}\mathfrak {N}\simeq MR^{-5}. \end{aligned}$$

(83)

Then, we see that this produces a fifth order equation in which finding its analytical solutions is not simple. However, in Fig. 20, we have shown that there are numerical solutions to (83) for some typical ranges of $\mathfrak {L}$ and $\mathfrak {N}$ parameters. This figure indicates that depending the parameter values, there are locations where the induced force, resulting from the radiation-accretion phenomena in the phantom background, is equal to the Newtonian gravitational force.

Fig. 19

The variation of $D_{s_1}$ (yellow plot) and $D_{s_2}$ (red plot) versus typical values of the M(u) and $N_p(u)$ parameters for the phantom background

Fig. 20

The variation of R versus typical values of the $\mathfrak {L}$ and $\mathfrak {N}$ parameters in (83) for the phantom background. We have set $M=1$ without loss of generality

8 Conclusion

In this work, we have studied the general surrounded Vaidya solution with the cosmological fields of dust, radiation, quintessence, cosmological constant-like and phantom, and investigated its nature describing the possibility of the formation of naked singularities or black holes. We have obtained the general equation describing the nature of the solution under a collapse, and have shown that depending on the parameter values, the formation of both naked singularity and black hole as the end state of the collapse are possible. We have given the corresponding analytical solutions as well as some plots indicating these possibilities. Then, motivated by the fact that real astrophysical black holes as non-stationary and non-isolated objects are living in non-empty backgrounds, we have focused on the black hole subclasses of the obtained general solution, namely the “surrounded Vaidya black hole”, describing a dynamical evaporating-accreting black holes in the mentioned dynamical cosmological backgrounds. In the following, we summarize some of our obtained results for this solution.

• Some of the subclasses of the obtained general solution for both the dynamical and stationary limits have been addressed. In particular, we have shown that the original Vaidya solution can be recovered by turning off the background field, and that the Kiselev static solution can be obtained in the stationary limit with an appropriate coordinate transformation. Also, the Schwarzschild solution can be obtained in the stationary limit with a turned off background.

• We have shown that for the background field possessing $\omega _s>0$, if $\dot{M}(u)$ and $\dot{N}_s(u)$ have a same order of magnitude, the surrounding background field contribution to the total density $\sigma (u,r)$ is dominant near the black hole while at far distances from the black hole it decreases faster than the contribution of the black hole mass changing term. In contrast, for the background field possessing $\omega _s<0$, the surrounding background field contribution is dominant at large distances while the black hole contribution is dominant near the black hole itself. Then, from astrophysical point of view, the detected amount of the radiation-accretion density by the observer not only depends on the distance from the black but also depends on the nature of background field.

• We have discussed that positive energy condition for the surrounding field is met by the constraint $\omega _s N_{s}(u)\le 0$, which determines the gravitational nature of the term associated to surrounding field in the metric function f(u, r). The positive energy condition for the radiation-accretion density is met by the constraints $\dot{N}_{s}(u)\le -\,2\,{r}^{3\omega _s }\,\dot{M}(u)$ and $\dot{N}_{s}(u)\ge -\,2\,{r}^{3\omega _s }\,\dot{M}(u)$ at any distance r for $\epsilon =-\ 1$ and $\epsilon =+ \ 1$, respectively. One astrophysical importance of such physical constraints is that the observer knows the dynamical range of the background field at any distance and then prior to any observation, he knows how to include or remove the background field contribution, if he is only interested in black hole’s contributions, or vice vera.

• We have addressed the solutions of the black hole in the dust ($\omega _s=0$), radiation ($\omega _s=\frac{1}{3}$), quintessence ($\omega _s=-\frac{2}{3}$), cosmological constant-like ($\omega _s=-1$) and phantom ($\omega _s=-\frac{4}{3}$) fields in detail. We have found that the effectively evaporating-accreting black hole in the dust background appears as an effectively evaporating-accreting black hole with an effective mass $M_{eff}(u) = 2M(u)+N_{d}(u)$. Then, the presence of new mass term changes the causal structure just up to a re-scaling in the original Vaidya solution. For the radiation background, the spacetime metric looks like the Bonnor–Vaidya and radiating dyon solutions with the dynamical charge $Q(u) = \sqrt{\mathcal {N}_r(u)}$. A similar effect in the causal structure of spacetime here happens when one adds charge to the static Schwarzschild black hole leading to Reissner–Nordström black hole. For the black hole in the dust and radiation backgrounds, the spacetime metrics are asymptotically flat while for the the quintessence, cosmological-like and phantom backgrounds, spacetime metrics are asymptotically non-flat quintessence, de Sitter-like and phantom, respectively. Consequently, the causal structure of these three latter spacetimes are quite different from the original Vaidya spacetime where the background is turned off.

• Regarding the obtained radiation-accretion density $\sigma (u,r)$, we have found that there are particular distances $r_*$ where $\sigma (u,r)$ vanishes, i.e $\sigma (u,r_*)=0$. These distances are given by $r_{*}(u)=\left( -\frac{\dot{N}_{s}(u)}{2\dot{M}(u)}\right) ^{\frac{1}{3\omega _{s}}}$ and $r_*=\infty$. Then, one realizes that in the first case, for (i) $-\frac{2}{3}<\omega _s<0$ with $|\dot{N}_s(u)|\ll 2|\dot{M}(u)|$ and for (ii) $\omega _s\ge 0$ with $2|\dot{M}(u)|\ll |\dot{N}_{s}(u)|$, we have $r_*\rightarrow \infty$. This means that for (i), the black hole evolves very faster than its background while for (ii), the black hole evolves very slow relative to its background. Also, for (iii) $\omega _s\le -\frac{2}{3}$ with $|\dot{N}_s(u)|\ll 2|\dot{M}(u)|$, representing a rapidly evolving black hole relative to its background, we have $r_*\rightarrow \infty$. Then, by satisfaction of these dynamical conditions to hold $r_*\rightarrow \infty$, the positive energy density is respected everywhere in the spacetime. The case (ii) includes the black hole surrounded by the rapidly evolving dust and radiation fields while the cases (i) and (iii) imply an evolving black hole in an almost static cosmological backgrounds (quintessence, cosmological constant-like or phantom fields) responsible for the accelerating expansion of the universe. In practice, an astrophysical observer detects a radiation-accretion density resulting from the interaction of the black hole with its surrounding field even at far distances, in which for (ii) the main contribution in the detected radiation-accretion density belongs to surrounding field while for (i) and (iii), it is the black hole which has the main contribution in the radiation-accretion density.

• In the case that there is a real, positive and finite value for $r_*$, the positive energy condition is violated in some regions of spacetime such that real features of those regions are hidden by the weak energy condition. Then, it is physically reasonable to do any astrophysical experiments in the regions respecting the energy condition. For the cases in which $r_*$ is not positive and real, the interpretation is as follows: the radiation-accretion density $\sigma (u,r)$ never and nowhere vanishes.

• We have classified the possible situations respecting or violating the energy condition for all the solutions of black hole in dust, radiation, quintessence, cosmological constant-like and phantom backgrounds in the Tables 1, 2, 3, 4, 5, 6 and 7. It is shown that there are cases for all the backgrounds in which the positive energy condition is respected in whole spacetime under the determined behaviors of black hole and its surrounding field. Also, we have given some plots for radiation-accretion density versus some typical values of black hole and its surrounding fields in Figs. 1, 2, 3, 4 and 5. Using these plots, one realizes for the dust and radiation backgrounds, although the radiation-accretion density is a decreasing function but is always positive, and consequently the positive energy condition is satisfied everywhere in spacetime. This is while for the quintessence, cosmological constant-like and phantom backgrounds, if the condition $|\dot{N}_{q,c,p}(u)|\ll |\dot{M}(u)|$ is not met, the positive energy condition is violated in some regions of spacetime. Comparing the plots with common values of the parameters, we observe that the radiation-accretion density $\sigma (r)$ for the radiation background is larger than the dust background at any distance r, i.e $\sigma _{r}(r)>\sigma _{d}(r)$. Similarly, for the quintessence, cosmological constant-like and phantom backgrounds, we have $\sigma _q(r)<\sigma _c(r)<\sigma _p(r)$ for the radiation-accretion density.

• We have analyzed the timelike geodesics associated with the obtained surrounded black holes and have found that two kinds of new correction terms arise relative to the case of Schwarzschild black hole. The first kind of corrections are due to the presence of the background fields which surround the Vaidya black hole. This corrections include two terms in which its first term is similar to the term of Newtonian gravitational potential, while its second term is similar to the relativistic correction of GR. For the various background fields, we have discussed that there are possibilities for the equality of Newtonian and GR correction terms to the corresponding background fields contributions. We have given some plots denoting these possibilities for each case. The second kind of corrections is also a non-Newtonian correction resulting from the dynamics of black hole and its surrounding field. We have shown that depending on the dynamical features of black hole and its background, there are also possibilities that dynamical correction terms can be equal to the Newtonian case. Some plots representing these situations are given for each case. Then, one realizes that for the more realistic non-empty and non-static backgrounds, the geodesic equation of any object depends strictly not only on the mass of the central object of the system and the angular momentum of the orbiting body, but also on the (i) background field type and (ii) black hole and its background field dynamics.

We have reported elsewhere on the causal structures and thermodynamical properties of our obtained solutions here [140]. We also aim to generalize this work to the case of charged black holes.