Constraint on hybrid stars with gravitational wave events

Motivated by the recent discoveries of compact objects from LIGO/Virgo observations, we study the possibility of identifying some of these objects as compact stars made of dark matter called dark stars, or the mix of dark and nuclear matters called hybrid stars. In particular, in GW190814, a new compact object with 2.6 $M_{\odot}$ is reported. This could be the lightest black hole, the heaviest neutron star, and a dark or hybrid star. In this work, we extend the discussion on the interpretations of the recent LIGO/Virgo events as hybrid stars made of various self-interacting dark matter (SIDM) in the isotropic limit. We pay particular attention to the saddle instability of the hybrid stars which will constrain the possible SIDM models.


Introduction
The LIGO/Virgo events GW170817 [1,2] and GW190425 [3] are in general considered as binaries of neutron stars (BNSs), as well as the possibility of being binaries of hybrid stars.
However, in [4], a new binary coalescence event GW190814 is reported, while one is a black hole with 23 M , and the other is a compact object with only 2.6 M (with a range between 2.50 and 2.67M ). Then, the question is what is the identity of the secondary compact object? There are several possibilities: It could be the lightest black hole, or the heaviest neutron star, or a dark star. In addition, there are also possibilities for this compact object to be a hybrid star made of neutron and dark matter [5]. Below, we will extend a few more discussions on the above possibilities, and then focus on elaborating more on the scenarios of dark and hybrid stars in the rest of the paper.
The scenario of black hole-The black holes observed by LIGO/Virgo are all known to have masses more than 5 M . However, the remnant of GW170817, which is identified to be the BNS by parameter estimation (PE), is estimated to have a mass of about 2.7 M , and is believed to be a black hole [2]. Thus, the black hole of mass less than 5 M can be formed via the coalescence of BNS. Since no information about tidal Love number is given in the released PE of GW190814, and no electromagnetic follow-up is reported, we cannot exclude the possibility that this 2.6 M object is the lightest black hole, which could likely be either formed astronomically via the coalescence of BNS, or formed in the early universe as a primordial black hole [6]. Alternatively, the light black hole may come from core-collapse supernova of a massive normal star [7].
The scenario of neutron star-Before this GW190814 event, masses of neutron stars are generally considered to be less than 2.3 M [8,9]. In fact, within our galaxy, the most massive known pulsar has 2.14 +0. 10 −0:09 M at 68.3% credible interval measured by the Shapiro time-delay of its white dwarf companion [10]. Most known equations of state (EoSs) for the dense nuclear matter used in the astronomical search of neutron stars cannot sustain such high mass as 2.3 M , although this can also be seen as a result of tuning the parameters of the EoSs to not violate the upper mass limit of the astronomical observations. On the other hand, rapid spinning of neutron star [11,12] can increase the maximal mass up to ∼20%, which will lead to ∼2.7 M . This case cannot be totally ruled out since there is no constraint on the spin of the secondary object in GW190814. However, such high spin NSBH system may hardly exist because of the expected subsequent collapse into black hole due to the loss of spin via gravitational wave radiation.
The scenario of dark star and hybrid star-Currently, dark stars have not yet been confirmed in observations so that it remains to pin down the accessible model space of dark matter and the associated EoSs by the future observations in the coming era of gravitational wave astronomy. In light of the varieties of dark matter models, it is easy for the resultant dark stars to cover a wide mass range, say 1 to 5 M , or even one order higher with proper parameters. This makes the dark stars or the hybrid stars of dark and nuclear matters to be highly possible candidates to explain the companion star of GW190814. In the remaining of this paper, we will explore this possibility by studying the mass-radius relations for the dark and hybrid stars based on various dark matter models inspired by the particle physics.
Specifically, we consider the massive bosonic field φ with the following self-interactions, φ 4 , φ 6 , their linear combinations and φ 10 .
One way to characterize a compact star is through its Tidal Love Number (TLN) 1 [13,14], which shows the star's tendency to deform under an external quadrupolar tidal field. It is known that the TLN of black holes in Einstein gravity is vanishing, and the overall TLN effect for a compact binary coalescence event is a weighted average of the individual TLNs. Thus, for GW190814, the overall tidal effect is insignificant in the resultant gravitational waveform due to the high mass ratio between black hole and the companion star. For example, even if the 2.6 M object has a large TLN such as 30,000, when combining with zero TLN of the companion black hole of 23 M , the overall TLNΛ 2 contributed to the gravitational waveform is about 43, which can hardly be observed. Since there is no available information on TLN from this event anyway, we need more future events with much smaller masses where the TLN judgement can be applied. For the above reason, we will not present the mass-TLN relation in this work.
The rest of the paper is organized as follows. In the next section, we will sketch how we extract the EoSs of the bosonic dark matter models in the isotropic limit. In Section 3, we discuss the stability of the dark and hybrid stars based on the famous Bardeen-Thorne-Meltzer (BTM) criteria, especially with the emphasis on the saddle instability when the mixed phase rule does not apply. The key result is exposed in Section 4 by plotting the massradius relation for various EoSs extracted from the respective dark matter models. Based on the mass-radius relation, we discuss the relevance to and interpretation of GW190814.
Finally, we conclude our paper in Section 5.

