1 Introduction

Dust is considered to be the simplest form of matter composed of pressureless, radiation and is abundant in galaxies, clusters and superclusters in a cosmological context. This form of matter has been widely studied in a variety of contexts such as scalar fields, rotational and irrotational distributions, forms of spatial symmetry, null geodesics, the modelling of gravitational collapse and black hole formation [1,2,3,4,5]. It has also been shown that stars are encompassed by these radiating pressureless particles which make up the atmosphere of the star. In a pioneering work, Vaidya [6] described the spacetime outside a spherically symmetric star emitting or absorbing null dust. Wang and Wu [7] generalised the Vaidya solution for a type II fluid. The gravitational collapse of spherically symmetric dust distributions can result in black holes or naked singularities depending on initial conditions. Of this nature, the Lemaître–Tolman–Bondi (LTB) metric is among one of the earliest solutions to the Einstein field equations that describe a spherically symmetric dust distribution undergoing collapse in a neutral setting [8,9,10]. In a similar treatment, Papapetrou [11] showed that the end result of the gravitational collapse of null dust in a Vaidya spacetime is a naked singularity. This Vaidya–Papapetrou model was then extended by [12] to investigate outgoing causal geodesics. Furthermore, the presence of the electromagnetic field has a significant effect on the gravitational behaviour of astrophysical objects. In particular, the gravitational collapse of a compact object such as a neutron star can be counteracted by the repulsive nature of the Coulombic force along with the pressure gradient. It has been noted that charge is one of the characterising features of a black hole along with its mass and angular momentum according to the ‘no hair’ theorem proposed by Misner et al. [13]. Significant treatments in this direction include the charged Vaidya solution which was used to study the thermodynamics of black holes [14] and evaporating black holes [15]. Breithaupt et al. [16] demonstrated that the rigid motion of a rotating disc of dust with constant charge results in a Kerr–Newman black hole in the ultra-relativistic limit. Therefore it remains important to study the simple form of dust matter in the presence of an electromagnetic field.

The dynamics of charged dust gravitating systems have been of great interest in general relativity and have been studied in the context of static distributions, rotating systems, and various astrophysical scenarios which include gravitational collapse and black hole physics. The earliest treatments in this direction can be attributed to Majumdar [17] and Papapetrou [18]. Following the result of Weyl [19], Majumdar demonstrated a functional relationship between the gravitational potential and the electrostatic potential, and a similar approach was considered by Papapetrou. In a nonstatic spherically spacetime, the general solution to the Einstein–Maxwell field equations was presented in implicit form by Vickers [20]. Raychaudhuri [21] considered the rigid motion of charged dust without specifying symmetry. Tikekar [22] obtained physically viable models for a spherically symmetric static sphere of charged dust in equilibrium. More recently, Hansraj et al. [23] presented exact solutions for static charged distributions and identified cases of the Finch-Skea geometry and the Schwarzchild metric. Brassel et al. [24] modelled the dynamics of a spherically symmetric radiating star with three concentric regions for which pressureless null radiation was considered and the geometry described by the Vaidya spacetime. Several families of solutions for various realistic equations of state were also obtained in this work. Other treatments of charged dust have been discussed extensively in [25,26,27,28,29]. In the context of black holes, Oppenheimer and Snyder [30] generated the first analytic model of black hole formation for a homogeneous dust cloud. Additionally, the final fate of gravitational collapse of a dust fluid with spherical symmetry has been considered in a number of works by Christodoulou [31], Newman [32] and Banerjee et al. [33].

In the context of modified gravity, fewer treatments of charged null dust models are known. It is important to investigate the effects of higher order corrections, higher dimensions and other possible modifications of the theories, on pressureless models. The gravitational collapse of null dust in f(R) gravity was analysed in detail by Ghosh and Maharaj [34] in higher dimensions. It was found that a naked singularity forms in the situation where a null dust is injected into an initially higher dimensional (anti) de-Sitter spacetime. A charged spherically symmetric model in modified f(G, T) gravity under the class I embedding condition was obtained by Maurya et al. [35]. An anisotropic compact stellar system in the presence of charged dust was studied by Mustafa et al. [36] in the framework of Rastall gravity in the phantom field regime. Generalised static solutions of an equilibrium state were found and as a result, strange stellar body configurations were identified in the physical features of the solution.

It is important to determine the gravitational dust dynamics in the presence of an electromagnetic field in the framework of Einstein–Gauss–Bonnet (EGB) gravity. Little information about such structures in a static spacetime are known in this framework of gravity. EGB gravity is one of the more promising theories of gravity which belongs to a special case (order two) of Lovelock gravity. In this theory, higher order curvature terms arise and appear as corrections to Einstein gravity. Moreover the field equations are second order and quasilinear which is a remarkable feature as similar properties are found in general relativity. As such numerous results have been conducted in this framework of gravity in a variety of astrophysical scenarios. For example, Kobayashi generalised the Vaidya type radiating solution in EGB gravity [37]. Maharaj et al. [38] obtained exact solutions for a static spherically symmetric distribution in five dimensions. Moreover, Brassel et al. [39] investigated the gravitational collapse of a null radiation shell in five dimensions which was then generalized to arbitrary dimensions in [40]. Some other pioneering results on black hole solutions, radiating stars, static systems and charged distributions in EGB gravity include those by [41,42,43,44]. Furthermore, seminal works on the nature of a neutral dust fluid distribution include those considered by Maeda [45] who modelled the gravitational collapse of a spherically symmetric neutral dust cloud in higher dimensions. Jhingan and Ghosh [46] presented an exact model of the gravitational collapse of an inhomogeneous dust distribution which represents a generalization of the LTB spacetimes to five dimensional EGB gravity. Ghosh and Maharaj [47] obtained exact nonstatic null dust solutions in the novel four dimensional limit of EGB gravity. Similar works in this direction were conducted by [48,49,50,51]. In addition, Brassel et al. [52] analysed the process of gravitational collapse of a null dust radiation shell in five dimensional EGB gravity in the presence of an electromagnetic field. It was found that the nature of collapse differs significantly to its uncharged counterpart. Hansraj [53] was the first to study static charged dust distributions in EGB gravity with spherical symmetry and obtained a number of exact solutions in closed form in dimensions \(N=5\) and \(N=6\). Brassel et al. [54] studied the formation of singularities in the context of the Cosmic Censorship Conjecture (CCC) for a charged radiating Boulware–Deser spacetime in five dimensions; the EGB analogue of the Vaidya spacetime with charge. In light of the above it is worthwhile to investigate the dynamics of charged dust in higher dimensional EGB gravity.

In this paper we consider the matter distribution to be a pressureless fluid in the presence of an electromagnetic field. We generate the EGB field equations for a spherically symmetric static spacetime in arbitrary dimensions. For such a fluid distribution, the pressure is zero and the governing equation that arises is an Abel differential equation of the second kind in one of the gravitational potentials. It is also a second order linear differential equation in the second potential. We demonstrate solutions to the governing equation by first treating it as a second order linear differential equation. Furthermore we show that the Abel differential equation can be transformed to a simpler canonical form which is still first order and nonlinear in nature. We generate several new classes of solutions to this equation by imposing restrictions on the canonical differential equation.

2 EGB gravity and Maxwell’s equations

The electromagnetic energy tensor \({{\varvec{E}}}\) in N dimensions is represented by

$$\begin{aligned} E_{a b}=\frac{1}{{\mathcal {A}}_{N-2}}\left( F_{a c}F_{b}{}^{{}c}-\frac{1}{4}g_{a b}F_{c d}F^{c d}\right) , \end{aligned}$$
(1)

which is comprised of the Faraday tensor F and the metric tensor g. The quantity \({\mathcal {A}}_{N-2}\) represents the surface area of the unit \((N-2)\)-sphere and is defined by

$$\begin{aligned} {\mathcal {A}}_{N-2}= & {} \frac{2 \pi ^{\frac{N-1}{2}}}{\Gamma \left( \frac{N-1}{2}\right) }, \end{aligned}$$
(2)

where \(\Gamma (\ldots )\) is the gamma function. Maxwell’s equations govern the electromagnetic field and these governing equations are expressed in a covariant manner as

$$\begin{aligned} F_{a b;c} + F_{bc;a} + F_{c a;b}= & {} 0, \end{aligned}$$
(3a)
$$\begin{aligned} F^{a b}{}_{{ };b}= & {} {\mathcal {A}}_{N-2}J^{a}, \end{aligned}$$
(3b)

where \({{\varvec{J}}}\) represents the current for a non-conducting fluid and is represented by

$$\begin{aligned} J^{a}= & {} \sigma u^{a}, \end{aligned}$$
(4)

where \(\sigma \) is the proper charge density and \(u^{a}\) is a unit and timelike N-velocity \((u^{a}u_{a}=-1)\).

In the case of a dust distribution, the energy momentum tensor is given by

$$\begin{aligned} T_{a b}= & {} \rho u_{a}u_{b}, \end{aligned}$$
(5)

where \(\rho \) is the energy density. Thus the total energy momentum tensor for a charged dust matter distribution \(\mathbf {{\mathcal {T}}}\) is then defined by

$$\begin{aligned} {\mathcal {T}}_{ab}= & {} T_{ab}+E_{ab}. \end{aligned}$$
(6)

The EGB field equations for a charged dust distribution are given by

$$\begin{aligned} G_{a b} - \frac{\alpha }{2} H_{ab} = \kappa _{N} \mathbf {{\mathcal {T}}}_{a b}. \end{aligned}$$
(7)

where \({{\varvec{G}}}\) is the Einstein tensor, \(\alpha \) is the Gauss–Bonnet constant and \(\kappa _{N}\) is the Einstein coupling constant denoted by

$$\begin{aligned} \kappa _{N}=\frac{2 \left( N-2\right) \pi ^{\frac{N-1}{2}} G}{c^{4}\left( N-3\right) \Gamma \left( \frac{N-1}{2}\right) }. \end{aligned}$$
(8)

