1 Introduction

The Newtonian theory of stellar structure is usually used to describe the stars with not extremely high densities [1, 2]. Even for high density compact objects like white dwarfs and neutron stars, applying Newtonian gravity leads to acceptable results, comparable to the results of the relativistic models [3].

The contents of the star is often modeled by a perfect fluid with an equation of state relating the mass density \(\rho \) and pressure P. Among all possible choices of this equation, the polytrope equation of state, \(P=K\rho ^{\gamma }\) is of great importance for many astrophysical situations [1, 2, 4] where K and \(\gamma =1+\frac{1}{n}\) are both constants called the polytropic constant and polytropic exponent, respectively (n is called the polytropic index). The polytrope equation is the simplest equation of state useful for a wide range of fluid densities. Moreover, it leads to the scale symmetry and therefore enables introducing homology invariants [2].

On the other hand there are some phenomena leading to the local anisotropy of pressure inside a self gravitating system. For low density objects, the Newtonian approximation is valid. Some of the physical mechanisms for anisotropy in this regime are [5]:

  • When there is an anisotropic velocity distribution in a collisionless gas [6, 7], the radial pressure satisfies the Jeans equation:

    $$\begin{aligned} \frac{\mathrm{d}P_\mathrm{r}}{\mathrm{d}r}=\rho \frac{\mathrm{d}\phi }{\mathrm{d}r}+\frac{2}{r}\Delta (r) \end{aligned}$$
    (1)

    where \(\phi \) is the Newtonian gravitational potential of the fluid and we have considered a static and spherically symmetric distribution of matter. \(\Delta =P_\mathrm{t}-P_\mathrm{r}\), where \(P_\mathrm{r}\) and \(P_\mathrm{t}\) are the radial and tangential pressure, respectively. In this case \(\Delta \) measures the anisotropy of the velocity distributions:

    $$\begin{aligned} \Delta =\rho (\langle V^2_\mathrm{r}\rangle -\langle V^2_\upphi \rangle ) \end{aligned}$$
    (2)

    where \(\langle V^2_\mathrm{r}\rangle \) and \(\langle V^2_\upphi \rangle \) are the radial and azimuthal velocity dispersions. (Note that according to the spherical symmetry assumption we have \(\langle V^2_\upphi \rangle =\langle V^2_\theta \rangle \)). Equation (1) describes the hydrostatic equilibrium of a self gravitating object. To produce an anisotropic velocity distribution, consider a galactic halo of fermionic dark matter. The conservation of angular momentum of neutrinos streaming into the halo leads to an anisotropic pressure [8].

  • For a slowly rotating system [2, 5], the equation of hydrostatic equilibrium in the first order is given by Eq. (1) in which \(\Delta =-\frac{1}{3}\rho \omega ^2r^2\) where \(\omega \) is the angular velocity.

  • For a mixture of two non-interacting perfect fluids, the energy-momentum tensor is that of an anisotropic fluid in which \(\rho \), \(P_\mathrm{r}\), and \(P_\mathrm{t}\) are some functions of the mass densities and pressure of each component [9].

  • For a low mass charged white dwarf, the repulsive electrostatic force can be regarded as a source of anisotropy [10].

According to the above cases, even in the Newtonian regime, generally we are dealing with two components of pressure.

Here we shall focus on the Newtonian anisotropic polytropes [4, 11]. In the next section, first we derive the anisotropic version of dimensionless Newtonian hydrostatic equation, Lane–Emden equation and its boundary conditions. It is written in terms of physical quantities like the anisotropy function, in contrast to the form used in [4]. We extend the homology theorem for anisotropic polytropes and then in Sect. 3 we derive a new set of analytical solutions. There are some exact solutions of the Lane–Emden equation for special values of barotropic index in arbitrary spatial dimension [2].

In this paper, we use the ansatz that the form of the Lane–Emden equation does not change after introducing the anisotropy factor. First, a numerical analysis of the mass–radius relation of an anisotropic star is presented. We find the anisotropic modified form of the existing isotropic solutions by the above ansatz in Sect. 3. The authors of [12, 13], propose another ansatz on \(\Delta \) in the modeling of relativistic stars. This is used in [4] to generate some anisotropic Newtonian polytrope solutions from isotropic solutions. Other various assumptions on \(\Delta \) can be found in [11]. In Sect. 4, we focus on some integral theorems on the physical quantities of a star in the hydrostatic equilibrium obtained by Chandrasekhar in his book [1]. These are derived from Newtonian equations directly without any special assumption on the internal structure of the star or the equation of state. These theorems contain some inequalities containing the central pressure, the gravitational potential and the mean gravitational acceleration of a star. Here we extend these theorems to the general anisotropic equilibrium configuration. At the end, we present a perturbative treatment to study the effects of anisotropy in a star.

2 Anisotropic polytropes

The hydrostatic equilibrium of a star is governed by Eq. (1). The Newtonian potential is related to the density of the fluid by Poisson’s equation:

$$\begin{aligned} \nabla ^2 \phi =-4\pi G\rho . \end{aligned}$$
(3)

Considering the fluid obeys the polytropic equation of state along with Eq. (1), one gets

$$\begin{aligned} \phi (r)-\phi _0=K(n+1)\left( \rho ^{\frac{1}{n}}-\rho _0^{\frac{1}{n}}\right) -\int _0^r \frac{2\Delta (x)}{x\rho (x)}\mathrm{d}x \end{aligned}$$
(4)

