{'fraction_non_alphanumeric': 0.05129379436684223, 'fraction_numerical': 0.02064727883367682, 'mean_word_length': 4.207273449920509, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Radial oscillations and stability of multiple-fluid compact stars

I derive a system of pulsation equations for compact stars made up of an arbitrary number of perfect fluids that can be used to study radial oscillations and stability with respect to small perturbations. I assume spherical symmetry and that the only inter-fluid interactions are gravitational. My derivation is in line with Chandrasekhar's original derivation for the pulsation equation of a single-fluid compact star and keeps the contributions from the individual fluids manifest. I illustrate solutions to the system of pulsations equations with one-, two-, and three-fluid examples.


I. INTRODUCTION
The perfect fluid model has been very successful in modeling the matter content of compact stars [1]. This is particularly the case with fermionic matter, because the perfect fluid model is able to capture properties of the quantization of fermions, such as Pauli exclusion.
Static compact star solutions using perfect fluids are found by solving the time-independent Einstein field equations with a perfect fluid energymomentum tensor. Recently there has been increased interest in compact stars made with multiple fluids. This has been motivated by the possibility that neutron stars might capture dark matter [2][3][4]. If dark matter cannot self-annihilate (as is the case, for example, with asymmetric dark matter [5,6]), it will accumulate in the neutron star. In such a scenario, dark matter can be modeled as an additional fluid in the star. There have been a number of recent studies of the properties of static two-fluid compact stars [7][8][9][10][11][12][13][14][15][16][17][18].
After a compact star solution is found, one can study the effect of small perturbations to the solution. This may be done by solving a pulsation equation, which is a differential equation describing the time-evolution of the perturbation. For spherically symmetric static solutions and perturbations, the solution to the pulsation equation gives the squared radial oscillation frequency. In addition to describing radial oscillations, the squared radial oscillation frequency tells us about the stability of the star with respect to small perturbations. If the squared oscillation frequency is negative, the compact star solution is unstable.
Chandrasekhar was the first to derive a pulsation equation for a compact star made of a single perfect fluid [19]. His pulsation equation was subsequently rewritten in various ways [20][21][22][23], in some cases to facilitate numerical solutions. More recently, a pulsation equation was derived for twofluid systems [24,25] and has been used to study the stability of two-fluid stars [9,10,12] (an alternative approach for determining the stability of two-component stars, which has been applied to bosonfermion stars [26,27], was studied in [28,29]). In this construction, there is a single pulsation equation. It is possible to solve the single pulsation equation in a two-fluid system, but without additional assumptions it can be difficult to do so.
One such assumption is to disallow nongravitational inter-fluid interactions. Indeed, this assumption is made in many of the studies of two-fluid compact stars [7][8][9][10][11][12][13][14][15][16]18] (a notable exception is [17]). This assumption allows the equations of motion to separate and leads to a system of coupled pulsation equations, with the same number of pulsation equations as fluids. This increase in the number of pulsation equations makes them easier to solve.
In this work, I present this system of pulsation equations. I shall assume spherical symmetry and, as mentioned, allow only gravitational interfluid interactions (I make no assumptions about selfinteractions). In doing this, I give a very different derivation than in [24,25], a derivation made possible by my assumption of there only being gravitational inter-fluid interactions and a derivation in line with Chandrasekhar's original derivation for a single fluid. The system of equations I arrive at keeps the contributions from the individual fluids manifest and allows them to be easily solved for, if desired. I illustrate the use of this system of pulsation equations with one-, two-, and three-fluid examples.
In the next section I derive the system of pulsation equations. In Sec. III, I present one-, two-, and three-fluid examples. I conclude in Sec. IV.

