Two-Fluid Schwarzschild Solution

The paper presents the interior two-fluid Schwarzschild solution, a model for compact (neutron) stars with admixed dark matter. The two constant energy densities generalize the interior Schwarzschild solution to the case of two fluids (representing standard model matter and dark matter, for instance), interacting mutually solely by the common gravitational field. In general, the two-fluid core is surrounded by a one-fluid corona (envelope). Only for special parameters (the model depends on energy densities e1, e2 and central pressures pc1, pc2 of the fluids), both fluids occupy the same region. Despite the sum of two energy densities in that region, the compactness bound of 8/9 is not exceeded. A conceivable object is one with a stark leakage of one fluid (e.g. dark matter) beyond the two-fluid core.

The paper is organized as follows: Section 2 presents the analytical solution of the two-fluid TOV equations. The impact of DM on the pressure profile is considered in Section 3. The special case of two fluids, uncovering the same spatial region, is dealt with in Section 4, to show that the Buchdahl limit is not exceeded. The general case of leaking out one of the fluids beyond the other one is sketched in Section 5. We summarize in Section 6.

In Section 2, the two-fluid TOV equations are given in geometric units:

m′1 = 4πr2e1, m′2 = 4πr2e2,

p′1 = −(e1 + p1)Φ, p′2 = −(e2 + p2)Φ,

Φ ≡ (m1 + m2 + 4πr3p1 + 4πr3p2) / (r2(1 − 2m1/r)(1 − 2m2/r)) = (4π/3)(e+ + 3p+)r / (1 − (8π/3)e+r2),

where e+ ≡ e1 + e2 and p+ ≡ p1 + p2, and e− ≡ e1 − e2 and p− ≡ p1 − p2.

The interior two-fluid Schwarzschild solution reads:

p+(r) = −(C − W(r))e+ / (3C − W(r)),

p−(r) = ((3C − 1)/(3C − W(r)))(e− + p−c) − e−,

where C ≡ (e+ + p+c)/(e+ + 3p+c), W(r) ≡ √(1 − (8π/3)e+r2), and the subscript "c" refers to the central pressure(s).

The pressure profiles follow as:

2p1(r) = p+ + p− = (−(C − W)e+ + (3C − 1)(e− + p−c) − (3C − W)e−) / (3C − W),

2p2(r) = p+ − p− = (−(C − W)e+ − (3C − 1)(e− + p−c) + (3C − W)e−) / (3C − W).

These pressure profiles are linked by p1(r) = −e1 + ((e1 + pc1)/(e2 + pc2))(e2 + p2(r)), that is, p1 is a linear function of p2.

If e1 pc2 = e2 pc1, then the zero of p2(r) coincides with zero of p1(r) at the same value of the radial coordinate r = R. Moreover, on this special hypersurface in parameter space (also called quadric cone) p1/p2 = (e1 + pc1)/(e2 + pc2) = pc1/pc2 = const; const = 1 is a further special case thereof.

If one is interested in the individual profiles p1,2(r) and the ordering of the pressure zeroes, then Equation (10) is useful too:

(i) p2(R) = 0 ⇒ p1(R) > 0 for e1 pc2 e2 pc1.

Equation (10) may be cast in the form p1(r)/pc1 = ((e1 + pc1)pc2)/((e2 + pc2)pc1) + (δ/((e2 + pc2)pc1))(p2(r)/pc2), where δ ≡ −e1 pc2 + e2 pc1.

As seen in Equation (7), W(r) ∈ [0, 1] must be respected since W(r = 0) = 1 and W(r = rH ≡ √(3/(8πe+))) = 0. Accordingly, p1(W = 1) = pc1 and p2(W = 1) = pc2. Otherwise,

p1(W = 0) = (e2(−e1 + 2pc1) − e1(e1 + pc1 + 3pc2)) / (3(e1 + e2 + pc1 + pc2)),

p2(W = 0) = (−e2(e1 + e2 + 3pc1 + pc2) + 2e1pc2) / (3(e1 + e2 + pc1 + pc2)).

Observe the symmetry 1 ↔ 2. At small values of e2, the leading orders are p1(W = 0) ∝ −e1(e1 + pc1 + 2pc2) 0 with common factors (3[e1 + pc1 + pc2])−1.