where \(\rho _0\) and \(\phi _0\) denote the central density and gravitational potential. We take \(\phi _0=0\). Substituting Eq. (4) into Eq. (3) leads to the following expression which is the fundamental equation of equilibrium:

$$\begin{aligned} K(n+1)\nabla ^2 \rho ^{\frac{1}{n}}-\frac{1}{r^{N-1}}\frac{\mathrm{d}}{\mathrm{d}r}\left( r^{N-1}\frac{2\Delta (r)}{r\rho (r)}\right) =-4\pi G\rho \qquad \end{aligned}$$
(5)

where N is the dimension of space.

Introducing the new dimensionless variable \(\xi =\left( \pm \frac{(n+1)K}{4\pi G {\rho _0}^{1-\frac{1}{n}}}\right) ^{-\frac{1}{2}}r\), and with the following definitions:

$$\begin{aligned} \rho ={\rho }_0 \theta ^n, P_\mathrm{r}=P_0 {\theta }^{n+1}, \end{aligned}$$
(6)

Equation (5) takes the following form:

$$\begin{aligned} \theta ''+ & {} \frac{N-1}{\xi }\theta '\nonumber \\- & {} \frac{2}{P_0 (n+1)\xi \theta ^n}\left[ \Delta '+\frac{N-2}{\xi }\Delta -n\frac{\theta '}{\theta }\Delta \right] =\pm \theta ^n \end{aligned}$$
(7)

where a prime denotes differentiation with respect to the new radial coordinate, \(\xi \), and the plus and minus signs correspond to \(-\infty <n<-1\), \(-1<n<+\infty \), respectively. Note that \(\theta =\theta (r)\) is a dimensionless function. Equation (7) is the anisotropic version of the well-known Lane–Emden equation.

For \(n=-1\), the polytrope equation of state gives \(P_\mathrm{r}=K\) if \(\rho \ne 0\). Thus according to Eq. (4), the term involving \(\Delta \) only contributes to the gravitational potential and the Poisson equation gives

$$\begin{aligned} \Delta =\frac{2\pi G\rho }{r^{N-2}}\int _0^r x^{N-1} \rho (x)\mathrm{d}x. \end{aligned}$$
(8)

This shows that the case with \(n=-1\) is a special case, which must be treated separately. Whereas for the isotropic case \((\Delta =0)\), one gets \(\rho =0\) and therefore, \(n=-1\) is excluded from the further study for isotropic fluids.

The other important case arises when \(n=\pm \infty \) and therefore the polytrope equation of state reduces to

$$\begin{aligned} P_\mathrm{r}=K\rho . \end{aligned}$$
(9)

Replacing this in (1), we obtain

$$\begin{aligned} \phi -\phi _0=K\ln \left( \frac{\rho }{\rho _0}\right) -\int ^r_0\frac{2\Delta (x)}{\rho (x) x}\mathrm{d}x. \end{aligned}$$
(10)

The combination of the above expression when \(\phi _0=0\) and (3) gives

$$\begin{aligned} \theta ''+\frac{N-1}{\xi }\theta '+\frac{2}{P_0\xi }e^{\theta }\left( \Delta '+\frac{N-2}{\xi }\Delta +\Delta \theta '\right) =e^{-\theta }\nonumber \\ \end{aligned}$$
(11)

where \(\rho =\rho _0 e^{-\theta }\), \(P_\mathrm{r}=P_0 e^{-\theta }\), and \(\xi =\left( \frac{K}{4\pi G\rho _0}\right) ^{-\frac{1}{2}}r\).

Either the Lane–Emden equation (7) (for \(n\ne -1,\pm \infty \)) or (11) (for \(n=\pm \infty \)), contains two unknown functions, \(\theta (\xi )\) and \(\Delta (\xi )\). Thus we need another equation to obtain a closed system of equations.

Before describing our procedure to introduce the other equation, let us assume that \(\Delta \) is a given function. To obtain a unique solution for \(\theta \) one has to specify two boundary conditions. The first is \(\theta (0)=1\) (for \(n\ne -1, \pm \infty \)) and \(\theta (0)=0\) (for \(n=\pm \infty \)). In order to obtain the other condition, let us integrate (3):

$$\begin{aligned} \frac{\mathrm{d}\Phi }{\mathrm{d}r}=-r^{1-N}{\int _0}^r 4\pi G\rho x^{N-1}\mathrm{d}x. \end{aligned}$$
(12)

Combination of the above relation with (1) and (6) produces

$$\begin{aligned} \frac{\mathrm{d}\theta }{\mathrm{d}\xi }|_{\xi \rightarrow 0}=-\frac{\xi }{N}|_{\xi \rightarrow 0}+\frac{2}{P_0(n+1)}\lim _{\xi \rightarrow 0}\left( \frac{\Delta (\xi )}{\xi }\right) . \end{aligned}$$
(13)

A similar calculation for \(n=\pm \infty \) leads to

$$\begin{aligned} \frac{\mathrm{d}\theta }{\mathrm{d}\xi }|_{\xi \rightarrow 0}=\frac{\xi }{N}|_{\xi \rightarrow 0}-\frac{2}{P_0}\lim _{\xi \rightarrow 0}\frac{\Delta (\xi )}{\xi }. \end{aligned}$$
(14)

Equations (13) and (14) lead to the extensions of the Chandrasekhar’s theorem [1] for the anisotropic polytrope case.

Theorem 1

The finite solutions of anisotropic Lane–Emden equation at the origin have to satisfy \(\frac{\mathrm{d}\theta }{\mathrm{d}\xi }|_{\xi =0}=\frac{2}{P_0(n+1)}\lim _{\xi \rightarrow 0}\left( \frac{\Delta (\xi )}{\xi }\right) \) if \(n\ne -1, \pm \infty \) and \(\frac{\mathrm{d}\theta }{\mathrm{d}\xi }|_{\xi =0}=\frac{2}{P_0}\lim _{\xi \rightarrow 0}\frac{\Delta (\xi )}{\xi }\) if \(n=\pm \infty \).

Another important theorem is a generalization of the homology theorem [2].

Theorem 2

If the anisotropic Lane–Emden equation is satisfied by \(\theta (\xi )\) and \(\Delta (\xi )\) then \(A^{\frac{2}{n-1}} \theta (A\xi )\) and \(A^{\frac{2(n+1)}{1-n}}\Delta (A\xi )\) satisfy Eq. (7), and \(\theta (A\xi )-\ln A^2\) and \(\frac{1}{A^2}\Delta (A\xi )\) satisfy Eq. (11), where A is a constant.

The proof is straightforward. It can easily be obtained by substituting these expressions in (7) and (11). Thus we obtain a set of solutions parametrized by A.

One usual approach to solve (7) (or 11) with boundary conditions obtained in theorem I, is to suppose a special form for \(\Delta \). Then putting it into (7) (or 11) yields a differential equation for \(\theta \). If this procedure leads to smooth and non-negative pressures and density for the fluid, then they can be regarded as the physical solutions. However, knowing the function \(\Delta \), Eq. (7) (or 11) is a non-linear differential equation which in general one does not expect to obtain an analytical solution. Here we introduce a heuristic approach to get some analytical solutions. In order to do this, we assume that the anisotropy factor has no influence on the form of isotropic Lane–Emden equation. This means that \(\Delta \) modifies only the coefficients of isotropic Lane–Emden equation. Hereafter, we adopt the minus sign on the right-hand side of Lane–Emden equation (7) (or 11), since the plus sign leads to an imaginary radial coordinate [2]. In the following we shall consider two cases.

Case 1

\(\Delta \) modifies the coefficient of \(\theta '\) in the isotropic Lane–Emden equation. In this case Eq. (7) can be considered as two separated equations as below

$$\begin{aligned}&\frac{-2}{P_0(n+1)\theta ^n}\left[ \Delta '+\frac{N-2}{\xi }\Delta -n\frac{\theta '}{\theta }\Delta \right] =N_1\theta ' \end{aligned}$$
(15)
$$\begin{aligned}&\theta ''+\frac{N_1+N-1}{\xi }\theta '=-\theta ^n \end{aligned}$$
(16)

where \(N_1\) is an arbitrary constant. Equation (15) is a first order differential equation with the solution

$$\begin{aligned} \Delta (\xi )=-\frac{(n+1)P_0N_1\theta ^n}{2\xi ^{N-2}}\int \theta '\xi ^{N-2}\mathrm{d}\xi . \end{aligned}$$
(17)

Equations (16) and (17) with the boundary conditions (13) and \(\theta (0)=1\) give the solutions of anisotropic polytropes.

Case 2

\(\Delta \) modifies the coefficient of \(\theta ^n\) in the right-hand side of the isotropic Lane–Emden equation,

$$\begin{aligned}&\frac{-2}{P_0(n+1)\theta ^n\xi }\left[ \Delta '+\frac{N-2}{\xi }\Delta -n\frac{\theta '}{\theta }\Delta \right] =N_2\theta ^{n} \end{aligned}$$
(18)
$$\begin{aligned}&\frac{\mathrm{d}^2\theta }{\mathrm{d}\eta ^2}+\frac{N-1}{\eta }\frac{\mathrm{d}\theta }{\mathrm{d}\eta }=-\theta ^n \end{aligned}$$
(19)

