{'fraction_non_alphanumeric': 0.0554399762236803, 'fraction_numerical': 0.022747536636108, 'mean_word_length': 4.03649971214738, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New Classes of Spherically Symmetric, Inhomogeneous Cosmological Models

We present two classes of inhomogeneous, spherically symmetric solutions of the Einstein-Maxwell-Perfect Fluid field equations with cosmological constant generalizing the Vaidya-Shah solution. Some special limits of our solution reduce to the known inhomogeneous charged perfect fluid solutions of the Einstein field equations and under some other limits we obtain new charged and uncharged solutions with cosmological constant. Uncharged solutions in particular represent cosmological models where the universe may undergo a topology change and in between is a mixture of two different Friedmann-Robertson-Walker universes with different spatial curvatures. We show that there exist some spacelike surfaces where the Ricci scalar and pressure of the fluid diverge but the mass density of the fluid distribution remains finite. Such spacelike surfaces are known as (sudden) cosmological singularities. We study the behavior of our new solutions in their general form as the radial distance goes to zero and infinity. Finally, we briefly address the null geodesics and apparent horizons associated to the obtained solutions.


I. Introduction
In the last two decades there is an increasing interest in studying and finding exact inhomogeneous cosmological solutions in general relativity. Observational effects of inhomogeneity in cosmology are discussed in several works. Among these we note that the collection of articles in [1]- [4] are worth mentioning.
There are many reasons to study inhomogeneous cosmological models in general relativity. Among these, the following three mentioned by Ellis [1] (see also the references therein) are important. Local inhomogeneity may effect the averaged large scale dynamics of the universe (see also [2] and the references therein ), local inhomogeneity may effect the photon propagation hence may change the cosmological observations and the inhomogeneity at Hubble scale with the violation of Copernicus principle may lead to acceleration of the universe (see also [3] and the references therein). In his book [5] Krasinski gives other reasons such as the formation of voids and interaction of the cosmic microwave background radiation with matter in the universe can be explained by exact solutions of the Einstein field equations in an inhomogeneous spacetime. For all these reasons it worths finding new inhomogeneous solutions of Einstein's field equations.
Spherically symmetric cosmological models were studied previously by many authors [4]- [18]. Historically, Lemaitre [19,20] and McVittie metrics [21] can be considered as the first inhomogeneous solutions of the Einstein-perfect fluid field equations. Recently, it has been shown that the McVittie solution represents a black hole in an expanding universe [22,23]. Charged version of the McVittie solution is known as the Vaidya-Shah metric [24]- [26] which is a spherically symmetric solution of the Einstein-Maxwell-perfect fluid field equations. This metric of this solution is given as follows ds 2 = −A 2 dt 2 + B 2 (dr 2 + r 2 dθ 2 + r 2 sin 2 (θ) dφ 2 ), where B = a(t) 1 + kr 2 where a(t) is any arbitrary function of time t, M and Q are constants representing the conserved quantities of mass and charge, and k is also a constant. Pressure, mass and charge densities are respectively given by The uncharged (Q = 0) Vaidya-Shah solution is more general than the McVittie solution. McVittie solution corresponds to Q = 0 and k = 0. In spite of this fact, Vaidya-Shah solution is sometimes named as the charged McVittie solution. The Vaidya-Shah metric reduces to the Reissner-Nordström metric when k = 0 and a(t) = 1 in isotropic coordinates. Note that the charge density (6) for Vaidya-Shah solution vanishes as k = 0 but the Maxwell's field F 01 remains nonzero. Vaidya-Shah solution [26] has been studied by several authors [6]- [9], [23], and like the McVittie solution it has been shown that it describes a charged black hole in an expanding universe. The charged and uncharged cosmological black holes were also discussed in the works [10]- [14].
In this work, we start with the spherically symmetric metric in the isotropic coordinates in four dimensions ds 2 = −a 2 dt 2 + b 2 dr 2 + r 2 dθ 2 + r 2 sin 2 (θ) dφ 2 , where a and b are differentiable functions of t and r. We first show that the Einstein Maxwell-Perfect Fluid field equations with cosmological constant reduce to a single nonlinear ordinary differential equation for the function b(t, r) (Theorem 1). Then we solve this differential equation as general as possible. We use the method of separation of variables and find two distinct classes of solutions (Theorem 2). For the charged case, we have the following distinct solutions: and class 2: where a(t) and β(t) are arbitrary functions of t, ν 0 (r) is an arbitrary function of r while b 0 , b 1 , c 0 , c 1 , c 2 , c 3 , δ and γ are arbitrary constants. For the uncharged case, the above solutions reduce to the following distinct solutions: and where b 2 is also an arbitrary constant. For all the above cases, we found a(t, r) = q(t)˙b b . We show that, in particular for the uncharged case, the first class of solutions exhibit a cosmological model describing a universe as a mixture of two different Friedmann-Robertson-Walker universes with different spatial curvatures. If the signs of the spatial curvatures are different then we show that there is a possibility of the change of topology of the universe. If the spatial curvatures turns out to be the same the spacetime becomes a single FRW universe. We then study the asymptotical properties of our solutions. We show that the six parameter solution which is the generalization of the Vaidya-Shah solution (1)- (6) is nonsingular as the radial distance goes to zero and to infinity (Theorem 3). The uncharged limit (Q = 0) of our solutions generalize the McVittie solution. We show that there are surfaces Σ 1 (b(t, r) = 0) and Σ 2 (a(t, r) = 0) where the Ricci scalar diverges (spacetime singularities). Σ 1 is a timelike but Σ 2 is a spacelike surface. Physical constraints eliminate the timelike surfaces Σ 1 and there remain only the spacelike singular surfaces Σ 2 . This surface is commonly named as the cosmological singularity [27] where the mass density is regular but the pressure diverges on this surface. This surface is also called a "sudden cosmological singularity" [28]- [30]. We also obtain the apparent horizons of our solutions which correspond to null (constant) areal distance surfaces. We give a plot of null geodesics, apparent horizons and singular surface Σ 2 for the N = 2 uncharged solution for particular values of the parameters of the solution. The layout of the paper is as follows. In Section II, we simplify and reduce the field equations into a single ordinary nonlinear differential equation. In Section III, we solve the resulting differential equation by the use of the method of separation of variables and obtain two different distinct solutions. In Section IV, we obtain the asymptotic behaviours of our solutions and show that the corresponding spacetimes are nonsingular with respect to the asymptotic values of r. In Section V, we study all possible special limits of our solutions. In Section VI, we study the uncharged versions of our solutions. In Section VII, we investigate the possible apparent horizons and null geodesics of the charged and uncharged solutions. In Appendix A, we write the differential equation obtained in Section II in a different form and in Appendixes B-F we give the long expressions obtained in Sections IV and V. In the last Appendix we give the mass densities when a → 0 and a → ∞ respectively.

II. Field Equations of the Charged Fluids in Four Dimensions
We consider the Einstein-(anti) de Sitter-Maxwell-Perfect Fluid field equations where where Λ, T µν , E µν , and F µν are the cosmological constant, energy-momentum tensor of the perfect fluid, Maxwell and Faraday tensors, respectively. To obtain our solutions, we consider the spherical symmetric metric where a and b are generic functions of both the time t and radial coordinate r, i.e a = a(t, r) and b = b(t, r).
Regarding the spherical symmetry in the spacetime metric (16), the only non-vanishing component of the antisymmetric electromagnetic Faraday tensor is where ψ = ψ(t, r). Using the non-zero source Maxwell equation (15) and the metric (16), we obtaiṅ and 4πσ = 1 where σ = σ(t, r) is the charge density and the dot and prime signs denote the derivatives with respect to time and radial coordinates, respectively. On the other hand, using (13), (14) and (17) and considering the perfect fluid velocity vector as u µ = a δ 0 µ , the 00 component of the Einstein-Maxwell-Perfect Fluid equations (12) gives The 01 components reads asḃ while the 11 and 22 (or 33) components lead to and respectively. We can integrate the equation (21) to obtain where q = q(t) andḃ(t, r) = 0. One notes that forḃ(t, r) = 0, the equation (21) disappears. Using (18), we arrive at where h = h(r) is an arbitrary functions of r. Using (24) and (25), the charge density σ in (19) takes the following form Then, the total charge Q T in a spherical region with radius R 0 (t constant, r constant regions) can be obtained as Hence, the total charge in this volume is given by Finally, the equations (22) and (23) reduce to One can integrate the differential equation (29) with respect time and obtain the following second order ordinary nonlinear differential equation for b where To summarize what we have till now, we introduce the following theorem.
Theorem 1: Einstein field equations of a charged perfect fluid with a cosmological constant of a spherically symmetric spacetime reduce to the following subclasses.
(i) Forḃ(t, r) = 0, the field equations reduce to a single ordinary nonlinear differential equation, Eq (30), with two arbitrary functions of r, h and h 1 functions. Then, the metric function a(t, r) and the charge density σ(t, r) are given by (24) and (26) respectively, and the energy density ρ(t, r) in (20) and the pressure p(t, r) in (22) (or (23)) respectively read as and (ii) Forḃ(t, r) = 0, there is no 01 component for the field equations, then the equation (21) and the relation between the metric functions as (24) disappears. For this case, the Maxwell equation (15) gives ψ(t, r) = h(r)a(t, r), and the equations (20), (22) reduce to where b(r) should satisfy the following equation Solving this single differential equation with three unknown functions h(r), a(t, r) and b(r) is not possible except by supposing relations between these functions. One possible ansatz can be considering a specific equation of state for the perfect fluid, leading to a relation between a(t, r) and b(r) functions.
In this work, we consider only the general dynamical case, i.eḃ(t, r) = 0, and then our aim is to solve the nonlinear ordinary differential equation (30) for the metric function b(t, r). In the next sections, we will solve this equation and determine all of our unknown functions a(t, r), ψ(t, r), σ(t, r), ρ(t, r) and p(t, r) accordingly.

III. Exact Solutions of the Field Equations
The main aim of this section is to find solutions of the equation (30). For this purpose, we use the method of separation of variables. Although the equation (30) is a nonlinear ordinary differential equation, we can use this method by equating the coefficients of the products of the time dependent functions to zero. Let where α n (r) and β n (t) are all independent functions of r and t, respectively such that n = 0, 1, 2, · · · , N . There are N +1 number of functions α n (r) depending on r in (36) and 2 arbitrary functions h 1 (r) and h(r) in the main equation (30). Then, totally we have N + 3 functions of r. The functions β n (t) (n = 0, 1, 2, · · · , N ) are left arbitrary but independent functions of t. The time independent term 2rh 2 in the main equation (30) forces us to choose one of the time dependent functions β n (t) (n = 0, 1, 2, · · · , N ) to be a constant. Thus, without losing any generality, we let β 0 = 1. Hence, we have By inserting (37) in (30), we obtain more than 2N + 1 equations. This means that when N > 2, the number of equations becomes more than the number of unknown functions (an overdetermined system). Hence, we use the ansatz (36) only for N = 2 and for N = 1, and we investigate these cases in detail in the sections III A and III B. Before to proceed, we refer the reader to Appendix A summarizing the method introduced in [31,32] for solving the equation (30) for the uncharged case, where some particular solutions are also addressed. To produce the most generic solutions including the charge, the approach in [31,32] seems not suitable for us and we will follow the method of separation of variables as discussed above.

A. Solutions for N = 2
Letting N = 2, we have b(t, r) = α 0 (r) + β 1 (t)α 1 (r) + β 2 (t)α 2 (r), where as mentioned before, α 0 , α 1 and α 2 are functions of r and β 1 and β 2 are functions of t. In Eq.(30), when the function b in (38) is inserted, the coefficients of the time dependent functions β 2 1 , β 2 2 , β 1 and β 2 are set to zero and functions α 0 , α 1 and α 2 satisfy the following equations The remaining equation depends how the functions β 1 and β 2 are related which reads as General solutions of (39) and (40) are given by where c 0 , c 1 , c 2 and c 3 are arbitrary constants. In (41) and (42), the function α 0 satisfies a second order linear differential equation. Multiplying (41) by α 2 and (42) by α 1 , and subtracting them, we obtain α 0 as where c 4 is an arbitrary constant, and we have the condition c 4 (c 3 c 0 − c 1 c 2 ) > 0. As we will see in the classification of the possible solutions, the negative sign of δ is not physical due to its identification relation to mass. The equation (42) can be considered as the definition of the function h 1 . Hence, we have solved all equations (39)-(42). There remains only h function to be determined. For determining function h, there are two possibilities as follows.
(i) If β 1 and β 2 have no relations (if β 1 β 2 = constant). For this case, using (43), we have Now, the equation (47) can be considered as the definition of the function h but the last equation (48) gives c 1 = c 3 = 0 which means that the function b depends only on t which is not our desired solution in general.
For this case, we have the following single differential equation This equation can be considered as the definition of the function h. Thus, by this consideration, we can solve the equation (30) completely. Then, the function b(t, r) takes the form of (50) B. Solutions for N = 1 where here ν 0 (r) and ν 1 (r) are functions of r and β(t) is a function of t. By inserting the function b(t, r) in (51) in the equation (30), ν 0 (r) and ν 1 (r) should satisfy the following equations General solution of (52) is where b 0 and b 1 are arbitrary constants. Equations (53) and (54) can be considered as the definitions of the functions h and h 1 . Hence, ν 0 (r) function is left arbitrary. Then, b(t, r) takes the following form Thus, the following theorem represents the summary of what is done till now.

Theorem 2:
The most general solutions of the ordinary nonlinear differential equation (30) by the method of separation of variables are given in two classes: The first one containing one arbitrary function of t and six arbitrary parameters is corresponding to N = 2 and the second one containing two arbitrary constants and two arbitrary functions where one depends on r and the other depends on t corresponding to N = 1.
A different approach is given in Appendix I for solving (30). Such an approach was introduced in [31] for the uncharged case (see also [4,32]).

IV. Properties of the Solutions to the Field Equations
In this section, we first investigate singular structure of the obtained spacetimes. There are surfaces Σ 1 and Σ 2 where the pressure p and mass density ρ diverge. Then, we explicitly check the properties of the general solutions for both the cases of N = 2 and N = 1 as r goes to zero and tends to infinity, in detail. Furthermore, we will address some specific subclasses of these general solutions and study their properties also in the next sections.

A. Singular Structure of the Solutions
Here, we assume that c 0 , c 1 , c 2 and c 3 are non negative constants. Regarding the field equations (12), the scalar curvature (Ricci scalar) is given by R = 8π(ρ − 3p) + 4Λ. Hence, if any one of the quantities p or ρ is singular on some surfaces then they are the spacetime singularities. Regarding (24), (31) and (32), if the functions b and a vanish on some surfaces then either p or ρ diverges. Hence, we will focus on the surfaces In the following, we will explore these singular surfaces.

Singular Surfaces for the Class of
Letting X = c0+c1r 2 c2+c3r 2 1 2 , then b(t, r) = 0 leads to the following equation When γ = 0 this equation has real solutions only when δ 2 − 4γ ≥ 0. Then, there are two different dynamical surfaces given by (depending on the sign of δ) When γ = 0 we have The normal vectors of these surfaces satisfy Thus, near the Σ 1 surface, it is clear that g µν ∂ µ b ∂ ν b ≥ 0. Hence, Σ 1 surfaces are timelike or null. The case δ 2 − 4γ = 0, representing only one singular dynamical surface, corresponds to the uncharged solutions which will be discussed in Section VI. For physical spacetimes both δ and γ are positive. Hence in such cases Σ 1 surface does not exist.
then Σ 2 is defined as The normal vector of this surface satisfies representing that Σ 2 is a spacelike surface or null, since a = 0 then g µν ∂ µ a∂ ν a ≤ 0 near Σ 2 . Such singularities are named as the cosmological singularities [27] or sudden cosmological singularities [28]- [30].

Singular Surfaces for the Class of N = 1
Regarding (56), the surface Σ 1 is given by which is a timelike or null surface. In this case, there exists no Σ 2 surface.

B. Properties of the Solution for N = 2
Our new solution (57) can be written as Using (41), (44), (46) and (49), the functions h and h 1 can be obtained as Then, using (28) and (68), the total charge Q T in a spherical region with the radius R 0 is given by In our solution b(t, r) in (67) there are six arbitrary constants. We can reduce this number to four by scaling. It is easy to show that the function b(t, r) is form invariant under the following scalings where m and n are arbitrary nonzero real numbers. Hence, out of 6 parameters only 4 of them can be considered generic. In the next sections, without loosing any generality we use the following two different parameterizations to represent our new solution.
where k 1 and k 2 are any real numbers. Here in the case of part A we will consider only the cases µ > 0, k 1 > 0, k 2 > 0. The reason for presenting the above two different representations of our solution is to show how it differs from the known exact solutions.
A. The First Representation: The case of µ = c2 c3 For c 3 = 0, we can consider the following identifications where c 3 > 0 and k is the spatial curvature constant corresponding to 0 for the flat and to ±1 for closed and open universes in general. Using the above identifications and the δ in (46), we can obtain our c 4 constant as We defined our constants c 0 , c 1 , c 2 , c 3 and c 4 in such a way that our solution (67) reduces to Vaidya-Shah solution (3) (for either c 0 = 0 or c 2 = 0), as we will see in Section V A 1. Then, our a(t, r), b(t, r), h(r), h 1 (r), σ(t, r), ρ(t, r), p(t, r) and F 01 (t, r) functions become where S(t, r), X(t, r) and Y (t, r) functions are given in Appendix B. Here, without losing any generality, we have set q(t)ȧ(t)/a(t) = 1. Hence, in our new solution (75) and (76), in addition to the mass M , charge Q and the spatial curvature constant k, we have a new parameter µ. When µ = 0, this solution reduces to the Vaidya-Shah solution (3), as we will see in Section V A 1. Our solution reduces to the Reissner-Nordström metric when µ = k = 0 and a(t) = 1 in isotropic coordinates. When µ = k = 0 and a(t) = e √ Λ 3 t then we obtain Schwarzschild-Reissner-Nordström-de Sitter metric with cosmological constant Λ.

Remark 1:
We point out that in contrast to the Vaidya-Shah solution, in our new solution, the current vector J µ (or the charge density σ) is non-zero for the flat spatial curvature constant, i.e. k = 0. On the other hand if µk = 1 where the charge density and the total charge in a volume of radius R 0 vanish our solution reduces to the FRW metric (see Remark 3).
For this solution, we have the following points.
(i) The surface Σ 1 is given as Hence, Σ 1 exists only when M < 0 and |Q| > M as we mentioned also in the general case in page 8. Then, we conclude that Σ 1 does not exist for physical cases.
(ii) The surface Σ 2 is given by the following equation which requires M 2 − Q 2 > 0.
(iv) At the spatial origin, i.e r → 0, the metric functions a(t, r) and b(t, r) as well as σ(t, r), ρ(t, r) and p(t, r) are nonsingular in general except for the cosmological models with a(t) → 0, see Appendix B for more details.
(v) At the spatial infinity, i.e r → ∞, the metric functions a(t, r) and b(t, r) as well as σ(t, r), ρ(t, r) and p(t, r) remains regular and the behavior of this model at the asymptotic region is different than the FRW solution, see Appendix B for more details. For c 2 = 0, one may also consider the following identifications where c 2 > 0 and k 1 and k 2 are two generally different spatial curvatures. Using the above identifications and the δ in (46), we can obtain c 4 constant as where k 1 = k 2 . For this case, the a(t, r), b(t, r), h(r), h 1 (r), σ(t, r), ρ(t, r), p(t, r) and F 01 (t, r) functions can be found as where S(t, r), X(t, r) and Y (t, r) functions are given in Appendix C.

Remark 2:
We point out that the case k 1 = k 2 = k reduces to a FRW metric with spatial curvature k (see Remark 3).
For this solution, one realizes the following points.
(i) Depending on the sign and values of M = 0 and Q = 0 parameters, we have Hence, we have exactly the same conclusion as the previous case that Σ 1 does not exist for physical cases.
(ii) The surface Σ 2 is given by the following equation which requires M 2 − Q 2 > 0.
(iv) At the spatial origin, i.e r → 0, the metric functions a(t, r) and b(t, r) as well as σ(t, r), ρ(t, r) and p(t, r) are nonsingular in general except for the cosmological models with a(t) → 0, see Appendix C for more details.
(v) At the spatial infinity, i.e r → ∞, the metric functions a(t, r) and b(t, r) as well as σ(t, r), ρ(t, r) and p(t, r) remains regular and the behavior of this model at the asymptotic region is different than the FRW solution, see Appendix C for more details.
We summarize this section with the following Theorem.
Theorem 3: The spacetime represented by our solution for N = 2 either in (75) and (76) or in (87) and (88) are nonsingular in the sense that all the functions a(t, r), b(t, r), p(t, r), ρ(t, r) and σ(t, r) either go to zero or to a finite value as r goes to zero or to infinity.

C. Properties of the Solution for N = 1
For N = 1, using (53) and (54), the functions h and h 1 can be obtained as which can be written also as Here, similar to the previous solutions and without losing any generality, we set β(t) = a(t), q(t)ȧ(t)/a(t) = 1, b 0 = 1 and b 1 = k. Then, we have a(t, r) = 1 Using (28) the total charge Q T in a bounded region with r = R 0 is given by For this solution, the behavior of the metric functions, charge and energy densities as well as the pressure at the spatial origin or asymptotic region in general depend on the explicit form of the arbitrary function ν 0 (r). Then, without its explicit form, we can not discuss accurately on the properties of the singular surfaces as well as the properties at the spatial origin or infinity. We will introduce some special subclasses of this general solution in the sections V B and VI B, and discuss briefly how these subclasses can be a reasonable physical solutions or not. For example, regarding (100), one can show that the second term in the pressure (106) diverges at r → ∞ for ν 0 (r) ∝ 1 r n with n < 2. One may argue that these types of solutions can not be reasonable physical solutions regarding their divergence at the asymptotic region. As our next work, we will classify various possible subclasses for N = 1 regarding the possible physical choices for the arbitrary ν 0 (r) function.

V. Special Subclasses of the General Solutions and Their Properties
In this section, we investigate some particular subclasses of our general solutions as well as their properties.
A. Subclasses of N = 2 1. The case of either c 0 = 0 but c 2 = 0, or c 2 = 0 but c 0 = 0 Both these cases correspond to the Vaidya-Shah solution (3). To show that, for example, we consider the case of c 2 = 0 but c 0 = 0 which leads to where with the identifications of c 0 = 1, Consequently, we can also find the metric function a(t, r) as where without losing any generality, we have set q(t)ȧ(t)/a(t) = 1. We see that our metric functions a(t, r) and b(t, r) exactly reduce to the Vaidya-Shah solution (3). Then, the Vaidya-Shah solution can be considered as one of the particular subclasses of our generalized solution (50). Also, the k parameter here is the spatial curvature constant which in general corresponds to zero for the flat, and to ±1 for the closed and open universes, respectively. One can realize the following points about this solution. (ii) The surface Σ 2 is given by the following equation which requires M 2 − Q 2 > 0.
(iv) As it is proved for the general solutions in the section , for r → 0 and r → ∞, charge density, mass density and pressure remain finite also for this subclass.
These cases correspond to the same spacetime geometry. Then, we discuss only the case of c 1 = 0 as follows. For this case, our solution (50), takes the following form Using (41), (44), (46) and (49), the functions h(r) and h 1 (r) can be obtained as Similar to the case of Vaidya-Shah solution, we can consider the following identifications where requires c 2 > 0. Then, the metric function b(t, r) takes the following form where k is the spatial curvature constant. Using the above identifications and the δ in (46), we obtain our c 4 constant as Then, our a(t, r), b(t, r), h(r), h 1 (r), σ(t, r), ρ(t, r), p(t, r) and F 01 (t, r) functions become b(t, r) = a(t) 1 + M a(t) where S(t, r), X(t, r) and Y (t, r) functions are given in Appendix D. Here, without losing any generality, we have set q(t)ȧ(t)/a(t) = 1. One can realize the following points about this solution.
(i) Depending on the sign and values of M = 0 and Q = 0 parameters, we have which exists only for unphysical cases, i.e M − |Q| < 0.
(ii) The surface Σ 2 is given by the following equation which requires M 2 − Q 2 > 0.
(iv) Again, as it is proved for the general solutions in the section IV B, for r → 0 and r → ∞, charge density, mass density and pressure remain finite for this subclass.

The case of c 4 = 0
Regarding (46), this case corresponds to δ = 0 and α 0 = 0. Then, the metric function b(t, r) in (50) takes the following form which similar to the sections IV A and IV B, can be demonstrated in both the (k, µ) and (k 1 , k 2 ) representations. Also, we find that h 1 = 0, and the function h(r) takes the following form where requires the condition γ < 0. One may consider 4γ = −c 2 Q 2 which reduces our solution here to the solutions with M = 0 in the sections IV A and IV B, i.e to the charged massless solutions. We consider the following identifications where here k 1 and k 2 are generally two different spatial curvatures. For this case, a(t, r), b(t, r), σ(t, r), ρ(t, r), p(t, r) and F 01 (t, r) functions read as where we have supposed c 2 > 0 and S(t, r), X(t, r) and Y (t, r) functions are given in Appendix E. Regarding (132) and (133), this solution is the generalization of Vaidya-Shah solution to the case of two spatial curvature with M = 0.
For this solution, one realizes the following points.
(i) The surface Σ 1 is given by the following equation where in contrast to the previous cases, it does exist as a physical case.
(ii) The surface Σ 2 is given by the following equation where can not exist as a physical case.
(iii) Similarly, as it is proved for the general solutions in the section IV B, for r → 0 and r → ∞, charge density, mass density and pressure remain finite for this subclass.

The case of γ = 0
Regarding the condition to obtain (50), i.e. β 1 (t)β 2 (t) = γ, this case corresponds to the situation where at least one of β 1 (t) and β 2 (t) in (38)  For this case, our metric functions take the following form as well as Then, one can find that for r → 0, all the quantities in this solution are regular while at the asymptotic region, i.e r → ∞, the pressure diverges by its second term in (148). Then, regarding this unusual asymptotic behavior, one may argue that this solution can not be a physical charged solution. However, in the section VI B 2, we will show that the solution for ν 0 = constant can be a physical uncharged solution for the flat universe (k = 0). Also, as we stated at the end of the section IV C, we will classify the possible choices by this kind of physical arguments in our next work.

VI. Uncharged Solutions and Their Properties
In this section, we explore the uncharged solutions and their properties for N = 2 and N = 1 in detail.

A. Uncharged Solutions for N = 2
To obtain the uncharged solutions for N = 2, regarding (25) and (68), we first assume that the constants c 0 , c 1 , c 2 and c 3 are nonzero. Then, we investigate special cases where some of these parameters vanish or they are related. Regarding (25) and (68), there are two main possibilities to obtain uncharged solutions.
By the coordinate transformations a 1 (t)dt = dT, r 1 + kr 2 = R, the new metric in the new coordinates T and R becomes whereā 2 (T ) = a 2 (t(T )). Hence, this special case is identical to the Friedmann-Robertson-Walker model.

The case of δ 2 = 4γ
For this case, the function h(r) in (68) and consequently the function F 01 = ψ(t, r) in (25) vanish and the uncharged case (σ(t, r) = 0) can be provided. Considering δ 2 = 4γ, the metric function b(t, r) in (50) takes the form of and the corresponding a(t, r) metric function will be a(t, r) = Similar to the previous solutions, one can set β(t) = a(t) and q(t)ȧ(t)/a(t) = 1. Here, we assume that a(t) is nonnegative for all t. Thus, we find We have the following subclasses of this general solution. Here also we define our c 0 , c 1 , c 2 and c 3 parameters in such a way that our general solution reduces to the McVittie solution as one of its particular subclasses.
For this case, the metric functions (156) and (157) take the following forms where by the identifications c 0 = 1, and √ c 3 = δ M , they read as This solution is the McVittie solution for the flat background universe (k = 0). The density and pressure profiles of this case can be obtained as Then, one finds the following points for this solution.
(a) Σ 1 surface exists only for M < 0 as (b) Σ 2 surface exists only for M > 0 given by Interestingly, one notes that regarding the sign of M parameter, the singular surfaces can be spacelike or timelike. However, for physical cases M > 0, the only existing singular surface is the spacelike surface Σ 2 . This singular surface corresponds to the big-bang singularity [22].
(2) The case of c 2 = 0 or c 0 = 0. For the case c 2 = 0, the metric functions (156) and (157) take the following forms where by the identifications c 0 = 1, c 1 = k and δ √ c3 = M , we have This solution is the generalization of the McVittie solution to a non-flat background universe (k = 0). This solution can be also identified to the Vaidya-Shah solution with Q = 0, where we have previously addressed its asymptotic behavior in section . Then, one finds the following points for this solution.
(a) Σ 1 surface is given by which exists only for the unphysical cases, i.e for M < 0.
For this case c 1 = 0, the metric functions (156) and (157) take the following forms where by the identifications c 0 = 1, c3 c2 = k and δ √ c2 = M , they read as This solution is also another generalization of the McVittie solution to a non-flat background universe (k = 0) with different identification set of our integration constant parameters. We have studied the charged generalization of this solution in section VA 2 with detail of its behavior at the spatial origin and infinity. Then, to avoid the repetition, one can set Q = 0 to realize the properties of this solution. Then, one finds the following points for this solution.
(a) Σ 1 surface exists only for M < 0 as (b) Σ 2 surface exists for M > 0 as (4) The case where none of the c i parameters are zero.
For this case, the metric functions (156) and (157) take the following forms where by identifications c 0 = 1, c 1 = k 1 , √ c 2 = δ M and k 2 = c3 c2 , they read as Remark 3: When the metric function b(t, r) takes the form we showed that, in Section VI A 1, that the spacetime reduces to the FRW universe with scale factor R(t) and spatial curvature 4k. Hence this suggests us that the metric function b(t, r) in (186) is a kind of nonlinear superposition of two different FRW b-functions If k 1 = k 2 each one describes different FRW universes. The function b 1 (t, r) belongs to a FRW universe with the scale factor a(t) and the spatial curvature 4k 1 and the function b 2 (t, r) belongs to another FRW universe with the scale factor M 2 4a(t) and the spatial curvature 4k 2 . Hence our uncharged solution is a kind of a nonlinear superposition of two different FRW metrics with different spatial curvatures. If initially a(t) → 0 then the function b 2 (t, r) is dominant in b(t, r) and the corresponding universe is initially a FRW universe with spatial curvature 4k 2 , see Appendix F. On the other hand if a(t) → ∞ as t → ∞ then the function b 1 (t, r) is dominant in the function b(t, r) and the universe is described by a FRW metric with the spatial curvature 4k 1 , see Appendix F. If k 1 k 2 ≤ 0 then we obtain an interesting result saying that the universe undergoes a kind of a topological change. In between, for t ∈ (0, ∞), the universe is a mixture of the above two FRW universes. If k 1 = k 2 = k then the two FRW universes collapse to a single one with the spatial curvature 4k.
One realizes the following points for this solution.
(i) Σ 1 surface exists only for M < 0 as (ii) Σ 2 surface exists for M > 0 as 3. The case of c 0 = c 2 = 0 For this case, δ, h and h 1 functions vanish and the metric function b(t, r) reads as where by definingã(t) = β(t) takes the following simple form Then, using suitable coordinate transformations, one can show that this solution can be identical to the flat FRW solution. Thus, the spatially flat FRW solution is one of the uncharged subclasses of our general solution (50) with the parameters of c 0 = c 2 = 0.
The case of c 1 = c 3 = 0 For this case, δ, h and h 1 functions vanish and the metric function b(t, r) will be only a time dependant function as B. Uncharged Solution for N = 1 1. Uncharged solution for ν 0 (r) One can find that both the functions h(r) and F 01 in (101) and (103), respectively, vanish for Thus, we have a(t, r) = 1 1 + c5 b(t, r) = a(t) 1 + k r 2 1 + Then, regarding the coordinate transformation in the section VI A 1, this solution gives also the FRW model.

Uncharged solution for ν 0 = constant.
For this case, one can find that both the functions h(r) and F 01 in (101) and (103), respectively, vanish only for k = 0 or ν 0 = 0. Then, the condition for having uncharged solution for ν 0 = constant = 0 is similar to the Vaidya-Shah solution, where the uncharged case is provided only for k = 0. For this case, we find b(t, r) = ν 0 + a(t).
Accordingly, one can show that this solution also is identical to the flat FRW solution.

VII. Apparent Horizons and Null Geodesics
The areal distance R is defined as R = rb(t, r). Among the constant R surfaces the null ones are called the apparent horizons. In our case there are two apparent horizons where a = qḃ/b. One can verify that the apparent horizons defined above reduces to the those given in [7] for the charged McVittie solution obtained by ∇ c R∇ c R = 0 where R is defined as the areal radius. Then, there are two possibilities as for the location of the apparent horizon H = H 1 ∪ H 2 .
According to the spacetime metric (16), ingoing and outgoing radial null geodesics x µ = (t, r(t), θ = θ 0 , φ = φ 0 ), where θ 0 and φ 0 are constants, are given by where " ± " signs represent the "outgoing" and "ingoing" geodesics, respectively. These null geodesics when entered in the apparent horizon H, they stay there. To see this, when the radial null geodesics lie in H, by taking the derivative of r(t) b(r(t), t) = c with respect to t, we obtain Eqs. (199) and (200) are consistent because the expressions in the right hand sides of these equations are equal due to the nullity condition (197) or (198) of the apparent horizon H.
To study the causal and global structures of the spacetime we have to maximally extend the existing coordinates {−∞ < t < ∞, r ≥ 0, 0 < φ < 2π, 0 < θ < π} to a coordinate system where the areal distance R is one of the coordinates as done in [22,23,33]. We postpone a detailed study of this case as our future work. However, just to give an idea how the the radial null geodesics (NG) behave we plot them in Figure 1.
In the same figure we give also apparent horizons H 1 , H 2 and singular surface Σ 2 of this uncharged solution given by H : a 2 (t, r) (b(t, r) + rb ′ (t, r)) 2 − r 2 b 2 (t, r)ḃ 2 (t, r) = 0, respectively, corresponding to the metric functions b(t, r) = a(t) 1 + kr 2 for typical values of M, µ and k parameters in a de Sitter background as in the terminology of FRW models.

VIII. Conclusion
We have found two classes of solutions of Einstein-Maxwell-Perfect Fluid field equations with a cosmological constant in a spherically symmetric spacetime. In particular the first class corresponding to the N = 2 case contains six parameters four of which are essential generalizes the Vaidya-Shah solution. The uncharged version of our solution generalizes the McVittie solution. We showed that there are some, depending on sign of the parameters, timelike and spacelike surfaces where the spacetime becomes singular. We then investigated some special limits of our solutions in both classes. The list of our new solutions is given in Tables 1 and 2.  Among all the solutions we found in this work there are new charged and uncharged solutions of the Einstein-Maxwell-Perfect Fluid equations with cosmological constant. For the uncharged case the solution corresponding to the N = 2 class is a model of a universe which is a mixture of two different FRW universes with different spatial curvatures. We will study in particular this solution in more detail in a forthcoming publication.

Class Parameters Solution
Mathematical and Physical properties of our solutions can be summarized as follows: We proved three theorems: (i) The first theorem is on the reduction of the Einstein field equations into a single ordinary nonlinear differential equations.
(ii) The second theorem is on the two classes of solutions corresponding to N = 1 and N = 2.
(iii) The third theorem is on the regularity of spacetime when the radial coordinate r goes to zero and to infinity.
Other properties are the following: (i) Our solution corresponding to N = 2 of Section IV when the two spatial curvatures are equal, i.e., k 1 = k 2 , reduces to FRW metric if the total charge in the universe vanishes.
(ii) There exits some spacelike surfaces where the pressure of the fluid diverges but the mass density of the fluid distribution remains finite. Such spacelike surfaces are known as (sudden) cosmological singularities.
(iii) Null geodesics crossing the above spacelike surfaces remain in these surfaces.
(iv) If µ = 0 in the first representation in Section IV, we obtain the Vaidya-Shah metric. When µ = 0 and spatial curvature k = 0 and the scale factor a(t) = 1 we obtain the Reissner-Nordström metric in isotropic coordinates. When µ = k = 0 and a(t) = e √ Λ 3 t then we obtain Schwarzschild-Reissner-Nordström-de Sitter metric with cosmological constant Λ. When the charge parameter Q vanishes we obtain a generalization of McVittie metric. If the charge parameter vanishes and µ = k = 0 we get the McVittie solution. Furthermore if a(t) = 1 we obtain the Schwarzschild metric in isotropic coordinates.
(v) In particular for the uncharged case our solution can be considered as a nonlinear superposition of two different FRW metrics with different scale factors and different spatial curvatures. Due to this effect in our model our universe may start with a FRW universe with spatial curvature k 2 and ends up with a FRW universe with a different spatial curvature k 1 , so that k 1 k 2 ≤ 0. This means that the universe may undergo a change of topology.
Note added in the proof.-Recently, we became aware of a paper by Mashhoon and Partovi [34], focusing on the gravitational collapse of charged fluid spheres. Although the problem studied by Mashhoon and Partovi differs from ours, they present an exact solution of Einstein field equations for inhomogeneous charged fluid distribution which corresponds our N = 2 solution, with zero cosmological constant, in Theorem 2.
Appendix A: Kustaanheimo-Qvist approach for the charged case Following the Kustaanheimo-Qvist [31] approach (see also [4,32]), one can use the change of variables L = b −1 and x = r 2 to transform the equation (30) for the uncharged case (h = 0) to the following ordinary differential equation 4x where using the identification F (x) = − h1 4xr (206) can be written in the form of For the case where the charge is also included we obtain where F 1 (x) = − h1 4xr and F 2 (x) = − h 2 2x . For the case where charge is zero, as represented in [32], there are three different approaches to finding solutions for (207). The first approach is based on an ad hoc ansatz for the function F (x) [35]. The second approach is based on the answering to the question of "for which functions F (x) the equation (207) admits one (or two) Lie point symmetries or Noether symmetries?" [31]. The third approach introduced by Wyman [36] is based on the solutions of (207) which have the Painleve property. All known solutions belong to this class. For the case of F = 0, the solution to (207) is 1 b = L = A(t)r 2 + B(t). Some other subclasses with F = (ax 2 + 2bx + c) − 5 2 where a, b and c are real constants are given in the following [32].
Our uncharged solution given in (154) as gives a new solution to the above equation (207) where the function L is given by then the function F (x) is found as For the charged case we will investigate all possible new exact solutions and the properties of the equation (208) in a later communication.
Appendix B: The case of µ = c 2 c 3 For this case, S(t, r), X(t, r) and Y (t, r) functions are given by the following forms.
(i) For a(t) → 0, regarding (186), we have b(t, r) → M 2 4a(t) The corresponding matter density ρ(t, r) can be read from (31) as where using q(t) = a(t) a(t) and R(t) = M 2 4a(t) reduces to the following standard Friedmann equation describing an FRW universe with the scale factor R(t) and spatial curvature 4k 2 .
(i) For the case of a(t) → ∞, from (186), we have b(t, r) → a(t) 1 + k 1 r 2 . (262) The corresponding matter density can be obtained from (31) as where it represents the matter density of an FRW universe with the scale factor R(t) = a(t) and spatial curvature 4k 1 .