In contrast to the one-fluid Schwarzschild pressure profile, which extends from pc to −e0/3, the profiles p1,2(r) extend from pc1,c2 to p1,2(W = 0) given in Equations (12, 13). Depending on the four parameters e1, e2, pc1, pc2 the latter values can be positive or negative. In the former case, this signals the possibility of an unstable configuration since the pressure zero is not reached prior to the horizon rH. For instance, let be the parameters e1, pc1, pc2 fixed and consider the dependence on e2. Then, W2 = 0 is reached by e2crit given by e2crit = (1/2)(−b + √(b2 − 4c)), with b = e1 + 2p+c + p−c, c = e1(−p+c + p−c).

The pressure zeroes are generically at x1,2 = 1 − W1,22 with W1,2 = (Ce+ ± e− ∓ (3C − 1)p−c) / (e+ ± e−) or at R1,2 = 3x1,2/(8πe+). The scaled radial coordinate is here x = (8π/3)e+r2. One has separately to check whether W1,2 ∈ [0, 1].

The pressure profiles p1 and p2 cross each other at W× = (e− − (3C − 1)p−c)/e−, where the corresponding radius r× follows from Equation (7). One has separately to check whether p1(r×) = p2(r×) ≥ 0. The special case W1 = W2 = W× = C with p1(r×) = p2(r×) = 0 due to δ = 0 has been mentioned above.

The special setting e2 = 0 and pc2 = 0, i.e. e+ = e− = e1 and p+c = p−c = pc1, leads to the Schwarzschild pressure profile, p1(r) = −e1(C − W(r))/(3C − W(r)) and p2(r) = 0.

In Section 3, the impact of DM on pressure profiles is considered. The four-dimensional parameter space is constrained to a three-dimensional parameter space by using e+ for a scale setting. Accordingly, p1,2/e+ depends on e+r2 and parameters p+c/e+, p−c/e+, e−/e+. Figure 1 exhibits the scaled pressure profiles for four sets of parameters. The case p−c = 0 facilitates p1(r = 0) = p2(r = 0), see blue curves. However, due to e− ≠ 0, the profiles split off. If additionally e− = 0, then p1(r) and p2(r) are on top of another, see magenta curve(s). Fat curves depict p1(r)/e+, while thin curves are for p2(r)/e+. As stated above, the sign of the combination −e1pc2 + e2pc1 ≡ δ determines the order of the pressure zeroes. In the present examples, p2 for cases (i, ii) has its zeroes first (δ = 0.065, δ = 0.1), while in case (iii), p1 drops first to zero (δ = −0.035).

Another glimpse on the impact of DM on the two-fluid configuration is offered by a plot of scaled pressures, pi/pci (i = 1, 2), as a function of r/rH. Note that rH depends on both, e1 and e2. Figure 2 exhibits such pressure profiles for (i) e1, e2, pc1, pc2 = 1, 0.1, 0.5, 0.02 (red curves), (ii) 1, 0.2, 0.5, 0.2 (cyan curves) and (iii) 1, 0.1, 0.5, 0.2 (blue curves). p1/pc1 is exhibited by solid curves, while dashed curves are for p2/pc2.

To understand the parameter dependence let us cast Equations (2, 3) in the form:

(1/pci)(dpi/dξ) = −(1/2)(pci/ei)−1(pi/pci + 1)(ξ/(1−ξ))(1 + 3(p1+p2)/(e1+e2)),

where ξ = r/rH = 0...1. For the considered parameters, the term (pci/ei)−1 ranges from 5 to 1 to 0.5 for (i) to (ii) to (iii). That is, the gradient of p2 in ξ is largest for (i). Accordingly, given the small value of pc2, this scaled profile drops most rapidly, see red dashed curve. In contrast, the blue dashed curve is flattest, owing the smallest value of (pci/ei)−1. The term p+/e+ ranges from 0.47 to 0.58 to 0.63 for (i) to (ii) to (iii). Accordingly, the scaled pressure gradients of p1 change moderately, i.e. p1/pc1 drops somewhat faster for (iii), see blue solid curve. Altogether, the deviations to the Schwarzschild case (depicted by the solid black curve) are small. Measured by ξ = r/rH, the radii R1/rH = 0 of p1(R1/rH) shrink somewhat by the impact of DM. For the selected parameters, the radii of the DM component change much more drastically.

The considerations above are not yet final, since they uncover positive and negative pressures. The negative pressure branches must be truncated appropriately which requires more dedicated analyses.

In Section 4, the special parameter case where both fluids occupy the same region is considered. Naively, one would expect that the compactness C = 2M/R increases if two the energy densities e1 and e2 occupy the same region enclosed in the radius R1 = R2 = R, since M = (4π/3)(e1 + e2)R3. However, both fluids are coupled by the TOV equations which set the limit 2m1(r) + 2m2(r) < r when considering configurations with pc1,2 < ∞ and smoothly dropping pressure in going outward.