where \(\eta =\sqrt{(N_2+1)}\,\ \xi \), \(N_2>-1\). The solution of Eq. (18) is

$$\begin{aligned} \Delta (\xi )=\frac{-(n+1)P_0N_2\theta ^n}{2(N_2+1)\xi ^{N-2}}\int \theta ^n \xi ^{N-1}\mathrm{d}\xi .\end{aligned}$$
(20)

Now let us to interpret physically Eq. (7). The quantities of physical interest, for definite values of the polytropic index, are the stellar radius and mass. These are

$$\begin{aligned} R=\left( \frac{K(n+1)}{4\pi G}\right) ^{\frac{1}{2}}\rho ^{\frac{1-n}{2n}}_0\xi _1, \end{aligned}$$
(21)
$$\begin{aligned} M={\int _0}^R 4\pi r^2 \rho \mathrm{d}r=4\pi \left( \frac{K(n+1)}{4\pi G}\right) ^{\frac{3}{2}}\rho _0 {\int _0}^{\xi _1} \theta ^n \xi ^2 \mathrm{d}\xi \nonumber \\ \end{aligned}$$
(22)

where \(\xi _1\) is the first root of \(\theta \) function and thus this point defines the surface of the star. Using the anisotropic Lane–Emden equation (7), we get

$$\begin{aligned} M= & {} -4\pi \left( \frac{K(n+1)}{4\pi G}\right) ^{\frac{3}{2}}\rho _0 {\int _0}^{\xi _1} \frac{\mathrm{d}}{\mathrm{d}\xi }\nonumber \\&\times \left( \xi ^2\frac{\mathrm{d}\theta }{\mathrm{d}\xi }-\frac{2}{P_0(n+1)}\Delta \xi \theta ^{-n}\right) \mathrm{d}\xi , \end{aligned}$$
(23)

