Hydrostatic equilibrium of stars without electroneutrality constraint
M. I. Krivoruchenko D. K. Nadyozhin A. V. Yudin







^1Institute for Theoretical and Experimental Physics, B. Cheremushkinskaya 25 
 117218 Moscow, Russia

^2Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear Research 
 141980 Dubna, Russia 

^3Kurchatov Institute, Akademika Kurchatova Pl. 1
 123182 Moscow, Russia 

^aE-mail: yudin@itep.ru


 

The general solution of hydrostatic equilibrium equations
for a two-component fluid of ions and electrons without a local electroneutrality constraint is found in the framework of Newtonian gravity theory.
In agreement with the Poincaré theorem on analyticity and
in the context of Dyson's argument, the general solution is demonstrated to possess a fixed (essential) singularity in the gravitational constant G at G = 0.
The regular component of the general solution can be determined by perturbation theory in G starting from a locally neutral solution.
The non-perturbative component obtained using the method of Wentzel, Kramers and Brillouin
is exponentially small in the inner layers of the star and grows rapidly in the outward direction.
Near the surface of the star, both components are comparable in magnitude, and their non-linear interplay
determines the properties of an electro- or ionosphere.
The stellar charge
varies within the limits of - 0.1 to 150 C per solar mass.
The properties of electro- and ionospheres are exponentially sensitive to variations of the fluid densities in the central regions of the star.
The general solutions of two exactly solvable stellar models without a local electroneutrality constraint are also presented.



97.10.-q,
97.10.Cv


Hydrostatic equilibrium of stars 
 without electroneutrality constraint
    M. I. Krivoruchenko,^1,2 D. K. Nadyozhin,^1,3 A. V. Yudin,^1,a
    December 30, 2023
=======================================================================













§ INTRODUCTION




Charged particles constitute an appreciable or dominant fraction of the matter in stars.
The electromagnetic interaction between particles is stronger than the gravitational
interaction by a factor of

    λG = e^2/G m_u^2≈ 1.25 × 10^36,

where e denotes the proton charge in CGS units, G is the gravitational constant,
and m_u = 931.5 MeV/c^2 is the atomic mass unit, with c being the speed of light.
Under the assumption of the complete screening of electrostatic interactions, gravitational forces should play the dominant role
in determining stellar structure.
The hydrostatic equilibrium (HE) equations
are usually supplemented with a local electroneutrality constraint (LEC);
however, it is not immediately clear how effective the screening is at compensating for the
enormous difference in the strengths of the electromagnetic and gravitational forces.

Lebowitz and Lieb <cit.> showed that the LEC is required 
for the existence of extensive thermodynamic potentials of isolated systems.
The stability of electrically neutral, isolated systems
of charged particles in the thermodynamic limit was proven by Dyson and Lenard <cit.>
and by Lieb and Thirring <cit.>.

The thermodynamic equilibrium
in an external potential field is determined by Gibbs' condition <cit.>:

    μ + V(x) = const,

where μ is the chemical potential and V (x) is the external potential.
A generalization of Eq. (<ref>) for a multicomponent fluid
within the framework of General Relativity Theory (GRT) was given by O. Klein <cit.>
and, in a broader context, by Kodama and Yamada <cit.> and Olson and Bailyn <cit.>.
In the presence of gravity or other external fields, the chemical potentials of a multicomponent fluid acquire coordinate dependence,
as already seen from Gibbs' condition, which is generally incompatible with the LEC.

As early as 1924, Rosseland <cit.> showed that in Newtonian gravity theory (NGT),
the unconstrained HE equations for an ionized self-gravitating gas
with a constant temperature admit solutions with similar distributions of ions and electrons.
The charge density of the gas was found to be positive and λG times smaller than the particle number density.
In the context of white dwarfs, the polarization of stars was discussed by Schatzman <cit.>.

The smallness of the uncompensated charge density justifies the use of the LEC
when describing effects not related to the electrostatic properties of stars.

Based on the variational energy principle of GRT, Olson and Bailyn <cit.>
derived a complete system of unconstrained HE equations for a charged multicomponent fluid.
In the stellar structure models discussed by Olson and Bailyn <cit.>, Bally and Harrison <cit.>, and Neslušan <cit.>,
the polarization effect leads to a stellar charge of
Qs∼ e Ms/(m_uλG) ∼ 100 Ms/M_⊙ C.

Rotondo et al. <cit.> and Belvedere et al. <cit.>
showed recently that the components of a fluid have different spatial distributions, leading to the appearance of an electrosphere,
the charge of which is sufficiently large to completely compensate for the positive charge of
the star.

The deviations from local electroneutrality and, accordingly, the appearance of
a large-scale electrostatic field in a star
affects the diffusion and separation of chemical elements in the Sun <cit.> and in neutron stars <cit.>.
The uncompensated stellar charge could influence possible mechanisms of supernova explosions <cit.>.

Electrically charged stars <cit.> – <cit.> are described by particular solutions of the unconstrained HE equations.
Some features of the general solution can be foreseen
on the basis of simple qualitative considerations, as discussed below.

Let us consider a neutron star composed of neutrons, electrons and protons.
Suppose that it is electrically neutral and in chemical equilibrium with respect to the weak interaction.
The gravitational field of the star creates a potential well,
in which a number of electrons ΔNe can be confined,
given that that the energy of their electrostatic repulsion ∼ e^2 ΔNe^2 /Rs
does not exceed the energy of their gravitational attraction to the star
∼ G MsΔNeme / Rs, where Rs and Ms are the radius and mass of the star, respectively, and
me is the electron mass.
This condition implies that
ΔNe≲ G Msme / e^2 ∼ 5 × 10^17Ms/M_⊙.
The charge of a neutron star may vary from ∼ -0.1 C to zero.
An excess of electrons, however, would violate the assumption of chemical equilibrium.
The relaxation of the neutron star to a new state through the inverse β decay
is accompanied by the emission of
energy in the form of neutrinos.
The star then turns into a lower-energy state with the same electric charge.
If gravitational stability is not lost during the relaxation process, then the above estimate remains valid.
Considering that the number of electrons ΔNe is many orders of magnitude smaller than
the number of neutrons M_⊙/m_u≈ 1.2× 10^57,
the neutron star remains gravitationally stable.
Similar arguments lead to the conclusion that the neutron star is capable
of confining ΔNp≲ G Msmp / e^2 ∼ 10^21Ms/M_⊙ additional protons
with a total charge not exceeding Qs∼ 150 Ms/M_⊙ C.
The electric charge of a neutron star can vary from -0.1 to 150 C per solar mass.
Similar considerations apply to white dwarfs.

The general solution of the unconstrained HE equations is expected to describe
gravitationally stable configurations
with a stellar charge of -0.1 to 150 C per solar mass.

A single-component self-gravitating fluid of hypothetical
particles of mass m_f and charge ef, with f ≪ 1,
is discussed by S. Ray et al. <cit.>.
A comparison of the gravitational force ∼ G m_fMs
and the electrostatic repulsion force ∼ (ef)^2 Ms/m_f acting on an elementary volume of the fluid
yields f ∼ m_f/(m_u√(λG)) and a maximum stellar charge of
Qs∼ efMs/m_f∼  10^20Ms/M_⊙ C.
In nature, stars are composed of multicomponent fluids and are subject to a √(λG) times more
stringent constraint on Qs because of charge separation in the
stellar electric fields.

Another fundamental property of the general solution
follows from the Poincaré theorem on analyticity (see, e.g., Ref. <cit.>) and Dyson's argument <cit.>.

The Poincaré theorem on analyticity is an intuitively clear theorem that in rough terms states that
analyticity regions in the parameter space of ordinary differential equations (ODEs) are normally inherited by solutions of these equations.
Analytic functions are determined by their singularities, so one can talk also about singularities instead of analyticity regions.
It is not difficult to see that the unconstrained HE equations are singular at
G = 0 or, equivalently, at λG = ∞.
According to the Poincaré theorem, the general solution
is thus expected to be singular at λG = ∞ as well.
One of the unconstrained HE equations becomes an LEC in the limit of λG = ∞.
It would be natural to consider the LEC solution as a zeroth-order approximation
of the general solution in the framework of perturbation theory in the small parameter 1/λG.
The Poincaré theorem suggests, however, that as long as λG = ∞ is a singular point, such an expansion fails
because the radius of convergence of the perturbation series is zero.
Dyson's argument is interpreted similarly, with the same conclusions for convergence of the perturbation series.

The general solution to the unconstrained HE equations is therefore expected to be nonperturbative,
which implies that the density distributions of the fluid components are singular at λG = ∞.
Such solutions are typical of singularly perturbed ODE systems (see, e.g., <cit.>)
whose order changes in the limit of a vanishingly small parameter.
The unconstrained HE equations belong to this class of singularly perturbed ODE systems,
with the small parameter being 1/ λG.

In this paper, the properties of the general solution to the unconstrained HE equations are investigated within the NGT framework.
The outline of the paper is as follows:  Section II presents the unconstrained HE equations
for a two-component fluid of ions and electrons described by polytropic equations of state (EoS).
A unique regular solution at λG = ∞,
which is a particular solution to these equations, is analyzed in detail.
A regular solution can be constructed using perturbation theory in the small parameter 1 / λG, starting from the LEC solution.
The solutions to the unconstrained HE equations that have previously been discussed in Refs. <cit.> – <cit.>
are identified as such regular solutions.
The regular solutions are investigated
(i) for equal polytropic indices of ions and electrons,
(ii) by expanding the solution in a power series in the stellar radius near the center, and
(iii) by considering an expansion in 1/λG around the LEC solution.
In Sect. III, the electron shell of a star is discussed.
In Sect. IV, we consider an exactly solvable model with unit polytropic indices of ions and electrons,
which perfectly illustrates the non-perturbative nature of the general solution
and clarifies the important role of the irregular part of the general solution in the formation of the electro- or ionospheres of stars.
The mass, radius and charge of a star are found as functions of the central baryon number and charge densities.
The exponential sensitivity of the electro- and ionospheric characteristics to the central charge density is established.
A useful example of an exactly solvable model of a
charged star consisting of a two-component incompressible fluid
is further discussed in Appendix A.
In Sect. V, we construct the general solution to the unconstrained HE equations for arbitrary polytropic indices.
In the central stellar region, the irregular, non-perturbative part of the general solution can be found using
perturbation theory of singular ODE systems (see, e.g., <cit.>).
Such problems arise in the semiclassical limit of the Schrödinger equation,
where the method of Wentzel, Kramers and Brillouin, commonly known as the WKB approximation,
is proven to be effective.
Perturbation theory applies everywhere, except in a thin subsurface layer of the star with a thickness of ∼ 1/√(λG).
To match the internal solution to the outer electron or ion shell, numerical methods can be used.
The conclusions section summarizes the results.