The condition for R1 = R2 was mentioned above, δ = −e1pc2 + e2pc1 = 0 or e−p+c = e+p−c. Then W1,2 = C from Equation (15). The compactness is accordingly C = (8π/3)(e1 + e2)R12 = 4X/(1 + 3X) with X ≡ (pc1 + pc2)/(e1 + e2). Whatever the values of e1,2 ≥ 0 and pc1,2 ≥ 0 are (supposed X ≥ 0), C runs from zero (at small X) to 8/9 (at large X), i.e. the Buchdahl limit as in the Schwarzschild case is recovered.

In Section 5, the corona surrounding a two-fluid core is considered. The pressure zeroes follow from Equation (15). If W1 = W2 due to δ = 0, then no corona is needed. This is case (1) dealt with in the previous section. If W1 ≠ W2, one has to distinguish two cases:

(2) p1(p2 = 0) > 0: δ > 0, and the corona consists of medium 1 with energy density e1.

(3) p1(p2 = 0) < 0: δ < 0, and the corona consists of medium 2 with energy density e2.

Let us attribute the parameters Rx, Mx, px to the core radius, mass and surface pressure. Then one integrates the one-fluid TOV equations, with energy density ex, outward up to vanishing pressure with Mx = (4π/3)(e1 + e2)Rx3 and

(2) Rx = R2, px = δ/(e2 + pc2), ex = e1,

(3) Rx = R1, px = −δ/(e1 + pc1), ex = e2.

In such a way the use of elliptic integrals mentioned below is avoided.

The reasoning behind that procedure is as follows. For a viable star model, one has to continue beyond the first pressure zero (e.g. of fluid 2, as in the above cases (i) and (ii) in Figure 1 or (i) in Figure 2) with respective zero energy density and pressure (e.g. e2 = 0, p2 = 0). Beyond the second pressure zero, also the other energy density and pressure are set to be zero, i.e. the exterior Schwarzschild solution is to be applied. In other words, the core consisting of fluids 1 and 2, both with p1,2 ≥ 0, is surrounded by a one-fluid shell, the envelope or corona, which may be either fluid 1 or fluid 2, depending on the chosen parameters, as exemplified in Figure 1. The special case p1(r) = p2(r) (this is case (iv) in Figure 1) does not require a corona. Instead, outside of p1(R) = p2(R) = 0, the outer Schwarzschild solution with M = m1(R) + m2(R) and e1,2 = 0, p1,2 = 0 applies at r > R as emphasized in Section 4.

Let us now assume, for definiteness, that p2(R2) = 0 and p1(R1) = 0 with R2 < R1, as in cases (i) and (ii) in Figure 1. The opposite case (iii) corresponds to the black oyster [15], to be obtained by swapping the labels 1 ↔ 2. We denote m1(R2) = M1 and p1(R2) = px, and analogously m2(R2) = M2. The one-fluid corona at r > R2 is determined by the TOV equations:

m′1 = 4πr2e1,

p′1 = −(e1 + p1)(M2 + (4π/3)e1r3 + 4πr3p1) / (r(r − 2M2 − (8π/3)e1r3)).

The terms with M2 destroy the integration of Equation (20) by separation of variables. Instead, it becomes a Bernoulli differential equation with the solution:

e1 + p1(r) = (e1 + px)ω / √(1 − 2m/r) (1 − (e1 + px)ωI)

with ω ≡ √(1 − (8π/3)e1R22), m(r) ≡ M2 + (4π/3)e1r3, and the integral I ≡ 4π ∫ dr r / (1 − 2m/r)(3/2) from R2 to r.

This integral can be expressed as a sophisticated combination of the elliptic integrals of the first and second kind. If one assumes M2 to be sufficiently small, a perturbative expansion in M2 leads to the approximation I ≈ I0 + M2I1 with I0 = (1/(2e1))(1 − (8π/3)e1r2)(-1/2), I1 = (8πr)/(2(1 − (8π/3)e1r2)(3/2)) − (3/2)(8π/3)e1r2/(1 − (8π/3)e1r2)(3/2). The integration constants are fixed by the requirement p1(R2) = px. The total radius Rtot is to be found numerically by p1(Rtot) = 0.