which with the boundary conditions \(\theta _0=1\), Eq. (13), and the definition of surface \(\xi _1\) gives

$$\begin{aligned} M=4\pi \left( \frac{K(n+1)}{4\pi G}\right) ^{\frac{3}{2}}\rho _0\left| \xi ^2\frac{\mathrm{d}\theta }{\mathrm{d}\xi }\right| _{\xi _1}. \end{aligned}$$
(24)

This has the well-known form of the mass relation for isotropic star. Eliminating the central mass density between (21) and (24) we have

$$\begin{aligned} M R^{\frac{3-n}{n-1}}=4\pi \left( \frac{K(n+1)}{4\pi G}\right) ^{\frac{n}{n-1}}{\xi _1}^{\frac{n+1}{n-1}}\left| \frac{\mathrm{d}\theta }{\mathrm{d}\xi }\right| _{\xi _1}; \end{aligned}$$
(25)

knowing the polytropic index, \(\left| \frac{\mathrm{d}\theta }{\mathrm{d}\xi }\right| _{\xi _1}\) can be obtained from (16) or (19) numerically and thus the constant on the right-hand side of Eq. (25) is found. Therefore the anisotropy factor does not appear explicitly in the mass or the mass–radius relation of a star. The effect of anisotropy is included in the \(\theta \) function through its equation (parameter \(N_1\) in (16) and variable \(\eta \) in (19)). Figures 1 and 2 show the mass–radius diagrams for low density white dwarfs with \(n=\frac{3}{2}\) obtained from Eqs. (16) and (19), respectively.

Fig. 1
figure 1

Plot of \(\frac{M}{M_\odot }\) as a function of \(\frac{R}{R_\odot }\) in the Case 1, for \(n=\frac{3}{2}\) with \(N_1=-1.5\) (dotted line), \(N_1=0\) (thin line, the isotropic case), \(N_1=0.5\) (dashed line), \(N_1=1.5\) (thick line), and \(N_1=2.5\) (dot-dashed line)

Fig. 2
figure 2

Plot of \(\frac{M}{M_\odot }\) as a function of \(\frac{R}{R_\odot }\) in the Case 2, for \(n=\frac{3}{2}\) with \(N_2=-0.8\) (dotted line), \(N_2=0\) (thin line, the isotropic case), \(N_2=0.5\) (dashed line), \(N_2=1.5\) (thick line), and \(N_2=2.5\) (dot-dashed line)

The influence of anisotropy on the mass ratio for a fix star radius is numerically plotted in Figs. 3 and 4. Note that in the Case 1 the parameter \(N_1\) scales the anisotropy, while for the Case 2, the parameter \(N_2\) plays that role.

Fig. 3
figure 3

Plot of \(\frac{M}{M_\odot }\) as a function of \(N_1\) in the Case 1, for \(n=\frac{3}{2}\) and \(\frac{R}{R_\odot }=10\)

Fig. 4
figure 4

Plot of \(\frac{M}{M_\odot }\) as a function of \(N_2\) in the Case 2, for \(n=\frac{3}{2}\) and \(\frac{R}{R_\odot }=10\)

3 Exact analytical solutions of anisotropic Lane–Emden equation

3.1 Case I

3.1.1 \(n=0\)

With \(n=0\), the density is constant, \(\rho =\rho _0\), and Eq. (16) takes the following form:

$$\begin{aligned} \theta ''+\frac{S-1}{\xi }\theta '+1=0, \end{aligned}$$
(26)

in which \(S=N_1+N\), where solution is