The Einstein coupling constant, \(\kappa _{N}\) and the surface area term \({\mathcal {A}}_{N-2}\) appear as a ratio in the field equations as \(\frac{\kappa _{N}}{{\mathcal {A}}_{N-2}}=\frac{N-2}{N-3}\) from (2) to (8). In EGB gravity, curvature terms that are quadratic in nature appear as corrections to Einstein gravity in the form of the EGB tensor \({{\varvec{H}}}\) which is represented by

$$\begin{aligned} H_{a b} = g_{a b}L_{G B}-4RR_{a b} + 8R_{a c}R^c{}_{{} b} + 8R^{c d}R_{a c b d}-4R_{a}{}^{{}c d e} R_{b c d e}. \end{aligned}$$
(9)

The Gauss–Bonnet term \(L_{G B}\) is expressed by

$$\begin{aligned} L_{G B} = R^{2} + R_{a b c d} R^{a b c d} - 4R_{c d}R^{c d}, \end{aligned}$$
(10)

which consists of a linear combination of the quadratic curvature terms. The field equations to be solved are (3) and (7).

3 Static charged dust

The spherically symmetric static metric in N dimensions is given by

$$\begin{aligned} ds^2=-e^{2\nu (r)} dt^2 + e^{2\lambda (r)} dr^2 + r^2 d \Omega ^{2}_{N-2}, \end{aligned}$$
(11)

where \(\nu (r)\) and \(\lambda (r)\) are arbitrary functions of r that represent the metric potentials. The metric for the unit \((N-2)\)-sphere is written as

$$\begin{aligned} d\Omega ^{2}_{N-2}= & {} \sum _{i=1}^{N-2}\left( \left[ \prod _{j=1}^{i-1} \sin ^2({\theta _{j}})\right] (d\theta _{i})^{2}\right) . \end{aligned}$$
(12)

We choose the electromagnetic potential to be of the form

$$\begin{aligned} A_{a}= & {} \left( \Phi (r), 0, 0,\ldots ,0\right) , \end{aligned}$$
(13)

and the nonzero Faraday tensor components are calculated as

$$\begin{aligned} F^{0 1}=-F^{1 0}= e^{-2(\nu +\lambda )}\Phi '(r)=e^{-\left( \nu +\lambda \right) }E(r), \end{aligned}$$
(14)

where E(r) is the electrostatic field intensity. The total charge within the hypersphere is given by

$$\begin{aligned} l(r)= & {} {\mathcal {A}}_{N-2}\int ^{r} \sigma e^{\lambda }r^{N-2} d{\tilde{r}}. \end{aligned}$$
(15)

For the charged spherically symmetric dust distribution we can obtain the EGB field equations in N dimensions. These are expressed as