EoS for Bosonic SIDM in the Isotropic Limit
Most of the dark matter models are motivated by particle physics, which have either weak or no interaction with the standard model particles. The former is called the WIMP, namely weak interacting massive particles, and the latter will also include the self-interactions to explain the core profile of dark halos well so that it is usually called SIDM, self-interacting dark matter. In this paper, we will mostly focus on SIDM. The simplest model of SIDM is the massive φ 4 bosonic field theory considered in [15]. Naively, one should solve the combined scalar-tensor field equations for the possible compact dark star configurations. 1 The TLN denoted by Λ is defined by Q ab = −M 5 Λ E ab , where M is the mass of the star, Q ab is the induced quadrupole moment, and E ab is the external gravitational tidal field strength.
where Mi and Λi with i = 1 or 2 is the mass and TLN for the i-th component compact object.
However, if we assume the scalar profile inside the star varies very slowly, we can neglect the spatial profile and obtain an isotropic dark star configurations. In the low energy regime, these kinds of isotropic configurations are more favored energetically than the nonisotropic ones. Thus, for simplicity, we will only consider such kind of dark and hybrid stars. One additional advantage for solving the isotropic dark stars is we can first extract the EoS in the isotropic limit, then solve the Tolman-Oppemheimer-Volkoff (TOV) equations for the mass-radius relation. This is numerically far easier than solving the scalar-tensor field equations by the shooting method 3 .
In the following, we sketch how to extract the EoS for the generic bosonic SIDM models by generalizing the procedure given in [15]. The Lagrangian of the bosonic SIDM model considered in this work is from which we can obtain the field equation where V := ∂ φ V .
We will consider the following metric ansatz with spherical symmetry for the dark or hybrid star This metric is sourced by the spherically symmetric scalar field configuration which should solve (2) in the space-time (3), i.e., On the other hand, the stress tensor associated with Lagrangian (1) is which satisfies the conservation law and also sources the Einstein equation with G µν the Einstein tensor and G N the Newton constant. The total configurations for a boson star specified by A(r), B(r), and Φ(r) should be determined by solving (2) and (8) together. In general, the stress tensor for the stationary configurations of (4) satisfying (5) in the space-time (3) takes the following form in the co-moving frame However, in the isotropic limit, it will further reduce to the form of a perfect fluid, i.e., p ⊥ = p. Now, we consider the following concrete example of bosonic dark matter model with the self-interactions as This can be thought as a model possessing of a UV Z n symmetry, which is however mildly broken at low energy by the mass term, that is, the Z n symmetry is approximately good at low energy if Φ 0 m. For this bosonic SIDM model, its stress tensor is given by which takes the form of (9) as expected.
To perform the isotropic limit, we first introduce the following dimensionless variables Here, M P l is the Planckian mass scale which is related to the Newton constant by G N = 1/M 2 P L . Then, the scalar field equation (5) is reduced into By further introducing a new dimensionless parameter Λ n defined by and the new scaled quantities then the energy density and pressures of of (11) in the form of (9) can be expressed as follows: and (13) is also further reduced to is large enough 4 so that the spatial kinetic term can be neglected in comparison to the other terms 5 . Thus, the resultant approximate energy density and isotropic pressure then Moreover, in this isotropic limit, both ρ and p depend only on σ * and not on its spatial derivative, and the first term of (19) is suppressed with respect to the second term 6 so that it can be solved for σ * to yield Using (22), the energy density and pressure can be expressed as a pure function of σ * , where the overall energy density scale is given by Note that (23) and (24) are already forming a parametric EoS for the dark matter model in the isotropic limit. However, we can further eliminate the parametric function σ * (r) to yield a compact form of EoS as following (in the units of c = 1 and = 1): One can then adopt this EoS to solve the TOV configurations for the dark or hybrid stars. Note that, for n = 4, the above EoS reduces to the known result given in [15], namely, For the other n's, the resultant EoSs are new and not considered before in the literature.
It is interesting to see that, as n goes higher, the EoS becomes stiffer, e.g., as n −→ ∞, the EoS goes to p = ρ, i.e., the causality limit.
Even though we have arrived at the EoS (27), it is written in the unit of ρ 0 . For the convenience of later use when solving TOV configurations, we will choose the following astrophysical units associated with the solar mass M : 4 Since the UV Zn symmetry is approximately good if Φ0 m, then the isotropic limit Λn = 1 is easier to achieve. It is noticed that a large derivative of σ * may balance the large Λn term, but here we only consider the simplest case. 5 Naively, the comparison should also take into account the spatial profile of A(x), B(x) and σ * (x). Here, we assume that Λn is large enough as argued in the previous footnote so that the variation of the spatial profiles of A, B and σ * will not affect the result. This of course can be justified after solving the TOV configurations. 6 The same assumption is adopted as in the previous footnote to justify the suppression.
Similarly, the EoS (26) can be turn into The above method of extracting the EoS in the isotropic limit from a given scalar field theory can be applied to the model with more general potentials other than (10). For example, for the model with the following potential, In this case, we have two scaling parameters Λ 4 and Λ 6 as defined in (14). However, it is easier to parameterize the EoS in the isotropic limit by the following two parameters: Λ := √ Λ 6 and β := Λ4 4Λ so that the extracted EoS takes the following form in terms of the astrophysical units of (28): where B := . Some comments about the above EoS are in order: (i) one can see B > 0 but β can be either positive or negative. When β = 0, it reduces to the case of pure φ 6 model with We should require both ρ ≥ 0 for positive energy condition and p ≥ 0 for not considering the gravastars, thus the corresponding physical range of σ should be chosen carefully. In particular, when β > 0, one should require σ ≥ √ 3β to keep p ≥ 0.
An example for this case is shown in Figure 1. On the other hand, both ρ and p are positive for all ranges of σ when β < 0.
Below, we will study the TOV configurations of dark and hybrid stars made of nuclear and dark matters based on the above EoSs of bosonic SIDM model in the isotropic limit.
In particular, we will focus more on the hybrid stars with saddle instability which causes more marginal parameter space to constrain the dark matter models by the observed LIGO/Virgo gravitational wave events. In this way, we demonstrate how to constrain the particle physics models of dark matter by the gravitational wave astronomical observations.