$$\begin{aligned} \theta =\left\{ \begin{array}{ll} -\frac{\xi ^2}{2S}+1, &{} S > 1, \\ -\frac{\xi ^2}{2S}+A\frac{\xi ^{2-S}}{2-S}+1 , &{} S\le 1, S\ne 0, \end{array} \right. \end{aligned}$$
(27)

in which A is a constant. In the above solution we have also implemented the boundary condition \(\theta (0)=1\) and only those terms that lead to the non-singular solutions have been written.

We can now read \(\Delta \) from (17) considering the boundary condition (13):

$$\begin{aligned} \Delta =\left\{ \begin{array}{ll} \frac{(S-N)P_0}{2SN}\xi ^2 &{} S>1 ,\\ -\frac{(S-N)P_0}{2}\left( -\frac{\xi ^2}{NS}+\frac{A}{N-S }\xi ^{2-S}\right) &{} S\le 1, S\ne 0, S\ne N \\ 0 &{} N=S=1 .\\ \end{array} \right. \nonumber \\ \end{aligned}$$
(28)

We see that, to satisfy the boundary conditions, we have to have \(\Delta (0)=0\). The resulting \(\Delta \) function is plotted in Fig. 5.

Fig. 5
figure 5

Plot of \(\Delta \) for \(S>1\) in three dimensions with \(S=4\) (thick line), \(S=5\) (dotted line), \(S=6\) (dashed line), \(S=7\) (dot-dashed line)

Also the radial and tangential pressures are determined from the expressions (6), (27), and (28).

3.1.2 \(n=1\)

In this case, the differential equation (16) reduces to

$$\begin{aligned} \theta ''+\frac{(S-1)}{\xi ^2} \theta ' + \theta =0. \end{aligned}$$
(29)

The solutions of (29) are given in terms of Bessel functions as follows:

$$\begin{aligned} \theta (\xi )=\xi ^{-\frac{S-2}{2}}\left[ C_1 J_{\frac{S-2}{2}}(\xi )+C_2Y_{\frac{S-2}{2}}(\xi )\right] \end{aligned}$$
(30)

where \(C_1\) and \(C_2\) are integration constants. Imposing the condition \(\theta (0)=1\), \(C_2\) vanishes and we get

$$\begin{aligned} \theta (\xi )=\Gamma \left( \frac{S}{2}\right) \left( \frac{\xi }{2}\right) ^{-\frac{S-2}{2}}J_{\frac{S-2}{2}}(\xi ),S>2, \end{aligned}$$
(31)

in which the constant \(C_1\) is expressed in terms of the Gamma function, \(\Gamma \left( \frac{S}{2}\right) \). Substituting (31) into (17) leads to

$$\begin{aligned} \Delta (\xi )&=\frac{2^{\frac{S}{2}}}{8}(S-N) P_{0}\xi ^{3-\frac{s}{2}}J_{\frac{S-2}{2}}(\xi )\Gamma ^2\left( \frac{S}{2}\right) \Gamma \left( \frac{N}{2}\right) F\nonumber \\&\quad \times \left( \frac{N}{2};1+\frac{S}{2},1+\frac{N}{2};-\frac{\xi ^2}{4}\right) \end{aligned}$$
(32)

where F(abcx) is the hypergeometric function. These results also satisfy the boundary condition (13) with \(\theta '(0)=\Delta (0)=0\). The function (32) is depicted in Fig. 6.

Fig. 6
figure 6

Plot of \(\Delta \) for \(n=1\) in three dimensions with \(S=4\) (thick line), \(S=5\) (dotted line), \(S=6\) (dashed line), \(S=7\) (dot-dashed line)

3.1.3 \(n=\frac{S+2}{S-2}\)

In this case, if one takes

$$\begin{aligned} \theta (\xi )=\left( \frac{4}{(n-1)^2}\right) ^{n-1}\xi ^{\frac{2}{1-n}}z(\xi ),\xi =e^{-t}, \end{aligned}$$
(33)

then Eq. (16) reduces to

$$\begin{aligned} \frac{\mathrm{d}^2z}{\mathrm{d}t^2}=4\frac{z\mp z^n}{(n-1)^2}. \end{aligned}$$
(34)

The above equation can be integrated as follows:

$$\begin{aligned} \frac{1}{2}\left( \frac{\mathrm{d}z}{\mathrm{d}t}\right) ^2=\frac{4}{(n-1)^2}\left[ \frac{z^2}{2}\mp \frac{z^{n+1}}{n+1}\right] +C. \end{aligned}$$
(35)

Considering \(n>1\) (\(S>2\)), turning back to the original variables and using Eqs. (33) and (35), a simple calculation shows that the initial conditions on \(\theta \) leads to \(z = 0\) and \(\mathrm{d}z/\mathrm{d}t=0\) at \(\xi =0\). This gives \(C=0\). Writing the above equation in terms of the original variables, we have

$$\begin{aligned} \frac{2\theta \theta '}{n-1}+\xi \frac{\theta '^2}{2}+\xi \frac{\theta ^{n+1}}{n+1}=0. \end{aligned}$$
(36)

We proceed with integrating the above equation

$$\begin{aligned} \xi \frac{\theta ^{\frac{n+1}{2}}}{\theta '}=D .\end{aligned}$$
(37)

The constant D might be determined taking into account (13) and \(\theta (0)=1\)

$$\begin{aligned} \frac{1}{D}=\lim _{\xi \rightarrow 0}\theta ^{-\frac{S}{S-2}}\frac{\theta '}{\xi }=-\frac{1}{N}+\frac{2}{P_0(n+1)}\lim _{\xi \rightarrow 0}\frac{\Delta (\xi )}{\xi ^2}. \end{aligned}$$
(38)

As we see, the constant D depends on the \(\Delta \) function near the origin. The differential equation (37) has the following solution:

$$\begin{aligned} \theta =\left[ 1-\frac{\xi ^2}{D(S-2)}\right] ^{-\frac{S-2}{2}}. \end{aligned}$$
(39)

Employing (17), we get the following expression for the difference between the radial and tangential pressures:

$$\begin{aligned} \Delta (\xi )= & {} -\frac{S P_0 (S-N)}{ND(S-2)}\xi ^2\left[ 1-\frac{\xi ^2}{D(S-2)}\right] ^{-\frac{S+2}{2}}\nonumber \\&\quad \times F\left( \frac{N}{2},\frac{S}{2},1+\frac{N}{2},\frac{\xi ^2}{D(S-2)}\right) \end{aligned}$$
(40)

Now, we are in the position to fix the constant D. Combining (38) and (40), one gets \(D=-S\). The resulting \(\Delta \) is shown in Fig. 7.

Fig. 7
figure 7

Plot of \(\Delta \) for \(n=\frac{S+2}{S-2}\) in three dimensions with \(S=4\) (thick line), \(S=5\) (dotted line), \(S=6\) (dashed line), \(S=7\) (dot-dashed line)

3.2 Case II

3.2.1 \(n=0\)

The solutions of differential equations (19) and (20) satisfying the boundary conditions are

$$\begin{aligned} \theta =\left\{ \begin{array}{ll} -\frac{N_2+1}{2N}\xi ^2+A\sqrt{N_2+1}\xi +1, &{} N=1 ,\\ -\frac{N_2+1}{2N}\xi ^2+1 , &{} N > 1, \end{array} \right. \end{aligned}$$
(41)
$$\begin{aligned} \Delta =\left\{ \begin{array}{ll} -\frac{P_0 N_2}{2(N_2+1)N}\xi ^2, &{} N\ge 2, \\ -\frac{P_0N_2}{2(N_2+1)}(\xi ^2-\frac{\sqrt{(N_2+1)^3}}{N_2}A\xi ) , &{} N =1. \end{array} \right. \end{aligned}$$
(42)

The behavior of \(\Delta \) is shown in Fig. 8.

Fig. 8
figure 8

Plot of \(\Delta \) for \(n=0\) in three dimensions with \(N_2=1\) (thick line), \(N_2=2\) (dotted line), \(N_2=3\) (dashed line), \(N_2=4\) (dot-dashed line)

3.2.2 \(n=1\)

The situation for this case is also similar to the previous one and the substitution of S and \(\xi \) by N and \(\eta \) in (31) produces the solution of (19) with \(n=1\). This outcome with the integration (20) gives the differences between the radial and tangential pressures:

$$\begin{aligned} \Delta (\xi )= & {} -2^{N-2}\frac{P_0 N_2}{(N_2+1)^{\frac{1}{2}(N+1)}}\Gamma ^2\left( \frac{N}{2}\right) \xi ^{3-N}J_{\frac{N-2}{2}}\nonumber \\&\times (\sqrt{N_2+1}\xi )J_{\frac{N}{2}}(\sqrt{N_2+1}\xi ), \end{aligned}$$
(43)

which are depicted in Fig. 9.

Fig. 9
figure 9

Plot of \(\Delta \) for \(n=1\) in three dimensions with \(N_2=1\) (thick line), \(N_2=2\) (dotted line), \(N_2=3\) (dashed line), \(N_2=4\) (dot-dashed line)

3.2.3 \(n=\frac{N+2}{N-2}\)

Similar to the two former cases, in this situation the solution to the (19) with \(n=\frac{N+2}{N-2}\) might be obtained by changing S and \(\xi \) to N and \(\eta \) in (39),

$$\begin{aligned} \theta =\left[ 1-\frac{N_2+1}{D(N-2)}\xi ^2\right] ^{-\frac{N-2}{2}} ,N> 2. \end{aligned}$$
(44)

Putting the above result in (20) gives

$$\begin{aligned} \Delta =-\frac{NP_0N_2}{(N_2+1)(N-2)N}\left[ 1-\frac{1+N_2}{D(N-2)}\xi ^2\right] ^{-(N+1)}\xi ^2.\nonumber \\ \end{aligned}$$
(45)

Combining (38) and (45), one gets \(D=-\frac{N(N_2+1)}{2N_2-1}\). The graph of (45) is plotted in Fig. 10 for various values of \(N_2\).

Fig. 10
figure 10

Plot of \(\Delta \) for \(n=\frac{N+2}{N-2}\) in three dimensions with \(N_2=1\) (thick line), \(N_2=2\) (dotted line), \(N_2=3\) (dashed line), \(N_2=4\) (dot-dashed line)

4 Integral theorems

Chandrasekhar in his famous book [1] has discussed some inequalities for the physical quantities describing a star in the Newtonian gravitational equilibrium. In this section we will extend some of his results for anisotropic stars. The results are general and are not restricted to the polytropic case.

Let m(r) be the mass contained inside the radius r, then

$$\begin{aligned} m'=4\pi \rho r^2. \end{aligned}$$
(46)

Equation (1) can now be written as

$$\begin{aligned} P_\mathrm{r}'=-\frac{\mathrm{Gm}\rho }{r^2}+\frac{2\Delta }{r} \end{aligned}$$
(47)

where we have used the definition of gravitational potential. From (46) and (47) one can get the following equation:

$$\begin{aligned} \frac{1}{r^2}\left[ \frac{r^2}{\rho }\left( P_\mathrm{r}'-\frac{2\Delta }{r}\right) \right] ' =-4\pi G\rho . \end{aligned}$$
(48)

These equations can be used to prove the following theorems.

Theorem 3

For any equilibrium configuration the function

$$\begin{aligned} I=P_\mathrm{r}+ \frac{\mathrm{Gm}^2}{8\pi r^4}-2\int _0^r\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}} \end{aligned}$$
(49)

does not increase outward.

Proof

Let us calculate

$$\begin{aligned} I'=P_\mathrm{r}'+\frac{\mathrm{Gmm}'}{4\pi r^4}-\frac{\mathrm{Gm}^2}{2\pi r^5}-\frac{2\Delta }{r}. \end{aligned}$$
(50)

Then by (47)

$$\begin{aligned} I'=-\frac{\mathrm{Gm}^2}{2\pi r^5}\le 0. \end{aligned}$$
(51)

\(\square \)

Corollary 1

For the central pressure, we have

$$\begin{aligned} P_0>P_\mathrm{r}+ \frac{\mathrm{Gm}^2}{8\pi r^4}-2\int _0^r\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}}>\frac{\mathrm{Gm}^2}{8\pi R^4}-2\int _0^R\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}},\nonumber \\ \end{aligned}$$
(52)

in which \(M=m(R)\) is the mass of the star. The last term in the right-hand side shows a lower bound on \(P_0\). Moreover, for any arbitrary radius, the above relation leads to

$$\begin{aligned} P_\mathrm{r}>\frac{G}{8\pi }\left( \frac{M^2}{R^4}- \frac{m^2}{r^4}\right) -2\int _r^R\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}}. \end{aligned}$$
(53)

Theorem 4

For any equilibrium configuration

$$\begin{aligned} I_\nu \equiv \int _0^R \mathrm{d}m \frac{\mathrm{Gm}}{\tilde{r}^\nu }=4\pi \int _0^R \mathrm{d}\tilde{r}\tilde{r}^{3-\nu }[(4-\nu )P_\mathrm{r}+2\Delta ] \end{aligned}$$
(54)

if \(\nu <4\).

Proof

From (47) we have

$$\begin{aligned} \frac{Gmm'}{4\pi r^4}=\frac{2\Delta }{r}-P_\mathrm{r}'. \end{aligned}$$
(55)

Multiplying this by \(r^{4-\nu }\) and then integrating from r to R gives

$$\begin{aligned} I_\nu =8\pi \int _0^R \mathrm{d}\tilde{r}\tilde{r}^{\nu -3}\Delta -4\pi \int _0^R \mathrm{d}P_\mathrm{r}\tilde{r}^{4-\nu }. \end{aligned}$$
(56)

Integrating by parts the first integral, leads to (54). \(\square \)

For \(\nu =4\), Eq. (56) gives

$$\begin{aligned} I_4=8\pi \int _0^R\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}}+4\pi P_0, \end{aligned}$$
(57)

which also has a lower bound equal to \(\frac{\mathrm{Gm}^2}{2R^4}\) according to (52).

The gravitational potential energy \(\Omega \) of the configuration is \(-I_1\). According to (56) it is given by

$$\begin{aligned} \Omega =-\int _0^R \mathrm{d}V (3P_\mathrm{r}+2\Delta ). \end{aligned}$$
(58)

Denoting the mean value of gravitational acceleration by \(\bar{g}\), we have

$$\begin{aligned} M\bar{g}=\int _0^R \mathrm{d}m\frac{\mathrm{Gm}}{r^2}=I_2 \end{aligned}$$
(59)

and hence by (56)

$$\begin{aligned} M\bar{g}=8\pi \int _0^R \mathrm{d}\tilde{r}\tilde{r}\left( P_\mathrm{r}+\Delta \right) =8\pi \int _0^R \mathrm{d}\tilde{r}\tilde{r}P_\mathrm{t}. \end{aligned}$$
(60)

This shows that for an anisotropic star only the tangential pressure contributes to finding the mean value of acceleration.

Theorem 5

For any equilibrium configuration

$$\begin{aligned}&\pi \nu P_0R^{4-\nu }+\frac{(4-\nu ) \mathrm{Gm}^2}{8 R^\nu }\nonumber \\&\quad +2\pi \nu (4-\nu )\int _0^R \mathrm{d}\tilde{r}\tilde{r}^{3-\nu }\int _0^R \mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}}\nonumber \\&\quad >2\pi \nu \int _0^R \mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}^{\nu -3}}+I_\nu >\frac{\mathrm{Gm}^2}{2R^\nu }. \end{aligned}$$
(61)