II. EQUATIONS
In this section I derive a number of equations. This includes the TolmanOppenheimerVolkoff (TOV) equations adapted to a multi-fluid system. Solutions to the TOV equations describe static compact stars and I will often refer to such solutions as equilibrium solutions. I also derive the system of pulsation equations, which describe radial oscillations about the equilibrium solutions. I consider only spherically symmetric spacetimes. In deriving equations-in particular, the pulsation equations-it is convenient to use a spherically symmetric metric of the form where dΩ 2 = dθ 2 + sin 2 θ dφ 2 . I note that later I will use a different form for the metric that is better suited for numerically solving the TOV and pulsation equations. I assume there are an arbitrary number of perfect fluids, whose individual energy-momentum tensors separate: where x labels the fluid, is the standard energy-momentum tensor for a perfect fluid, u µ x is the four-velocity of the fluid, and x and p x can be thought of as the fluid's energy density and pressure. T µν tot in Eq. (2) is the (total) energymomentum tensor of the system and is what goes on the right hand side of the Einstein field equations: In addition to its energy-momentum tensor, the matter sector is defined by equations of state, which I also assume separate: i.e. p x depends only on its associated energy density x and not on any y =x . The assumption that T µν x and p x only depend on the variables of fluid x means that there are only gravitational inter-fluid interactions. This assumption greatly facilitates solving both the equilibrium and pulsation equations, in part because it means that the individual T µν x are conserved, in addition to the requisite ∇ µ T µν tot = 0. In a spherically symmetric system, fluids cannot flow in the θ-or φ-directions and thus u θ = u φ = 0. Following Chandrasekhar [19], I define v x ≡ e ν/2 u r x . Using that g µν u µ x u ν x = −1, I have u t x = exp(−ν/2) 1 + exp(λ − ν)v 2 x . It is now straightforward to write the components of T µν x in Eq. (3) in terms of v x . Doing so is not particularly illuminating, nor will I be using T µν x in such a form. Instead, I write the system variables as perturbations about their equilibrium values: The equilibrium values are denoted by a subscripted 0 and describe equilibrium, or static, solutions and hence are time-independent. I note that the equilibrium value of v x vanishes and v x is itself at the order of a perturbation. My interest is in using the perturbations to describe radial oscillations about the equilibrium solutions. I therefore consider equations only to first order in the perturbations. To this order, the nonzero components of T µν Summing these components over x gives (T tot ) µ ν . Setting the perturbations to zero, the equilibrium energy-momentum tensor is Summing these components over x gives (T tot 0 ) µ ν . For the equilibrium solutions, we see that x0 and p x0 are fluid x's contributions to the (total) energy density and pressure of the system: The Einstein field equations in (4) lead to a number of equations that are written in terms of the components of the energy-momentum tensor and determine the metric functions ν and λ. Three such equations are [19,30] where a prime denotes an r-derivative and a dot denotes a t-derivative. The equilibrium version of these equations are straightforward to write down: where the equilibrium version of Eq. (19) vanishes identically. The perturbed version of Eqs. (17)- (19), after using Eqs. (20) and (21) to cancel the equilibrium parts, are There are two more equations for the metric functions that I will need. The first is obtained by combining Eqs. (20) and (22): The second is less commonly used than Eqs. (17)- (19), but like them follows from the Einstein field equations. I only need its equilibrium version [19,30]: (26) Lastly, I need the equations of motion. I will obtain these from ∇ µ T µν x = 0, i.e. from conservation of the individual T µν x . Using the metric in Eq. (1) to evaluate the divergence, I find where the first equation is for ν = t and the second is for ν = r. Moving to equilibrium variables, which requires treating everything as time-independent, Eq. (27) vanishes identically and Eq. (28) becomes The perturbed versions of Eqs. (27) and (28) are