BTM Criteria and Saddle Instability for Hybrid Stars
Neutron stars and some other compact stars are relativistic objects that their structure should be analyzed using general relativity. TOV equations [16][17][18] can be applied for these kinds of calculations, which are derived from the Einstein equations and the conservation of the energy-stress tensor, assuming zero space velocity, spherical symmetry, and an ideal fluid model.
In this paper, we will mainly consider the hybrid stars by solving the following TOV configurations with multiple component fluids inside the star [19,20]: where the index I labels the fluid components, and the total pressure and energy density are given by p := I p I and ρ := I ρ I , respectively; and m(r) = I m I (r) is the mass profile inside the star, and the Newton potential φ := 1 2 ln(−g tt ) with g tt the tt-component of the metric. In this paper, we will consider the cases with only two-component hybrid stars made of nuclear matter and dark matter.
The resulting configurations of hybrid stars depend on if there are strong interactions between component fluids. If there are strong interactions, it is expected to form a domain wall between phases of different fluids. Following our previous consideration in [5], we refer the hybrid stars with neutron core and dark matter crust as scenario I, and the ones with dark matter core and neutron crust as scenario II. On the other hand, if there is almost no interaction among component fluids, then all components prevail inside the star and mix together right from the core. We refer to this case as scenario III [5]. To solve the TOV equations, we need to provide the core pressure as the initial condition, and then, by changing the initial pressure, we will obtain different TOV configurations to plot the mass-radius relation. For scenarios I and II, one only needs to provide one initial condition for the core pressure since only one kind of fluid dominating the core, but needs two for scenario III.
Not every TOV configuration is stable, and one needs to judge its stability by either solving the Sturm-Liouville eigenmodes of radial oscillation (another method of judging the equilibrium configuration is used in [21][22][23]) or by some equivalent empirical criteria. One set of such criteria for single component fluid is the so-called BTM (Bardeen-Thorne-Meltzer) criteria [24], which states that, in the direction of increasing core pressure along the massradius curve, whenever an extremum is passed, one stable mode becomes unstable if the curve bends counterclockwise. Contrarily, one unstable mode becomes stable if it rotates clockwise.
These criteria were originally formulated by starting from the stable planet configuration with low enough core pressure. It, however, will cause some trouble if one does not start from such kind of planet configuration when solving the TOV equations. Thus, in practicality, it is useful to formulate the reverse BTM criterion by traveling the mass-radius curve in the direction of decreasing core pressure, and to require the consistency with the BTM ones.
This then leads to following (Reverse) BTM stability criteria: Traveling the mass-radius curve along either directions of increasing or decreasing core pressure, whenever passing an extremum, one stable mode becomes unstable if the cure bends counterclockwise; otherwise, one unstable mode becomes stable.
By applying the (reverse) BTM criteria, one can make sure that, for the unstable regions on the mass-radius curve, not even one starts from a known stable region. However, one can only pin down the stable regions if one starts from a stable region, or knows how many    since both branches are unstable. X 3 and X 4 show the saddle instability that one of the two branches is unstable. In some cases, the two branches may cross twice like points Y and Z; then, at least one of them is a fake intersection. Importantly, the above criteria is only checked rigorously for the stars with single component fluid. One should be careful when applying the (reverse) BTM criteria to judge the stability of hybrid stars. For the hybrid stars of scenarios I and II, there is only one dominant fluid component in each region, one will then expect that the (reverse) BTM criteria should still work without the need of much modification; see [25] for the recent discussion.
On the other hand, for the hybrid stars of scenario III, one needs to set the core pressures for both fluid components, say one for nuclear matter and one for dark matter. Thus, one needs to plot two kinds of mass-radius curves, namely, one by fixing the neutron core pressure but tuning the dark matter core pressure, and the other by fixing the dark matter core pressure but tuning the neutron core pressure. Let us call the former the dark branch and the latter the neutron branch. In the left panel of Figure 2, we show some representatives of both branches, e.g., the curve o a b c is the neutron branch, and the three curves in the same figure are the dark branches 7 . For any intersection point of dark and neutron branches, e.g., point X 1 to X 4 on the right panel of Figure 2, there could be saddle instability 8 if one of the branches is unstable at this intersection point. The difficulty is how to judge the stability for each branch. Without rigorous study of the Sturm-Liouville eigenmodes of radial oscillation, it is hard to answer this question. In this work, we assume that the (reverse) BTM criteria still work for each dark or neutron branch, and apply the criteria to find out the saddle unstable regions. Furthermore, we shall also assume that the cusp points such as B and B (B is not even an extremum) on the left panel of Figure 2 will not induce the change of stability of any radial oscillation eigenmode. Otherwise, there could be more complications [26].
Even with the above assumptions holding when traveling along the mass-radius curve, to judge the stability with (reverse) BTM criteria, it is better to start from some known stable region. One such region is the stable part of the pure neutron star curve, and the other is the stable part of the pure dark star curve. We then need to further assume that the small doping with the other component will not change the stability; it then implies that the nearby regions of the stable part of single component fluid stars are also stable. 8 For a single-component star, the stability is solely determined by altering the core pressure. When the core pressure has some perturbation, if the radial oscillation is stable, then the star is stable. However, for a two-component star, the core pressure is a sum of two partial pressures of each component. Then, a necessary condition for such a star being stable is that it is stable when changing either of the partial pressure while keeping the other fixed. It might happen that a star is stable when changing the neutron (or dark matter) pressure but unstable when changing the other partial pressure. This will still lead to an unstable star and we call it the saddle instability.
A being stable is a necessary condition for A to be stable. In addition, there are some fake intersection points such as point X for the case when one curve intersects with the other curve more than once. Then, only the end point/branch point is real intersection and others are fake. We have summarized these situations in the right panel of Figure 2. We will apply the lessons from the above discussions to the case studies in the next section to judge the stability of hybrid stars, especially the saddle instability of scenario III.