Proof

By Theorem 3

$$\begin{aligned}&\frac{\mathrm{Gm}^2}{8\pi R^4}-\frac{\mathrm{Gm}^2}{8\pi r^4}-2\int _0^R\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}}+2\int _0^r\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}}<P_\mathrm{r}<P_0\nonumber \\&\quad -\frac{\mathrm{Gm}^2}{8\pi r^4}+2\int _0^r\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}} \end{aligned}$$
(62)

and by Theorem 4

$$\begin{aligned}&4\pi (4-\nu )\int _0^R\mathrm{d}\tilde{r}\tilde{r}^{3-\nu }\left[ \frac{\mathrm{Gm}^2}{8\pi R^4}-\frac{\mathrm{Gm}^2}{8\pi r^4}-2\int _0^R\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}}\right. \\&\quad \left. +2\int _0^r\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}}\right] \end{aligned}$$
$$\begin{aligned}&\quad <I_\nu +8\pi \int _0^R \mathrm{d}\tilde{r}\tilde{r}^{3-\nu }\Delta <4\pi (4-\nu )\int _0^R\mathrm{d}\tilde{r}\tilde{r}^{3-\nu }\nonumber \\&\quad \times \left[ P_0-\frac{\mathrm{Gm}^2}{8\pi r^4}+2\int _0^r\mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}}\right] . \end{aligned}$$
(63)