$$\begin{aligned} \kappa _{N}\left( \rho +\frac{E^{2}}{2 {\mathcal {A}}_{N-2}}\right)= & {} \frac{\left( N-2\right) }{r^{4}e^{4\lambda }}\left. \bigg [r^{3}e^{2\lambda } \lambda '+ \frac{\left( N-3\right) r^{2}e^{4\lambda }}{2}-\frac{\left( N-3\right) r^{2}e^{2\lambda }}{2}\right. \nonumber \\{} & {} \quad \left. +{\hat{\alpha }}\left( e^{2\lambda }-1\right) \left( 2r\lambda '+ \frac{\left( N-5\right) \left( e^{2\lambda }-1\right) }{2}\right) \right] , \end{aligned}$$
(16a)
$$\begin{aligned} \frac{\kappa _{N}}{2{\mathcal {A}}_{N-2}}E^{2}= & {} -\frac{\left( N-2\right) }{r^{4}e^{4\lambda }}\left[ r^{3}e^{2\lambda }\nu '+\frac{\left( N-3\right) r^{2}e^{2\lambda }}{2}-\frac{\left( N-3\right) r^{2}e^{4\lambda }}{2}\right. \nonumber \\{} & {} \quad \left. +{\hat{\alpha }}\left( e^{2\lambda }-1\right) \left( 2r\nu '- \frac{\left( N-5\right) \left( e^{2\lambda }-1\right) }{2}\right) \right] , \end{aligned}$$
(16b)
$$\begin{aligned} \frac{\kappa _{N}}{2{\mathcal {A}}_{N-2}}E^{2}= & {} \frac{1}{r^{2}e^{2 \lambda }}\left[ \frac{\left( N-3\right) \left( N-4\right) }{2}+ r^{2}\left( \nu ''+\nu '^{2}-\nu '\lambda '\right) \right. \nonumber \\ {}{} & {} \quad \left. + \left( N-3\right) r\left( \nu '-\lambda '\right) +2{\hat{\alpha }}\left( \nu ''+\nu '^{2}-\nu '\lambda '\right) \bigg ]\right. \nonumber \\ {}{} & {} \quad + \frac{1}{r^{2}e^{4\lambda }}\bigg [6{\hat{\alpha }} \nu '\lambda '-2{\hat{\alpha }}\left( \nu ''+\nu '^{2}\right) -\frac{\left( N-3\right) \left( N-4\right) e^{4\lambda }}{2}\bigg ]\nonumber \\ {}{} & {} \quad + \frac{2 {\hat{\alpha }}\left( N-5\right) }{r^{3}e^{4\lambda }}\left[ \left( e^{2\lambda }-1\right) \left( \nu '- \lambda '\right) \right] \nonumber \\{} & {} \quad -\frac{{\hat{\alpha }}\left( N-5\right) \left( N-6\right) \left( e^{2\lambda }-1\right) ^{2}}{2r^{4}e^{4\lambda }}, \end{aligned}$$
(16c)
$$\begin{aligned} \sigma= & {} \frac{e^{-\lambda } \left( r^{\left( N-2\right) }E\right) ' }{r^{\left( N-2\right) }{\mathcal {A}}_{N-2}}, \end{aligned}$$
(16d)

in terms of canonical coordinates and \({\hat{\alpha }}=\alpha \left( N-3\right) \left( N-4\right) \). Note that primes represent differentiation with respect to the spacetime coordinate r. System (16) describes the gravitational behaviour of a charged dust distribution in higher dimensional EGB gravity. Note that we can eliminate \(E^{2}\) from (16b) and (16c) leading to an integrability condition involving only the potentials \(\nu \) and \(\lambda \).

The conservation equations \({\mathcal {T}}^{ab}{}_{{};b}=0\) lead to the result

$$\begin{aligned} \rho \nu '= & {} \frac{E}{{\mathcal {A}}_{N-2}r^{(N-2)}}\left( r^{\left( N-2\right) }E\right) ', \end{aligned}$$
(17)

also called the continuity equation. Equation (17) can also be obtained directly from the field equations (16). The continuity equation (17) is useful in generating exact solutions.

We now make use of a change of coordinates in order to simplify system (16). We let

$$\begin{aligned} e^{2\nu (r)} =y^{2}(x),\hspace{10pt} e^{-2\lambda (r)} = Z(x),\hspace{10pt} x = r^2, \end{aligned}$$
(18)

which was first implemented by Durgapal and Bannerji [55] in general relativity for a neutron star. In terms of the new variables the line element (11) becomes

$$\begin{aligned} ds^2=-y^{2}(x) dt^2 + \frac{1}{4x Z(x)} dx^2 + x d \Omega ^{2}_{N-2}. \end{aligned}$$
(19)

Then the charged dust EGB field equations (16) can be written in terms of the new variables y and Z. We obtain the matter quantities

$$\begin{aligned} \kappa _{N}\rho= & {} \left( N-2\right) \left[ \frac{\left( N-3\right) \left( 1-Z\right) -2x{\dot{Z}}}{2x}\right. \nonumber \\ {}{} & {} \left. + \frac{{\hat{\alpha }}\left( 1-Z\right) }{2 x^{2}}\left( -4x{\dot{Z}}+ \left( N-5\right) \left( 1-Z\right) \right) \right] \nonumber \\ {}{} & {} -\frac{\left( N-2\right) }{2\left( N-3\right) }E^{2}, \end{aligned}$$
(20a)
$$\begin{aligned} \sigma ^{2}= & {} \frac{ Z \left[ 2 x^{\frac{(N-1)}{2}}{\dot{E}} + \left( N-2\right) x^{\frac{(N-3)}{2}}E\right] ^{2}}{({\mathcal {A}}_{N-2})^2 \, x^{\left( N-2\right) }}, \end{aligned}$$
(20b)
$$\begin{aligned} \frac{\left( N-2\right) }{\left( N-3\right) }E^{2}= & {} \frac{2}{y}\left[ 2 x Z \ddot{y}+ x {\dot{Z}} {\dot{y}}+ \left( N-2\right) {\dot{y}} Z\right] \nonumber \\ {}{} & {} + \left( N-3\right) \left[ {\dot{Z}}+ \frac{\left( N-4\right) \left( Z-1\right) }{2x}\right] \nonumber \\ {}{} & {} +{\hat{\alpha }}\left[ \frac{4{\dot{Z}}{\dot{y}} \left( 1-3Z\right) }{y}+ \frac{8Z\left( 1-Z\right) \ddot{y}}{y}\right. \nonumber \\ {}{} & {} \left. +\frac{4\left( N-4\right) Z\left( 1-Z\right) {\dot{y}}}{ x y} + \frac{2\left( N-5\right) {\dot{Z}}\left( 1-Z\right) }{x}\right. \nonumber \\{} & {} \left. -\frac{\left( N-5\right) \left( N-6\right) \left( 1-Z\right) ^{2}}{2 x^{2}}\right] \nonumber \\{} & {} -\left( N-2\right) \left[ \frac{2Z{\dot{y}}}{y}+ \frac{\left( N-3\right) \left( Z-1\right) }{2x}\right. \nonumber \\ {}{} & {} \left. + {\hat{\alpha }} \left( 1-Z\right) \left( \frac{4Z{\dot{y}}}{x y}-\frac{\left( N-5\right) \left( 1-Z\right) }{2 x^2}\right) \right] , \end{aligned}$$
(20c)

where dots represent differentiation with respect to x, subject to the integrability condition

$$\begin{aligned}{} & {} 4x^2Z\left[ x+2{\hat{\alpha }} \left( 1-Z\right) \right] \ddot{y}+\left[ 2x^3 {\dot{Z}} + 4{\hat{\alpha }}x^2 {\dot{Z}}-12{\hat{\alpha }}x^2 Z {\dot{Z}}\right. \nonumber \\{} & {} \quad \left. -8{\hat{\alpha }}\left( N-3\right) x Z^2 +4\left( N-2\right) x^2 Z+ 8{\hat{\alpha }} \left( N-3\right) xZ\right] {\dot{y}}\nonumber \\{} & {} \quad +\left[ \left( N-3\right) x^2 {\dot{Z}}+ \left( N-3\right) ^2 x\left( Z-1\right) + 2{\hat{\alpha }}\left( N-5\right) x\left( 1-Z\right) {\dot{Z}} \right. \nonumber \\{} & {} \quad \left. -{\hat{\alpha }} \left( N-4\right) \left( N-5\right) Z^2 + 2{\hat{\alpha }}\left( N-4\right) \left( N-5\right) Z -{\hat{\alpha }}\left( N-4\right) \left( N-5\right) \right] y=0. \nonumber \\ \end{aligned}$$
(21)

For gravitating charged dust, Eq. (21) must be solved to obtain forms for y and Z. It is a nonlinear differential equation, however it may be viewed as a second order linear differential equation in y if Z is specified. Then \(\rho \), E and \(\sigma \) follow from the system (20). It is interesting to observe that the dynamical behaviour of charged gravitating dust is influenced by the Gauss–Bonnet coupling constant \(\alpha \) and the spacetime dimension N. When \(\alpha =0\) then we obtain N dimensional general relativity. The EGB case \(N=5\) is special leading to the simpler integrability condition

$$\begin{aligned}{} & {} 2xZ\left[ x+4\alpha \left( 1-Z\right) \right] \ddot{y}+\left[ x^2 {\dot{Z}} + 4\alpha x \left( 1-3Z\right) {\dot{Z}}-16\alpha Z^2 +6x Z+ 16\alpha Z\right] {\dot{y}}\nonumber \\{} & {} \quad +\left[ x {\dot{Z}} +2\left( Z-1\right) \right] y=0, \end{aligned}$$
(22)

where \({\hat{\alpha }}=2\alpha \), with dynamics that is clearly different for \(N\ge 6\) in (21).

We find exact solutions to (21) by generating an Abel differential equation and a second order linear differential equation.

3.1 Abel differential equation: \(y={\tilde{a}}\)

The first case we consider are a class models for which the metric potential y is given by

$$\begin{aligned} y= & {} {\tilde{a}}, \end{aligned}$$
(23)

where \({\tilde{a}}\) is some constant. The assumption (23) affects the spacetime geometry and ensures that the timelike congruences are acceleration-free. In EGB gravity stellar solutions have been found using the condition (23), and the resulting model has been discussed by Chilambwe et al. [56]. In the present investigation we show that dust distributions are also possible in EGB gravity. Here the Gauss–Bonnet contributons from the curvature nullifies the effect of the electric field to generate dust. We observe that static dust models also rise in f(R) gravity models and they display stable behaviour [57]. Equation (21) then becomes

$$\begin{aligned}{} & {} x\left[ \left( N-3\right) x + 2{\hat{\alpha }}\left( N-5\right) \left( 1-Z\right) \right] {\dot{Z}} -{\hat{\alpha }}\left( N-4\right) \left( N-5\right) Z^{2} \nonumber \\{} & {} \quad + \left[ \left( N-3\right) ^{2}x + 2{\hat{\alpha }} \left( N-4\right) \left( N-5\right) \right] Z \nonumber \\{} & {} \quad - \left( N-3\right) ^{2}x- {\hat{\alpha }}\left( N-4\right) \left( N-5\right) =0, \end{aligned}$$
(24)

which is an Abel differential equation of the second kind in Z. Equations of this type are difficult to solve. We observe that it is a first order differential equation which can be written as

$$\begin{aligned} S(x,Z)dZ+M(x,Z)dx= & {} 0, \end{aligned}$$

where

$$\begin{aligned} S(x,Z)= & {} x\left[ \left( N-3\right) x + 2{\hat{\alpha }}\left( N-5\right) \left( 1-Z\right) \right] ,\\ M(x,Z)= & {} -{\hat{\alpha }}\left( N-4\right) \left( N-5\right) Z^{2} - \left( N-3\right) ^{2}x - {\hat{\alpha }}\left( N-4\right) \left( N-5\right) \\{} & {} + \left[ \left( N-3\right) ^{2}x + 2{\hat{\alpha }} \left( N-4\right) \left( N-5\right) \right] Z. \end{aligned}$$

It can be observed that \(\frac{\partial M(x,Z)}{\partial Z}\ne \frac{\partial S(x,Z)}{\partial x}\). We now require an integrating factor \({\mathcal {R}}(x)\) such that \({\mathcal {R}}(x)S(x,Z)dZ+{\mathcal {R}}(x)M(x,Z)dx=0\) is an exact differential equation. This implies that

$$\begin{aligned} \frac{\partial \left( {\mathcal {R}}(x) M(x,Z)\right) }{\partial Z}= & {} \frac{\partial \left( {\mathcal {R}}(x) S(x,Z)\right) }{\partial x}. \end{aligned}$$

The above equation simplifies to

$$\begin{aligned} \frac{d{\mathcal {R}}(x)}{dx}= & {} \frac{\left( N-5\right) }{x}{\mathcal {R}}(x). \end{aligned}$$

Integrating the above yields \({\mathcal {R}}(x)=x^{N-5}\). We now multiply (24) by this form of \({\mathcal {R}}(x)\) to generate the general solution

$$\begin{aligned} {\hat{\alpha }}\left( N-5\right) x^{\left( N-4\right) }\left( Z-1\right) ^2-\left( N-3\right) x^{N-3}\left( Z-1\right) +C=0, \end{aligned}$$
(25)

where C is an integration constant. Observe that for this class of models, \(E^{2} \sim \frac{1}{x^{3}}\) (\(N=5\)) and \(E^{2} \sim \frac{1}{x^{N-2}}\) (\(N>5\)). In both cases, upon substitution of \(E^{2}\) in (20b), we obtain for the proper charge density, \(\sigma =0\), and as a result there is no current J. We do not pursue this case further.

3.2 Linear equation: \(Z=k\)

We choose the potential Z to be of the form

$$\begin{aligned} Z= & {} k, \end{aligned}$$
(26)

where k is an arbitrary constant. Then (21) becomes

$$\begin{aligned}{} & {} 4kx^2\left[ x+2{\hat{\alpha }}\left( 1-k\right) \right] \ddot{y}+4kx\left[ 2{\hat{\alpha }}\left( N-3\right) \left( 1-k\right) +\left( N-2\right) x \right] {\dot{y}}\nonumber \\{} & {} \quad +\left[ {\hat{\alpha }} \left( N-4\right) \left( N-5\right) \left( -k^{2}+2k-1\right) + \left( N-3\right) ^2 x\left( k-1\right) \right] y=0. \end{aligned}$$
(27)

This equation is a second order linear differential equation in the variable y. It has a solution in terms of hypergeometric functions. The explicit form of y can easily be found in terms of these special functions using the software package Maple. Some examples with specific forms of y are studied in Hansraj [53]. In general the form of the solution in terms of hypergeometric functions is complicated and not useful in applications. It is therefore necessary to consider solutions in which \({\dot{Z}}\ne 0\).

In the general relativistic limit, \({\hat{\alpha }}=0\), equation (27) takes the form

$$\begin{aligned}{} & {} 4kx^2\ddot{y}+4\left( N-2\right) kx{\dot{y}}+\left[ \left( N-3\right) ^2 \left( k-1\right) \right] y=0. \end{aligned}$$
(28)

This equation is of Euler–Cauchy type which can easily be solved to obtain the solution

$$\begin{aligned} y=A_{1}x^{\frac{\left( 1-\sqrt{k}\right) \left( N-3\right) }{2\sqrt{k}}}+A_{2}x^{-\frac{\left( 1+\sqrt{k}\right) \left( N-3\right) }{2\sqrt{k}}}. \end{aligned}$$
(29)

In the above, \(A_{1}\) and \(A_{2}\) represent integration constants. Equation (29) represents charged dust models in N dimensions in the general relativistic case. When \(N=4\), we then obtain

$$\begin{aligned} y=A_{1}x^{\frac{\left( 1-\sqrt{k}\right) }{2\sqrt{k}}}+A_{2}x^{-\frac{\left( 1+\sqrt{k}\right) }{2\sqrt{k}}}, \end{aligned}$$
(30)

which regains the result presented in Hansraj et al. [23]. We have therefore shown that the EGB dust models do regain known dust models in the \(N=4\) general relativity limit \({\hat{\alpha }}=0\). When \({\hat{\alpha }}\ne 0\) the potentials are expressed in terms of hypergeometric functions as mentioned above.

3.3 Linear equation: \(Z=1+{\tilde{b}}x\)

To find solutions when \({\dot{Z}}\ne 0\), we can try

$$\begin{aligned} Z= & {} {\tilde{a}}+{\tilde{b}}x, \end{aligned}$$
(31)

where \({\tilde{a}}\) and \({\tilde{b}}\) are arbitrary constants. Then Eq. (21) leads to solutions in terms of special functions, namely Heun functions. This can be easily verified with the software package Maple. Again the resulting analytic form of the solution is complicated and cannot easily be applied to study the physical features of the model. If we restrict the values of the parameters \({\tilde{a}}\) and \({\tilde{b}}\) then substantial simplification is possible.

The choice \({\tilde{a}}=1\) gives

$$\begin{aligned} Z= & {} 1+{\tilde{b}}x, \end{aligned}$$
(32)

which leads to the differential equation

$$\begin{aligned}{} & {} 4x\left( 1+{\tilde{b}}x\right) \left[ 1-2{\hat{\alpha }}{\tilde{b}}\right] \ddot{y}+\left[ 4\left( N-2\right) \left( 1-2{\hat{\alpha }}{\tilde{b}}\right) +2\left( 2N-3\right) {\tilde{b}}\left( 1-2{\hat{\alpha }}{\tilde{b}}\right) x \right] {\dot{y}}\nonumber \\{} & {} \quad +\left( N-2\right) {\tilde{b}}\left[ \left( N-3\right) -{\hat{\alpha }}\left( N-5\right) {\tilde{b}}\right] y=0. \end{aligned}$$
(33)

This equation has a solution in terms of hypergeometric functions. We can also show that (33) has solution in terms of simple power series and polynomials.

As \(x=0\) is a regular singular point of (33) we can show that (33) has the general solution represented by

$$\begin{aligned} y= & {} Ay_{1}+By_{2}, \end{aligned}$$
(34)

where A and B are constants. The functions \(y_{1}\) and \(y_{2}\) are given by

$$\begin{aligned} y_{1}= & {} c_{0}+\sum _{n=1}^{\infty }\frac{\left( -1\right) ^{n}}{4^{n}\left( 1-2{\hat{\alpha }}{\tilde{b}}\right) ^{n}n!\left( n-1+N-2\right) !} \nonumber \\ {}{} & {} \times \left. \bigg [2{\tilde{b}}\left( 1-2{\hat{\alpha }}{\tilde{b}}\right) \left( n-1\right) \left[ 2\left( n-2\right) +2N-3\right] \right. \nonumber \\ {}{} & {} \left. +\left( N-2\right) {\tilde{b}}\left( \left( N-3\right) -{\hat{\alpha }}\left( N-5\right) {\tilde{b}}\right) \bigg ]\right. ! \, c_{0}x^{n}, \end{aligned}$$
(35)

and

$$\begin{aligned} y_{2}= & {} {\tilde{Q}}y_{1}\ln x +\sum _{n=0}^{\infty }d_{n}x^{n+3-N}. \end{aligned}$$
(36)

The relevant details on the calculations for \(y_{1}\) and \(y_{2}\) are contained in Appendix A. Hence exact solutions in terms of simple power series exist.

We demonstrate that polynomial solutions are possible. For this to be true there must be a fixed natural number \(n=m\) such that the last bracketed term in (35) vanishes. This gives the condition

$$\begin{aligned}{} & {} \left[ -4{\hat{\alpha }}\left( m-1\right) \left( 2\left( m-2\right) +2N-3\right) -{\hat{\alpha }}\left( N-2\right) \left( N-5\right) \right] {\tilde{b}}^{2}\nonumber \\{} & {} \quad + \left[ 2\left( 2\left( m-2\right) +2N-3\right) \left( m-1\right) +\left( N-2\right) \left( N-3\right) \right] {\tilde{b}}=0, \end{aligned}$$
(37)

which is quadratic in \({\tilde{b}}\). This gives the root

$$\begin{aligned} {\tilde{b}}= & {} \frac{\left[ 2\left( m-1\right) \left( 2\left( m-2\right) +2N-3\right) +\left( N-2\right) \left( N-3\right) \right] }{{\hat{\alpha }}\left[ 4\left( m-1\right) \left( 2\left( m-2\right) +2N-3\right) +\left( N-2\right) \left( N-5\right) \right] }. \end{aligned}$$
(38)

For this value of \({\tilde{b}}\) we obtain

$$\begin{aligned} y_{1}= & {} c_{0}+\sum _{n=1}^{m-1}c_{n}x^{n}, \end{aligned}$$
(39)

for the first solution which is a polynomial of degree \((m-1)\). The elementary function (39), written in terms of a polynomial, will be helpful when discussing the physical features of the model (as we show in Sect. 6).

4 Canonical form

We can integrate (21) by writing it in canonical form. An equivalent form of (21) is given by

$$\begin{aligned}{} & {} \left[ 2x^3{\dot{y}}+ \left( N-3\right) x^2y+4{\hat{\alpha }}x^2{\dot{y}}-12{\hat{\alpha }}x^2{\dot{y}}Z + 2{\hat{\alpha }}\left( N-5\right) xy\left( 1-Z\right) \right] {\dot{Z}}\nonumber \\{} & {} \quad -{\hat{\alpha }} \left[ 8\left( N-3\right) x{\dot{y}}+8x^2\ddot{y}+ \left( N-4\right) \left( N-5\right) y\right] Z^{2}\nonumber \\{} & {} \quad + \left[ 4x^3\ddot{y}+ 4\left( N-2\right) x^2{\dot{y}}+ 8{\hat{\alpha }} \left( N-3\right) x{\dot{y}} \right. \nonumber \\{} & {} \quad \left. + 8{\hat{\alpha }}x^2 \ddot{y}+ \left( N-3\right) ^{2}xy + 2{\hat{\alpha }}\left( N-4\right) \left( N-5\right) y\right] Z\nonumber \\{} & {} \quad - \left( N-3\right) ^{2}xy -{\hat{\alpha }}\left( N-4\right) \left( N-5\right) y=0. \end{aligned}$$
(40)

This equation is a first order nonlinear ordinary differential equation in Z and is further identified as an Abel differential equation of the second kind. It can be simplified by applying a transformation given in Polyanin and Zaitsev [58]. We introduce the new variable

$$\begin{aligned} u=\left( Z-\frac{\left[ \left[ 2{\hat{\alpha }}\left( N-5\right) +\left( N-3\right) x\right] y+ 2x\left( x+2{\hat{\alpha }}\right) {\dot{y}} \right] }{2 {\hat{\alpha }}\left[ 6x{\dot{y}}+\left( N-5\right) y\right] }\right) U, \end{aligned}$$
(41)

where

$$\begin{aligned} U= & {} \exp \left( -\int \frac{\left[ -\left( N-4\right) \left( N-5\right) y-8\left( N-3\right) x{\dot{y}} -8x^{2}\ddot{y}\right] }{2x\left[ 6x{\dot{y}}+\left( N-5\right) y\right] } \ dx\right) , \end{aligned}$$
(42)

and \(N\ge 5\). Note that the transformation (41) holds when \({\hat{\alpha }}\ne 0\) and \(6x{\dot{y}}+\left( N-5\right) y\ne 0\).

Equation (40) then reduces to

$$\begin{aligned} u {\dot{u}}= & {} u {\mathscr {F}} + {\mathscr {G}}, \end{aligned}$$
(43)

where \(u= u(x)\) and we have introduced the new functions \({\mathscr {F}}\) and \({\mathscr {G}}\) which depend on the potential y and its derivatives. These are given by

$$\begin{aligned} {\mathscr {F}}= & {} -\frac{x\left[ y-2\left( x+2{\hat{\alpha }}\right) {\dot{y}}\right] \left[ \left( N-1\right) {\dot{y}}+2x\ddot{y}\right] U}{{\hat{\alpha }}\left[ 6x{\dot{y}}+\left( N-5\right) y\right] ^2}, \end{aligned}$$
(44)
$$\begin{aligned} {\mathscr {G}}= & {} \bigg (-\frac{{\hat{\alpha }}\left( N-4\right) \left( N-5\right) y+\left( N-3\right) ^{2}xy}{2{\hat{\alpha }}x\left[ 6x{\dot{y}}+\left( N-5\right) y\right] }\nonumber \\{} & {} +\frac{\left[ y\left( 2{\hat{\alpha }}\left( N-5\right) +\left( N-3\right) x\right) +2x{\dot{y}}\left( x+2{\hat{\alpha }}\right) \right] ^{2}}{8{\hat{\alpha }}^{2} x\left[ 6x{\dot{y}}+\left( N-5\right) y\right] ^{3}}\nonumber \\{} & {} \times \left[ -\left( N-4\right) \left( N-5\right) y-8x\left( \left( N-3\right) {\dot{y}}+x\ddot{y}\right) \right] \nonumber \\{} & {} +\frac{\big [y\left( 2{\hat{\alpha }}\left( N-5\right) +\left( N-3\right) x\right) +2x{\dot{y}}\left( x+2{\hat{\alpha }}\right) \big ]}{4{\hat{\alpha }}^{2} x\left[ 6x{\dot{y}}+\left( N-5\right) y\right] ^{2}}\nonumber \\ {}{} & {} \times \bigg [y\big (2{\hat{\alpha }}\left( N-4\right) \left( N-5\right) +\left( N-3\right) ^{2}x\big )+ 4x\big (\left[ \left( N-2\right) x\right. \nonumber \\{} & {} \left. +2{\hat{\alpha }}\left( N-3\right) \right] {\dot{y}} + x\ddot{y}\big [x+2{\hat{\alpha }}\big ]\big )\bigg ]\bigg )U^{2}. \end{aligned}$$
(45)

In order to find a solution for \(u = u(x)\), we must integrate (43). Equation (43) is in canonical form and is simpler than (40), but it remains an Abelian differential equation. Integration is a nontrivial task since \({\mathscr {F}}\) and \({\mathscr {G}}\) both depend on the metric function y and its derivatives in a complex manner. However we can demonstrate solutions by restricting the functions \({\mathscr {F}}\) and \({\mathscr {G}}\) in specific cases where a solution for u(x) is possible. These cases are provided below.

4.1 Case I: \({\mathscr {G}}=0\)

In this case we place the constraint

$$\begin{aligned} {\mathscr {G}}=0. \end{aligned}$$
(46)

With condition (46), the canonical form (43) reduces to

$$\begin{aligned} {\dot{u}}= & {} {\mathscr {F}}, \end{aligned}$$
(47)

which is a separable equation that can be integrated to obtain the potential

$$\begin{aligned} Z= & {} \left( \int -\frac{x\left[ y-2\left( x+2{\hat{\alpha }}\right) {\dot{y}}\right] \left[ \left( N-1\right) {\dot{y}}+2x\ddot{y}\right] U}{{\hat{\alpha }}\left[ 6x{\dot{y}}+\left( N-5\right) y\right] ^2} \, dx +C\right) \frac{1}{U}\nonumber \\{} & {} +\frac{ \left[ \left( N-3\right) x+2{\hat{\alpha }}\left( N-5\right) \right] y+ 2x\left( x+2{\hat{\alpha }}\right) {\dot{y}} }{2{\hat{\alpha }}\left[ 6x{\dot{y}}+\left( N-5\right) y\right] }. \end{aligned}$$
(48)

Therefore the metric potential Z is given explicitly and a form for y has to be chosen that satisfies the constraint (46).

4.2 Case II: \({\mathscr {F}}=0\)

We let

$$\begin{aligned} {\mathscr {F}}=0, \end{aligned}$$
(49)

which is equivalent to the constraint equation

$$\begin{aligned} \left[ y-2\left( x+2{\hat{\alpha }}\right) {\dot{y}}\right] \left[ \left( N-1\right) {\dot{y}}+2x\ddot{y}\right] =0. \end{aligned}$$
(50)

The above equation admits two solutions given by

$$\begin{aligned} y= & {} {\tilde{C}}\sqrt{x+2{\hat{\alpha }}}, \end{aligned}$$
(51)

and

$$\begin{aligned} y= & {} -\frac{2{\mathcal {C}}_{1}x^{-\frac{N-3}{2}}}{N-3}+{\mathcal {C}}_{2}, \end{aligned}$$
(52)

where \({\tilde{C}}\), \({\mathcal {C}}_{1}\) and \({\mathcal {C}}_{2}\) represent integration constants.

With condition (49), the canonical form (43) reduces to

$$\begin{aligned} u{\dot{u}}= & {} {\mathscr {G}}, \end{aligned}$$
(53)

which is a separable differential equation. Integrating we have the potential

$$\begin{aligned} Z= & {} \pm \left[ 2 \int \bigg (-\frac{{\hat{\alpha }}\left( N-4\right) \left( N-5\right) y+\left( N-3\right) ^{2}xy}{2{\hat{\alpha }}x\left[ 6x{\dot{y}}+\left( N-5\right) y\right] }\right. \nonumber \\{} & {} \left. +\frac{\left[ y\left( 2{\hat{\alpha }}\left( N-5\right) +\left( N-3\right) x\right) +2x{\dot{y}}\left( x+2{\hat{\alpha }}\right) \right] ^{2}}{8{\hat{\alpha }}^{2} x\left[ 6x{\dot{y}}+\left( N-5\right) y\right] ^{3}}\right. \nonumber \\{} & {} \left. \times \left[ -\left( N-4\right) \left( N-5\right) y-8x\left( \left( N-3\right) {\dot{y}}+x\ddot{y}\right) \right] \right. \nonumber \\{} & {} \left. +\frac{\big [y\left( 2{\hat{\alpha }}\left( N-5\right) +\left( N-3\right) x\right) +2x{\dot{y}}\left( x+2{\hat{\alpha }}\right) \big ]}{4{\hat{\alpha }}^{2} x\left[ 6x{\dot{y}}+\left( N-5\right) y\right] ^{2}}\right. \nonumber \\{} & {} \left. \times \bigg [y\big (2{\hat{\alpha }}\left( N-4\right) \left( N-5\right) +\left( N-3\right) ^{2}x\big )+ 4x\big (\left[ \left( N-2\right) x\right. \right. \nonumber \\{} & {} \left. \left. +2{\hat{\alpha }}\left( N-3\right) \right] {\dot{y}} + x\ddot{y}\big [x+2{\hat{\alpha }}\big ]\big )\bigg ]\bigg )U^{2} \, \, dx + 2C\right] ^{\frac{1}{2}} \frac{1}{U}\nonumber \\{} & {} +\frac{\left[ \left[ 2{\hat{\alpha }}\left( N-5\right) +\left( N-3\right) x\right] y+ 2x\left( x+2{\hat{\alpha }}\right) {\dot{y}} \right] }{2{\hat{\alpha }}\left[ 6x{\dot{y}}+\left( N-5\right) y\right] }. \end{aligned}$$
(54)

The integration in (54) can be performed explicitly. Firstly for the form in (51) we obtain the metric function

$$\begin{aligned} Z= & {} \pm \left. \bigg [2\left. \bigg [\left( x+2{\hat{\alpha }}\right) ^{1/3} \left( \left( N-2\right) x+2{\hat{\alpha }}\left( N-5\right) \right) ^{2/3} \left( N-1\right) x^{N-2}\right. \right. \nonumber \\{} & {} \left. \left. \times \left( \text{ Appell }F_{1}\left[ N-2, -\frac{1}{3}, -\frac{2}{3}, N-1, -\frac{x}{2{\hat{\alpha }}}, -\frac{\left( N-2\right) x}{2{\hat{\alpha }}\left( N-5\right) } \right] \right) \right] \right. \nonumber \\{} & {} \left. \times \left. \bigg [4{\hat{\alpha }}^2 \left[ 2+\frac{x}{2{\hat{\alpha }}}\right] ^{1/3} \left[ 2+\frac{\left( N-2\right) x}{2{\hat{\alpha }}\left( N-5\right) }\right] ^{2/3} \right] ^{-1} +2C\bigg ]\right. ^{\frac{1}{2}}\nonumber \\{} & {} \times \left( \frac{x^{\frac{N-4}{2}} \left( \left( N-2\right) x +2{\hat{\alpha }}\left( N-5\right) \right) ^{5/6}}{\left( x+2{\hat{\alpha }}\right) ^{1/3}}\right) ^{-1}+1+\frac{x}{2{\hat{\alpha }}}. \end{aligned}$$
(55)

As a result, the potential Z is defined in terms of special functions, namely Appell functions and elementary functions of x. Secondly for the form (52) we find the metric function

$$\begin{aligned} Z= & {} \pm \left. \bigg [2\left[ -\frac{256{\hat{\alpha }}}{\left( N-3\right) } {\mathcal {C}}_{1}{}^{ {3}} \left( N-1\right) \left( N-2\right) ^2 x^4 -256{\hat{\alpha }} \left( N-2\right) \left( N-5\right) {\mathcal {C}}_{1}{}^{ {2}}{\mathcal {C}}_{2}x^{\frac{N+5}{2}}\right. \right. \nonumber \\{} & {} \left. \left. +64{\hat{\alpha }}^2 {\mathcal {C}}_{1}{}^{{}2}{\mathcal {C}}_{2} \left( N-5\right) ^2x^{\frac{N+3}{2}}+{\mathcal {C}}_{2}{}^{{}3} (N-3)^3 (N-5)^2 x^{\frac{3N+1}{2}} \right. \right. \nonumber \\{} & {} \left. \left. -32{\hat{\alpha }} {\mathcal {C}}_{1}{\mathcal {C}}_{2}{}^{{}2}\left( N-3\right) \left( N-5\right) ^2 x^{N+1}\bigg ]\right. \left[ \frac{8{\hat{\alpha }}^{2}}{\left( N-3\right) } {\mathcal {C}}_{2} \left( N-5\right) ^2 x^{\frac{5}{2}} \right. \right. \nonumber \\ {}{} & {} \left. \left. \times \left( 4{\mathcal {C}}_{1}\left( N-2\right) x^{\frac{3}{2}}+\left( N-3\right) \left( N-5\right) {\mathcal {C}}_{2}x^{\frac{N}{2}}\right) \bigg ]\right. ^{-1} +2C\right] ^{\frac{1}{2}}\nonumber \\{} & {} \times \frac{x^{-\frac{N-4}{2}}}{ \sqrt{4\left( N-2\right) {\mathcal {C}}_{1}x^{\frac{N-3}{2}} + {\mathcal {C}}_{2}\left( N-3\right) \left( N-5\right) }}\nonumber \\{} & {} +\frac{\frac{8{\hat{\alpha }}}{\left( N-3\right) } {\mathcal {C}}_{1}x^{\frac{3}{2}}+{\mathcal {C}}_{2} x^{\frac{N}{2}}\left( \left( N-3\right) x+2{\hat{\alpha }}\left( N-5\right) \right) }{2\frac{{\hat{\alpha }}}{\left( N-3\right) }\left( 4{\mathcal {C}}_{1}\left( N-2\right) x^{\frac{3}{2}}+\left( N-3\right) \left( N-5\right) {\mathcal {C}}_{2}x^{\frac{N}{2}}\right) }. \end{aligned}$$
(56)

In this class of models we obtain two analytic forms for the function y in terms of elementary functions of x. The condition (50) is satisfied and as result two functional forms for the potential Z have been found.

4.3 Case III: \({\mathscr {F}}=K{\mathscr {G}}\)

We now consider a case where the functions \({\mathscr {G}}\) and \({\mathscr {F}}\) are proportional to each other. This yields the following constraint

$$\begin{aligned} {\mathscr {F}}=K{\mathscr {G}}, \end{aligned}$$
(57)

where K is some constant. With condition (57), the canonical form (43) becomes

$$\begin{aligned} u {\dot{u}}= & {} {\mathscr {G}}\left( K\hspace{2pt} u + 1\right) , \end{aligned}$$
(58)

which is a separable equation. We can integrate to obtain an equation containing the potential Z:

$$\begin{aligned}{} & {} \frac{1}{K}\left( Z-\frac{\left[ \left[ 2{\hat{\alpha }}\left( N-5\right) +\left( N-3\right) x\right] y+ 2x\left( x+2{\hat{\alpha }}\right) {\dot{y}} \right] }{2{\hat{\alpha }}\left[ 6x{\dot{y}}+\left( N-5\right) y\right] }\right) U\nonumber \\{} & {} \qquad - \frac{1}{K^2}\ln \left( 1 + K \hspace{2pt} \left( Z-\frac{\left[ \left[ 2{\hat{\alpha }}\left( N-5\right) +\left( N-3\right) x\right] y+ 2x\left( x+2{\hat{\alpha }}\right) {\dot{y}} \right] }{2{\hat{\alpha }}\left[ 6x{\dot{y}}+\left( N-5\right) y\right] }\right) U \right) \nonumber \\{} & {} \quad = \int \bigg (-\frac{{\hat{\alpha }}\left( N-4\right) \left( N-5\right) y+\left( N-3\right) ^{2}xy}{2{\hat{\alpha }}x\left[ 6x{\dot{y}}+\left( N-5\right) y\right] }\nonumber \\{} & {} \qquad +\frac{\left[ y\left( 2{\hat{\alpha }}\left( N-5\right) +\left( N-3\right) x\right) +2x{\dot{y}}\left( x+2{\hat{\alpha }}\right) \right] ^{2}}{8{\hat{\alpha }}^{2} x\left[ 6x{\dot{y}}+\left( N-5\right) y\right] ^{3}} \nonumber \\{} & {} \qquad \times \left[ -\left( N-4\right) \left( N-5\right) y-8x\left( \left( N-3\right) {\dot{y}}+x\ddot{y}\right) \right] \nonumber \\{} & {} \qquad +\frac{\big [y\left( 2{\hat{\alpha }}\left( N-5\right) +\left( N-3\right) x\right) +2x{\dot{y}}\left( x+2{\hat{\alpha }}\right) \big ]}{4{\hat{\alpha }}^{2} x\left[ 6x{\dot{y}}+\left( N-5\right) y\right] ^{2}} \nonumber \\{} & {} \qquad \times \bigg [y\big (2{\hat{\alpha }}\left( N-4\right) \left( N-5\right) +\left( N-3\right) ^{2}x\big )+ 4x\big (\left[ \left( N-2\right) x\right. \nonumber \\{} & {} \qquad \left. +2{\hat{\alpha }}\left( N-3\right) \right] {\dot{y}} + x\ddot{y}\big [x+2{\hat{\alpha }}\big ]\big )\bigg ]\bigg )U^{2}\,dx +C. \end{aligned}$$
(59)

Therefore the potential Z has an implicit form and a functional form for y has to be chosen that satisfies the constraint (57).

5 Exceptional cases

We now consider the case when the canonical form cannot be applied, that is, when \(\alpha =0\) (the Einstein limit) and \(6x{\dot{y}}+ (N-5)y=0\). Note that the term \(6x{\dot{y}}+(N-5)y\) appears in the denominator in (41) and therefore cannot be zero for the transformation to hold.

We first consider \(\alpha =0\), which is the Einstein case. Then the governing equation (40) reduces to the form

$$\begin{aligned}{} & {} \left[ 2x^2{\dot{y}}+ \left( N-3\right) xy\right] {\dot{Z}} \nonumber \\{} & {} \quad + \left[ 4x^2\ddot{y}+ 4\left( N-2\right) x{\dot{y}}+ \left( N-3\right) ^{2}y \right] Z- \left( N-3\right) ^{2}y=0, \end{aligned}$$
(60)

which is a first order linear differential equation in the variable Z if y is known. Several solutions to this equation were illustrated in four dimensions by Hansraj et al. [23]. We can show that the particular form \(y=\sqrt{x}\) in Eq. (60) yields the solution

$$\begin{aligned} Z= & {} \frac{\left( N-3\right) ^2}{\left( N-2\right) ^2}+ {\mathscr {C}}x^{2-N}, \end{aligned}$$
(61)

where \({\mathscr {C}}\) is an integration constant.

Secondly the case \(6x{\dot{y}}+ (N-5)y=0\) yields

$$\begin{aligned} y=Cx^{\frac{5-N}{6}}, \end{aligned}$$
(62)

where C is an integration constant. Then (40) becomes

$$\begin{aligned}{} & {} \left[ 6\left( N-2\right) x^2 + 12{\hat{\alpha }}\left( N-5\right) x\right] {\dot{Z}} + {\hat{\alpha }}\left( N-2\right) \left( N-5\right) Z^{2}\nonumber \\{} & {} \quad + \left[ \left( N+1\right) \left( N-5\right) x + 2{\hat{\alpha }}\left( N-5\right) \left( 4N-17\right) + 3\left( N^2-4N+7\right) x\right] Z \nonumber \\{} & {} \quad -9\left( N-3\right) ^2x -9{\hat{\alpha }}\left( N-4\right) \left( N-5\right) =0. \end{aligned}$$
(63)

Equation (63) is a Riccati equation which is difficult to solve in general.

For the particular spacetime dimension \(N=5\), expression (63) is a first order linear differential equation. It has the solution

$$\begin{aligned} Z= & {} \frac{{\tilde{C}}_{1}}{x^2}+1, \end{aligned}$$
(64)

where \({\tilde{C}}_{1}\) is an integration constant. When \(N\ne 5\) linearity is lost and we have a Riccati equation in Z. The solution is then given in terms of hypergeometric functions as can be verified using the software package Maple. We do not pursue this case further as special functions complicate a physical analysis.

6 Physical features

The simpler form of the exact solutions obtained, in the particular series solution, allows us to consider the physical viability of the models generated. For physical accessibility, we let \(B=0\) and \(A=1\) in (34). Equation (35) implies the form

$$\begin{aligned} y= & {} y_{1}=c_{0}\left[ 1+{\tilde{c}}x+{\tilde{d}}x^2+O(x^3)\right] , \end{aligned}$$
(65)

and

$$\begin{aligned} \frac{1}{y}= & {} \frac{1}{c_{0}}\left[ 1-{\tilde{c}}x+\left( {\tilde{c}}^2 -{\tilde{d}}\right) x^2+O(x^3)\right] , \end{aligned}$$
(66)

where the constants \({\tilde{c}}\) and \({\tilde{d}}\) are presented in Appendix B.

The forms \(Z=1+{\tilde{b}}x\) and y in (65) will allow physically viable solutions for the density and electromagnetic field variables in (20). These are expressed by

$$\begin{aligned} \kappa _{N}\rho= & {} B_{0}-\frac{1}{2}B_{1}x-\frac{1}{2}B_{2}x^2+O(x^3), \end{aligned}$$
(67a)
$$\begin{aligned} \kappa _{N}\left( \frac{E^{2}}{ {\mathcal {A}}_{N-2}}\right)= & {} B_{1}x+B_{2}x^2 +O(x^3), \end{aligned}$$
(67b)
$$\begin{aligned} \sigma ^2= & {} \frac{4}{({\mathcal {A}}_{N-2})^2}\left( D_{0}+D_{1}x + D_{2}x^2 +O(x^3)\right) , \end{aligned}$$
(67c)

where the coefficients \(B_{0}\), \(B_{1}\), \(B_{2}\), \(D_{0}\), \(D_{1}\) and \(D_{2}\) are provided in Appendix B.

We can observe that the energy density (67a), the electromagnetic field (67b) and the charge density (67c) are expressed explicitly in terms of a series in x, and they are given to second order in x. We observe that \(\rho \), \(E^{2}\) and \(\sigma ^{2}\) are finite at the centre. These quantities are regular and well behaved away from the centre. Varela et al. [59] indicate that the electric field E(x) should vanish at the centre to assure regularity of the proper charge density. This desirable feature of the vanishing of the electric field at the centre is satisfied in our model.

Charged dust spheres have been generated in general relativity. The dust sphere extends to a finite radius and at this point the interior is matched to the exterior Reissner-Nordström metric. A physically acceptable model of a charged dust sphere, satisfying the Einstein field equations and the boundary conditions, is given by Tikekar [22]. We demonstrate in this section that charged dust spheres arise in EGB gravity with interiors described by the exact solutions of this paper. At the matching surface, the exterior spacetime is the Boulware–Deser–Wiltshire metric.

The matching conditions across a comoving surface \(\Sigma \) are important for modelling relativistic stars as this will provide insights on the evolution and stellar structure of these astrophysical bodies. In EGB gravity, the boundary conditions on \(\Sigma \) have been recently derived in [60] extending the earlier work of Davis [61]. These are provided by

$$\begin{aligned}{} & {} \left( ds^{2}_{-}\right) _{\Sigma }=\left( ds^{2}_{+}\right) _{\Sigma }, \end{aligned}$$
(68a)
$$\begin{aligned}{} & {} \left[ K_{ab}-Kh_{ab}\right] ^{\pm }+2\alpha \left[ 3J_{ab}-Jh_{ab}+2{\hat{P}}_{a b cd}K^{bc}\right] ^{\pm }=0, \end{aligned}$$
(68b)

where \(K_{ab}\) is the extrinsic curvature. In the above we have

$$\begin{aligned} {\hat{P}}_{abcd}= & {} {\hat{R}}_{a b c d}+2{\hat{R}}_{b[c}h_{d]a}-2{\hat{R}}_{a[c}h_{d]b}+{\hat{R}}h_{a[c}h_{d]b}, \end{aligned}$$
(69)

where the caret “\({^{\hat{\,}}}\)” indicates quantities associated with the induced metric and \(P_{abcd}\) is the divergence-free part of the Riemann tensor. The tensor \(J_{ab}\) is defined by

$$\begin{aligned} J_{ab}= & {} \frac{1}{3}\left( 2KK_{ac}K^c{}_{ {}b}+K_{cd}K^{cd}K_{ab}\right. \nonumber \\{} & {} \left. -2K_{ac}K^{cd}K_{db}-K^{2}K_{ab}\right) , \end{aligned}$$
(70)

and J is its trace. As shown in [60] the equations in system (68) are satisfied if the following conditions

$$\begin{aligned} \left( ds^{2}_{-}\right) _{\Sigma }= & {} \left( ds^{2}_{+}\right) _{\Sigma }, \end{aligned}$$
(71a)
$$\begin{aligned} \left( K^{-}_{ab}\right) _{\Sigma }= & {} \left( K^{+}_{ab}\right) _{\Sigma }, \end{aligned}$$
(71b)

are valid. Note that (71) represent the matching conditions in general relativity. Consequently the Israel-Darmois conditions (71) imply that the Davis conditions (68) are satisfied.

We now consider a charged dust static distribution contained within a finite radius \({\mathcal {R}}\) for a particular solution found in this paper and utilize the boundary conditions for a dust matter distribution in EGB gravity. The higher dimensional interior static spacetime in terms of the variables y, Z and x is provided by

$$\begin{aligned} ds^{2}_{-}=-y^{2}(x) dt^2 + \frac{1}{4x Z(x)} dx^2 + x d \Omega ^{2}_{N-2}. \end{aligned}$$
(72)

For the exterior spacetime, we make use of the N-dimensional Boulware–Deser–Wiltshire exterior static metric. The Boulware–Deser–Wiltshire metric represents the charged counterpart of the Boulware–Deser metric was first studied by Wiltshire [62]. When \(\alpha \rightarrow 0\), the Reissner-Nordström spacetime is regained. Note that the Reissner-Nordström metric cannot be used to model the exterior spacetime in EGB gravity as it does not satisfy the EGB field equations (16). Brassel et al. [60] have demonstrated that the exterior spacetime has to be the Boulware–Deser–Wiltshire metric in a detailed study of the matching surface in EGB gravity. The Boulware–Deser–Wiltshire metric has the form

$$\begin{aligned} ds^{2}_{+}= & {} -F(r)dt^{2}+\frac{1}{F(r)}dr^{2}+r^{2}d\Omega ^{2}_{N-2}, \end{aligned}$$
(73)

where

$$\begin{aligned} F(r)= & {} 1+\frac{r^{2}}{2{\hat{\alpha }}}\left. \bigg (1 -\sqrt{1+4\alpha \left( N-4\right) \left( \frac{2M}{r^{N-1}}-\frac{Q^{2}}{\left( N-3\right) r^{2N-4}}\right) } \right) . \end{aligned}$$
(74)

In the above, M is the gravitational mass of the hypersphere and Q represents the charge contribution.

We match the particular interior charged dust solution given by (52) and (56) to the Boulware–Deser–Wiltshire metric (74) at the boundary \(r={\mathcal {R}}\). This gives the following expressions

$$\begin{aligned}{} & {} 1+\frac{{\mathcal {R}}^{2}}{2{\hat{\alpha }}}\left( 1 -\sqrt{1+4\alpha \left( N-4\right) \left( \frac{2M}{{\mathcal {R}}^{N-1}}-\frac{Q^{2}}{\left( N-3\right) {\mathcal {R}}^{2N-4}}\right) }\right) \nonumber \\{} & {} \quad =\left[ -\frac{2{\mathcal {C}}_{1}{\mathcal {R}}^{-\left( N-3\right) }}{N-3}+{\mathcal {C}}_{2}\right] ^{2}, \end{aligned}$$
(75a)
$$\begin{aligned}{} & {} \frac{1}{4{\mathcal {R}}^{2}}\left( 1+\frac{{\mathcal {R}}^{2}}{2{\hat{\alpha }}}\left( 1 -\sqrt{1+4\alpha \left( N-4\right) \left( \frac{2M}{{\mathcal {R}}^{N-1}}-\frac{Q^{2}}{\left( N-3\right) {\mathcal {R}}^{2N-4}}\right) }\right) \right) \nonumber \\{} & {} \quad =\left. \bigg [2\left[ -\frac{256{\hat{\alpha }}}{\left( N-3\right) } {\mathcal {C}}_{1}{}^{ {3}} \left( N-1\right) \left( N-2\right) ^2 {\mathcal {R}}^{8} -256{\hat{\alpha }} \left( N-2\right) \left( N-5\right) {\mathcal {C}}_{1}{}^{ {2}}{\mathcal {C}}_{2}{\mathcal {R}}^{N+5}\right. \right. \nonumber \\{} & {} \qquad \left. \left. +64{\hat{\alpha }}^2 {\mathcal {C}}_{1}{}^{{}2}{\mathcal {C}}_{2} \left( N-5\right) ^2{\mathcal {R}}^{N+3}+{\mathcal {C}}_{2}{}^{{}3} (N-3)^3 (N-5)^2 {\mathcal {R}}^{3N+1} \right. \right. \nonumber \\{} & {} \qquad \left. \left. -32{\hat{\alpha }} {\mathcal {C}}_{1}{\mathcal {C}}_{2}{}^{{}2}\left( N-3\right) \left( N-5\right) ^2 {\mathcal {R}}^{2(N+1)}\bigg ]\right. \left[ \frac{8{\hat{\alpha }}^{2}}{\left( N-3\right) } {\mathcal {C}}_{2} \left( N-5\right) ^2 {\mathcal {R}}^{5} \right. \right. \nonumber \\{} & {} \qquad \left. \left. \times \left( 4{\mathcal {C}}_{1}\left( N-2\right) {\mathcal {R}}^{3}+\left( N-3\right) \left( N-5\right) {\mathcal {C}}_{2}{\mathcal {R}}^{N}\right) \bigg ]\right. ^{-1} +2C\right] ^{\frac{1}{2}}\nonumber \\{} & {} \qquad \times \frac{{\mathcal {R}}^{-(N-4)}}{ \sqrt{4\left( N-2\right) {\mathcal {C}}_{1}{\mathcal {R}}^{N-3} + {\mathcal {C}}_{2}\left( N-3\right) \left( N-5\right) }}\nonumber \\{} & {} \qquad +\frac{\frac{8{\hat{\alpha }}}{\left( N-3\right) } {\mathcal {C}}_{1}{\mathcal {R}}^{3}+{\mathcal {C}}_{2} {\mathcal {R}}^{N}\left( \left( N-3\right) {\mathcal {R}}+2{\hat{\alpha }}\left( N-5\right) \right) }{2\frac{{\hat{\alpha }}}{\left( N-3\right) }\left( 4{\mathcal {C}}_{1}\left( N-2\right) {\mathcal {R}}^{3}+\left( N-3\right) \left( N-5\right) {\mathcal {C}}_{2}{\mathcal {R}}^{N}\right) }. \end{aligned}$$
(75b)

The free parameters in (75a)–(75b) are \({\hat{\alpha }}\), C, \({\mathcal {C}}_{1}\), \({\mathcal {C}}_{2}\) and \({\mathcal {R}}\). As these are algebraic equations there are sufficient free parameters for a real solution; this is a system of two equations in five unknowns. We observe that the dimension N appears explicitly in the boundary conditions and these are affected by Gauss–Bonnet coupling constant \({\hat{\alpha }}\). Hence the first boundary condition (68a) is satisfied since (75a)–(75b) hold. The second boundary condition (68b) is identically satisfied since the matter distribution is dust and the pressure is zero.

The total gravitational mass M at a finite radius \({\mathcal {R}}\) for the dust sphere then reduces to

$$\begin{aligned} M= & {} \left[ \frac{\left( N-3\right) }{2}\left. \bigg ({\mathcal {R}}^{N-3}\left( 1-Z({\mathcal {R}}^2)\right) +\frac{Q^{2}}{\left( N-3\right) ^{2}{\mathcal {R}}^{N-3}}\right. \right. \nonumber \\ {}{} & {} \left. \left. +\,{\hat{\alpha }}{\mathcal {R}}^{N-5}\left( 1-Z({\mathcal {R}}^2)\right) ^{2}\bigg )\right. \bigg ]\right. . \end{aligned}$$
(76)

The term \({\hat{\alpha }}{\mathcal {R}}^{N-5}\left( 1-Z({\mathcal {R}}^2)\right) ^{2}\) arises from the higher order curvature corrections in EGB gravity. The mass M is also clearly affected by the dimension N. The charge density at \(r={\mathcal {R}}\) is expressed by

$$\begin{aligned} \sigma= & {} \frac{ \sqrt{Z({\mathcal {R}}^{2})} \left[ {\mathcal {R}}^{N-2}E'({\mathcal {R}}^{2}) + \left( N-2\right) {\mathcal {R}}^{N-3}E({\mathcal {R}}^{2})\right] }{({\mathcal {A}}_{N-2}) \, {\mathcal {R}}^{\left( N-2\right) }}. \end{aligned}$$
(77)

The total charge Q within a radius \({\mathcal {R}}\) of the hypersphere is given by

$$\begin{aligned} Q= & {} \int _{0}^{{\mathcal {R}}} \left[ r^{N-2}E'(r^{2}) + \left( N-2\right) r^{N-3}E(r^{2})\right] dr, \end{aligned}$$
(78)

as measured by an external observer at infinity. In order to perform the integration in (78) we require the explicit expression for the electric field intensity E. This electric field intensity can be obtained by using (20c) and substituting for the potentials y, Z from (52) and (56) respectively. We observe that the electric field intensity depends on the potentials y, Z and its derivatives. The expression for the gravitational potential Z in (56) is given in a closed form but the integration in (78) is nontrivial. If necessary the integration can be performed numerically. The total charge Q remains finite. There is some simplification when \(N=5\). The charge Q is influenced by both the spacetime dimension N and the Gauss–Bonnet coupling constant \({\hat{\alpha }}\).

7 Feasible regions and graphical plots

A detailed analysis of the physical features of particular charged dust models is difficult for two reasons. Firstly the solutions in terms of special functions, such as hypergeometric and Appell functions in addition to elementary functions, arise in several cases which complicates a physical analysis. Secondly feasible regions for the existence of the metric functions y and Z may not be found which makes the existence of physically reasonable models questionable. To overcome these issues we follow the novel approach of Hansraj [53] to generate feasible spatial potential regions and plot the thermodynamical variables for various spacetime dimensions N.

We observe that the energy density can be written as \(\kappa _{N}\rho =\frac{\left( N-2\right) }{xy}\left[ 2Z{\dot{y}}-{\dot{Z}}y\right] \left[ x+\beta \left( 1-Z\right) \right] \) where \(\beta =2{\hat{\alpha }}\). If we assume \(x+\beta \left( 1-Z\right) >0\) the positivity of \(\rho \) requires

$$\begin{aligned} \frac{{\dot{y}}}{y}\ge \frac{{\dot{Z}}}{2Z}. \end{aligned}$$
(79)

In addition positivity of \(E^{2}\) provides the condition

$$\begin{aligned} \frac{{\dot{y}}}{y}\le \frac{\left( 1-Z\right) \left[ \left( N-3\right) +\frac{\beta \left( N-5\right) \left( 1-Z\right) }{2x}\right] }{4Z\left[ x+\beta \left( 1-Z\right) \right] }. \end{aligned}$$
(80)

We observe that (79) and (80) implies

$$\begin{aligned} \frac{{\dot{Z}}}{2Z}\le \frac{\left( 1-Z\right) \left[ \left( N-3\right) +\frac{\beta \left( N-5\right) \left( 1-Z\right) }{x}\right] }{4Z\left[ x+\beta \left( 1-Z\right) \right] }. \end{aligned}$$
(81)

The above can be rewritten as

$$\begin{aligned}{} & {} 2\left[ x+\beta \left( 1-Z\right) \right] x{\dot{Z}}-\frac{\beta }{2}\left( N-5\right) Z^{2}+\left[ \left( N-3\right) x+\beta \left( N-5\right) \right] Z\nonumber \\{} & {} \quad -\frac{\beta }{2}\left( N-5\right) -\left( N-3\right) x\le 0. \end{aligned}$$
(82)

In the case of equality (82) is an Abel differential equation of the second kind in Z which can be integrated to obtain the solution

$$\begin{aligned} Z= & {} \frac{x+\beta \pm \sqrt{x^{2}+\left( \beta ^{2}+2\beta C\right) x^{-\frac{(N-5)}{2}}}}{\beta }. \end{aligned}$$
(83)

This solution regains the result obtained in [53] when \(N=5\). For \(\rho >0\) and \(E^{2}\ge 0\), the potential Z must satisfy

$$\begin{aligned} 1+\frac{x}{\beta }\left( 1-\sqrt{1+\frac{\beta ^{2}+2\beta C}{x^{\frac{\left( N-1\right) }{2}}}}\right) \le Z\le 1+\frac{x}{\beta }\left( 1+\sqrt{1+\frac{\beta ^{2}+2\beta C}{x^{\frac{\left( N-1\right) }{2}}}}\right) . \end{aligned}$$
(84)

The above shows that the potential Z is constrained between the negative and positive branches of the N dimensional Boulware–Deser spacetime for uncharged fluids. Furthermore, the case \(Z=1+\frac{x}{\beta }\) represents the higher dimensional Schwarzschild interior metric.

We first consider the case when \(Z<1+\frac{x}{\beta }\). Then (84) becomes

$$\begin{aligned} 1+\frac{x}{\beta }\left( 1-\sqrt{1+\frac{\beta ^{2}+2\beta C}{x^{\frac{\left( N-1\right) }{2}}}}\right)< Z<1+\frac{x}{\beta }. \end{aligned}$$
(85)

This implies that for potentials found below the Schwarzschild interior metric, the lower bound is the negative branch of the N dimensional Boulware–Deser spacetime. Therefore the potentials that lie in this region yield physically acceptable charged dust models. We demonstrate these feasible regions of Z for (85) for the particular dimensions \(N=5,6,7\), in Figs. 1, 2 and 3 for \(C=0\) and \(\beta =1\). Then the case \(C=2\) and \(\beta =1\), the feasible regions are illustrated in Figs. 4, 5, and 6. Similarly we present the feasible regions for \(C=-2\) and \(\beta =1\) in Figs. 7, 8 and 9. The parameter values for C and \(\beta \) are chosen to coincide with the analysis in Hansraj [53]. Note that in the figures we have used the graphical illustrations

  • Dotted line..............: Higher dimensional Boulware–Deser \(+\) spacetime.

  • Solid line —————: Higher dimensional Schwarzschild spacetime.

  • Dashed lines \(- - - - -\): Higher dimensional Boulware–Deser − spacetime.

to describe the relevant spacetimes.

Next we consider the case when \(Z>1+\frac{x}{\beta }\), then the positivity of the energy density results in the condition positivity of the energy density requires

$$\begin{aligned} \frac{{\dot{y}}}{y}\le \frac{{\dot{Z}}}{2Z}, \end{aligned}$$
(86)

and \(E^{2}\ge 0\) gives

$$\begin{aligned} \frac{{\dot{Z}}}{2Z}\le \frac{\left( 1-Z\right) \left[ \left( N-3\right) +\frac{\beta \left( N-5\right) \left( 1-Z\right) }{x}\right] }{4Z\left[ x+\beta \left( 1-Z\right) \right] }. \end{aligned}$$
(87)

As in the above and due to \(Z>1+\frac{x}{\beta }\), we then obtain the constraint

$$\begin{aligned} 1+\frac{x}{\beta }< Z< 1+\frac{x}{\beta }\left( 1+\sqrt{1+\frac{\beta ^{2}+2\beta C}{x^{\frac{\left( N-1\right) }{2}}}}\right) . \end{aligned}$$
(88)

The potential Z must now lie between the Schwarzschild interior metric and the positive branch of the Boulware–Deser spacetime to achieve physically viable charged dust models. The upper bound for the metric potentials is then the higher dimensional positive branch of the Boulware–Deser metric. We illustrate the feasible spatial regions for (88) with \(N=5,6,7\) in Figs. 10, 11 and 12 for \(C=0\) and \(\beta =1\). The case \(C=2\) and \(\beta =1\), the feasible regions are illustrated in Figs. 13, 14 and 15. Similarly we present the feasible regions for \(C=-2\) and \(\beta =1\) in Figs. 16, 17 and 18. Again the values for C and \(\beta \) are selected to be consistent with the analysis in Hansraj [53].

Fig. 1
figure 1

5D Spatial region with \(C=0\) (\(\beta =1\))

Full size image
Fig. 2
figure 2

6D Spatial region with \(C=0\) (\(\beta =1\))

Full size image
Fig. 3
figure 3

7D Spatial region with \(C=0\) (\(\beta =1\))

Full size image
Fig. 4
figure 4

5D Spatial region with \(C=2\) (\(\beta =1\))

Full size image
Fig. 5
figure 5

6D Spatial region with \(C=2\) (\(\beta =1\))

Full size image
Fig. 6
figure 6

7D Spatial region with \(C=2\) (\(\beta =1\))

Full size image
Fig. 7
figure 7

5D Spatial region with \(C=-2\) (\(\beta =1\))

Full size image
Fig. 8
figure 8

6D Spatial region with \(C=-2\) (\(\beta =1\))

Full size image
Fig. 9
figure 9

7D Spatial region with \(C=-2\) (\(\beta =1\))

Full size image
Fig. 10
figure 10

5D Spatial region with \(C=0\) (\(\beta =1\))

Full size image
Fig. 11
figure 11

6D Spatial region with \(C=0\) (\(\beta =1\))

Full size image
Fig. 12
figure 12

7D Spatial region with \(C=0\) (\(\beta =1\))

Full size image
Fig. 13
figure 13

5D Spatial region with \(C=2\) (\(\beta =1\))

Full size image
Fig. 14
figure 14

6D Spatial region with \(C=2\) (\(\beta =1\))

Full size image
Fig. 15
figure 15

7D Spatial region with \(C=2\) (\(\beta =1\))

Full size image
Fig. 16
figure 16

5D Spatial region with \(C=-2\) (\(\beta =1\))

Full size image
Fig. 17
figure 17

6D Spatial region with \(C=-2\) (\(\beta =1\))

Full size image
Fig. 18
figure 18

7D Spatial region with \(C=-2\) (\(\beta =1\))

Full size image

We now consider the graphical behaviour of \(\rho \), \(E^{2}\) and \(\sigma ^2\). Again the earlier analysis of Hansraj [53] proves to be valuable in the investigation. Following [53] we choose the potential \(Z=\frac{x+\beta \pm \sqrt{x^{2}+\left( \beta ^{2}+2\beta C\right) x^{-\frac{(N-5)}{2}}}}{\beta }\) and \(y=\frac{1}{x^4}\). In our case Z depends on the spacetime dimension N. We plot the energy density, electric field and charge density for these potentials for various dimensions \(N=5,6,7,8\). We consider the parameter values \(C=1\) and \(\alpha =1/4\). These behaviours are illustrated in Figs. 19, 20 and 21. We observe that the energy density is well defined and positive. It is initially a maximum and decreases gradually, eventually approaching a stable finite value. The energy density is a steeper profile as the dimension N increases. The electric field intensity \(E^{2}\) is plotted in Fig. 20 for various dimensions. It is positive and well defined through out the fluid distribution. It is initially a maximum value and begins decreasing thereafter. A similar behaviour is observed for the charge density \(\sigma ^2\) plotted in Fig. 21.

Fig. 19
figure 19

Energy density plot for various dimensions N

Full size image
Fig. 20
figure 20

Electric field plot for various dimensions N

Full size image
Fig. 21
figure 21

Charge density plot for various dimensions N

Full size image

8 Discussion

We investigated the dynamics of a dust matter distribution in the presence of an electromagnetic field in a higher dimensional EGB gravity setting. The EGB field equations were generated for such a fluid distribution. The governing equation is an Abel differential equation of the second kind in the gravitational potential Z, a nonlinear first order differential equation. The master equation may also be considered as a second order linear differential equation in the metric potential y. The Gauss Bonnet constant \(\alpha \), the dimension N and the charge E have a significant impact on the gravitational behaviour of charged dust. In particular when \(\alpha =0\), we regain N dimensional general relativity, here the governing equation for charged dust is a linear differential equation in Z. We found exact solutions to the EGB field equations with the restrictions: \(y={\tilde{a}}\) and \(Z=1+{\tilde{b}}x\). In the latter case a series solution was demonstrated which contains polynomials for a specific value of \({\tilde{b}}\). In general, solutions that are found can be given in terms of special functions and elementary functions. The governing dynamical equation can be transformed to the canonical form \(u{\dot{u}}=u{\mathscr {F}}+{\mathscr {G}}\). We show that three classes of solutions are possible in this canonical case. We analysed the physical features close to the centre. The density, electric field and the proper charge density are finite at the centre and regular in a region containing the centre. This suggests that the model is well behaved. We have also shown the existence of charged dust spheres in EGB gravity in analogue to the dust spheres in general relativity. Our treatment shows that a systematic study of charged dust is possible in EGB gravity. There are some immediate extensions that can follow this treatment. Firstly we have focussed on spherical symmetry in our treatment in this paper. Following earlier treatments in general relativity we can consider other spacetimes geometries in EGB gravity such as axially symmetric models. New physical features will definitely arise if the spacetime geometry is not spherical. Furthermore it will be interesting to study charged dust in the general class of Lovelock gravity theories which extend EGB. The Lovelock tensor will be different from EGB gravity and we will obtain new and novel features arising from the additional curvature corrections.