§ REGULAR SOLUTION TO THE UNCONSTRAINED HYDROSTATIC EQUILIBRIUM EQUATIONS





We consider the problem of the polarization of a substance in a spherically symmetric star.
The thermodynamic equilibrium equations for the chemical potentials of the ions, μi, and the electrons, μe,
coincide with Gibbs' condition (<ref>) and can be written as

    μi + miφG + Ze φE   =   const, 
    μe + meφG - e φE   =   const,

where φG and φE are the gravitational and electrostatic potentials,
eZ and mi = A m_u are the charge and mass of the ions.
Bearing in mind possible applications to physics of white dwarfs,
we set temperature equal to zero, and suppose the dominance of the contribution of electrons in thermodynamic functions.
The chemical equilibrium with respect to the weak interactions is neglected.
We also assume the ion and electron EoS are polytropic with P ∼  n^γ
and attempt to reduce the problem to one known in the theory of polytropes (see, e.g., Ref. <cit.>).




 §.§ Two-fluid model



The gravitational and electrostatic potentials 
satisfy the Poisson equations

    φG   =     4π Gρm,
    φE   =    -4πρe,

where the mass and charge densities are, respectively,

    ρm   =   mini + mene,
    ρe   =    Zeni-e ne.

The radial derivative term of the Laplacian  is equal to

    _r = 1/r^2d/dr(r^2d/dr).


By applying  to Eqs. (<ref>) and (<ref>) and constructing appropriate linear combinations of the chemical potentials,
we obtain

    _r(Z μe+ μi)   =   -4π G (mi+Zme)ρm,                   
    _r( miμe-meμi)   =   -4π e (mi+Zme)ρe.


According to Eqs. (<ref>) and (<ref>), the gradient -∇ (Zμe+μi) 
is equal to the gravitational force acting on an ion and Z electrons.
The thermodynamic relationship dP = n d μ allows the gradients of the chemical potentials to be presented in terms of pressure gradients.
As a special case, the hydrostatic limit of the Euler equation is recovered for an electroneutral fluid with ne = Z ni:

    1/ρm∇(Pi+Pe) + ∇φG = 0.

As mentioned above, we assume that Pi≪Pe holds everywhere.

First, we define the standard relations between pressures and concentrations from the theory of polytropes:

    Pk = Kknk^1 + 1/k,   k = (i,e),

where Kk and k are fixed constants and k is called the polytropic index.
The ion and electron components interact with each other only through the gravitational and electrostatic forces.
No fundamental difficulties arise when considering the influence of the interactions of the components (see, e.g., <cit.>), although this consideration requires special attention.
The dimensionless functions θk are defined as follows

    nk≡nk0θk^k,

where the nk0 are the concentrations at r = 0.
Excluding the rest mass, the chemical potentials are equal to μk = μk0k, where μk0 = Kk(1 + k)nk0^1/k.
The functions k =  k (x) depend on a dimensionless coordinate x, which is related to r as follows:

    r=r_0 x,   r_0^2 = Z μe0/4π G mi(mi + Zme)ni0.

These expressions allow Eqs. (<ref>) and (<ref>) to be rewritten in the following dimensionless form:

    _x(e+ii)    =    - (i^i+mee^e),
    _x(e-mii)    =    -G (i^i-ee^e),

with the initial conditions

    k(0)=1,     θ^'_k(0)=0.

The Laplacian _x that appears in Eqs. (<ref>) and (<ref>) acts on the coordinate x.
The dimensionless Λ parameters are defined below.

The ratio between the densities ne0 and ni0 determines the parameter

    e=ne0/Zni0≈ 1,

which characterizes the deviation from electroneutrality in the center of the star.
In the limit of G = ∞, Eqs. (<ref>) and (<ref>)
require e = 1, which is the LEC at r=0.

i characterizes the smallness of the ion pressure:

    i=μi0/Z μe0≈(1 + i)Pi0/(1 + e) Pe0≪ 1,

where Pk 0 is the pressure of the components at the center of the star.
The last equality in (<ref>)
is valid to the same accuracy as the approximation e≈ 1.
For e≠i,
the parameters e and i are independent.

The parameter m is a measure of the smallness of the electron mass relative to the ion mass:

    m=Zme/mi=Zme/A mu≈ 5.5 × 10^-4Z/A.


The main parameter G is expressed as follows:

    G=Z^2e^2/Gmi^2 = (Z/A)^2 λG .

A high value of ΛG
plays the predominant role in specifying the nature of the problem.

Notably, e and i depend on the EoS and the stellar structure through nk0 and Kk,
whereas m and G depend only on the chemical composition of the fluid (A and Z).

Equations (<ref>) – (<ref>) constitute the system of unconstrained HE equations in dimensionless form.
This system defines a Cauchy problem with the initial conditions given in Eq. (<ref>).
In Eq. (<ref>), the small parameter 1/G appears as a multiplier of the Laplacian; consequently,
Eqs. (<ref>) and (<ref>) belong to the class of singularly perturbed ODEs.






 §.§ Two-fluid model with equal polytropic indices



The key aspect of the problem is best illustrated by first considering
the case of identical polytropic indices i = e≡η.
A particular solution to the system defined in Eqs. (<ref>) and (<ref>)
can be represented as
i (x) = e(x) ≡θ (x̃), with

    x̃=x√(1+me/1+i),

where the function θ(x̃) is a solution to the Lane-Emden equation

    _x̃θ(x̃) = -θ^η (x̃)

with the standard initial conditions θ(0)=1 and θ'(0)=0.
In the case of i = e, the parameters e and i are dependent:

    i=Ki/Kee^1/η Z^1 + 1/η.

The self-consistency of Eqs. (<ref>) and (<ref>) requires e = e^reg, where

    e^reg=G(1+i)-(1-mi)/G(1+i)+m(1-mi).

Equations (<ref>) and (<ref>) define i and e
as functions of G.

In the model considered here, the ion and electron densities appear to be similar, as in the Rosseland model <cit.>.
The parameter e differs from unity, so the LEC is violated at O(1/G).
The solution describes a uniformly polarized star with a positive electric charge.

The concentrations ni0 and ne0 are considered as input parameters of the model.
One would expect the general solution to represent a two-parameter set of functions.
The constraint e = e^reg restricts the freedom to a one-parameter set.
This seeming paradox is explained in Sects. IV and V. Here, we merely remark that an additional
degree of freedom due to variations of e
arises from the singular component of the general
solution to Eqs. (<ref>) and (<ref>).




 §.§ A power-series expansion near the center



A series expansion in powers of x is effective for solving Eqs. (<ref>) and (<ref>) for arbitrary polytropic indices,
and it sheds further light on the origin of the constraint e = e^reg.
Considering the initial conditions, we search for a solution in the following form:

    k=1+∑_p=1^∞β_k px^2p,   k=(i,e).

Terms with odd powers of x are absent because of the even parity of Δ_x.
By substituting these expressions into Eqs. (<ref>) and (<ref>) and equating the coefficients at the same powers of x,
we obtain an infinite system of algebraic equations.

The zeroth-order terms yield

    6(βe1+iβi1)   =   -(1+me),
    
    6(βe1-miβi1)   =   -G(1-e).

The left-hand side of Eq. (<ref>) is of O(1),
so we require e = 1 - αe / G with αe = O(1).
βe1 and βi1 can be found from Eqs. (<ref>) and (<ref>), whereas αe remains a free parameter.

The ∼ x^2 terms yield the equations

    20( βe2+iβi2)   =   -(βi1i+meβe1e),
    
    20( βe2-miβi2)   =   -G(βi1i-eβe1e).

Again the boundedness of the expansion coefficients in last of the above equations yields

    βi1i1=βe1e1 + γ_1/G,

where γ_1  = O(1) is an unknown parameter.
This equality, with the use of Eqs. (<ref>) and (<ref>), allows one to determine the lowest-order parameters,
including α_e, to an accuracy of O(1/ Λ_G):

    βe1/i   =   βi1/e = -1 + m/6(i + ie), 
    
    αe   =   (i-mie)1 + m/i + ie.

In the limit of i = e, αe is consistent with Eq. (<ref>).
The higher-order terms in 1/ G are determined, in turn, from the condition of the boundedness of the expansion coefficients at x^4, etc.
The regularity condition at G = ∞ thus imposes a constraint on the concentrations of ions and electrons in the center of the star.
We are left with a set of one-parameter solutions.

The results of this subsection can be generalized to GRT <cit.>.




 §.§ A power-series expansion in G



The occurrence of a large parameter G in Eqs. (<ref>) and (<ref>) makes
it possible to search for solutions throughout the entire x range using a power-series expansion:

    k   =   k0+k1G^-1  +O(G^-2), 
    e   =   e0+e1G^-1+O(G^-2).

To a zeroth-order approximation, the regularity of the right-hand side of Eq. (<ref>) yields
i0^i (x) = e0e0^e (x).
The initial conditions e0 (0) = i0 (0) = 1 imply that

    i0^i(x)≡e0^e(x),  e0=1.

Equation (<ref>) yields

    _x(e0+ii0)   =   
    -(1+m)e0^e
       =   -(1+m)i0^i,

which is the generalized polytropic equation. We solve it for the initial conditions (<ref>)
in the interval from x = 0 up to the stellar surface at x = xb, where
e0 (xb) = i0 (xb) = 0.

The first-order terms in Eqs. (<ref>) and (<ref>) yield


    _x(e1+ii1)= 
                  -e0^e(ii1/i0 + mee1/e0+me1),
    
          _x(e0-mii0)= 
                   -e0^e(ii1/i0 - ee1/e0-e1),


with the initial conditions k1 (0) = θ^'_k1 (0) = 0.
Let us consider Eq. (<ref>) for x → 0. Considering the asymptotic behavior
of k (Eq. (<ref>) for βk1 and αe),
we obtain the identity to zeroth order. Considering e1 = - αe,
Eq. (<ref>) can be written as follows:

    (1+m)i/i+ie_x(ee0-ii0)=-e0^e(ii1/i0-
    ee1/e0).