A. Equations for equilibrium solutions
In this subsection I put together the equations whose solutions describe equilibrium, or static, compact stars. In doing this, I will write the equations using a different parametrization of the spherically symmetric metric than written in Eq. (1). I will use instead where which is better suited for numerical solutions. Dropping the time-dependence, so that the new equilibrium metric variables are σ 0 (r) and m tot 0 (r), the equilibrium equations in (20), (21), and (29) become where Equations (35)-(37) are the Tolman-Oppenheimer-Volkoff (TOV) equations adapted to a multi-fluid system [31] that has only gravitational inter-fluid interactions. The individual m x0 (r) can be interpreted as giving the total mass inside a radius r for fluid x.
To solve the TOV equations, and thus find static compact star solutions, one specifies, say, the central values of the energy densities, x0 (0), from which the central values of the pressures, p x0 (0), are obtained from the equations of state, p x0 ( x0 ). Equations (35)-(37) can then be integrated outward from some small r once appropriate inner and outer boundary conditions are determined. Inner boundary conditions can be determined by subbing into Eqs. (35)-(37) power law expansions of the variables, which gives σ 0 (0) = σ c +O(r 2 ) and m x0 (0) = O(r 3 ), where σ c is an as-yet-undetermined constant. For outer boundary conditions, let R x be the smallest value of r for which p x0 (r) = 0, which I will take to define the edge of fluid x. For r > R x , I set p x0 (r) = 0. The edge of the star, R, is taken as the largest R x . At the edge of the star there is no longer any matter and the spacetime must match up to a Schwarzschild spacetime, and thus σ 0 (R) = 1.
I have so far that the inner boundary condition σ 0 (0) = σ c is as-yet-unknown, while the outer boundary condition σ(R) = 1 is known. It is convenient to flip these by defininĝ In usingσ 0 , the only equation that changes is Eq. (35), which becomeŝ The inner boundary conditions are nowσ 0 (0) = 1 and m x0 (0) = 0, and the TOV equations in (36), (37), and (40) can be integrated outward from some small r to the edge of the star at r = R, where the total mass of the star is given by and σ c is given by σ c = 1/σ 0 (R). In this way equilibrium/static compact star solutions can be found. I shall do this in Sec. III.