Figure 3 exhibits a few examples of the pressure profile p1(x) with dimensionless measure x̂ ≡ (8π/3)e1r2 of the radial coordinate r. Various parameter combinations of e1,2 and pc1,c2 are chosen. Those sections displayed by solid curves have p2 ≥ 0. In cases where p2 drops below zero at p1 > 0, e2 is put to zero in the TOV equations, and the profile of p1 is continued by the dashed sections. For comparison, the Schwarzschild solution is exhibited. The cyan and red curves in the left panel have δ 0 all the way of p1 dropping from pc1 to 0. The fluid 2 leaks out fluid 1 in such a case. This is the black oyster case. The radius R2 of the oyster is to be determined by p2(R2) = 0. Changing e2 from 0.4 to 0.6 changes the sign of δ. The zero of the pressure p2 is reached at p1 > 0 (displayed as solid section). At this zero one has to put e2 = 0 and continue the profile p1 up to its zero (dashed section). In such case, the DM admixed core is inside the whole configuration.

The boundary case δ = 0 is achieved by e2 = 0.5, where p1 = 0 and p2 are reached simultaneously (not displayed). This is why we selected the nearby values e2 = 0.4 and 0.6. The right panel keeps pc2 = pc1/2, but changes to e1 = 0.1. The value e2 = 0.01 again corresponds to the black oyster, i.e. p1 dropped to zero, while p2 > 0 leaks out far away. Increasing values of e2 make the one-fluid corona larger, contrary to expectation that a larger energy density makes the pressure gradient steeper.

Figure 4 exhibits the pressure profile of p2 as a function of x̂ ≡ (8π/3)e1r2 analog to Figure 3. Note that, for the sake of comparison, the abscissa measure x̂ is here used too. The side conditions p1 ≥ 0 and p2 ≥ 0 are imposed. Solid curves are for sections where p1 ≥ 0. The dashed curves continue p2 with one-fluid Schwarzschild solution, where e1 = 0 and p1 = 0.

The masses follow simply from M1 = (4π/3)e1R13 and M2 = (4π/3)e2R23. The values of R1 can be read off the figure, while R2 needs the information on p2(x) given in Figure 4. For the sake of comparison with Figure 3, we employ here also the radial coordinate measure x̂ = (8π/3)e1r2. In those cases where δ > 0 (left panel with e2 = 0.6 and 2, and right panel with e2 = 0.1, 0.5 and 1), the pressure p2 drops to zero within the region where p1 > 0, i.e. the corona consists of fluid 1. Conversely, for those cases where δ < 0 (left panel with e2 = 0.4 and 0.1, and right panel with e2 = 0.01), the pressure p2 drops to zero outside the fluid 1. That is, the corona consists of fluid 2, thus forming the black oyster. Particularly striking is the wide leakage of the corona like an extended halo for small values of e2. If one uses as measures of the radial coordinate x̄ ≡ (8π/3)e2r2, then the zero of p2 is seen at x̄ < 1. For parameter choices constrained to pc2 < e2, the corona is not so extremely extended.

In Section 6, the summary is given. The famous interior Schwarzschild solution of Einstein's equations for a static fluid sphere with energy density e0 and central pressure pc delivers the pressure profile p(r) = −e0(C0c − W0(r))/(3C0c − W0(r)) with C0c ≡ (e0 + pc)/(e0 + 3pc) and W0(r) ≡ (1 − (8π/3)e0r2)(1/2). Considering the energy density e0 ≥ 0, the pressure runs from pc = p(r = 0) > 0 to negative values with increasing values of the circumferential radial coordinate r. Defining the radius R0 by the pressure zero, p(R0) = 0, the Misner-Sharp mass becomes M = (4π/3)e0R03 and serves as parameter of the exterior Schwarzschild solution. The present paper extends this one-fluid Schwarschild solution to the two-fluid Schwarzschild solution. That is, two fluids 1 and 2 (with constant energy densities e1,2 and central pressures pc1,c2) interact solely via their common gravitational field. We present the profiles p1,2(r) and their parametric dependencies. In general, we face a four-dimensional parameter space.

Our focus is on the two-fluid core. In special cases, the pressure profiles of both components agree and the radius is determined by their agreeing pressure zeroes. This is the special case of a quadric cone in the parameter space, where both fluids occupy the same region, surrounded by vacuum. In general, however, the pressure zeroes of both fluids differ. Then, one defines the core radius by that pressure zero with smaller radius, i.e. the first zero. One has to continue the solution outside the core by a one-fluid zone (termed corona) up to the second pressure zero. Depending on the parameters, the corona may consist of fluid 1 or fluid 2. In the outer space region the famous exterior vacuum Schwarzschild solution applies.

It is surprising to see many studies of the impact of dark matter on properties of compact (neutron) stars. What was missed in the literature is a comprehensive analytic reference dealing with the two-fluid Schwarzschild solution. We fill this gap and emphasize the possible leakage of one fluid far beyond the two-fluid core. Low energy density and high central pressure of dark matter can facilitate such an object (black oyster).