For i = e, this equation is satisfied by the similar distribution i1 = e1.
Equation (<ref>), when combined with Eq. (<ref>), can also be
written as follows:

    _x(e1+ii1+m[e0-mii0])=                  
    
    -(1+m)e0^eii1/i0.


The zeroth- and first-order approximations are defined for all values of x.
The region near the surface of the star should be treated separately,
because for k0→ 0, the validity of the expansion (<ref>) can be violated.




 §.§ Global stellar parameters



The stellar radius (more precisely, the boundary of its baryon component) is defined as Rs = r_0 xb,
where r_0 is given by Eq. (<ref>) and xb is fixed by the condition θ (xb) = 0
from the solution to Eqs. (<ref>) and (<ref>). The mass of the star takes the form

    Ms   =   ∫_0^Rs 4π r^2(nimi + neme) dr 
       =   -4π r_0^3 ni0mi x^2 d/dx(e+ii)|_x=xb,

where Eq. (<ref>) is used to proceed to the second line.
The total stellar charge (at the boundary of the baryon component) is equal to

    Qs   =   ∫_0^Rs 4π r^2(Zeni-ene)dr 
       =   -4π r_0^3 ni0Ze/G x^2 d/dx(e-mii)|_x=xb,

where we have used Eq. (<ref>). The stellar charge can be estimated as

    Qs∼ZeMs/A m_uG≃ 10^21 e A/ZMs/M_⊙,

where M_⊙ is the mass of the Sun.

The total uncompensated electric charge of a star of one solar mass, on the boundary of its baryon component,
is on the order of 100 C, in agreement with the estimates of Pikel'ner <cit.>, Bally and Harrison <cit.> and Neslušan <cit.>.
The maximum charges of neutron stars and white dwarfs are estimated to be 50 C and 500 C, respectively <cit.>.
For comparison, one mole of ^12C contains 6 N_A protons, whose total charge is 6 × 10^5 C.
The total negative charge of the Earth at its surface is estimated to be (4 - 5.7) × 10^5 C <cit.>.
In absolute value, the solar charge is less than (0.4 - 1) × 10^18 C <cit.>.





§ ELECTROSPHERE





In Sect. II, we have seen that
the self-consistent regular solution to the unconstrained HE equations 
describes a star with a positive charge. To move an electron from the surface of such a star
to infinity, it is necessary to perform work against the forces of gravity and Coulomb attraction.
In fact, a star can confine more electrons than it should be able to according to the regular solution.
This possibility leads to the formation of an electrosphere that is quite independent of
the internal structure of the star. We therefore discuss the electrosphere separately.

Under normal conditions, a metal is surrounded by a double charged layer, where an excess of positive charge at the subsurface
is created by the atomic lattice and a negative charge outside is created by electrons.
The thickness of the layer is equal to several Bohr radii, and the work function is several eV. An electrosphere surrounds
the surface of hypothetical strange stars <cit.>. Within the strange star, the positive charge of the baryon component
is almost completely neutralized by electrons, whereas it is only partly neutralized in the subsurface.
As a result, an electron shell forms above the surface, with a thickness exceeding hundreds of fm, and
the chemical potential of electrons near the surface is several tens of MeV.
The electrons also give rise to an electrostatic potential jump across the surface separating different
phases of nuclear matter <cit.>.

Earlier, we presented arguments for the existence of stars with a total charge of -0.1 to 150 C per solar mass,
which are able to support either an electrosphere or an ionosphere.
Here, we restrict ourselves to a discussion of the electrosphere. 
Section IV describes the properties of the ionosphere.



 §.§ A polytropic model


The regular solution explicitly yields a positive charge excess at r=0. 
Furthermore, the balance of the gravitational and electrostatic forces ensures that the charge excess must be positive everywhere
from the central region up to the surface.
The electron component extends somewhat beyond the ion component and compensates for the stellar charge.

Applying the Laplacian to Eq. (<ref>), we obtain

    _rμe=4π e^2ne(1-G me^2 /e^2).

At the surface (r = Rs), an equilibrium condition holds:

    -1/ned Pe/dr=eQs/Rs^2+
    GmeMs/Rs^2,

where Ms and Qs are given by Eqs. (<ref>) and (<ref>).

Proceeding in the standard way, we define a dimensionless polytropic function based on the equation
ne = na0a^e, where na0 is the electron density at the boundary
of the ion component.
The thickness of the envelope is assumed to be much smaller than Rs.
This assumption will be checked at the end of the calculation. The Laplacian simplifies to _r→d^2 / dr^2.
For convenience, we introduce the dimensionless parameter

    λm = me/m_u.

The dimensionless variable y is defined by

    r=Rs + ray,  ra^2 = μa0/4π e^2 (1 -λm^2/λG)na0,

where μa0 is the chemical potential of electrons at the surface. Thus, Eq. (<ref>) becomes

    d^2a/dy^2=a^e(y).

At the surface, y = 0 and a(0) = 1.
Once the stellar structure equations are solved and the stellar mass, radius, and charge are determined,
the initial condition for a^' (0) follows from Eq. (<ref>).
Considering Eq. (<ref>), we write the charge in the form

    Qs=Ms/m_ueqs/λG,

where qs = O(1) can be found by solving the unconstrained HE equations. We remark that M_⊙/m_u∼ (λG/α)^3/2,
where α = e^2/(ħ c) ≈ 1/137 is the fine structure constant,
and therefore, Qs∼ e √(λG)/α^3/2.
The second initial condition can be written as

    da/dy|_y=0=-(qs+λm)
    √(Pun/(1 + e)Pa0(λG - λm^2)),

where Pa0=Kena0^1 + 1/e is the pressure of the electrons at the boundary
and Pun is the natural unit of pressure in the star:

    Pun = G Ms^2/4 πRs^4.

Integrating Eq. (<ref>) yields

    da/dy=-√(2(W_0 - V(a))),

where 
    W_0   =   a^' 2(0)/2 - 1/1+e, 
    
    V(a)    =    - a^1 + e/1+e,

and a'(0) is given by Eq. (<ref>).
A second integration yields

    y = ya   -   a/√(2(W_0 - V(a)))
       ×   _2F_1 (1/2,1;2+e/1+e; - V(a)/W_0 - V(a)),

where _2 F_1 is the Gauss hypergeometric function. Equation (<ref>)
implicitly specifies a = a (y).
With the initial condition a (0) = 1,
the boundary of the electrosphere can be found from a (ya) = 0 to give

    ya= 1/|a^'(0)| _2 F_1(1/2,1;2+e/1+e;
    2/(1+e)a^' 2(0)).


Figure 1 shows the behavior of ya as a function of W_0. With a decrease in W_0, the radial thickness of the envelope increases,
and for W_0 → 0, the radial thickness goes to infinity. However, the function ya is still well defined at zero:

    ya(W_0 = 0) = - √(2(e+1))/e-1 < 0,

because of the irregularity of 
ya for W_0 = 0. In comparison with the nonrelativistic case,
the shell of relativistic electrons extends higher.









Equation (<ref>) can be interpreted as the equation of motion of a point particle of unit mass in a potential V(a),
where y plays the role of a time variable. Equation (<ref>) follows from the conservation of the total energy:

    W = a^' 2/2 + V(a) .

There are three types of solutions. If the initial momentum is small, the particle moves to the left relative to its initial position a (0) = 1,
stops at some point 0 < a < 1, and then accelerates in the positive direction and goes to infinity.
The total energy of the particle is negative: W = W_0 <0. In the Thomas-Fermi theory of multielectron atoms, solutions of this type are called compressed;
they terminate at a^' = 0 for a finite electron density,
for which the global electroneutrality of the atom holds <cit.>.
In our problem, however, such solutions are not physical since the stars are, of course, isolated.
The second type of solution corresponds to the case of zero energy: the particle arrives at the point a = 0 and stops there, such that a^' = 0.
In the Thomas-Fermi theory, such solutions describe neutral atoms. In our case, they describe stars with a maximum negative
total charge, i.e., Q_tot^max < 0.
The third type of solution is obtained for W_0> 0. The particle arrives at the point a = 0 with a finite velocity of a^' <0,
and then the solution terminates. In the Thomas-Fermi theory, solutions of this type correspond to ionized atoms with a total positive charge.
In our case, such a solution describes an electrosphere of a finite thickness y_a, and the total charge of the star can be negative, zero, or positive.

The existence of solutions, therefore, depends on the derivative of the chemical potential of the electrons at the boundary of the ion component.
The absolute value of the derivative is constrained from below by the condition W_0 ≥ 0, which imposes an upper limit on the number
of electrons in the shell.

At the boundary of the ion component, the total positive charge of the star is equal to Qs, and the charge of the electrosphere Qe either partially, completely,
or excessively compensates for Qs. Considering Eqs. (<ref>) and (<ref>),
the charge of the electrosphere can be written in the form

    Qe   =    - 4πRs^2∫_0^ya e nera dy 
       =    4πRs^2 e na0ra( a^'(0) + √(a^' 2(0) - 2/1 + e)).

Let κ
be a parameter characterizing the screening, such that Qe = - κQs, where 0 ≤κ≤κ_max.
In dimensionless units (see Eq. (<ref>)), the charge of the electrosphere is equal to
qe = - κqs. In particular, for κ = 1, the screening is complete, and the total stellar charge vanishes.
From Eq. (<ref>), one can find

    Pa0=κqsPun/λG[qs + λm-
    κqs/2(1-λm^2/λG)],

implying that the pressure in the envelope is λG times smaller than Pun.
It is now convenient to rewrite Eq. (<ref>) in the form


    da/dy|_y=0=-
    qs+λm/√(κqs (1 + e)(1 -λm^2/λG)(q_s + λm -
    κqs/2(1 -λm^2/λG)))
    .



The maximum negative charge of the envelope corresponds to W_0 = 0. From Eq. (<ref>),
one can find the critical value of the compensation parameter:

    κ_max = 1 + λm/qs/1 - λm^2/λG > 1.