This inequality can be written as

$$\begin{aligned}&4\pi P_0 R^{4-\nu }+8\pi \nu (4-\nu )\int _0^R \mathrm{d}\tilde{r}\tilde{r}^{3-\nu }\int _0^r \mathrm{d}\tilde{r}\frac{\Delta }{\tilde{r}} \end{aligned}$$
$$\begin{aligned}&>I_\nu +8\pi \int _0^R \mathrm{d}\tilde{r}\tilde{r}^{3-\nu }\Delta +\frac{(4-\nu )}{2}\int _0^R \mathrm{d}\tilde{r} \frac{\mathrm{Gm}^2}{\tilde{r}^{\nu +1}}\nonumber \\&>\frac{\mathrm{Gm}^2}{2R^\nu }. \end{aligned}$$
(64)

Inserting

$$\begin{aligned} \int _0^R \mathrm{d}\tilde{r} \frac{\mathrm{Gm}^2}{\tilde{r}^{\nu +1}}=\frac{\mathrm{Gm}^2}{R^\nu }-2I_\nu \end{aligned}$$
(65)

in (64) and simplifying it, we get (61). \(\square \)

Corollary 2

Setting \(\nu =1\), we have

$$\begin{aligned}&\left( \pi P_0 R^3+\frac{3\mathrm{Gm}^2}{8R}+6\pi \int _0^R d\tilde{r}\tilde{r}^2\int _0^r d\tilde{r}\frac{\Delta }{\tilde{r}}\right) \nonumber \\&\quad >2\pi \int _0^R d\tilde{r}\tilde{r}^2\Delta -\Omega >\frac{\mathrm{Gm}^2}{2R}, \end{aligned}$$
(66)