B. Pulsation equations
I now take up one of the main parts of this paper, which is the derivation of a system of pulsation equations whose solution describes radial oscillations for a system of perfect fluids with only gravitational inter-fluid interactions. Further, the derivation is in line with Chandrasekhar's original derivation of a pulsation equation for a single fluid. One convenience of the resulting system of pulsation equations is that the contributions from individual fluids are manifest and easily solved for, if desired. I begin first by putting together a number of equations, which I will then use to derive the pulsation equations. The first step is to definė Plugging this into Eq. (24) and integrating gives which can be combined with the metric equilibrium equations (20) and (21) and Eq. (25) to obtain The equation of motion in (30), when written in terms of ξ x using Eq. (42), can immediately be integrated to give Replacing δλ with Eq. (43) and then combining the result with Eqs. (20) and (21), I find The equations of state, p x = p x ( x ), after being perturbed and canceling their equilibrium terms, becomes δp x = (∂p x0 /∂ x0 )δ x . Plugging in Eq. (46), this can be written as where is the adiabatic index for fluid x.
Having established the above results, I will now outline the derivation of the pulsation equations, which follow from the equation of motion in (31). Introducing a harmonic time-dependence for all fields: δν(t, r) = δν(r)e iωt , δλ(t, r) = δλ(r)e iωt , δ x (t, r) = δ x (r)e iωt , and δp x (t, r) = δp x (r)e iωt , where ω is the radial oscillation frequency, Eq. (31) becomes where I defined Consider first the first term in Eq. (50). Plugging in δp x from Eq. (47) and then using Eqs. (26) and (29), I find For the second term in Eq. (50), I have, after using Eqs. (20), (21), (29), (44), and (47), Finally, using Eqs. (29), (46), and (53) the second and third terms in Eq. (50) are Equations (52) and (54) can be plugged into Eq. (50). The point is to construct an equation in which the only perturbation present is ζ x and its derivatives. In putting everything together, I will move to the metric functionsσ 0 (r), m x0 (r), and N 0 (r) that were introduced in Sec. II A and which are more convenient for numerical solutions. The system of pulsation equations is wherê Note that all functions with a hat are functions that have been scaled by powers of σ c = σ 0 (0), as explained in Sec. II A. Importantly, we see that the only perturbation present in Eqs. (55) and (56) isζ x and its derivatives. Before discussing boundary conditions, a couple comments are in order. The right hand side of Eq. (55) vanishes for a single fluid, in which case the left hand side of Eq. (55) is equivalent to Chandrasekhar's pulsation equation [19]. Though equivalent, the left hand side of Eq. (55) is not written in an identical form to Chandrasekhar's because, in the presence of multiple fluids, terms cannot cancel and combine in the same way. When solving Eq. (55) I have found it best to do so in conjunction with solving the equilibrium TOV equations in (36), (37), and (40). Specifically, I have found that solving them all simultaneously is faster than first solving the TOV equations and then interpolating between points in the equilibrium solution for use in solving the pulsation equations.
To actually solve the pulsation equation I define η x ≡ Π xζ x and solve the system of first order differential equations made up of Eq. (55), but written in terms of η x , andζ Doing this requires inner and outer boundary conditions. Inner boundary conditions are found by plugging power series expansions of the variables into the equations, which reproduce the inner boundary conditions for the equilibrium variables given in Sec. II A and ζ x = O(r 3 ) and η x = η x c + O(r 2 ), where the η x c are as-yet-unknown constants. A look at the equations shows that the set ofζ x and η x that solve them can all be divided by the same constant and still be solutions. I can use this to set i.e. η x=1 c = 1. I still need to determine the η x =1 c . My method for doing this is explained in the next paragraph. The outer boundary conditions for the equilibrium variable were discussed in Sec. II A. The outer boundary conditions for the perturbations are η x (R x ) = 0. Intuitively this makes sense, as outside the fluid the perturbation ζ x should stop changing. It can be derived by requiring the Lagrangian pressure for fluid x to be zero at r = R x [20][21][22].
The as-yet-undetermined constants are η x =1 c and the squared radial oscillation frequency ω 2 . I determine them using the shooting method. That is, I choose values for ω 2 and the η x =1 c . The remaining inner boundary conditions areσ 0 (0) = η 1 c = 1 and m x0 (0) = ζ x (0) = 0. I can then integrate the TOV equations (36), (37), and (40) and the pulsation equations (55) and (57) outward from some small r. I vary ω 2 and the η x =1 c until the smallest value of r for which η x (r) = 0 is exactly r = R x , and hence the outer boundary conditions are satisfied. I will have then found a solution. If ω 2 < 0, the solution is unstable to small radial perturbations.

III. EXAMPLES
In this section I present one-, two-, and three-fluid examples in which I solve for the squared radial oscillation frequency, ω 2 , and determine stability. In all cases, for simplicity, I take each fluid to be a free Fermi gas. The well-known equation of state and number density for a free Fermi gas are [1,32] where m xf is the fermion mass and k xF is the Fermi momentum for fluid x. It is convenient to define so that the equation of state and number density can be written in the parametric form The equilibrium versions of these equations are given by simply replacing x , p x , n x and t x with x0 , p x0 , n x0 , and t x0 .
To find static compact star solutions, I can numerically solve the TOV equations in (36), (37), and (40), as explained in Sec. II A. However, instead of evolving the pressure p x0 using Eq. (37), it is much simpler to evolve t x0 using where p x0 is given by Eq. (37) and which follows from the equilibrium version of Eq. (65). The general procedure for finding static compact star solutions and their squared radial oscillation frequencies is as described in Secs. II A and II B. The only changes specific to the free Fermi gas that I make is in the use of the variable t x0 : To find static compact star solutions, I use the central values t x0 (0), which are related to the central values for the energy density, pressure, and number density through Eqs. (64)-(66). I integrate outward from some small r using the TOV equations in (36), (40), and (67). I take the edge of fluid x, whose location I label as r = R x , to be the smallest value of r such that t x0 (r) = 0. Defining the edge of the fluid this way is equivalent to defining it as p x0 (r) = 0. For r > R x , I set t x0 = 0. The edge of the star is given by the largest R x .
In presenting results in the following subsections, I will make use of the dimensionless variables where m P = 1/ √ G is the Planck mass. I have 1f /m P )R x ,m tot = xm x , and M tot = x M x . Regardless of the number of fluids, I choose to scale variables using m 1f , the first fluid's fermion mass.