Solutions with a total negative charge, Q_tot =  Qs+ Qe= Qs (1 - κ) <0, are admissible
since an electrically neutral star has a gravitational field capable of binding electrons under the condition
that the energy of their Coulomb repulsion does not exceed the energy of their gravitational attraction to the star.
An estimate of Qs∼ 100 C was obtained above. The maximum negative charge Q_tot^max
is proportional to the parameter λm≪ 1; thus, Q_tot^max∼ - 0.05 C per solar mass.

The gravitational collapse of a charged single-component spherically symmetric fluid is discussed
in Refs. <cit.>.

The maximum charge of a neutron star is determined in Ref. <cit.> on the basis of the requirement of gravitational stability.
We showed that the maximum charge can instead be determined from the requirement for the existence of a solution of finite extent
for the electro- or ionosphere of the star.
The latter requirement yields a limit that is √(λG) times stronger than that obtained from the requirement of gravitational stability.

The positive charge excess of strange stars at the quark surface is found in Refs. <cit.> to be about 10^20 C.
The electrosphere compensates this charge, although the detailed mechanism has not yet been
studied. The particular values of the bag constant B used in Refs. <cit.>
rule out the existence of massive neutron stars metastable against conversion to strange stars <cit.>.




 §.§ Electrically neutral stars



For W_0 = 0, Eqs. (<ref>) - (<ref>) yield

    a^max(y)=(1 + y/y_0)^-2/e-1,

where y_0 = - ya(W_0 = 0) > 0.
This solution extends to infinity, which means that an electron can be moved from the center of the star to infinity without doing work.
In the case of a relativistic electron gas, e = 3, the solution decays rather slowly: a^max∼ 1/y.
At low densities, the electrons become nonrelativistic, with e = 3/2, and the solution
approaches a^max∼ 1/y^4.
It is not difficult to directly verify that a^max (y) is a solution to Eq. (<ref>) with the proper boundary conditions.










In view of the approximate equality κ_max≈ 1, the electrosphere of a neutral star with κ = 1 is described by a solution close to
a^max  (y).
We return to Eq. (<ref>) and set κ = 1. From this equation,
by neglecting the small terms ∼ 1/λG and ∼λm,
the derivative of a  (y) at y=0 can be found to be

    dθ/dy|_y=0≈-√(2/1+e)(1 + ξ),

where

    ξ=1/2(λm/qs)^2 ≪ 1.

We are looking for a solution to Eq. (<ref>) in the form a = a^max + σ, where σ is a small correction.
The initial conditions are of the form σ (0) = 0 and σ^'(0) = - ξ√(2/(1+e)).
To the lowest order, σ satisfies the equation

    d^2σ/dy^2=eσ/(1+y/y_0)^2.

The solution to this equation is given by

    σ(y)=2ξ/3e+1[(1+y/y_0)^-e+1/e-1-(1+y/y_0)^2e/e-1].











The dimensionless chemical potential of the electrosphere of a neutral star
as a function of the distance to the surface is shown
in Fig. <ref> for e = 3/2 and
in Fig. <ref> for e = 3.
The equilibrium conditions (<ref>) and (<ref>)
are valid for nonrelativistic particles. For completeness, we also consider the relativistic case η_e = 3,
which can be interesting for qualitative estimates. 
The total stellar charge is assumed to be zero (κ = 1),
whereas the parameter ξ = 1.5× 10^-7 corresponds to qs = 1.
The dashed curves show the asymptotic solutions corresponding to Eq. (<ref>),
the solid curves labeled E show the exact solutions corresponding to Eq. (<ref>), and the solid curves labeled A show the approximations a^max + σ.
The approximations reproduce the exact formula in Eq. (<ref>) with reasonable accuracy.
The thickness of the electrosphere can estimated to be

    ya/y_0≈(3e+1/2ξ)^e-1/2(e+1)-1.


An order-of-magnitude estimate of the thickness of the electrosphere yara can be found from Eq. (<ref>).

For nonrelativistic electrons with e = 3/2, the pressure and density are related by Pe∼ħ^2ne^5/3 / me.
Considering that Pa0∼Pun / λG, we obtain

    ra/Rs   ∼   [(aB/Rs)^3 m_u/MsλG]^1/5
       ≈    2.2× 10^-16(R_⊙/Rs)^3/5(M_⊙/Ms)^1/5,

where aB = ħ^2 / (e^2 me) is the Bohr radius.
In comparison with the stellar radius, the thickness of the electrosphere is quite small,
although it is not small on an atomic scale:
yara∼ 6 × 10^4aB for an electrically neutral star of mass Ms = M_⊙ and radius Rs = R_⊙.

In the case of e = 3, 
the substitution of the relativistic relationship
Pe∼ħ c ne^4/3 into Eq. (<ref>) yields

    ra/Rs∼√(λG/α^3/2m_u/Ms)≈ 1.3×10^-9(M_⊙/Ms)^1/2.

The thickness of the electrosphere of an electrically neutral neutron star with a radius of Rs = 10 km and a mass of 1.4   M_⊙
can be estimated to be yara∼ 1 mm.





§ EXACTLY SOLVABLE MODEL





As seen from the previous discussion,
the regular solution describes a polarized star with a nonzero positive charge as measured at the ion surface.
The properties of the electrosphere are uniquely determined by the internal structure of the star.
The electrosphere ensures that the total stellar charge is generally not equal to zero.
This feature is particularly evident in the model with equal polytropic indices considered in Sect. II.B.
An electrosphere occurs provided that θe (xb) ≠ 0, whereas for identical polytropic indices,
the distributions vanish at the same point, x = xb.

Earlier, we presented simple physical arguments in favor of the existence of stars with electric charges in the range of - 0.1 to 150 C per solar mass.
A regular solution corresponds to a single representative of this class.
The regular solutions described in Refs. <cit.> – <cit.> and in Sect. II above
do not exhaust all physical possibilities, nor are they general solutions to the unconstrained HE equations.

Fourth-order systems of ODEs, similar to Eqs. (<ref>) and (<ref>), depend on four initial conditions.
In our case, two of the initial conditions constrain the derivatives of the chemical potentials (cf. Eq. (<ref>)).
These conditions ensure smooth behavior of the solutions at r=0 and cannot be modified.
The remaining two initial conditions fix μi and μe and thereby define a two-parameter set of stellar configurations.
As shown in Sect. II, the requirement of regularity imposes a constraint on the value of Λe,
such that only a one-parameter set of stellar configurations survives.

Some of the solutions to Eqs. (<ref>) and (<ref>) appear to be lost
with the imposition of a requirement for regularity in the neighborhood of
ΛG = ∞.




 §.§ General properties of solutions



The analytic properties of solutions to ODE systems are determined by the Poincaré theorem on analyticity (see, e.g., <cit.>),
which asserts, roughly speaking, that solutions to ODE systems, when they exist,
are analytic functions of the initial coordinates and parameters in the region of analyticity of the ODEs.
The origin of non-regular solutions can also be understood on the basis of Dyson's argument,
which provides an effective qualitative criterion for the non-analyticity of observables in terms of the system parameters <cit.>.
From these perspectives, we discuss the analytic properties in the parameter ΛG
of the solutions to Eqs. (<ref>) and (<ref>).




  §.§.§ Poincaré theorem on analyticity



Using the independent variables πk = k^',
Eqs. (<ref>) and (<ref>) can be reduced to a system of four first-order ODEs.
In terms of Φ = (πe, e, πi, i),
these equations
take the form Φ^' = 𝔉 (x,  Φ, Λ),
where 𝔉 is a vector function of x, Φ,
and Λ = (e, i, m, G).
The initial conditions are Φ(x = 0) = (0,1,0,1).

The components of 𝔉, as can be easily seen, are analytic functions of Φ and Λ.
The domain of analyticity of 𝔉 is a Cartesian product of the complex planes of Φ and Λ
without Φ = ∞ and Λ = ∞.
Given that the polytropic indices are not integers,
k = 0 and k = ∞ are branch points of 𝔉.

According to the Poincaré theorem, the solutions, when they exist, depend analytically on Φ (0) and Λ
in the domain of analyticity of 𝔉 (x, Φ, Λ).
From Eq. (<ref>),
we observe that the dimensionless chemical potentials at x=0 are frozen,
so that the initial values Φ(0) are not of interest.
The essential information is contained in e
and the dimensional quantities r_0 and ni0.
The analyticity theorem can be interpreted to mean that throughout the entire domain of their existence, solutions to Eqs. (<ref>) and (<ref>)
are analytic in Λ for Λ≠∞.
At infinity, the solutions are either regular or irregular.

The dependence on G is of particular interest because G = ∞ corresponds to the LEC.
It would be natural to use the LEC solution
as an initial approximation of the general solution to
Eqs. (<ref>) and (<ref>)
and to account for deviations from local electroneutrality
by means of perturbation theory in 1 / G≪ 1.
The analyticity theorem, however, asserts that the general solution is typically not regular at G = ∞.
If this is so, the perturbation series in 1/G diverges
because the radius of convergence, which is defined by the closest singularity on the complex G plane, is equal to zero.

The convergence of the perturbation series for solutions that are regular at G = ∞
can be determined by, e.g., the singularities of e and i as functions of G,
as in the model of Sect. II.B.
According to Eqs. (<ref>) and (<ref>),
e, i, and the corresponding system of ODEs
are not analytic for |G| ∼ 1.
If there are no singularities close to infinity,
the point G = ∞ lies within the circle of convergence of
the perturbation series centered at the physical value (<ref>).
In the model of Sect. II.B, one could obtain then the LEC solution
in the form of a convergent series starting from the regular solution.
The reciprocal statement is also true: starting from the LEC solution,
one could obtain the regular solution in the form of a convergent series.




  §.§.§ Dyson's argument



A physical argument leading to similar conclusions has been presented by Dyson <cit.>.
In an extended context, the argument 
concerns the relationship
between singular points in the parameter space and
qualitative changes in the behavior of the system.
As applied to quantum electrodynamics (QED), Dyson's argument leads to the well-known
conclusions that α = 0 is a singular point of physical observables, QED series in powers
of the fine structure constant have radii of convergence equal to zero and therefore are asymptotic.

In our case, G = ∞ corresponds to G = 0.
The physical picture for G < 0 undergoes obvious qualitative changes; namely,
the attraction between particles is replaced by repulsion, and as a result, 
bound states (stars) no longer exist.

We also remark that the many-body Hamilton function for G <0 is nonnegative,
whereas for positive G, it is not bounded from below.