Subsequent investigations should envisage the tidal deformability, stability properties and the formal aspects of the sign of the energie densities e1,2 and the central pressures Pc1,c2. Relations to metric quantities need to be considered. Interesting are additionally the impact of energy conditions and the occurrence conditions of photon spheres. Also multi-layered structures should be considered to bridge over to realistic equations of state.

Improvements to AI Systems Based on This Paper:

1. Analytic Two-Fluid Stellar Structure Solver

• Implement the closed-form two-fluid TOV solution (Equations 2–3) as a fast, exact module for neutron star models with admixed dark matter, replacing purely numerical integrators.

• The AI can instantly compute pressure profiles p1(r), p2(r), core radii, and corona boundaries for arbitrary e1, e2, pc1, pc2 without iterative shooting methods.

2. Phase-Space Classifier for Core–Corona Configurations

• Use the sign of δ = −e1pc2 + e2pc1 to automatically classify solutions into: (a) both fluids co-located (δ=0), (b) corona of fluid 1 (δ>0), (c) corona of fluid 2 / black oyster (δ<0).

• The AI can predict the topological structure of the star (core, corona, vacuum) from parameter inputs alone, enabling rapid parameter-space surveys.

3. Compactness Bound Validator

• Encode the proof that the Buchdahl limit (C ≤ 8/9) is never exceeded even with two energy densities, using the derived compactness formula C = 4X/(1+3X) with X = (pc1+pc2)/(e1+e2).

• The AI can automatically flag any proposed two-fluid configuration that would violate this bound, serving as a safety check for astrophysical models.

4. Corona Continuation Algorithm

• Implement the Bernoulli-equation solution (Equation 20) for the one-fluid corona, including the elliptic-integral expression and the perturbative M2 expansion.

• The AI can generate complete pressure profiles from core to surface without numerical ODE solvers, with controllable accuracy via the perturbation order.

5. Parameter-Space Explorer for Dark Matter Effects

• Use the scaling relations (pi/e+ as functions of e+r2, p+c/e+, p−c/e+, e−/e+) to reduce dimensionality from 4 to 3.

• The AI can rapidly map how dark matter energy density and central pressure alter pressure gradients, zero-crossing radii, and corona extent, identifying regimes of extreme leakage (e.g., low e2, high pc2 → extended halo).

6. Automated Stability Pre-Screening

• Detect configurations where p1(W=0) or p2(W=0) remain positive (indicating no pressure zero before the horizon), signaling potential instability.

• The AI can pre-filter parameter sets to exclude unphysical or unstable models before detailed dynamical analysis.

7. Black Oyster Detection Module

• Identify parameter regimes (δ < 0, small e2, high pc2) that produce a black oyster — a dark matter halo leaking far beyond the ordinary matter core.

• The AI can generate synthetic observational signatures (mass–radius relations, tidal deformability) for such exotic objects, aiding in dark matter searches.

8. Multi-Layer Extension Framework

• Use the corona-matching procedure (core + one-fluid shell + vacuum) as a building block for multi-layered stars with piecewise constant densities.

• The AI can iteratively stack solutions to approximate realistic equations of state, bridging the gap between this analytic model and numerical stellar structure codes.

9. Cross-Checking Tool for Numerical Relativity Codes

• Provide exact analytic solutions as benchmark tests for general-relativistic hydrodynamics codes.

• The AI can automatically generate test cases with known pressure profiles, masses, and radii to validate numerical implementations of two-fluid TOV equations.

10. Educational and Explanatory Interface

• Convert the paper’s derivations into an interactive AI assistant that explains the physics, derives the key equations step-by-step, and visualizes pressure profiles for user-chosen parameters.

• The AI can answer what-if questions (e.g., How does increasing dark matter pressure affect the core radius?) using the analytic formulas, providing immediate intuition.

We consider the interior two-fluid Schwarzschild solution. That is a model of compact (neutron) stars with admixed dark matter. The two constant energy densities generalize the interior Schwarschild solution to the case of two fluids (representing standard model matter and dark matter, for instance), interacting mutually solely by the common gravitational field. In general, the two-fluid core is surrounded by a one-fluid corona (envelope). Only for special parameters (the model depends on energy densities e 1,e 2 and central pressures p c1,p c2 of the fluids), both fluids occupy the same region. Despite the sum of two energy densities in that region, the compactness bound of 8/9 is not exceeded. A conceivable object is one with a stark leakage of one fluid (e.g. dark matter) beyond the two-fluid core.