A. One fluid
With only a single fluid I can drop the subscripted xs on all variables. Figure 1(a) shows the well-known (total) mass versus radius curve for the static/equilibrium solutions for a free Fermi gas [33]. Each point on the curve corresponds to a different choice for the central value t 0 (0) (or, equivalently, different choices for the central values 0 (0), p 0 (0), or n 0 (0)). An important question is whether each point on the curve corresponds to an equilibrium solution that is stable or unstable with respect to small radial perturbations.
Stability may be determined by computing the squared radial oscillation frequency, ω 2 , for the particular equilibrium solution, which is done by solving the pulsation equation. An example of a solution to the pulsation equation is shown in Fig. 1(b). The top solid curve plots t 0 (r) for the equilibrium solution indicated by the dot in Fig. 1(a). The bottom  Fig.  1(a). The squared oscillation frequency is seen to transition from positive to negative, and hence the equilibrium solutions transition from stable to unstable, at the equilibrium solution with the largest mass. The curves in Fig. 1(b) correspond to the solutions located at the black dots.
dashed curve plots η(r). As explained in Sec. II B, η(0) = 1 and a solution to the pulsation equation is found when the smallest value ofr for which η(r) = 0 occurs precisely at the edge of the star atr = R. The advantage in plotting t 0 in Fig. 1(b) instead of, say, p 0 , is that t 0 heads toward t 0 = 0 relatively abruptly at the edge of the star, making the edge of the star easily discernible in plots such as Fig. 1(b).
For a single fluid, it is well-known that the pulsation equation can be written in Sturm-Liouville form [20,23]. This tells us that a single equilibrium solution has an infinite number of radial oscillation frequencies, corresponding to an infinite number of solutions to the pulsation equation [32]. When the solution to the pulsation equation is found by having r = R be the smallest r such that η(r) = 0, the oscillation frequency found is the fundamental, which is the smallest of the oscillation frequencies.
As the oscillation frequency is increased, η(r) develops nodes, i.e. additional locations where η(r) = 0. A solution is still found when η(R) = 0, so that the outer boundary condition is satisfied. In this paper I consider only fundamental frequencies, and thus the smallest frequencies. Figure 2(a) plots the squared radial oscillation frequencies as a function of the central energy densitȳ 0 (0) when they are positive, and hence for the stable equilibrium solutions. Figure 2(b) plots the (total) mass of the equilibrium solutions also as a function of the central energy density¯ 0 (0) and is an alternative to Fig. 1(a) as a way of presenting the space of equilibrium solutions. It is well-known that for a single fluid, the transition from stable equilibrium solutions to unstable ones occurs at the equilibrium solution with the largest mass [1,32]. This is exactly what is seen in Fig. 2, where the squared frequencies transition from positive to negative at precisely the location of the largest mass equilibrium solution. The dots in Fig. 2 indicate the solutions shown in Fig. 1(b), just as the dot in Fig. 1(a) does.

