Examples of Numerical Solutions

Let us now describe in detail our strategy for the numerical solution of Equation (37). We will consider this to be a system of two coupled first-order differential equations for the functions $\gamma (\xi )$ and $\pi (\xi )$, leading to a pair of numerical development equations, which can be easily obtained, and that are given by

$$\delta\gamma = \pi\,\delta\xi,$$ (70)

$$\delta\pi = \pi\, \frac {4(1+1/n)F\xi-[1+(5+4/n)F]\gamma-F(1+F)\xi\pi} {2(1+1/n)F\xi(\xi-\gamma)}\, \delta\xi,$$ (71)

where $F$ is the value of the auxiliary function $F(\xi,\pi)$, where $\delta\xi$ is a small increment of the variable $\xi$ and where $\delta\gamma$ and $\delta\pi$ are the corresponding increments of the functions $\gamma (\xi )$ and $\pi (\xi )$. This system can be integrated by the use of the Runge-Kutta fourth-order algorithm in a straightforward way, so long as one has a sensible starting point for the integration process.

We must now discuss where to start the integration process. There are two alternatives that we have used. The first alternative is to start at the inflection point $\xi_{e}$ of the function $\gamma (\xi )$. As we saw in Section 3, Equation (37) itself provides us with a fixed relation between $\gamma(\xi_{e})$ and $\pi(\xi_{e})$ at this point. This means that, given a value of $\pi(\xi_{e})$ at this point, which will have the role of a free parameter of the system, we have the corresponding value of $\gamma(\xi_{e})$, which then allows us to start the integration process. The second alternative is to start at the root $\xi_{z}$ of $\gamma (\xi )$, which we know necessarily exists and is necessarily located within the matter region. Once again the value of $\pi(\xi_{z})$ at this point will have the role of a free parameter of the system, and in this case we just start with $\gamma(\xi_{z})=0$. The first alternative seems to work better for the larger initial values of $\pi(\xi_{e})$, which corresponds to the more dense objects, and the second alternative seems to work better for the smaller initial values of $\pi(\xi_{z})$, which corresponds to the less dense objects.

In either case, we must start with a choice for the initial value of $\xi$. Taking the position $\xi_{e}=r_{e}/r_{0}$ of the inflection point as the example of a starting point for this discussion, we note that a choice of value for $\xi_{e}$ is tantamount to a choice of a relation between $r_{e}$ and $r_{0}$, without implying the choice of a definite value for either one. Therefore the choice of the initial value of $\xi$ can be quite arbitrary, without leading to any loss of generality. For example, if we start with $\xi_{e}=1$, this only means that we choose the arbitrary parameter $r_{0}=r_{e}$ to be the position of the inflection point, without actually choosing a value for $r_{e}$. Therefore, a solution found in this way is not just a single solution, but a one-parameter family of solutions, parametrized by the values of $r_{0}$. A similar argument can be made in the case in which we start at the position $\xi_{z}=r_{z}/r_{0}$ of the zero of $\gamma (\xi )$.

Note that the physical constants $K$, $\kappa$ and $r_{0}$ do not appear individually in either the numerical propagation equations or the initial values. They appear only as the combination seen in the definition of the dimensionless constant $C$, which in turn appears only in the expression for $F(\xi,\pi)$,

$$F(\xi,\pi) = C\left(\frac{\pi}{\xi^{2}}\right)^{1/n},$$ (72)

$$C = \frac{K}{\left(\kappa r_{0}^{2}\right)^{1/n}}.$$ (73)

We will therefore adopt $C$ as one of our free parameters for the numerical work, the others being $n$ and the initial value of $\pi (\xi )$, which will then result in a certain value of $\xi_{M}$. Note that the parameters $\xi_{M}$ that contains the mass $M$ also does not appear in either the equation or the initial values. It appears only as the asymptotic boundary condition for $\gamma (\xi )$, which is obtained only at the end of the integration process in the outward direction.

Figure 1: Graph of the functions $\gamma (\xi )$, $\pi (\xi )$, $\pi '(\xi )$ and $\rho (\xi )$ in a typical case, with the parameters $C=0.05$, $n=3/2$ and $\pi _{e}=1.0$, resulting in $\xi _{M}\approx 0.5392$ and $\xi _{\mu }\approx 0.1668$. The vertical dotted lines mark the positions of the matter-vacuum interfaces.

Figure 2: Graph of the functions $\nu (\xi )$ and $\lambda (\xi )$ in a typical case, with the parameters $C=0.05$, $n=3/2$ and $\pi _{e}=1.0$, resulting in $\xi _{M}\approx 0.5392$ and $\xi _{\mu }\approx 0.1668$. The vertical dotted lines mark the positions of the matter-vacuum interfaces.

Figure 3: Logarithmic graph depicting $\pi (\xi )$ in the approach to the inner interface from within the matter region. That approach proceeds to the left in the graph. Small error bars are barely visible at the left end of the graph. The line seen was obtained by linear regression from the data points shown.

Figure 4: Graph of the functions $\gamma (\xi )$, $\pi (\xi )$ and $\rho (\xi )$ in a high-density case, with the parameters $C=0.01$, $n=3$ and $\pi _{e}=100.0$, resulting in $\xi _{M}\approx 1.0096$ and $\xi _{\mu }\approx 15.0838$. The vertical dotted lines mark the positions of the matter-vacuum interfaces.

Figure 5: Graph of the functions $\nu (\xi )$ and $\lambda (\xi )$ in a high-density case, with the parameters $C=0.01$, $n=3$ and $\pi _{e}=100.0$, resulting in $\xi _{M}\approx 1.0096$ and $\xi _{\mu }\approx 15.0838$. The vertical dotted lines mark the positions of the matter-vacuum interfaces.