It can be observed that the characteristic length given in Eq. (<ref>) has branch points at G = 0 and ∞.
The stellar mass (<ref>) and charge (<ref>) inherit these singularities,
which are further superimposed with the singularities of k.

Dyson's argument indicates that
the general solution to Eqs. (<ref>) and (<ref>) should be sought in the class of functions
that depend on G in an irregular manner.




 §.§ Two-fluid model with unit polytropic indices



The above statements can be illustrated and precisely formulated
using the model of Sect. II.B with unit polytropic indices, which admits an explicit analytic solution.
The adiabatic index of the model, γ = 1 + 1/η = 2, is sufficiently close to that of an ideal gas model with γ = 5/3.
In terms of the functions φk≡ x k (x), Eqs. (<ref>) and (<ref>)
are reduced to a linear ODE system:

    φe^'' + iφi^''   =    -(φi+meφe),
    φe^''-miφi^''   =    -G(φi - eφe),

where the prime denotes differentiation in x. We look for solutions of the form φk = αk e^β x.
Equations (<ref>) and (<ref>) yield

    (αe+iαi)β^2           =   -(αi+meαe),
    
    (αe-miαi)β^2   =   -G(αi-eαe).

The algebraic system formed by Eqs. (<ref>) and (<ref>) has nontrivial solutions provided that the determinant is equal to zero.
This condition gives rise to a quadratic equation in β^2, whose solutions are


    β^2_±   =   1/2i(1+m){G(1+ie)-(1+m^2ie) 
       ±   √(G^2(1+ie)^2-
    2G[(1-mie)^2-ie(1+m)^2]+
    (1+m^2ie)^2)}.


To the leading order in 1/G, the values of β are given by

    β^2_+   =   G1+ie/i(1 + m) + O(1),  
    β^2_-   =    -e(1+m)/1+ie+ O(G^-1).

The values ±√(β^2_+) are real and large, whereas the values ±√(β^2_-) are imaginary.
We adopt the notation β_1 = |β_+| > 0 and β_2 = |β_- |> 0.
The general solution has the form

    φk(x)   =   αk1exp(-β_1 x)+αk2exp(β_1 x)
       +   αk3sin(β_2 x)+αk4cos(β_2 x).

The initial conditions for k lead to
φk (0) = 0, φk^'(0) = 1, and φk^''(0) = 0.
The general solution becomes

    φk(x)=αk1sinh(β_1 x)+αk2sin(β_2 x),

with the constraint

    αk1β_1+αk2β_2=1.

By substituting β^2 = β_1^2 and - β_2^2 into Eq. (<ref>) for αk1 and αk2, respectively,
and using Eq. (<ref>), we can uniquely determine αk1 and αk2 and find the solution.




  §.§.§ Closer inspection of analyticity in G



The analytic properties of φk (x) in the complex G-plane can be clarified by
transforming the system defined by Eqs. (<ref>) and (<ref>) into normal form with the
independent variables φk and πk = φk^'.
In terms of Φ≡ ( πi, φi, πe, φe ),
Eqs. (<ref>) and (<ref>)
take the form

    Φ^' = 𝔄Φ,

where 𝔄 is a real 4 × 4 matrix that
depends linearly on the parameters
Λ = (e, i, m, G). The solution takes the form

    Φ (x) = exp (𝔄 x) Φ (0).

The series expansion of the exponential matrix converges,
so the evolution operator 𝔇 = exp (𝔄 x) is well defined.
The matrix elements of 𝔇 are analytic functions in Λ for Λ≠∞.
A simple pole of 𝔄 at Λ = ∞ is transformed
into an essential fixed singularity of 𝔇, which does not depend on the initial values.
The functions Φ (x) are, therefore, analytic functions
of Λ in a Cartesian product of the complex planes of e, i, m, and G.
Equation (<ref>) also shows that Φ (x) are analytic on the initial values Φ (0),
in agreement with the Poincaré theorem on analyticity.

Let 𝔞_n denote the eigenvalues of 𝔄 that coincide with ±β_±∈ℂ^1,
and let 𝔯_n and 𝔩_n be their corresponding right and left eigenvectors
satisfying 𝔄𝔯_n = 𝔞_n𝔯_n and
𝔩_n^T 𝔄 = 𝔩_n^T𝔞_n. The eigenvectors are normalized according to
𝔩_n^T𝔯_m = δ_nm, where n = 1,2,3,4.
The evolution operator admits the representation

    𝔇 = ∑_n𝔯_n𝔩_n^Texp (𝔞_n x).

The eigenvalues 𝔞_n are roots of the characteristic equation,
so they have additional singularities compared with 𝔄.
The analyticity of 𝔇 implies that additional singularities cancel out with
singularities of the projection operators 𝔯_n𝔩_n^T.

A more explicit proof of this claim uses the Frobenius representation

    𝔯_n𝔩_n^T = ∏_m ≠ n𝔄 - 𝔞_m/𝔞_n - 𝔞_m.

We first note that the radicand in Eq. (<ref>) vanishes for a complex-conjugate
pair of G, where β_± have square-root branch points.
The point G = ∞ is a simple pole of β_+^2,
a square-root branch point of β_+, and an essentially singular point of φ_k (x).
G = 0 is a simple zero of β_+^2 and a square-root branch point of
β_+.
The eigenvalues 𝔞_n (±β_±) accordingly have four square-root branch points in the complex G plane.
By passing around one of these points, we obtain simply a permutation of 𝔞_n.
The representation of 𝔇 using Eq. (<ref>) is explicitly invariant under these permutations,
so the matrix elements of 𝔇 do not
inherit additional singularities of 𝔞_n.

Consequently, the eigenvalues ±β_± and the projection operators of 𝔄 have four singular points,
whereas φ_k (x) has a unique singular point G = ∞.

Given that φk(x) is analytic in G∈ℂ^1 ∖∞,
φk(x) if not a constant must be singular at G = ∞. 
The existence of solutions independent of G is excluded by Eq. (<ref>).
Regular solutions at G = ∞ can be constructed by replacing the parameters, e.g., e and i,
with functions of G, as in the model of Sect. II.B.
Solutions inherit the singularities of these functions,
which makes it possible to construct a regular solution.
The dependence of the parameters on G is a necessary 
condition for the existence of a regular solution in the neighborhood of G = ∞.

A qualitative analysis of Eqs. (<ref>) and (<ref>), developed above, reveals
the mathematical reason why the parameters e and i
are not independent in the regular solution.
The analyticity in G follows from the representation given in Eq. (<ref>).
The requirement of regularity at G = ∞
would imply regularity in the extended complex G plane.
In such a case, φk(x) must be a constant.
By introducing a dependence of the parameters on G,
particular, non-trivial, regular solutions at G = ∞ can be constructed,
which can further be used to approximate the general solution.




  §.§.§ Electro- and ionospheres



The first term in Eq. (<ref>) depends exponentially on β_1 ∼√(G).
For x ≠ 0, sinh (β_1 x) has no zeros. For the chemical potential to vanish at the boundary,
the solution must contain the component ∼sin (β_2 x). The first term in Eq. (<ref>)
is maximal on the boundary of the star, where it should not exceed α_2. The ratio α_1 / α_2
is therefore exponentially small.

The solutions are recovered with exponential accuracy already for α_2 β_2 = 1.
Substituting this expression into Eqs. (<ref>) and (<ref>), we obtain

    (1+i)β_2^2           =   1+me,
    
    (1-mi)β_2^2   =   G(1-e).

The self-consistency of the system leads to e = e^reg,
where e^reg is given by Eq. (<ref>).
The corresponding β parameters (with an additional index 0) become

    β_10^2   =   G(1+i)-(1-mi)/i(1+m),
    β_20^2   =   G(1+m)/G(1+i)+m(1-mi).

Equations (<ref>), (<ref>) and (<ref>) are exact to all orders in 1 / G.
i is a function of G, according to Eq. (<ref>).
In the approximation α_1=0, with exponential accuracy, the function θ (x) is given by

    θ(x)=φ(x)/x≈sin(β_20 x)/β_20 x,

in agreement with the analysis of Sect. II.B for unit polytropic indices.
The function θ (x) is regular for G = ∞ but inherits the singularities of β_ 20.

A variation in the electron density at r=0 leads to a variation e  = e - e^reg≠ 0,
corresponding variations of i and β_± and, ultimately, to αk 1≠ 0.
We return to the general solution (<ref>) and write the function k in the form

    k(x)=αk 1sinh(β_1 x)/x + αk 2sin(β_2 x)/x,

where, with account of Eq. (<ref>),

    αk 2 = 1-αk1β_1/β_2.

The regular solution (<ref>) vanishes for x_0 = π / β_ 20.
For αk1≠ 0, one of the functions k (x) vanishes
already at x_ b = x_0 +  x, whereas the second function is still finite.
We write this condition in the form

    i(x_b)   =   0,
    e(x_b)   =   eb,

where eb is the value of the electron function at the boundary. We assume first that an electrosphere is formed.

As follows from Sect. III, the electron pressure at the boundary of the star is Peb∼Pun / G,
where Pun is the characteristic pressure scale. The relation Pe∼ne^1 + 1/ e∼e^1+ e
implies that eb∼G^- 1 / (1 + e)≪ 1. For the case under consideration, e = 1, and accordingly,
eb∼ 1 / √(G). The typical scale of the electrosphere is ra / r_0∼eb∼G^-1 / (1 + e),
in agreement with the estimates given in Eqs. (<ref>) and (<ref>). Using Eq. (<ref>) and the inequality W_0 ≥ 0, we obtain


    eb≤λG^-1/(1+e)[ ( 1+e) (qs+λm)^2/2e(1+m)(1-λm^2/λG)( Ms/
    4π r_0^3ni0mi) ^2( 
    r_0/Rs) ^4] ^1/(1+e).



Equations (<ref>) and (<ref>) allow to express the unknown coefficients αk 1 in terms of eb.
Given the relation in Eq. (<ref>), the functions k (x) are completely defined, but the parameters x and eb
on which these functions depend must be found. We first substitute αk 1 into the first of Eqs. (<ref>) and (<ref>), setting β^2 = β_1^2, and then substitute αk 2, setting β^2 = - β_2^2. The two resulting equations are linear in eb. By excluding eb, we find

    (1+i)β_2^2 - 1 - me/β_2^2  - mesinh( β_1 xb)/β_1 xb
       =    (1 + i)β_1^2 +1+me/β_1^2 + mesin( β_2 xb )/β_2 xb.