B. Two fluids
For two fluids, the first set of results I show is for when the fermion mass of the second fluid is seven times larger than the fermion mass of the first fluid, m 2f = 7. I fix the central energy density of the second fluid,¯ 2,0 (0), and scan through values of the central energy density of the first fluid,¯ 1,0 (0). The (total) mass versus radius curves for the resulting equilibrium/static solutions are shown in Fig. 3(a). Each curve is for a different fixed value of¯ 2,0 (0). I have also included in this plot, for comparison, the single fluid curve from Fig. 1(a) as the dashed line.
To get a sense of the equilibrium solutions, in Figs. 3(b)-(e) I show the maximum (total) mass equilibrium solutions from each of the two-fluid curves in Fig. 3(a). Moving from Fig. 3(b) to Fig. 3(e) means increasing¯ 2,0 (0). The solid curves in Figs. 3(b)-(e) display¯ 1,0 (r) and¯ 2,0 (r) and the dashed curves displaym 1,0 (r) andm 2,0 (r), which gives the mass of fluid x inside a radiusr. In Figs. 3(b)-(e) we can can see that fluid 2 dominates the core of the star, while fluid 1 dictates the total mass (M tot =m 1,0 (R) +m 2,0 (R)) and radius of the star. Figure 4(a) shows the squared radial oscillation frequencies for the equilibrium solutions shown in Fig. 4(b) (which are the same equilibrium solutions shown in Fig. 3(a)). In both figures, the curves are given as functions of the first fluid's central energy density,¯ 1,0 (0). I have only plotted the oscillation frequencies and equilibrium solutions for when fluid 1 affects M tot . This is the reason why the oscillation frequency curves in Fig. 4(a) begin at the same value of¯ 1,0 (0) as the equilibrium solution curves do in Fig. 4(b). Just as in the single-fluid case, we see that the squared frequencies transition from positive to negative, and hence the equilibrium solutions transition from stable to unstable, at the equilibrium solution with the largest mass.
I now present a different set of results. This time I fix¯ 2,0 (0)/¯ 1,0 (0) = 10 and find solutions  Fig. 3(a). Just as in the single fluid case (see Fig. 2), we see here for two fluids that the squared oscillation frequencies transition from positive to negative, and hence the equilibrium solutions transition from stable to unstable, at the equilibrium solution with the largest M tot. form 2f = 3, 4, 5, and 6. Figure 5(a) displays the (total) mass versus radius curves and Figs. 5(b)-(e) displays the maximum (total) mass equilibrium solutions, analogously to Fig. 3. Figure 6 displays the squared radial oscillation frequencies (in the top plots) for the same equilibrium solutions (shown in the bottom plots) as displayed in 5(a). In the bottom plots, the vertical lines mark the maximum mass solutions and the crosses mark the transition from stable to unstable. One can see that the transition from stable to unstable does not occur at the equilibrium solution with the largest mass. This is not unexpected, and, in general, occurs in multiple-component stars (see, for example, [28,29]).

C. Three fluids
With three fluids there are three parameters to be determined using the shooting method:ω 2 , η 2 (0), and η 3 (0). It is significantly more time consuming to determine three parameters than it is to determine two parameters and I present results for only a single equilibrium solution in Fig. 7. Figure 7 is analogous to Fig. 1(b), with the solid curves plotting the equilibrium variables t x0 (r) and the dashed curves plotting the pulsation variables η x (r).

IV. CONCLUSION
In this work I derived a system of pulsation equations for spherically symmetric, multiple-fluid compact stars, under the assumption that the only interfluid interactions are gravitational. The solution to the pulsation equations is the squared radial oscillation frequency for the star. I solved the system of pulsation equations in one-, two-, and three-fluid examples and used the squared radial oscillation frequency to determine stability with respect to small   Fig. 1(b). The fermion masses for the fluids arem 2f = 1.5 andm 3f = 2. The solid curves plot the equilibrium variables tx0(r), with (from top to bottom) t1,0(0) = 1.5 (black), t2,0(0) = 1.2 (blue), and t3,0(0) = 0.8 (purple). The dashed curves plot the pulsation variables ηx(r). This solution has M tot = 0.1072, R = 5.0756, andω 2 = 0.003091 and is stable.
perturbations. As far as I am aware, this is the first time that a pulsation equation has been solved for a three-fluid compact star.