Dark Star and and Hybrid Star Interpretations
In this section, we consider the dark stars and hybrid stars of all three scenarios based on the EoSs extracted from the bosonic SIDM models discussed in Section 2 in the isotropic  We find that this dark star model can cover any mass range by adjusting the free parameter B 4 . In particular, for B 4 ≤ 0.047, the maximal mass exceeds 2.6 M . The maximal mass grows without an upper limit as the value of B 4 drops, and, for example, when B 4 = 0.035, we have 3.5 M . Since for scenarios I and II the 2.6 M hybrid stars can always be achieved as long as the pure dark star has a maximal mass higher than 2.6 M , in the following, we only concentrate on the more nontrivial scenario III, which is more realistic since it is assumed that there is no interaction between SIDM and baryonic matter.
For scenario III, by applying the (reverse) BTM stability criteria to the resultant massradius relation shown in Figure 4, no stable stars around 2.6 M can be formed. The pure SLy4 neutron stars are marked in red as the reference configurations. For the hybrid stars, B 4 = 0.045 case is denoted by the brown lines, and 0.035 case by the green lines. Now, we need the core pressures for both dark and nuclear matters to solve the TOV configurations.
By tuning one of the core pressures and fixing the other, we can obtain the so-called dark branch and neutron branch as discussed in Section 3. In Figure 4, the dark branches are denoted by the dash-dotted lines and the neutron branches by the solid lines.
To apply the method in Section 3 to check the saddle (in)stability, we can apply the (reverse) BTM criteria to both neutron and dark branches and then determine the saddle stability for their intersection points. As discussed in Section 3, we also need to assume that the nearby regions of the pure stable star configurations are stable when applying (reverse) BTM criteria. After checking this way for most of the intersections in Figure 4 one, which are all much shorter than the 25 km of φ 4 model.
In Figure 5, it shows the mass-radius relations for the hybrid stars of scenario III made of nuclear matter of SLy4 EoS and dark matter of φ 6 model's EoS with B 6 = 0.011 (brown lines) and 0.008 (green lines). The pure SLy4 neutron stars are marked in red as before. Similarly, the neutron and dark branches are denoted by solid and dash-dotted lines, respectively; and then we apply the (reverse) BTM criteria to judge the saddle (in-)stabilities as before. We find that the BDA line marked in Figure 5 bends higher than that in Figure 4, but still point  The unit of mass is M , and the unit of radius is km.  Bosonic φ 4 + φ 6 model. The EoS for this model is given in (32) and (33). Naively, it seems to be the intermediate model between φ 4 and φ 6 models. However, as we have seen in Section 2, it contains two tuning parameters B and β which lead to some novelty: there are some regions with p < 0 or even ρ < 0. We will only consider the region for positive p and ρ.
In Figure 7, it shows the mass-radius relations for the hybrid stars of scenario III made  Figure 4. We then run through the same check of saddle (in-)stabilities as before to pin down the stable and unstable regions by applying the (reverse) BTM criteria. As a result, we find that the maximal mass is mainly affected by the value of B, with lower B leads to higher mass, similar to the behavior of B n . On the other hand, β affects the compactness, namely, the radius becomes smaller as β grows from negative values to the positive ones. However, the change of radius significantly slows down when β ≥ 0.1. Though the β = 0 case is not included in Figure 7, this case is equivalent to φ 6 model with B 6 = (3B 2 ) 1/3 as discussed in Section 2. For example, B = 0.0006 corresponds to B 6 = 0.010, which is comparable to the B 6 = 0.011 case (brown lines) in Figure 5 .
As shown in Figure 7, although the EoS becomes considerably stiffer when β increases, we still cannot have stable stars around 2.6 M because there is a local minimum at D so that A is a saddle point. The blue closed path roughly encloses the stable region for