This equation defines xb as a function of the parameters e, i, m, and G.
Equation (<ref>) is not used because for β^2 = β_±^2, it is not independent. Both parts of Eq. (<ref>)
define eb≥ 0; therefore, a solution can exist for β_2^2 ≥β_20^2 and x ≤ 0.
As noted above, eb∼ 1 /β_1∼ 1 / √(G), and so
x ∼ra / r_0∼ 1 / √(G).
Considering x ∼ 1 / β_ 1≪ 1, the right-hand side of Eq. (<ref>) can be expanded into a series
in the neighborhood of xb = x_0; as a result, we find

    x = - 1/β_1𝒲( eexp(πβ_1/β_2)/2(1 + i)).

The displacement x is expressed in terms of the Lambert W-function, which gives a solution to the equation 𝒲 (x) exp (𝒲 (x)) = x. Since the argument of the function is of order unity, the correction to e is exponentially small:

    e∼exp(-O(√(G))).

For eb, we obtain

    eb = -   x /x_0(1 + i).

In Eqs. (<ref>) and (<ref>), the ∼exp (- πβ_ 1 / β_2) terms and
higher terms in 1 / √(G) are omitted. At this level of accuracy, we can replace β_1 and β_2
with the expressions in Eqs. (<ref>) and (<ref>) everywhere, including exponents. The positive definiteness of eb implies that
x <0 and e > 0.

For asymptotically large values of the argument of 𝒲 (x), the expansion 𝒲 (x) = ln (x) - ln (ln (x)) + … is valid.
A correction e∼ 1 / √(G) could lead
to the estimate x ∼eb∼ 1, which contradicts the upper limit of Eq. (<ref>).
The set of admissible solutions, therefore, consists only of solutions with exponentially small deviations e.

Equations (<ref>) and (<ref>) also permit solutions corresponding to an ionosphere.
Instead of Eqs. (<ref>) and (<ref>), we require

    i(x_b)   =   ib, 
    e(x_b)   =   0.

The equation for determining x is as follows:

    (1+i)β_2^2 - 1 - me/β_2^2 i - 1sinh( β_1 xb)/β_1 xb
       =    (1 + i)β_1^2 +1+me/β_1^2 i + 1sin( β_2 xb )/β_2 xb.

Both parts are equal to ib> 0. The relationship between x, e and ib
is established by the relations (<ref>) and (<ref>) with the substitutions
eb→ib, e→ - e and i→ 1 / i, while the other parameters remain unchanged.
The condition ib> 0 implies that x <0 and e <0.
An upper limit on ib is found from the inequality W_0 ≥ 0
and the balance between the pressure and the forces of gravitational attraction and electrostatic repulsion acting on the ions at the boundary of the electron component (cf. (<ref>)):


    ib≤λG^-1/(1+i)[ (1+i)(- qs+ λm/m)^2/2i(1+m)(1-1/G)( Ms/4π r_0^3ni0mi) ^2( r_0/Rs) ^4] ^1/(1+i).


An ionosphere exists provided that qs < λm / m = A / Z.

For positive values of e, the substance in the center of the star
is closer to the electrically neutral state than in the regular solution.
On the surface of the baryon component, the bulk charge is positive, and the electrosphere compensates for this charge either partially, completely, or excessively.
Suppression of the electron fraction increases the deviation from electroneutrality and corresponds to e < 0.
Under the condition qs < A / Z, the electron density vanishes when the ion density is still finite.
On the surface of the electron component, the stellar charge is positive (by continuity) and increases with | e |.
The ionosphere makes an additional contribution to this charge.
| e | can be increased until it violates (<ref>) in the case of an electrosphere or
(<ref>) in the case of an ionosphere.




  §.§.§ Charge-mass-radius relation



The general properties associated with a deviation from local electroneutrality are expected to be inherent in realistic models. We thus consider the charge-mass-radius relation for stars. The mass and radius are expressed as follows:

    Ms/4π r_0^3 n_i0m_i   =   π (1 + i)/β_2, 
    Rs/r_0   =   π/β_2.

The violation of the LEC affects the mass and radius at higher orders of the expansion in 1 / √(G).
The correction to the radius is ΔRs = r_0  x ∼ r_0 / √(G),
so the mass-radius relation remains unchanged to within ∼ 1 / √(G).
The bulk charge of the star, qs, and the electrosphere charge, qe, are as follows:


    Z/Aqs   =   1-mi/1+i-
    (1+m)i/(1+i)^2πΘeb, 
    Z/Aqe   =    -1+m/1+i[ -√(( 1-i/1+iπΘeb) ^2-i/
    1+iπ ^2Θeb^2)+( 1-i/
    1+iπΘeb) ],


where Θeb = β_1 eb / β_2 ∼ 1.
Both formulae are valid with an accuracy of ∼ 1 / √(G). The existence of an electron shell
requires an excess of electrons relative to the solution eb = ib = 0.
Accordingly, the second term in qs is negative, and as Θeb increases, its contribution reduces
the stellar charge. By equating the root expression to zero, we can obtain the maximum value of Θ_eb:

    Θeb^max = (√(1 + 1/i) - 1)(1 + i)/π.

This value corresponds to the upper limit in Eq. (<ref>). For Θeb = Θeb^max, the shell overcompensates
for the positive charge of the star: (qs + qe)^min = - λm <0.

An ionosphere forms for a lower concentration of electrons.
The bulk charge, qs, measured at the boundary of the electron component and the ion envelope charge, qe, are given by the expressions


    Z/Aqs   =   1-Λ _mΛ _i/1+Λ _i
    +(1+Λ _m)Λ _i/(1+Λ _i)^2πΘib, 
    Z/Aqe   =   (1+Λ _m)Λ _i/1+Λ _i[ -√(( 1-1/1+Λ _iπΘib) ^2-1/
    1+Λ _iπ ^2Θ _ib^2)+( 1-1/
    1+Λ _iπΘib)],


where Θib = β_1 ib / β_2 ∼ 1.
The charge qs increases with Θib. The radicand vanishes
when Θib takes its maximum value

    Θib^max = (√(1 + i) - 1)(1 + 1/i)/π.

One can check that for Θib = Θib^max, the inequality in Eq. (<ref>)
becomes an equality. The maximum total charge is equal to (qs + qe)^max = A / Z.
The condition qs < A / Z constrains the value of Θib from above by
(1 + i) / π. However, this limit is weaker than Θib≤Θib^max.

Figure <ref> plots the bulk charge of the star, qs, and the total charge including the envelope, qs + qe,
as functions of the difference Θeb - Θib
in the interval (- Θib^max, Θeb^max) for a mixture of electrons and protons, i.e., for
Z = A = 1 and i = m. The region Θeb - Θib∼ 0 is presented separately in Fig. <ref> on an enlarged scale.

In summary, in the two-fluid model with unit polytropic indices,
it is possible to explicitly construct the general solution describing charged stars with electro- and ionospheres.
An electrosphere forms if the ion chemical potential vanishes at the boundary, i (xb) = 0,
whereas the electron chemical potential at the boundary is finite, e (xb)> 0, according to Eqs. (<ref>) and (<ref>).
An ionosphere forms in the opposite case; see Eqs. (<ref>) and (<ref>).
The solution for the internal stellar structure is unambiguously matched with the solution for the electro- or ionosphere.
The guiding parameter of the problem is e = e - e^reg.
The general solution also describes stars with zero total charge, but such stars are not locally neutral.
In an exponentially small neighborhood of e^reg, the electron and ion densities
can be varied without restrictions. The stellar envelope is sensitive to exponentially small deviations of e
from e^reg.
The same fact can be interpreted as indicating remarkable overall stellar stability: changes in the envelope, and, accordingly, in the total stellar charge (within acceptable limits) have an exponentially weak influence on the internal structure of a star.






















§ GENERAL SOLUTION TO THE UNCONSTRAINED HYDROSTATIC EQUILIBRIUM EQUATIONS





From the analysis of the exactly solvable model, the specific nature of the problem becomes rather transparent:
the general solution to Eqs. (<ref>) and (<ref>) contains a component that is irregular at G = ∞.
In the model of Sect. IV, the regular component dominates up to the surface of the star,
whereas the irregular component is exponentially small everywhere except for the subsurface layer
Δ r ∼ra∼ r_0 /√(G), where
the contributions of the two components are comparable. The exponential smallness of the irregular component makes perturbation
theory appropriate for describing the inner layers of the star.

There are four independent parameters, e, i, m, and G,
on which Eqs. (<ref>) and (<ref>) and their solutions depend.
The regular solution is a three-parameter solution with the independent parameters i, m, and G.
The function e = e^reg is defined in Sect. II.C to O(1/G ).
The general solution and the regular solution satisfy the equations with different values of e.
We denote the variance by e = e - e^reg, |e|  ≪ 1.
The parameter i is not modified since for the general case ei considered here,
e and i are independent.

Let k0 denote the regular solution. The general solution is assumed to be of the form

    k = k0 + χk,

where χk is a small correction. This representation is valid everywhere except within a thin subsurface layer.
By substituting k and k0
into Eqs. (<ref>) and (<ref>) and linearizing these equations in χk and e, we obtain

    _x(χe+iχi)   =   -(ηii0^i-1χi
          + me0ee0^e-1χe +mee0^e), 
    _x(χe-miχi)   =   -G(ii0^i-1χi
          -e0ee0^e-1χe -ee0^e).

The initial conditions

    χk(0)=χk^'(0)=0

show that e determines the overall scale of χk.




 §.§ Irregular component in the WKB approximation



When both sides of Eq. (<ref>) are divided by G,
the Laplacian acquires the small coefficient 1 / G.
A similar case occurs in the Schrödinger wave equation, in which the kinetic term with the Laplacian contains ħ^2.
In the semiclassical limit, ħ is small, which makes the similarity of Eqs. (<ref>) and (<ref>)
to the semiclassical limit of the Schrödinger wave equation quite obvious. In quantum mechanics,
the Wentzel, Kramers, and Brillouin approximation is used to describe the semiclassical limit; this approach is known as the WKB method (see, e.g., <cit.>).

Following the general scheme, we look for a solution of the form

    χk^irr(x) = gk(x) exp(√(G) S(x)).