which shows the upper and lower bounds of the potential energy, \(-\Omega \).

5 Concluding remarks

In this paper we discussed how the anisotropy factor modifies the Lane–Emden equation and homology theorem. We obtained some theorems governing the characteristic functions of anisotropic star. These are the extension of Chandrasekhar’s theorems. We performed a procedure to find the anisotropy factor, the radial pressure and the density functions exactly satisfying the Lane–Emden equation. We had two cases. In the first case the effect of anisotropy is to change the dimension of space from N to \(N+N_1\) (which may be a non-integer number). For the second case, anisotropy shows itself as a rescaling of the radial coordinate in the function \(\theta \).

Moreover, it is straightforward to find some approximate analytical solutions of the polytropic stars with anisotropic pressure. For example, let us assume that the including anisotropy factor is equivalent to a slight modification of the polytropic index from its value in the absence of anisotropy, \(n_0\). This means that if an exact solution \(\theta _0\) of the isotropic Lane–Emden equation is known:

$$\begin{aligned} \theta _0''+\frac{N-1}{\xi }\theta '_0=-\theta _0^{n_0}. \end{aligned}$$
(67)

Inclusion of the anisotropy factor leads to

$$\begin{aligned}&\theta ''+\frac{N-1}{\xi }\theta '-\frac{2}{P_0 (n_0+1)\xi \theta ^{n_0}}\nonumber \\&\quad \times \left[ \Delta '+\frac{1}{\xi }\Delta -n_0\frac{\theta '}{\theta }\Delta \right] =-\theta ^{n_0}, \end{aligned}$$
(68)

which is equivalent to perturbing \(\theta _0\) and \(n_0\) in Eq. (67) in such a way that \(\theta =\theta _0+\epsilon \theta _*\) and \(n=n_0+\epsilon \) satisfy (67)

$$\begin{aligned} \theta ''+\frac{N-1}{\xi }\theta '=-\theta ^n, \end{aligned}$$
(69)

where \(\vert \epsilon \vert \ll 1\). Thus we have the following equation for \(\theta _*\) and \(\Delta \):

$$\begin{aligned} \theta _*''+\frac{N-1}{\xi }\theta _* '=-\theta _0^{n_0}\left( \ln {\theta _0}+\frac{n_0\theta _*}{\theta _0}\right) , \end{aligned}$$
(70)
$$\begin{aligned} \frac{2}{P_0 (n_0+1)\xi \theta ^{n_0}}\left[ \Delta '+\frac{1}{\xi }\Delta -n_0\frac{\theta '}{\theta }\Delta \right] =-\theta ^{n}+\theta ^{n_0}.\nonumber \\ \end{aligned}$$
(71)

Expanding the right-hand side of (71) and assuming \(\Delta =\epsilon \Delta _*\), the following equation for \(\Delta _*\) is obtained:

$$\begin{aligned}&\frac{2}{P_0 (n_0+1)\xi \theta _0^{n_0}}\left[ \Delta _*'+\frac{1}{\xi }\Delta _*-n_0\frac{\theta _0'}{\theta _0}\Delta _*\right] \nonumber \\&\quad =-\theta _0^{n_0}(\ln {\theta _0}+\frac{n_0\theta _*}{\theta _0}). \end{aligned}$$
(72)

Having \(n_0\), the functions \(\theta _0\) and thus \(\theta _*\) can be determined from (67) and (70) and then \(\Delta _*\) is obtained from (72). With the initial conditions \(\theta _0(0) = 1\), \(\theta _0'(0)=\theta _*(0)=\theta _*'(0) = 0\) and for the case \(n_ 0 = 0\) and \(N=3\), we get \(\theta _ 0 = 1 - \xi ^2 /6\). Inserting \(\theta _ 0\) into (69), we find that [14]

$$\begin{aligned} \theta _* = (3 - \frac{\xi ^2}{6}) \ln {(1 - \frac{\xi ^2}{6})} + \frac{2\sqrt{6}}{\xi } \ln {\frac{\sqrt{6} + \xi }{\sqrt{6} - \xi }} + \frac{5}{18}\xi ^ 2 - 4\nonumber \\ \end{aligned}$$
(73)

where \(\xi \le \sqrt{6}\). Thus:

$$\begin{aligned} \Delta= & {} -\frac{P_0 \epsilon }{2}\left[ \frac{\xi ^2}{9} [3\ln {(1 - \frac{\xi ^2}{6})-2}]\right. \nonumber \\&\left. + \frac{4\sqrt{6}}{\xi } \arctan {\frac{\xi }{\sqrt{6}}} - 4\right] . \end{aligned}$$
(74)

The plot of \(\Delta \) is shown in Fig. 11.

Fig. 11
figure 11

Plot of \(\Delta \) for the analytical solution close to the exact solution (\(n_0=0\), \(N=3\)) as explained in the text

One can find other solutions of (67), (70), and (72) numerically for arbitrary values of the barotropic index. For example, for \(n_0=1\), the graph of \(\Delta \) is plotted in Fig. 12.

Fig. 12
figure 12

Plot of \(\Delta \) for the numerical solution close to the exact solution (\(n_0=1\), \(N=3\)) as explained in the text