(1) must be fully neutralized by the negative charge of electrons. However, since electrons are light and interact only via the electromagnetic force, they will penetrate through the boundary and generate a local charge disbalance around the star surface. The induced electrostatic potential φ ( z ) is determined from the Poisson equation (see e.g. [14]) where e = √$_{4}$$_{πα}$$_{=0.3028 is the proton charge and}$ α is the fine-structure constant[26]. For a given proton distribution, Eq. (1), the electron charge distribution ρ$_{e}$ ( z ) should be determined self-consistently. For this purpose we use the Thomas-Fermi approximation [15, 16], which should work well for an extended object like a heavy nucleus or star. The relativistic version of this method was considered e.g. in refs. [17, 18]. In a semi-classical approximation the electron energy at point z can be written as where k is its 3-momentum and − V ( z ) = − eφ ( z ) is the potential energy. At zero tem- perature all electronic states with k ≤ k$_{F}$ ( z ) are occupied, where k$_{F}$ ( z ) is the local Fermi momentum. It is determined from the condition that gives The local electron density is found by integration over k : By inserting this expression into Eq. (2) we obtain a non-linear differential equation for φ ( z ). To find its solutions we need to specify boundary conditions. According to our assumption, the boundary is in direct contact with the vacuum and therefore the electron density as well as the electric potential must vanish at z → ∞ . This condition can be fulfilled if µ ≤ m$_{e}$ . For µ < m$_{e}$ the electrons are bound to the surface, i. e. a finite energy is needed to extract an electron from the star. This energy is well known in ordinary metals as the exit work. In this case the electron density should vanish at a certain d $^{2}$φ dz 2 = − e [ ρ$_{p}$ ( z ) − ρ$_{e}$ ( z )] ≡ − eρ$_{ch}$ ( z ) , (2) ϵ ( k , z ) = √ k 2 + m 2 e − V ( z ) , (3) ϵ ( k$_{F}$ ( z ) , z ) = √ k 2 $_{F}$( z ) + m 2 e − V ( z ) = µ = const , (4) k$_{F}$ ( z ) = √ [ µ + V ( z )] 2 − m 2 e . (5) ρ$_{e}$ ( z ) = 2 k$_{F}$ ( z ) ∫ 0 d $^{3}$k (2 π ) 3 = [ ( µ + V ( z )) 2 − m 2 e ] 3 / 2 3 π 2 . (6) 3