In the neighborhood of G = ∞, the functions gk admit a series expansion in 1 / √(G):

    gk(x) = gk 0(x) + gk 1(x)/√(G) + gk 2(x)/G + … .

The exponential factor causes the point G = ∞ to be an essentially singular point of χk^irr (x).
Applying _x to χk^irr (x) yields

    _xχk^irr   =   exp(√(G) S )
    
       ×   [Ggk (S^' )^2  + √(G)(x^2 gk^2 S^')^'/x^2 gk
    +(x^2 gk^')^'/x^2].

The coefficients of powers of √(G) are as follows:

    G   :     gk0(S^' )^2,    
    √(G)   :     gk1(S^' )^2 + (x^2 gk0^2 S^')^'/x^2 gk0,  
    
    1   :     gk2(S^' )^2 + (x^2 gk1^2 S^')^'/x^2 gk1+(x^2 gk0^')^'/x^2,

where the common exponential factor is suppressed.

The G term in Eq. (<ref>) yields the relation

    ge0+igi0=0.

By substituting this relation into Eq. (<ref>), we obtain with the same accuracy

    S^' = ±√(ii0^i - 1  + ie^regee0^e - 1/i (1 + m) ).

The exponent can be found by integrating Eq. (<ref>). In the model with
unit polytropic indices,
this expression gives a simple linear dependence,  S ∼ x.

The √(G) term in Eq. (<ref>) leads to the relation

    ge1+igi1=0.

Using this equality in Eqs. (<ref>) and (<ref>),
we find that the second differential term in Eq. (<ref>) vanishes. We thus obtain

    x^2 ge0^2 S^'=const.

Near the center, gk0 behaves as gk0∼ 1/x.
Equation (<ref>) leads to

    S_±(x) = ±∫_0^x√(ii0^i - 1(x^')+
    ie^regee0^e - 1(x^')/i (1 + m)) dx^'.

The terms with S_± appear in the expression for χk0^irr (x ) in a certain combination.
The singularity at x = 0 can be eliminated by taking the difference
exp(√(G) S_+ (x) ) - exp(√(G) S_-(x) ).
Finally,

    χk0^irr(x) = Ck0sinh(√(G)S_+(x))/x√( S_+^'(x))  + O(G^-1/2),

where the Ck 0 are arbitrary constants.

The higher-order terms of the expansion
can be found in a similar manner.

The functions χk0^irr (x) do not vanish at x = 0. The normalization condition k(0) = 1
is restored when considering a particular solution of the inhomogeneous ODE system defined by Eqs. (<ref>) and (<ref>).




 §.§ Correction to the irregular component


A particular solution to Eqs. (<ref>) and (<ref>) can be found
from the original ODE system defined by Eqs. (<ref>) and (<ref>).
We can construct the regular solution to this system, starting from the LEC solution,
in the form of a series expansion in 1 / G.
The functions k0 (x) and e^reg that are obtained in this way
depend on the parameters i, m, and G.

Suppose that we are looking for the regular solution to Eqs. (<ref>) and (<ref>)
with the initial conditions k (0) ≠ 1 and k^' (0) = 0
for |k (0) - 1| ≪ 1.
Let θk0 (x) denote the corresponding solutions,
which depend on the parameters i, m, and G.
The same parameters determine e≡e^reg'.

It is not difficult to see that by rescaling x and θk0(x),
the problem can be reduced to the initial problem with the standard initial conditions, the rescaled variable x,

    x=x√(i^i(0)/e(0)),

and the new parameters e and i:

    Λe=e^reg'e^e(0)/i^i(0),  Λi=ii(0)/e(0).

By contrast, m is not modified. Obviously, Λe is a known function of Λi,
namely, Λe = e^reg (Λi); thus,

    e^reg'(i) = i^i(0)/e^e(0)e^reg(ii(0)/e(0)).

Under the assumption of the smallness of k = k (0) - 1, we obtain

    e = e^reg' - e^reg
       =   i(ie^reg+
    i∂e^reg/∂i) 
       -   e(ee^reg+
    i∂e^reg/∂i).

The solution to Eqs. (<ref>) and (<ref>) with the non-standard boundary conditions can be written as follows:

    θk0(x)=k(0)k0(x√(i^i(0)/e(0)),
    ii(0)/e(0)).

The function

    χk0 = θk0-k0

satisfies the linearized equations (<ref>) and (<ref>) for e = e^reg' - e^reg and the initial conditions
χk0(0) = k.
For small k, we can write

    χk0 (x)    =   k0(x)k + x ∂k0(x)/∂ xii-e/2
       +   i∂k0(x)/∂i (i-e).




 §.§ General solution


The irregular component of the general solution is given by the first term of the expansion in Eq. (<ref>).
Equation (<ref>) is a solution to the homogeneous equation,
whereas Eq. (<ref>) is a solution to the non-homogeneous equation.
The general solution is defined by Eq. (<ref>) with χk = χk^irr + χk0.
Among these functions, k0 and χk0 are regular at G = ∞, whereas χk^irr is singular at G = ∞.
The initial conditions (<ref>) lead to the relation

    Ck0√(GS_+^'(0)))+k=0,

where the Ck0 are the constants that appear in Eq. (<ref>).
The derivative of χk0 at x = 0 vanishes by construction.
In the general case, S_± (x) = x S_±^' (0) + O (x^3);
consequently, the derivative χk^irr also vanishes automatically at x = 0.
As a consequence of Eq. (<ref>), the coefficients Ck0 satisfy

    Ce0 + iCi0=0.

By expressing the unknown coefficients in terms of Ci0, we can write


    χi(x)    =   Ci0/x√(S_+^'(x))sinh(√(G) S_+ (x))
    - Ci0√(GS_+^'(0))[i0(x) + x ∂i0(x)/∂ xi+i/2+i∂i0(x)/∂i (1 +i)], 
    χe(x)    =   -Ci0i/x√(S_+^'(x))sinh(√(G)S_+(x))
    + Ci0√(GS_+^'(0))[ie0(x) -
     ∂e0(x)/∂ xi+i/2-i∂e0(x)/∂i (1+i)].


In each of these expressions, the first term increases exponentially with x, dominates the second (regular) term for x ≳ 1/√(G),
and becomes comparable to the regular component k0 near the surface in the layer | x - xb | ∼eb,
 
where perturbation theory in χk is no longer applicable.
The constant Ci0 is exponentially small.
This feature has been noted earlier in the analysis of the exactly solvable model of Sect. IV.
The exponential smallness of the irregular solution in the inner layers of the star is obviously a common property.
The functions χi (x) and χe (x) differ in sign, and the electron component is suppressed by approximately i.

The matching of k and k^' with the electrosphere solution requires

    χe(xb) ∼eb,   χe^'(xb)∼ 1.

The first of the equations allows to fix Ci0.
For i(xb) = 0, Eq. (<ref>) leads to S_+^'(xb) ∼eb^(e - 1)/2;
the second equation then follows automatically.
One can verify that for ηe = ηi = 1, χk^irr
is analytic in G∈ℂ^1 ∖∞, with G = ∞ being a unique singular point of χk^irr,
in agreement with the analysis of Sect. IV.B.

Standard numerical methods for solving ODE systems do not allow the separation of the regular solution from the irregular one
due to the exponential smallness of the latter. A numerical solution to Eqs. (<ref>) and (<ref>) for x ≲ 1
will inevitably contain an admixture of the irregular component, which is interpreted by the program as a machine zero.
In the subsurface layer | x - xb | ∼eb, due to the exponential growth of χk^irr (x),
this machine zero is converted into a finite but random value.
The stellar properties that are sensitive to e cannot be determined using numerical methods.
A change in the integration, e.g., of the grid step, is thus expected to result in uncontrolled changes in the stellar charge, the thickness of the electro- or ionosphere, and other quantities that depend on e. In view of these remarks,
the numerical solutions of Refs. <cit.> are likely to be unstable.

At the same time, as the example of Eqs. (<ref>) and (<ref>) shows,
relationships which do not depend on e are numerically stable.
Each integration with a varying grid step
leads to random values of qs and eb, but all points on the plane (qs, eb) must lie on a smooth curve.
In the exactly solvable model of Sect. IV, the quantities qs and eb are related by a linear dependence.





§ CONCLUSIONS





In this paper, we considered the unconstrained HE equations for stars in the absence of local electroneutrality.
The condition of local electroneutrality is usually imposed on the EoS phenomenologically as an additional constraint (LEC).
In our approach, the electrostatic interactions are considered ab initio;
as a result, the LEC is satisfied approximately, with the uncompensated charge density being 1/λG∼ 10^-36
of the particle number density.
The smallness of the parameter 1/λG makes it possible to use the LEC to describe with a high accuracy properties that are not related to stellar electrostatics.

We found that the mass-radius relationship and the condition for gravitational stability for charged stars
are not modified to an accuracy of ∼λG^-1/(1 + e).
The admissible states of charged stars are determined by the condition that
the electro- and ionospheres do not extend to infinity.
This  requirement provides lower and upper bounds on the total stellar charge of -0.1 and 150 C per solar mass, respectively.

The main system of unconstrained HE equations, Eqs. (II.12) and (II.13), belongs to the class of singularly perturbed ODEs.
In such a system, the perturbation parameter (in our case, 1/λG)
appears as a coefficient of the higher derivative and thus plays the main role in specifying the nature of the problem.
The use of conventional (regular) methods for solving ODEs via series expansion in a small parameter
leads to the loss of certain solutions in this case.
In the limit of λG→∞, Eq. (II.13) becomes a constraint on the particle concentrations,
which we identify as an LEC.

The Poincaré theorem on analyticity and Dyson's argument indicate that the general solution to Eqs. (II.12) - (II.13)
is singular for λG = ∞.
A particular, regular solution can exist provided that one of the parameters, e.g., e,
is replaced by a function of λG.
Such a solution can be constructed as a formal power series in G (∼ 1/λG) starting from a locally neutral solution.
Only these two types of solutions, namely, regular ones and locally neutral ones, have been discussed in the literature to date.

In this paper, the general solution to the HE problem was first constructed in NGT.

In Sect. IV, we presented an exactly solvable model of a two-component mixture of ions and electrons with unit polytropic indices.
The general solution to the unconstrained HE equations appears to be irregular in the neighborhood
of λG = ∞. The irregular component is exponentially small in the inner layers of the star.

In the general case, the WKB method can be used to construct the irregular component.
Its magnitude is comparable to that of the regular component in a subsurface layer of thickness ∼λG^-1/(1 + e),
where these components interact non-linearly and determine the structure of the electro- or ionosphere.
The properties of these envelopes are exponentially sensitive to variations in the charge density in the central regions of the star.

In general, stars are characterized by their mass, radius and charge.
Nuclear reactions in the interior of a star cause variations in the parameters e, i, and m
and thereby affect the bulk stellar charge,
which, due to the conservation of electric charge, leads to a charge exchange between the envelope and the inner region.

The structure of the general solution to the unconstrained HE equations
can be described as follows:

Starting from a locally neutral solution, one can construct a regular solution  (k0 in Sect. V) in the form of a series in powers of G.
On the basis of this solution, an irregular component can be constructed using the WKB method (χk^irr in Sect. V).
This component is singular for λG = ∞.
A correction to this irregular component is needed to satisfy the initial conditions (χk0 in Sect. V).
The sum of these three functions yields the general solution, excluding only a thin subsurface layer of thickness ∼λG^-1/(1 + e),
where the magnitude of the irregular component becomes comparable to that of the regular component and the WKB method is no longer applicable.
Numerical methods should then be used to match, within the subsurface layer, the general solution
with the known solution for the electro- or ionosphere.

Extension of the formalism developed here to GRT
is left as an open problem.
An important direction of study is generalization to
three or more component fluids in chemical equilibrium.
Yet another interesting application is to find the general solution of the unconstrained HE equations for strange stars.



The authors are grateful to S. I. Blinnikov and I. L. Iosilevsky for stimulating discussions.
M.I.K. wishes to acknowledge M. D. Malykh for discussion of analytic properties of the solutions.
The work is supported in part by the RFBR Grants Nos. 16-02-00228 and 18-02-00733.






§ INCOMPRESSIBLE TWO-FLUID MODEL
 



As a second example of an exactly solvable, albeit somewhat exotic, model,
we consider the case of ηi = ηe = 0, which corresponds to an incompressible fluid
(n = n_0 θ^η∼ const). 
Although the speed of sound in such a substance is infinite,
this case is of methodological interest with regard to the restrictions on the parameter e in the regular solution,
especially as S_± (x) → 0 for ηi = ηe = 0.
For an incompressible fluid, Eqs. (<ref>) and (<ref>) become

    _x(e+ii)   =   -(1+me),
    _x(e-mii)   =   -G(1-e).

Their solutions for the initial conditions k (0) = 1 and k^'(0) = 0 have the form

    i   =   1-(1+me)-G(1-e)/6i(1+m) x^2, 
    e   =   1-m(1+me)+G(1-e)/6(1+m) x^2.

In principle, these solutions describe a charged stellar configuration free from constraints on e.
This freedom is due to the irregular nature of the solution for G = ∞.
The requirement of regularity leads to the LEC e = 1, given that e is a parameter.
If e is considered to be a function of G, then e = 1 + O(1/G).
The solutions (<ref>) and (<ref>) are valid in the inner region of the star.
The outer region requires special consideration.

According to Eq. (<ref>), the ion concentration vanishes for

    xi=√(6i(1+m)/(1+me)-G(1-e)).

If we require the electron concentration to vanish at the same point, xe = xi,
then we obtain e = e^reg, where e^reg is given by Eq. (<ref>).
In this case, the star still has a charge. Suppose now that we do not require the concentrations to vanish simultaneously.
Again, let us consider the outer part of the star, the electrosphere, where the ion concentration is zero.
The equation for e reduces to

    (1 + m)_x(e)=e(G-m^2).

At the boundary x = xi, the function e and its first derivative are continuous.
The solution to Eq. (<ref>) takes the form

    ea=1-xi^2(G+m)/2(1+m)(1-2xi/3x)
    +e(G-m^2)/6(1+m)x^2.

For an electrically neutral star, the boundary x = xe, where ea (xe) = 0,
determines the total charge. The bulk charge is proportional to
(cf. Eq. (<ref>))

    Qs∼ -x^2(e'-mii')|^xi_0,

whereas the charge of the envelope is proportional to

    Qe∼ -G(1+m)/G-m^2x^2θea^'|^xe_xi.

The condition Qs + Qe = 0 leads to the relation

    xe=xi/√(e),

whereas the condition θea (xe) = 0 leads to the following equation:

    e(G+m)+i√(e)(3G+2m-m^2) 
                      =
    G-1+3i(G+m),

whose solution, e, coincides with e^reg to an accuracy of 1/G.

99

Lebowitz1969 J. L. Lebowitz and E. H. Lieb, Phys. Rev. Lett. 22, 631 (1969).
Dyson1967 F. J. Dyson, A. J. Lenard, Math. Phys. 8, 423 (1967).
Lieb1975 E. H. Lieb, W. Thirring, Phys. Rev. Lett. 35, 1116 (1975).
Landau1976 L. D. Landau, E. M. Lifshitz, Statistical Physics. Vol. 5, 3rd ed. (Butterworth-Heinemann, 1980).
Klein1949 O. Klein, Rev. Mod. Phys. 21, 531 (1949).
Kodama1972 T. Kodama, M. Yamada, Prog. Theor. Phys. 47, 444 (1972).
Olson1975 E. Olson and M. Bailyn, Phys. Rev. D 12, 3030 (1975).
Rosseland1924 S. Rosseland, Mon. Not. Roy. Astron. Soc. 84, 720 (1924).
Schatzman1958 E. L. Schatzman, White Dwarfs, (North-Holland Pub. Co., NY, Int. Publ., 1958).
Olson1978 E. Olson and M. Bailyn, Phys. Rev. D 18, 2175 (1978).
Pikelner1964 S. B. Pikel'ner, Fundamentals of cosmic electrodynamics, (Fizmatgiz, Moscow, 1961; NASA technical translation. TT F-175, 1964).
Bally1978 J. Bally, E. R. Harrison, Astrophys. J. 220, 743 (1978).
Neslusan2001 L. Neslušan, Astron. Astrophys. 372, 913 (2001).
Rotondo2011 M. Rotondo, J. A. Rueda, R. Ruffini, S.-S. Xue, Phys. Lett. B 701, 667 (2011).
Belvedere2012 R. Belvedere, D. Pugliese, J. A. Rueda, R. Ruffini, S.-S. Xue, Nucl. Phys. A 883, 1 (2012).
Belvedere2014 R. Belvedere, J. A. Rueda and R. Ruffini, J. Korean Phys. Soc. 65, 897 (2014).
Gorshkov2008 A. B. Gorshkov, V. A. Baturin, Astron. Rep. 52, 760 (2008).
Beznogov2013 M. V. Beznogov, D. G. Yakovlev, Phys. Rev. Lett. 111, 161101 (2013).
Nadyozhin2005 D. K. Nadyozhin, V. S. Imshennik, Int. J. Mod. Phys. A 20, 6597 (2005).
Krivoruchenko2011 M. I. Krivoruchenko, D. K. Nadyozhin, T. L. Razinkova, Yu. A. Simonov, M. A. Trusov, A. V. Yudin, Phys. Atom. Nucl. 74, 371 (2011).
Nadyozhin2016 D. K. Nadyozhin, Mechanisms of supernova explosion: modern status, in: Proceeings of the Conference
Quark Phase Transition in Compact Objects and Multimessenger Astronomy: Neutrino Signals, Supernovae and Gamma-Ray Bursts, 7-14 October, 2015, the Special Astrophysical Observatory of Russian Academy of Sciences (SAO RAS) and Baksan Neutrino Observatory of Institute for Nuclear Research of Russian Academy of Sciences (BNO INR RAS). Eds. V. V. Sokolov et al. (Sneg, Pyatigorsk 2016), pp. 69 - 75.
Ray2003 S. Ray, A. L. Espíndola, M. Malheiro, J. P. S. Lemos and V. T. Zanchin, Phys. Rev. D 68, 084004 (2003).
Ilyashenko2008 Y. Ilyashenko, S. Yakovenko, Lectures on Analytic Differential Equations, (American Mathematical Society, 2008), Theorem 1.1.
Dyson1952 F. J. Dyson, Phys. Rev. 85, 631 (1952).
Malley1990 R. E. O'Malley, Jr., Singular Perturbation Methods for Ordinary Differential Equations, (Springer-Verlag, New York 1991).
Horedt2004 G. P. Horedt, Polytropes. Applications in astrophysics and related fields, (Kluwer Academic Publishers, Boston 2004).
Yoselevsky2009 I. L. Iosilevskiy, J. Phys. A 42, 214008 (2009).
Yudin2018 A. V. Yudin et al., in preparation (2018).
Chavalier2007 G. Chavalier, The Earth's Electrical Surface Potential. A Summary of Present Understanding.
California Institute for Human Science, Encinitas (2007).
Iorio2012 L. Iorio, Gen. Relativ. Gravit. 44, 1753 (2012).
Witten1984 E. Witten, Phys. Rev. D 30, 272 (1984).
Hu2002 J. Hu and R. X. Xu, Astron. Astrophys. 387, 710 (2002).
Krivoruchenko1995 M. I. Krivoruchenko, B. V. Martemyanov, Yad. Fiz. 58, 536 (1995) [Phys. Atom. Nucl. 58, 484 (1995)].
LandauV3 L. D. Landau, E. M. Lifshitz, Quantum Mechanics: Non-Relativistic Theory. Vol. 3, 3rd ed. (Pergamon Press, 1977).
Bekenstein1971 J.D. Bekenstein, Phys. Rev. D 4, 2185 (1971).
Lasky2007 P. D. Lasky and A. W. C. Lun, Phys. Rev. D 75, 104010 (2007).
Ghezzi2005 C. R. Ghezzi, Phys. Rev. D 72, 104017 (2005).
Negreiros2009 R. P. Negreiros, F. Weber, M. Malheiro, and V. Usov, Phys. Rev. D 80, 083006 (2009).
Arbanil2015 J. D. V. Arbanil and M. Malheiro, Phys. Rev. D 92, 084009 (2015).
Krivoruchenko1987 M. I. Krivoruchenko, JETP Lett. 46, 3 (1987) [Pisma Zh. Eksp. Teor. Fiz. 46, 5 (1987)].
Krivoruchenko1991 M. I. Krivoruchenko, B. V. Martemyanov, Astrophys. J. 378, 628 (1991).