Conclusions
In this paper, we have extended the usual study of dark and hybrid stars for φ 4 SIDM to more general types of bosonic SIDM models by extracting their EoSs in the isotropic limit so that we can have more access to the complete mass-radius relations for the TOV configurations. Among them, we are especially interested in the φ n models which can have a stiff EoS when n is large enough, say n = 10. These kinds of models can be motivated by the UV Z n flavor symmetry, which may be a natural symmetry for some higher theories. It is fascinating to further explore this connection between particle physics and the resultant astrophysical compact objects in the future.
In general, it is easy to tune the parameters in SIDM to yield compact objects with masses comparable to or higher than the ones of neutron stars. Therefore, it is interesting to check the dark star possibilities by the future observations of the compact objects via the gravitational wave observations. Similar conclusions can be reached for the hybrid stars of scenarios I and II made of nuclear and bosonic SIDMs. In particular, many such kinds of dark and hybrid stars can have masses more than 2.6 M , and thus can be adopted to explain the recent GW190814 in which one companion compact object with such a mass has been identified by parameter estimation, which is hard to be explained by the usual EoSs for neutron stars.
On the other hand, the hybrid stars of scenario III are subjected to the saddle instability, and it is difficult for such kind of stars to have higher mass such as 2.6 M . However, in this paper, we do find such a stable massive configuration when the n of the φ n model rises up to 10 for which the saddle instability around 2.6 M is lifted. Although in practice it is hard to tell if a compact object can be the hybrid stars of scenario III in the near future, we may hope that the obstacle will be overcome in the long run to have precise measurements on the structure of the stars via gravitational wave detection to pin down the scenarios.
Despite that, theoretically, it is interesting to see the existence of some mass bound by the saddle instability. In this paper, we simply assume that the BTM criteria still work for the mixing fluids without the mutual interaction, applying it to judge the saddle instability. It will be illuminating to study the saddle instability by directly examining the eigenspectrum of the radial oscillation modes.
Overall, our work demonstrates the intriguing interplay between the particle physics models of dark matter and the astrophysical observations via the gravitational wave detection. We hope this will encourage more works to explore the dark matter physics via the study of dark and hybrid stars in the new era of gravitational astronomy.