Figure 6: Graph of the functions $\gamma (\xi )$ and $\pi (\xi )$ in a low-density case, with the parameters $C=1/3$, $n=3/2$ and $\pi _{z}=10^{-3}$, resulting in $\xi _{M}\approx 0.1848$ and $\xi _{\mu }\approx 3.3039\times 10^{-4}$. The vertical dotted lines mark the positions of the matter-vacuum interfaces.

Figure 7: Graph of the functions $\nu (\xi )$ and $\lambda (\xi )$ in a low-density case, with the parameters $C=1/3$, $n=3/2$ and $\pi _{e}=10^{-3}$, resulting in $\xi _{M}\approx 0.1848$ and $\xi _{\mu }\approx 3.3039\times 10^{-4}$. The vertical dotted line marks the position of the outer matter-vacuum interface.

Figure 8: Graph of the function $\rho (\xi )$ in a low-density case, with the parameters $C=1/3$, $n=3/2$ and $\pi _{e}=10^{-3}$, resulting in $\xi _{M}\approx 0.1848$ and $\xi _{\mu }\approx 3.3039\times 10^{-4}$. The vertical dotted lines mark the positions of the matter-vacuum interfaces.

Having chosen an arbitrary initial value for $\xi$, and depending on whether we are starting at the zero of $\gamma (\xi )$ or at its inflection point, for example with either $\xi_{z}=1$ or $\xi_{e}=1$, we then choose a value for either $\pi(\xi_{z})$ or $\pi(\xi_{e})$, and we use either the value $\gamma(\xi_{z})=0$ or the value of $\gamma(\xi_{e})$ given in Equation (53), as implied by $\pi_{e}=\pi(\xi_{e})$, which is given by

  $$\gamma(\xi_{e}) = \xi_{e} F(\xi_{e},\pi_{e})\, \frac {(4+4/n)-\left[1+F(\xi_{e},\pi_{e})\right]\pi_{e}} {1+(5+4/n)F(\xi_{e},\pi_{e})}.$$ (74)

Departing from the chosen initial point, we then integrate out in both directions, of increasing and decreasing $\xi$, until $\pi (\xi )$ approaches zero in each case, with a high degree of numerical precision. Subsequently we concatenate the results into a single function from $\xi_{1}$ to $\xi_{2}$. The values obtained for $\gamma (\xi )$ at the two boundary points are then used to generate analytically the correct exact solutions in the inner and outer vacuum regions, which are plotted alongside the corresponding numerical solution, from some point to the left of $\xi_{1}$ to some point to the right of $\xi_{2}$. In order to do this, at the outer boundary we just use the fact that

  $$\gamma(\xi_{2}) = \xi_{M},$$ (75)

as given in Equation (51), in order to generate the correct exterior Schwarzschild solution as given in Equations (41) and (42). In order to obtain the dimensionfull physical parameters, we must recall that

  $$\xi_{M} = \frac{r_{M}}{r_{0}}.$$ (76)

We can then simply choose any positive value that we want for $r_{M}$, and thus obtain the corresponding value of the parameter $r_{0}=r_{M}/\xi_{M}$. This value of $r_{0}$ can then be used to obtain the values of all the other dimensionfull physical parameters of the system. In the case of the inner boundary we use the fact that

  $$\gamma(\xi_{1}) = -\xi_{\mu},$$ (77)

as given in Equation (50), as well as the fact that we have for the integration constant $A$

  $$A = \frac{1}{2}\, \ln\! \left( \frac{\xi_{1}}{\xi_{2}}\, \frac{\xi_{2}-\xi_{M}}{\xi_{1}+\xi_{\mu}} \right),$$ (78)

in order to generate the correct interior vacuum solution as given in Equations (39) and (40), with

  $$\xi_{\mu} = \frac{r_{\mu}}{r_{0}}.$$ (79)

Note that the integration constant for the function $\nu (\xi )$ within the matter region must be chosen so that at the outer interface we have $\nu(\xi_{2})=\nu_{s}(\xi_{2})$ has the correct value, given by Equation (42). This is easily accomplished by just correcting the values of $\nu (\xi )$ afterwards, by simply adding a constant to them. From Equation (32) we can see that, since at the outer interface we have that $F(\xi_{2},\pi_{2})=0$, where $\pi_{2}=\pi(\xi_{2})$, it follows that if we just ignore the integration constant we get $\nu(\xi_{2})=0$. Therefore, all that we have to do is to add $\nu_{s}(\xi_{2})$ to $\nu (\xi )$ afterwards, for which we may take advantage of the fact that $\nu_{s}(\xi_{2})=-\lambda_{s}(\xi_{2})$ is already known.

We have used three sets of values of the parameters in order to generate the data seen on the graphs, the first two starting from the inflection point of the function $\gamma (\xi )$, and the last one starting from the zero of that function. These sets of parameters are as follows.

$$\begin{array}{cllll} \mbox{Set 1:} & C=0.05, & n=1.5, ... ...5, & \xi_{z}=1.0, & \pi_{z}=1.0\times 10^{-3}. \end{array}$$

The first set uses some arbitrary mid-range values of the parameters, that have the property of just displaying in a simple and clear way all the main characteristics of the solutions. The second set of parameters corresponds qualitatively to the configuration of a high-density object such as a white dwarf or neutron star. In these two cases the integration was started from the inflection point $\xi_{e}$ of $\gamma (\xi )$. The third set of parameters corresponds qualitatively to the configuration of a low-density object, such as a normal main sequence star. In this case the integration was started from the zero $\xi_{z}$ of $\gamma (\xi )$.