Probing elastic quark phases in hybrid stars with radius measurements

The internal composition of neutron stars is currently largely unknown. Due to the possibility of phase transitions in quantum chromodynamics, stars could be hybrid and have quark cores. We investigate some imprints of elastic quark phases (only when perturbed) on the dynamical stability of hybrid stars. We show that they increase the dynamical stability window of hybrid stars in the sense that the onset of instabilities happen at larger central densities than the ones for maximum masses. In particular, when the shear modulus of a crystalline quark phase is taken at face value, relative radius differences between elastic and perfect-fluid hybrid stars with null radial frequencies (onset of instability) would be up to $1\%-2\%$. Roughly, this would imply a maximum relative radius dispersion (on top of the perfect-fluid predictions) of $2\%-4\%$ for stars in a given mass range exclusively due to the elasticity of the quark phase. In the more agnostic approach where the estimates for the quark shear modulus only suggest its possible order of magnitude (due to the many approximations taken in its calculation), the relative radius dispersion uniquely due to a quark phase elasticity might be as large as $5\%-10\%$. Therefore, the above dispersion of radii suggest that it might be possible to constrain the elasticity of the quark phase in neutron stars with electromagnetic missions such as NICER, eXTP and ATHENA.


INTRODUCTION
Neutron stars (NSs) are very compact remnants of stellar evolution that offer ways to probe particle and dense-matter physics aspects not possible in terrestrial laboratories. So far, the constitution and internal structure of an NS is not known in detail. The observational evidence for NSs with masses around 2 M (Demorest et al. 2010;Antoniadis et al. 2013;Fonseca et al. 2016;Cromartie et al. 2020) has helped constrain some microphysical models, but there are still many viable possibilities. The advent of the multi-messenger astronomy, thanks to direct detections of gravitational waves (GWs) by the Advanced LIGO (Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2015) detectors' network, promises further progress. Specifically, the observation of the GW170817 event, undoubtedly an inspiral of a binary NS system by the LIGO-Virgo collaboration (Abbott et al. 2017(Abbott et al. , 2018(Abbott et al. , 2019, has provided the first measurement of the late inspiral mutual tidal deformability of the components (in the form of the component mass-weighted sum of individual tidal deformabilities). This value, albeit burdened by large measurement errors, already constrains some microphysical models for the dense matter equation of state (EOS) with and without phase transitions (see, e.g., Chatziioannou (2020); Morawski and Bejger (2020); Blacker et al. (2020); Miao et al. (2020); Ferreira et al. (2020) and references therein). However, this single observation is not sufficient to definitely answer the question whether the GW170817 NS components were one-phase (hadronic) stars or exhibited phase transitions (De et al. 2018 Essick et al. 2020;. Additional inspiral and merger GW observations Chatziioannou and Han 2020), as well as the post-merger (Bauswein et al. 2019;Most et al. 2019;Weih et al. 2020) observations might resolve this dichotomy. It is of general interest to combine these future measurements with other observable quantities which could differentiate purely hadronic stars from stars with quark cores (hybrid stars). The same goal may be achieved with radius measurements by means of electromagnetic observations; note that GW measurements of binary NS inspirals do not provide direct radius measurements. This might now be possible with the NICER mission (Özel et al. 2016), which has already constrained an NS radius with an uncertainty smaller than 10% (at 1sigma confidence level) (see Raaijmakers et al. 2019;Bogdanov et al. 2019a,b;. Future missions such as the eXTP (Zhang et al. 2016) and ATHENA (Majczyna et al. 2020) promise even smaller uncertainties to NS radii (a few percent). ¶ An important question is whether radius differences of stars with and without phase transitions could be larger than NICER and eXTP uncertainties, and if there are suitable observational candidates presenting such features. When compared to purely hadronic stars, some models for the "third family" of NSs (Alvarez-Castillo and Blaschke 2017; Alford and Sedrakian 2017;Christian and Schaffner-Bielich 2020;Maslov et al. 2019) represent such a theoretical possibility. They are hybrid stars exhibiting substantial matter softening (resulting in a local minimum in the mass-radius sequence), for instance due to a first-order phase transition, where the density at the quark-hadron interface is discontinuous (Zdunik et al. 2008). So far-due to the lack of a quark ¶ Such constraints are not direct but rather byproducts of NS emission models, which attempt to explain for instance their light-curves and pulse profiles (Özel et al. 2016;de Lima et al. 2020). phase EOS-, there is not a definite threshold mass for hybrid stars, although there are statistical suggestions they might be large (Annala et al. 2020). Any possibility could be observationally covered since there is evidence that the NS mass function supports a bimodal distribution centered at approximately 1.4 M and 1.8 M , with spread around 0.2 M (Alsing et al. 2018). In summary, one has interesting candidates and detectors to probe aspects pertaining to hybrid stars.
However, in spite of the above optimistic observational scenario, ambiguities still remain concerning other aspects of their quark phases. It has been recently shown that the quark phase relevant to NSs could also be crystalline (Rajagopal and Sharma 2006), with a shear modulus possibly up to thousands of times larger than the one of usual crusts (Mannarelli et al. 2007). The question arises how it would be possible to differentiate between an elastic and a perfect-fluid model for a hybrid star when the quark phase is concerned. * * This would be of huge interest because it might allow for constraints of the quark phase. Tidal deformations seem excellent observables (Lau et al. 2017(Lau et al. , 2019 in this regard, although only third generation GW detectors, many combined GW measurements of the Advanced detectors, or an (unlikely) nearby high signalto-noise binary NS merger might reach the required precision; see, e.g., ) and references therein.
Electromagnetic observations may provide a similar and complementary information in the near future if models of stable elastic and perfect-fluid stars exhibit radii differences larger than the NICER and eXTP uncertainties. For a given EOS, radial stability analyses could easily determine this. However, one would not expect that null eigenfrequencies in elastic hybrid stars should satisfy the usual dynamical stability rule (Harrison et al. 1965;Friedman et al. 1988;Takami et al. 2011) ∂M/∂ρ c = 0, where M is the background mass of the star and ρ c is its central density, though this is the case for elastic stars without phase transitions (Karlovini et al. 2004). Indeed, ∂M/∂ρ c ≥ 0 for real radial eigenfrequencies arises in the context of cold, catalyzed matter constituting a perfect fluid (Harrison et al. 1965;Friedman et al. 1988;Takami et al. 2011), or a fully elastic star (Karlovini et al. 2004). For hybrid stars with elastic and liquid components, perturbation analyses are important in order to come up with the minimum radii of stable systems.
The article is arranged as follows. In Sec. 2, we lay out the models for elastic hybrid stars we will work with. Section 3 is devoted to the deduction of the field equations for radial perturbations in elastic hybrid stars and the appropriate boundary conditions to be taken in this case. In Sec. 4 we present our main results; we mainly show that an elastic quark phase increases the stability window of a hybrid star. Finally, we discuss all points raised in Sec. 5. Unless otherwise stated, we work with geometric units and our metric signature convention is (−, +, +, +).

ELASTIC HYBRID STARS
In this work we assume that a hybrid star is constituted of a quark phase and a hadronic phase and that they are split by a sharp phase-transition surface. By this we mean that we do not consider here a mixed phase and that the density at the phase-splitting surface is in general discontinuous.
2.1. Models for hybrid stars * * The elasticity of the crust of an NS is expected to change negligibly its dynamical stability because it contributes to a mass around 1% of the stellar mass.

Chiral effective model + MIT bag-like
In recent years, the chiral effective field theory (cEFT) has provided a framework that allows a systematic expansion of nuclear forces at low energies based on the symmetries of quantum chromodynamics (QCD) (Epelbaum et al. 2009;Hammer et al. 2013). The cEFT interactions can be employed in microscopic many-body frameworks in order to derive an EOS for neutron-rich matter. Since cEFT is an effective lowenergy theory, it contains a breakdown scale that imposes an upper density limit in the calculation of the EOS. A conservative choice of 1.1 ρ sat (ρ sat = 2.7 × 10 14 g cm −3 ) was adopted in Hebeler et al. (2013). Above this density, the EOS can be extended by means of piecewise polytropic EOSs which take into account the validity of causality and the consistency with observed 2M pulsars. Here we adopt the model of Hebeler et al. (2013) which employs a set of three polytropes valid in three consecutive density regions. This procedure leads to a very large number of EOSs, which verify the physical and observational constraints mentioned before. For the use in astrophysical applications, Hebeler et al. (2013) provides detailed numerical tables for three representative EOS labeled as soft, intermediate and stiff. For densities below 0.5 ρ sat the BPS crust EOS is used (Baym et al. 1971;. For the quark (innermost) phase, we take a phenomenological microscopic model put forth by Alford et al. (2005) which leads to the following EOS (see Pereira et al. (2018) for further details): where (B, a 4 , a 2 ) are free parameters. Here, B emulates the quantum chromodynamics confinement; a 4 accounts for the quark strong interactions and a 2 encompasses aspects such as the strange quark mass and color superconductivity (Alford et al. 2005). Further microscopic details on this model can be found in Pereira et al. (2018). These EOSs are glued together through Gibbs criteria for phases in mechanical, thermal and chemical equilibrium (oftentimes called "the Maxwell construction") (Shapiro and Teukolsky 1986). They lead to a unique phase transition pressure p pt and density jump at the phase-splitting surface [ρ] + − .

SLy4 + effective polytropic + simple MIT bag
The SLy4 EOS is based on Skyrme-type effective interactions in nuclear matter, originated in the ideas presented in the seminal works of Skyrme (Skyrme 1956(Skyrme , 1958. The Skyrme-type effective interactions were first successfuly applied, within the density functional theory (DFT), to atomic nuclei (Vautherin and Brink 1972) as well as to the nuclei in NS crusts . The effective interaction contains two-body terms and terms resulting from the averaging of the three-body interaction. The SLy effective hamiltonians (Chabanat et al. 1997(Chabanat et al. , 1998) contain a number of parameters that have been adjusted to reproduce experimental data on selected neutron rich atomic nuclei. A specific SLy4 effective hamiltonian is consistent with many-body calculations of the EOS of pure neutron matter with realistic two-nucleon and three-nucleon force (UV14+UVII model of Wiringa et al. (1988)), and so is suitable to describe the NS core. A unified SLy4 EOS for NS was obtained, based on SLy4 effective interaction model. It describes in a physically consistent way (i.e. starting from the same nuclear effective interaction) the structure and EOS of the crust and the core, including the transition between them (Douchin and Haensel 2001).
In order to construct equations of state for the hadronic part, we connect the SLy4 EOS to a relativistic polytrope p = κ ef n γ b [ρ = p/(γ − 1) + n b m b ] at a baryon density n 0 and extend it up to a higher baryon density n 1 . In the above expressions, n b is the baryon density, (n 0 , n 1 , γ) are free parameters, and m b (baryon mass) and κ ef are adjustable parameters assuring the continuity of the pressure and the chemical potential at n 0 . The simple MIT bag EOS p = c 2 s (ρ − ρ ) describing the quark phase is connected to the above hadronic EOS in particular. In this work, we choose γ = 4.5 and n 1 = 0.335 fm −3 (above 2ρ sat ), and take c 2 s = 1 in order to also probe aspects of stiff quark matter. The baryon density jump at the quark-hadron interface is a free parameter and chemical and mechanical equilibrium there imply a unique density jump and ρ . For further details on this model, see Sieniawska et al. ( Here we take into account crystalline aspects expected to the crust of NSs (Chamel and Haensel 2008), obtained from known constraints to nuclear matter. We work with a simple linear model for the shear modulus (Penner et al. 2011), namelyμ whereμ 0 and κ are tunable constants, and p is the pressure. We takeμ 0 0 and κ = 0.015 (Chamel and Haensel 2008;. We assume that the onset of the crust elasticity is at 2 × 10 14 g cm −3 . It is already known that such shear modulus leads to negligible tidal deformation changes in ordinary NSs (Penner et al. 2011;, described in terms of zero-frequency nonradial perturbations. Therefore, one would expect that for other types of perturbations they would also lead to negligible changes when compared to the perfect-fluid results. For completeness, we nevertheless include them in our analysis.

Shear modulus for the quark phase
Following the model of Mannarelli et al. (2007), we investigate crystalline aspects of color superconducting quark matter (the LOFF phase) (Alford et al. 2008), whose shear modulus can be hundreds of times larger than the shear modulus of the crust for some crystal structures (Mannarelli et al. 2007). Under certain assumptions, it has been shown that the shear modulus of the LOFF phase is (Mannarelli et al. 2007) where µ q is the quark chemical potential and ∆ is the crystalline pairing gap parameter. Conservative values for ∆ are in the range 5 − 25 MeV; µ q typically varies from 350 to 500 MeV (Mannarelli et al. 2007). Therefore, taken at face value, the shear modulus of the crystalline quark phase should be in the interval However, the pre-factor of Eq. (3) might vary considerably due to the many approximations and crystalline structures assumed for its calculation (Anglani et al. 2014;). This means that the above range ofμ q might also change appreciably. Due to this, we take the agnostic point of view that the estimates of Mannarelli et al. (2007) only point to the correct order of magnitude of the crystalline quark shear modulus. Thus, we allow the shear modulus of the quark phase (taken at the fiducial chemical potential 400 MeV) to be as large as 50 − 60 MeV fm −3 . In order to encompass the above uncertainty, hereafter we replace Eq. (3) by the phenomenological shear modulus where 1 α 4. In our numerical analysis, we take µ q as directly coming from the assumed microphysics.

RADIAL PERTURBATIONS IN ELASTIC NEUTRON STARS
The issue of linear radial perturbations in elastic stars has been addressed for the first time in Karlovini et al. (2004). They have shown that for one-phase stars the ordinary dynamical stability rules are not affected by elastic aspects of matter. However, to the best of our knowledge, this issue has not been addressed for hybrid stars with elastic and liquid parts in the presence of phase conversions.

Field equations for radial perturbations in perfect fluids
The equations describing the dynamics of perturbations have been obtained in the 1960s by Chandrasekhar (Chandrasekhar 1964a,b), and they lead to a Sturm-Liouville problem. For numerical analyses, though, it is more convenient to use a set of coupled first order differential equations derived in Gondek et al. (1997), written in terms of the Lagrangian displacements of the pressure and the radial coordinate.
We start with the equations for the perfect fluid case and then generalize them for elastic systems. For perfect fluids, they are (Gondek et al. 1997) with whereξ ≡ ∆r/r(∝ e iωt ), Γ is the adiabatic index of a star and ∆p(∝ e iωt ) is the Lagrangian displacement of the pressure; for hybrid stars Γ is in general a discontinuous distribution. The functions V(r), X(r) and Y(r) are related to background quantities, satisfying the Tolman-Oppenheimer-Volkoff (TOV) system of equations, namely, where p is the pressure, ρ is the mass-energy density and m the gravitational mass at r. The background metric is given by the Ansatz with ν given by whereas the function λ(r) is related with m(r) by means of 3.2. Adiabatic index and reaction timescales It is important to note that the adiabatic index, defined in general as where n b is the baryon number density, need not be taken as the equilibrium value, Γ eq (in the presence of perturbations). Its actual value depends upon weak reactions in the NS matter (Haensel et al. 1989;Gourgoulhon et al. 1995;Haensel et al. 2002). In particular, when weak reaction timescales, τ weak , are much smaller than the one of perturbations, τ pert , matter is in weak equilibrium at all times, and one is allowed to take Γ = Γ eq (Gourgoulhon et al. 1995). However, if τ weak τ pert , then one must take Γ = Γ frozen (Gourgoulhon et al. 1995), where "frozen" stands for a situation where matter is not in full thermodynamic equilibrium, and its composition does not change with perturbations. For practical purposes, it remains in a metastable equilibrium state. Clearly, the full thermodynamic equilibrium and the frozen composition are idealized cases. Real systems would have adiabatic indices between Γ eq and Γ frozen (Haensel et al. 2002).
Let us consider first the issue of weak reactions in the hadronic phase. Processes such as Urca and modified Urca lead to τ weak > τ pert (Haensel et al. 2002), and therefore are too slow to get Γ = Γ eq . Notwithstanding, it has been shown that in the regime of densities smaller than the threshold for the presence of hyperons-where the hadronic phase would usually live for first order phase transitions in hybrid stars-(Γ frozen −Γ eq )/Γ eq 15% (Haensel et al. 2002). The difference obtained in Haensel et al. (2002) is on the high side. For the SLy4 model of the hadronic phase, used in the present paper, (Γ frozen − Γ eq )/Γ eq 5% (see Fig. 3 of (Douchin and Haensel 2001)).
In order to estimate (Γ frozen − Γ eq )/Γ eq for the quark core, we have used the approach applied in Haensel and Zdunik (2006) within the MIT bag model for uds quark matter. At ρ ∼ 10 15 g cm −3 , the dependence of (Γ frozen − Γ eq )/Γ eq on the QCD coupling constant and B was very weak, while the dependence on the quark strange mass, m s , was strong (∝ m 6 s ). As a reference, for the quark mass m s = 200 MeV/c 2 , we get (Γ frozen − Γ eq )/Γ eq 0.5%. Therefore, the equilibrium adiabatic index is a good approximation for the adiabatic index of the quark phase in the presence of perturbations.
Finally, the assumption Γ = Γ eq is expected to result in the minimum violation of the usual dynamical stability rules due to the quark phase elasticity. (We already know that the stability rules are violated when Γ = Γ frozen for a liquid hadronic phase (Gourgoulhon et al. 1995).) Thus, this choice is important for back-of-the-envelope estimates which give the order of magnitude of the effect and single out relevant cases for future analysis. Due to this and the simplicity in obtaining Γ eq , we follow the above-mentioned route in this work.

Elastic stresses
We advance now to the case where some star phases might be elastic. We assume that this could just be the case in the presence of perturbations. In other words, the background phases of stars are perfect fluids and their aspects are obtained normally via the TOV system of equations. In order to investigate the dynamical stability regions in the M(R) or M(ρ c ) relations, we restrict our analysis to radial perturbations. In this case, shear forces manifest themselves as restoring forces due to the radial Lagrangian displacements of volume elements from their equilibrium positions. We follow the shear description given by Penner et al. (2011), which we rewrite to the radial case (for more details on the formalism, see Andersson et al. (2019) and references therein). In the presence of perturbations, the stress-energy momentum tensor gains a (Eulerian perturbed) component due to shear stresses and it is given by (Penner et al. 2011;Andersson et al. 2019) whereμ is the shear modulus, P ab ≡ g ab +u a u b is the projector onto the orthogonal direction of the four-velocity u a , and In the above equation use has been made of the shortcut ∂ c ≡ ∂/∂x c and Γ a cd are the usual Christoffel symbols (connection coefficients) (Landau and Lifshitz 1975). Note here that ξ a ≡ ∆x a = (0, ξ, 0, 0) and ξ b = g ab ξ a . In addition, from Eq. (14), where δλ and δν are functions to be fixed by the perturbed (first order) Einstein equations and would generalize the quantities valid for perfect fluids (see, e.g., Misner et al. (1973); Chandrasekhar (1964b)). For the case of radial oscillations, it is also simple to show that from the definition of u a and u a u a = −1, where we have defined the "dot" operation as the time derivative (Ȧ ≡ ∂ t A). By using the Christoffel symbols for the metric (14) and keeping terms up to first order in ξ, it follows that the only nonzero components of δΠ b a are the diagonal ones (F ≡ dF/dr), and δΠ θ θ = δΠ φ φ = −δΠ r r /2. Just for completeness, where the first term is the perturbation of the energy momentum tensor of the perfect-fluid background.
3.4. Generalized perturbation equations Now we proceed with the generalization of the radial perturbation equations for elastic stars. Note first that δΠ b t = 0 implies that the [rt] component of the Einstein equations does not change with respect to a perfect fluid. Thus (Misner et al. 1973;Chandrasekhar 1964b), where the subscript "0" has been used for the background quantities. In the second equality of the above equation, use has been made of the background field equations. From the [tt] components of the Einstein equations, it also follows that (Misner et al. 1973;Chandrasekhar 1964b) From the perturbed [rr] components of the Einstein equations and Eq. (24), one has that which could also be further simplified by the background equations, but we choose not to do so here.
One could also work with other Einstein equations for obtaining the final pulsation equation, but it turns out to be more efficient to work with T a r ; a = 0. In the spherically symmetric case, it implies that In the case of radial perturbations, it is easy to show that (Chandrasekhar 1964b;Misner et al. 1973) From Eulerian perturbations, One already has all ingredients to calculate δT a b . From the perfect-fluid background and Eq. (23) (We note that under shear stresses a star becomes anisotropic, and our model could be seen as particular realization of the more general analysis of Raposo et al. (2019). Indeed, as we shall show, more compact stars could emerge in this case.) Therefore, when Eqs. (28), (29) and (30) are replaced into Eq. (27), one has (δΠ θ θ = −δΠ r r /2) Note that the above equation is the general equation for radial perturbations in elastic stars once δp is found, given that all other terms are already known. Actually, for adiabatic processes, δp can be easily obtained. Indeed, from Eq. (17), one learns that Given the perturbations, one can always find ∆n b /n b geometrically by means of (Penner et al. 2011;Andersson et al. 2019) In the spherically symmetric case, where, in the last equality of the above equation, we have used Eq. (24). Formally, Eq. (34) is the same as its perfect-fluid counterpart.
As in the case of perfect fluids, for elastic stars it is also numerically convenient to solve equations for ξ and ∆p as if they were independent variables. This can be easily found with the help of Eqs. (31), (32) and (34). From Eqs. (32) and (34), when the variableξ ≡ ξ/r is inserted, it follows that [(ξ, ∆p) ∝ e iωt is assumed from now on] which is exactly Eq. (6). For obtaining the above equation, use has been made of Eq. (15). The second equation for the coupled system of equations involvingξ and ∆p can be easily found from Eq. (31) when the first equality of Eqs. (32) and (26) are taken into account. It can be schematically written as Naturally, in addition to the perfect-fluid part, several other terms also appear due to the elasticity of the star. In order to see their influence on d∆p/dr, we rewrite δΠ r r based on Eq. (35). With the use of Eq. (24), Eq. (22) simplifies to 3.5. Boundary conditions to elastic hybrid stars We now set out the boundary conditions to the problem of radial perturbations in elastic hybrid stars. It turns out not all of them are the same as in the perfect-fluid case (Pereira et al. 2018). Let us explicitly show this. Equation (35) is only well defined at the origin if which is always meaningful whenξ(0) is finite. From Eq. (36), we note that one has also to impose that δΠ r r (0) = 0.
However, this is automatically satisfied through Eq. (38). At the surface of the star (r = R), one should impose that with the obvious condition thatξ(R) is finite.
In hybrid stars, one should also provide boundary conditions toξ and ∆p at the phase-splitting surfaces.
Here we choose to work with nucleation processes around these interfaces. Essentially, they could be classified into two asymptotic regimes: when the timescales of nucleation (τ nucl ) are either much larger or much smaller than the ones of perturbations (τ pert ). When τ nucl τ pert , volume elements are instantaneously converted from one phase to the other and the phase conversion is said to be "rapid". For the case τ nucl τ pert , nucleation processes occur at a very low speed and for all purposes volume elements just squash or stretch when perturbed. For this reason we call them "slow phase conversions".
Just for definiteness, assume that a phase-splitting surface in a hybrid star is at r = R pt ("pt" stands for "phase transition") and jumps of physical quantities are with respect to it. It has been shown (Pereira et al. 2018;Haensel et al. 1989) that for slow (s) phase conversion processes, the extra-boundary condition to be taken for ξ is while for rapid (r) phase transitions, Naturally, Eq. (40) is not only related to slow nucleation processes. It is a standard boundary condition to be taken into account when there is no mass fluxes from one phase to the other. Slow conversions happen to be a special case where these conditions are fulfilled.
We turn now to the jump condition for ∆p at R pt . This is obtained as a byproduct of the perturbation equations, and it could be more easily obtained from Eq. (36). When it is promoted to a distribution, the Dirac delta terms in Eq. (36) imply that [∆p − (∆p) perf + δΠ r r ] + − = 0. Since for perfect fluids [(∆p) perf ] + − = 0 (Pereira et al. 2018; Tonetto and Lugones 2020), it trivially follows that which is a consequence of the continuity of the radial traction at R pt   (41) might also have important implications for the dynamical stability of stars with elastic parts. We recall that it is already known that when ∆p is continuous, Eq. (40) results in ω s = 0 for ∂M 0 /∂ρ c < 0 and Eq. (41) leads to the presence of the reaction mode (Haensel et al. 1989;Pereira et al. 2018) (an extra-mode which does not have an one-phase counterpart), though for rapid phase conversions ∂M/∂ρ c = 0 when ω r = 0. At least one would expect quantitative changes of these results in the presence of elasticity.
Put in the above way, the problem of radial perturbations in stars with elastic parts and phase transitions is a Sturm-Liouville problem. Thus, the set of characteristic eigenfrequencies ω is discrete and hierarchic (ω 2 0 < ω 2 1 < ω 2 2 ...., where ω 2 i is the ith-mode). For dynamical stability purposes, it suffices to investigate the eigenfrequencies of the fundamental mode.
The appropriate quark-hadron boundary condition critically depends on the physical processes and scenarios one wants to investigate. When a hadronic fluid element in the neighbourhood of the sharp interface is compressed beyond the transition pressure, a rapid conversion to LOFF quark matter is expected to be strongly suppressed at low enough temperatures. To understand this, let us focus first on unpaired quark matter. Typical timescales for nucleation driven by quantum fluctuations are much larger than the age of the Universe for temperatures below 10 MeV (Bombaci et al. 2016). If thermal nucleation is considered, the timescale is also larger than the age of the Universe for T < 5 MeV but decreases to ∼ 1 s for T ≈ 10MeV (Bombaci et al. 2016). These numbers are significantly larger than the typical period of radial oscillations strongly suggesting that the assumption of slow phase conversions is appropriate for cold enough matter. In the case of quark matter in a LOFF crystalline phase, this conclusion is expected to remain valid because an additional time (with respect to unpaired matter) would be needed to form the specific three-dimensional lattice structure needed for LOFF formation (Rajagopal and Sharma 2006). As stressed previously, nucleation timescales can decrease significantly for hot NSs. It has been suggested that they may be smaller than typical perturbation timescales for temperatures around 20 MeV (Tonetto and Lugones 2020). In this case, conversions around a phase-splitting surface should be rapid. Examples of stars fulfilling these conditions would be those in the aftermath of binary NS mergers. However, a natural caveat here would be the need to work with hot equations of state.
Given that we are mainly interested in cold stars (both due to the possibility of probing the crystalline colorsuperconducting phase and the fact that such stars are the main targets of missions such as NICER and eXTP), we mostly focus our analysis on slow conversions. Nonetheless, for completeness, we also comment on rapid conversions.

RESULTS
In the forthcoming analysis we do not make a systematic study of the space of parameters of the quark phase, but just choose some representative parameters. We choose (B, a 4 , a 2 ) in such a way to respect the constraints already measured for NSs such as tidal deformabilities and maximum masses, as well as assume that the hadronic phase generally starts at densities larger than the nuclear saturation density (ρ sat = 2.7 × 10 14 gcm −3 ). For the hadronic phase, we work either with chiral models or the SLy4 + polytropic equation of state, as explained in Sec. 2.1. Table 1 summarizes the main aspects of the models used in this work and Fig. 1 shows some associated M − R relations for hybrid stars. Just to check the consistency of numerical calculations, we also make use of the (book-keeping) polytropic EOS p = Kρ 2 , with K = 100 km 2 for the hadronic phase. Details about the connection (through the Maxwell construction) of this EOS with the MIT bag-like model are given in . We compute the tidal deformability (Λ) based on perfect fluids (we follow the formalism of Hinderer (2008); Damour and Nagar (2009);Hinderer et al. (2010)), even though we are assessing elastic stars. The reason is because an elastic star should have a smaller Love number than a perfect-fluid star (Penner et al. 2011; and if the latter already satisfies the GW constraints, so will the former. In what follows, we focus more on the case of slow conversions given that they would be expected for cold stars. We use the phenomenological shear modulus given by Eq. (5) in order to reflect the uncertainties in the calculations of Mannarelli et al. (2007). TABLE 1 Main aspects of the hybrid star models used in our work. Our notation is such that ρ q is the density at the top of the quark phase and ρ h is the density at the bottom of the hadronic phase. Therefore, η is directly related to the density jump at the quark-hadron phase-splitting surface   Table 1. We mostly focus on hybrid models with the SLy4 EOS for different density jumps (H3 models), but also consider the chiral effective theory (H1 and H4 models). From bottom to top of the M(R) thick curves (from gray to red) we have chosen n 0 = (0.235, 0.21, 0.185, 0.16) fm −3 . The dashed lines have the same n 0 as the thick yellow line (n 0 = 0.21 fm −3 ), but different density jumps. For all H3 models, we take n 1 = 0.335 fm −3 (∼ 2.1ρ sat , the density marking the end of the hadronic phase). The upper horizontal band corresponds to the observed mass of the pulsar PSR J0348+0432, while the lower one relates to PSR J1614-2230 (Arzoumanian et al. 2018;Antoniadis et al. 2013). The diamonds on each EOS mark the terminal masses (ω 2 = 0) for perfect-fluid hybrid stars under slow conversions. Stable stars are on the right of the markers. For elastic hybrid stars, ω 2 = 0 is always more to the left of the diamonds and their precise locations depend on the shear modulus of the quark phase.
For cold, catalyzed stars with no phase transitions (onephase stars), powerful theorems deduced in Harrison et al. (1965) state that for dynamically stable stars, ∂M/∂ρ c ≥ 0. Here M is the background (gravitational) mass of a star with a given EOS and a central density ρ c . It is already known that one-phase stars under slow weak reactions (slightly) violate this condition (Gourgoulhon et al. 1995). However, appreciable violations might happen for hybrid stars under slow phase conversions (Pereira et al. 2018), another reason why we choose to investigate them. Table 2 summarizes some stability aspects of hybrid stars with elastic quark phases for α = 1 and ∆ = 25 MeV (μ (p) q,400 ∼ 15 MeV fm −3 ). (Table 3 gives the phase transition parameters for the H1, H3 and H4 EOSs.) The numbers are obtained via the numerical solution of the perturbation equations from Sec. 3 with different boundary conditions (slow and rapid phase conversions). One sees that relative radius changes for the zero eigenfrequency for both conversions are 0.1 − 2%. (Slow conversions lead to relative radius changes smaller than rapid ones.) Such differences are with respect to the radii of null eigenfrequencies of perfect-fluid stars. Observationally speaking, this would translate in the following way. Due to elasticity and slow conversions, a given observed mass between M elas ω=0 and M perf ω=0 could now be on either side of the maximum of the M − R relation. Thus, the above changes roughly mean that the expected relative radius dispersion exclusively due to a quark phase elasticity for stars with masses in the above mass range would be up to 0.2 − 4%. Larger radius differences are also possible for larger values of α. Table 4 shows the associated changes for α = 4 and slow conversions. Stiff EOSs suggest that the maximum increase of relative radius differences solely due to elastic quark aspects will surpass the 5% threshold forμ (p) q,400 ∼ 50 MeV fm −3 . Soft EOSs-such as some representatives of the H3 modelsuggest that the largest radius difference increases due to the elasticity of the quark phase are around 2−3%. Possible ways to connect the above-mentioned radius changes with observations and associated difficulties are discussed in the next section.
Tables 2 and 4 also suggest a nontrivial consequence of slow conversions for EOS models similar to the H4 model: masses that result in ω 2 s = 0 may also be very different from the maximum one. † † However, for the H4 model in particular, the terminal radius and mass of elastic NSs only change by approximately 0.5% with respect to their perfect-fluid counterparts. (When stiff EOSs for the chiral model are taken into account, maximum relative changes for the terminal mass and radius of elastic hybrid stars are around 1.5%.) We discuss more about the above class of EOSs in the next section.
One also notes from Tabs. 2 and 4 that the relative radius change depends on the density jump at the quark-hadron phase transition. However, for slow conversions, the relationship is not trivial. They depend, for instance, on the stiffness of the EOS, the values the chemical potential take in the quark phase, the relative size of the crust (when compared to the core), and where in the M − R diagram the terminal radius is reached for a perfect-fluid star. ‡ ‡ In addition, it has a de- † † It should be mentioned that the physics around the maximum mass for the H4 model is different from the other models we have considered. It resembles more the situation around the phase transition masses (small quark cores) for the H1, H2 and H3 models in the case of large energy jumps. Slow conversions render non-negligible regions of their M−R relations meta-stable when ∂M/∂ρ c < 0, similarly to what happens for the H4 model. This is the reason why for the H4 model the terminal mass for a perfect-fluid star is around 15% smaller than its maximum mass, and radius differences for a mass within this extended branch could be up to 2 km. ‡ ‡ For rapid conversions, there is a simple trend for the relative radius change (due to elasticity) and the quark-hadron density jump-mostly owing pendency on the EOS of the quark phase. For the H3 models, changes are small even for stiff EOSs, while for the H2 model, the radius differences are more pronounced for stiffer EOSs.
For completeness, we plot in Fig. 2 the square of the fundamental eigenfrequencies of the H3 and H1 models (∆ = 25 MeV) for slow conversions with and without taking into account the elasticity of the quark core. For α = 1, one sees that the imprint of elasticity on the eigenfrequencies is mostly relevant for large masses: relative changes would be larger than 10% (∼ 20% for the frequency squared) for NS masses above 2.0 M (for stiff EOSs). Naturally, if the value of the shear modulus is larger (larger α), significant differences between the eigenfrequencies of a liquid and an elastic quark phase take place at smaller masses. For instance, taking α = 3.5 and the H1 model, one finds that relative frequency changes above 10% occur for masses larger than 1.9 M . We note that the masses (central densities) where ω = 0 are always larger than the ones where relative frequency changes become relevant. Thus, if direct or indirect observations of NS frequencies are possible, quark matter aspects for a given EOS could also be constrained with them. We return to this issue in the next section too.

TABLE 2
Hybrid star models (see Table 1) and some of their properties for ∆ = 25 MeV and α = 1 [μ (p) q,400 ∼ 15 MeV fm −3 ; see Eq. (5)] in the case of slow and rapid phase conversions. Here, "elas" stands for elastic, "perf" for perfect fluid, "s" for slow and "r" for rapid. In addition, Λ perf 1.4M is the tidal deformability of a 1.4M perfect-fluid star and ρ c is the central density of an NS. All parameters have been chosen so that their maximum masses are above 2M and Λ 1.4M 660, in agreement with gravitational and electromagnetic wave observations from NSs. We note that maximum masses in the M − ρ c relation coincide with null frequencies for rapid phase conversions in the perfect-fluid case.

DISCUSSION AND CONCLUSIONS
Elastic aspects to hybrid stars make them distinctly different from their perfect-fluid counterparts because they have dissimilar restoring forces upon perturbations. This naturally influences the dynamics of perturbations. In particular, the points in the M(R) or M(ρ c ) relations where eigenfrequencies are null are not their critical points anymore. This is expected because the classical rules of dynamical stability of to the fact that for perfect fluids the terminal radius is always related to the maximum mass. Around such a mass, a larger energy jump would imply a smaller hadronic phase and hence stars "on the way" of being one-phased; in this one-phased case, it is known that the usual dynamical stability rule holds (Karlovini et al. 2004 q,400 ∼ 60 MeV fm −3 ) in the case of slow conversions. The maximum relative radius change for a given mass range exclusively due to quark elasticity would be approximately twice as large for all the models. The largest differences relate to stiff EOSs.
and references therein), differently from rapid conversions. The nondotted vertical lines mark the densities where the eigenfrequencies of hybrid star sequences are null (we respect the style and color of each curve for clarity). They are all larger than the densities respectively associated with their maximum masses (ρ c = 7.142 ρ sat for the H1 EOS and ρ c = 10.584 ρ sat for the H3 EOS). At face value, the maximum shear modulus already leads to non-negligible changes to the eigenfrequencies at smaller densities than the ones for their maximum masses. Indeed, for the H1(H3) model with α = 1, relative frequency differences of 10% take place at ρ c = 5.78 ρ sat (ρ c = 10.17 ρ sat ), which corresponds to M = 1.99 M (M = 2.166 M ). When α = 3.5, instead, we have that for the H1 model (stiff) 10% differences happen at ρ c = 4.41 ρ sat (M = 1.92 M ). In addition, ω 2 s is zero at ρ c = 9.02 ρ sat (M = 1.99 M , R = 11.60 km). In this case, when compared to their perfect-fluid star counterpart (see Tab. 2), the largest increase of relative radius differences would be around 5.2 %. minal central density (where ω = 0) of an elastic star is larger than those of a perfect fluid is also reasonable, because elastic shear stresses increase the pressure of the system, allowing for denser objects. This is in agreement with the more general analysis of Raposo et al. (2019), who showed that even ultracompact stars are possible when larger anisotropic pressures take place. An interesting consequence for systems with ∂M/∂ρ c < 0 is that they could be more compact, which could lead to a richer phenomenology that may be probed with electromagnetic missions and relativistic ray-tracing models (see, e.g., Vincent et al. (2018); Raaijmakers et al. (2019)).
We have shown that the radial stability of hybrid stars strongly depends on the value of the shear modulus of the elastic quark phase, which at present is not yet known. The key point is whether or not stability changes due to it could be detected. In terms of GW measurements, the hope for very precise radius constraints (inferred from the component masses and tidal deformability measurements) lies basically in the third generation detectors, or in collecting many detections which will lead to analyses unraveling the NS interior through features visible in the population of measurements, or in a lucky nearby GW170817-like strong signal, etc. Therefore, it is an issue in principle for the future.
The alternative relies on electromagnetic observations. The NICER mission, which is currently collecting data from some compact systems, has already delivered NS radius constraints with 10% accuracy. For stars which are brighter, have favourable geometries (e.g., particular inclination angles with respect to the observer and colatitudes of the polar caps), etc., it might be possible that even smaller uncertainties are reached, for instance breaking the 5% threshold (Özel et al. 2016). This means that according to our analysis in the most promising scenario, NICER might be able to set upper limits to the shear modulus of the quark phase for a given hybrid star EOS. Our estimates suggest that it might be 50 MeV fm −3 for stiff EOSs. Future missions, such as the eXTP and ATHENA, might be more promising in constraining elastic quark parameters for a given equation of state because they are supposed to measure NS radii more accurately. However, one should also bear in mind that the maximum increase of radius differences for stable stars due to elasticity is EOS-dependent. Naturally, this brings ambiguities to the shear modulus extraction. A possible way to observationally address this problem is to have several radius measurements of stars with similar masses. The required number of measurements depends on the parameters of the EOS and the shear modulus. The mass value on which to focus should be, broadly speaking, between M elas,s ω=0 and M max (for a given EOS and shear modulus). Since the minimum value for M elas,s ω=0 is not known, no mass should be disregarded in principle; however, our analysis roughly suggests that stars with masses larger than approximately (1.7 − 1.8) M might be interesting candidates for probing the elasticity of the quark phase. It is still to be checked if such stars could be feasible candidates for the NICER, eXTP and ATHENA missions.
The elasticity of the quark phase could also leave an important imprint on the stellar eigenfrequencies. We have found that the relative changes (with respect to the perfect-fluid case) are larger for stiff EOSs, and they could be larger than 10% (∼ 20% for the squared frequencies-related to the energies of the modes) for the massive stars to be probed. Hence, if NS radial modes could be observed, the quark elasticity might in principle be probed. A possible way of doing that would relate to modulations of the NS lightcurve (Chirenti et al. 2019). However, this might work only for large enough amplitudes of oscillation, which are generally unknown. Another possibility concerns direct GW observations due to the coupling of radial modes with rotation (Chau 1967). However, this seems a possibility only for postmerger stars (Dai 2019), where our analysis is not reliable due to the use of cold EOS candidates and the unknown impact rapid rotation could have on the terminal mass. In this case, though, boundary conditions related to rapid conversions seem more appropriate; we leave these analyses for future work.
Finally, we point out that depending on the nature of the EOS, a large extended branch of (metastable) NSs could appear if conversions are slow. In this case in particular, the terminal mass (where ω = 0) could be significantly different from the maximum mass. As a result, it might reach the range of masses more commonly observed for NSs. One question in this context is whether this is relevant for probing an elastic quark phase. Our analysis suggests that it seems unlikely to probe the elasticity of the quark phase electromagnetically if the EOS is soft or intermediate. The possibility might exist for stiff EOSs, but satellite missions such as eXTP and ATHENA might be required.
Summing up, in this work we have laid out the problem of radial stability for hybrid stars with elastic quark phases and have shown that their mass-radius region of dynamical stability is extended with respect to their perfect-fluid counterpart. Our analysis suggests that significant increases of radius differences between elastic and perfect-fluid stars with null frequencies might take place when the shear modulus of the quark phase is roughly larger than 50 MeV fm −3 .

ACKNOWLEDGEMENTS
We thank Andreas Schmitt and Nils Andersson for helpful discussions in an earlier version of this work. J.P.P. and M.B. acknowledge the financial support from the Polish National Science Centre grant No. 2016/22/E/ST9/00037. J.P.P. is also thankful for partial support at an earlier stage of this paper given by