© The Authors 2025

1. Introduction

Neutron stars (NSs) can sustain a variety of oscillation modes due to their stratified internal structure and composition. These normal (or quasi-normal, when the frequency is complex) oscillation modes include fundamental (f), pressure (p), and gravitational (g) modes, among others, each characterised by distinct frequencies and damping times (Thorne & Campolattaro 1967; Reisenegger & Goldreich 1992; Andersson et al. 1996; Kokkotas & Schmidt 1999). A detection of gravitational waves (GWs) emitted by quasi-normal oscillations would allow the direct observation of the dominant mode frequencies, enabling a new way of probing NSs’ internal properties and dynamical processes (Andersson 2021; Jones 2022; Andersson 2019). For example, it has been suggested that the p₁ mode carries information that can be used to distinguish between nucleonic, hybrid, and strange stars (Vásquez Flores & Lugones 2014). To date, forthcoming runs of LIGO, Virgo, Kagra gravitational wave interferometers, and the planned Einstein Telescope and Cosmic Explorer represent a promising avenue for detecting these oscillations (Andersson et al. 2011; Piccinni 2022; Jones 2022). This holds the appealing prospect of integrating such observations with other data – such as results from NICER (Özel et al. 2016) and the planned Athena spacecraft (Majczyna et al. 2020) – to constrain the equation of state of dense matter.

For cold NSs, the subject of this study, non-radial f modes are expected to be excited during magnetar flares (Levin & van Hoven 2011; Ball et al. 2024) and pulsar glitch events (van Eysden & Melatos 2008; Bennett et al. 2010; Ho et al. 2020); we invite the reader to consult (Antonelli et al. 2022; Haskell & Jones 2024) for a recent review and Yim et al. (2024) for an analysis of glitching pulsar candidates as priority targets of future observations. This is an attractive possibility, as future detections of these modes may be used to discriminate between neutron and quark stars (Wilson & Ho 2024; Sotani et al. 2011).

In addition to probing the internal structure of NSs, oscillation modes can be used to disentangle macroscopic characteristics such as mass and radius when used in tandem with other observations. In fact, Andersson & Kokkotas (1998) found a set of quasi-universal (QU) relations – in the sense that they are almost EoS-independent relations (Yagi & Yunes 2017) – between the normal mode frequencies and the average density or the compactness. To date, there are numerous studies presenting different QU relations for mode frequencies, usually tested with a small sample of EoSs (Tsui & Leung 2005; Benhar et al. 2004; Pradhan et al. 2022; Sotani 2021) or a large set of purely barotropic (namely, zero temperature and β-equilibrated) agnostic matter models of the kind used in, for example, Lindblom (2010), Breu & Rezzolla (2016), Fasano et al. (2019), Moustakidis et al. (2017), Yao et al. (2024). We will perform a systematic study of proposed QU relations for nucleonic NSs’ oscillation modes by using a large set of EoS models that are compatible with the latest astrophysical observations and nuclear physics constraints. This is done by using the phenomenological metamodel technique (Margueron et al. 2018), which allows us to explore the parameter space of cold npeμ EoSs and, at the same time, to include the constraints from the chiral effective theory, experimental nuclear physics, and astrophysical observations via a Bayesian framework (Zhang et al. 2018; Carreau et al. 2019a; Güven et al. 2020; Dinh Thi et al. 2021a). Furthermore, the metamodel is able to reproduce existing realistic nucleonic models and interpolate between them (Mondal & Gulminelli 2022; Davis et al. 2024).

A downside of exploring a wide parameter space for the metamodel representation of the EoS is that we have to find the f- and p₁-mode frequencies for a large set of stellar structures, making it impractical to calculate the frequencies in full general relativity. Therefore, we chose to work within the Cowling approximation, which greatly speeds up the computation of the frequencies. In doing so, we also tested the impact of assuming two opposite idealised limits¹ for matter undergoing time-dependent compression (Haensel et al. 2002; Andersson 2019), which are described below.

1. Frozen regime – In this limit, the local relaxation processes that bring back npeμ matter back to β equilibrium do not have time to occur, as the compression-expansion cycle imparted by the oscillation is faster than the typical reactions mediated by the weak interaction. This limit is characterised by the local conservation of the chemical fractions, meaning that fractions are purely advected by the fluid motion. This is the limit expected to hold in cold NSs.

2. Equilibrium regime – In this limit, the relaxation processes are so fast that each fluid element has a negligible departure from β equilibrium, so the matter model reduces to that of a perfect barotropic fluid. Given that the relaxation processes are mediated by the weak interaction, this limit might be expected to hold only in high temperature processes, such as proto-NSs and post-merger oscillations.

Similar to the approach taken for non-compact stars (Hansen & Kawaler 1994), evaluating the mode frequencies requires knowledge of the adiabatic index, which determines how pressure responds to changes in local baryonic density (Thorne & Campolattaro 1967; Shapiro & Teukolsky 1983). The choice of one of the two limits can significantly impact the local value of the adiabatic index and, consequently, the pressure response of npeμ matter (Haensel et al. 2002; Andersson 2019). In particular, calculating mode frequencies using a barotropic equation of state and a consistent adiabatic index inherently assumes an equilibrium regime.

To investigate the frozen regime, we need an EoS model that is not purely barotropic, allowing the pressure (or adiabatic index) to be calculated at fixed chemical fractions. The metamodel representation of the cold (neutrinoless) npeμ EoS provides this possibility, enabling a systematic comparison of mode frequencies derived from a purely barotropic EoS versus those that account for the effects of a frozen composition.

In this work, we extended the type of Bayesian analysis performed in previous studies (Zhang et al. 2018; Carreau et al. 2019a; Güven et al. 2020; Dinh Thi et al. 2021a; Davis et al. 2024) by solving, for a large set of metamodel instances, the perturbation equations in the Cowling approximation in the two idealised – that is, frozen and equilibrated – regimes, testing possible deviations from the proposed QU relations. In Sect. 2, we recall the relevant properties the metamodel representation of the energy of cold npeμ matter. Sect. 3 outlines the Bayesian technique developed for our inference: a large number of metamodel instances are assigned with a likelihood depending on how they satisfy astrophysical and nuclear constraints. Then, in Sect. 4, we summarise how the mode frequencies are obtained for each metamodel instance. Finally, the resulting mode frequencies and their posterior distributions – which may be interpreted as a possible frequency range for a future detection – are given in Sect. 5.

2. Metamodel representation of the equation of state and internal composition

The metamodel representation of the nucleonic EoS of an NS was introduced in Margueron et al. (2018). The fundamental assumption is that an NS’s core consists of npeμ matter in weak equilibrium, disregarding the possibility of having other degrees of freedom; albeit, it is possible to modify it to account for phase transitions to quark matter (Mondal et al. 2023). The EoS for the uniform npeμ matter in the core is then consistently prolonged to the lower density layers of the solid crust thanks to the compressible liquid-drop model approach described in Carreau et al. (2019b) and Dinh Thi et al. (2021b). Although not as microscopic as other approaches, this method reproduces results that are consistent with extended Thomas-Fermi calculations at both zero (Grams et al. 2022) and finite temperature (Carreau et al. 2020). Furthermore, it enables the quantitative estimation of a unified EoS for both the core and the crust at a relatively low computing cost.

Within the metamodel technique, each unified² EoS model is represented by ten independent empirical parameters that correspond to the coefficients of a fourth-order Taylor expansion of the uniform matter binding energy in the isoscalar and isovector channels around saturation density. For non-homogeneous matter, they are supplemented by five further surface and curvature parameters (Dinh Thi et al. 2021b), which are selected by fitting the experimental atomic-mass evaluation nuclear-mass table (Huang et al. 2021) for each set of the ten aforementioned parameters. The density dependence of the symmetry energy and the energy in symmetric matter are characterised by these parameters, and over a wide range of nuclear data, their prior distribution is in agreement with current empirical information (Margueron et al. 2018). Three more parameters are needed: two to account for the density dependence of the effective mass and the effective mass splitting; and one that enforces the correct behaviour at zero density, for a total of 13 independent parameters.

As far as this study is concerned, the metamodel can be thought of as a procedure, denoted as ℳ,

$$\begin{array}{c}\hfill \mathcal{M}:\mathbf{X}\to \{ϵ(n_B),P(n_B),\delta (n_B),v_{\beta }(n_B),v_{\mathrm{FR}}(n_B),\dots \},\end{array}$$(1)

which takes the values of 13 nuclear matter parameters X as input and outputs a β-equilibrated equation of state (EoS) and the composition of the entire star, including the crust. In practice, ℳ provides the β-equilibrated total energy density ϵ, pressure P, electron and muon fractions, and nuclear asymmetry δ (i.e. δ = 1 − 2x_p, where x_p is the proton fraction), all as functions of the baryon number density n_B. We refer the reader to Margueron et al. (2018), Mondal & Gulminelli (2022), and Davis et al. (2024) for an extensive presentation of the nuclear metamodel and its astrophysical applications. Here, we only note that we added the equilibrated v_β and frozen v_FR sound speeds to the metamodel output, which will be important in Sect. 4. Given non-informative priors on the 13 nuclear matter parameters X, the resulting metamodel realisation³ ℳ(X) undergoes a Bayesian filtering process that assigns a likelihood ℒ(X), which is detailed in the next section.

3. Likelihood of metamodel realisations

The metamodel instances ℳ(X) are not all equally realistic, in the sense that some give rise to, say, an EoS that is inconsistent with astrophysical observations, or they are not able to reproduce some experimental nuclear phenomenology. Therefore, we assigned a likelihood ℒ(X) to each ℳ(X) via a sequence of Bayesian filters, similar to the ones detailed in Dinh Thi et al. (2021a) and Davis et al. (2024). The reader may also consult Scurto et al. (2024), Char et al. (2023) and Malik et al. (2024) for a similar approach with the relativistic mean field. This filters are set out as follows.

1. The nuclear model ℳ(X) must be consistent with the energy per nucleon of pure neutron matter obtained by ab initio calculations employing chiral effective interactions (χ-EFT) and renormalisation group methods. The conflation of results in the literature obtained from different many-body methods results in an energy band (Huth et al. 2021), which is used to build an informed prior.

2. The nuclear model ℳ(X) must reproduce the nuclear mass measurements in the AME2020 mass table (Huang et al. 2021).

3. The β-equilibrated EoS obtained from ℳ(X) must support a maximum TOV mass greater than that of PSR J0348+0432, as measured by Antoniadis et al. (2013). Additionally, β-equilibrated matter must be stable and causal at least up to the central density of the star with the maximum TOV mass.

4. Similarly, we implement the constraint on the tidal deformability from the GW170817 event (Abbott et al. 2019); see Appendix B.

5. Finally, we implement constraints from the mass-radius X-ray pulse-profile estimates of the masses and radii of PSR J0030+0451, PSR J0437−4715, and PSR J0740+6620 (Vinciguerra et al. 2024; Choudhury et al. 2024; Salmi et al. 2024).

Compared to the previous works on this matter, the main difference lies in how we implement the causality constraint, which is part of point iii and is discussed later, and how we handle the information from χ-EFT calculations. In fact, we take care of point i by constructing an χ-EFT-informed prior via a Metropolis-Hastings sampling, as discussed in Sect. 3.1. Then, we randomly extract 10⁵ models ℳ(X) from this informed prior and pass them trough the sequence of Bayesian filters (ii-v), each of which is assigned a partial likelihood ℒ_i(X). The total likelihood of each metamodel instance ℳ(X) is

$$\begin{array}{c}\hfill \mathcal{L}(\mathbf{X})={\displaystyle \prod _j}{\mathcal{L}}_j(\mathbf{X})={\displaystyle \prod _j}p\left(D_j|\mathcal{M}(\mathbf{X})\right),\end{array}$$(2)

where p(D_j|ℳ(X)) is the conditional probability of reproducing the data, D_j, assuming the metamodel instance, M(X), and the index, j, runs over all the aforementioned constraints. Clearly, ℒ(X) is automatically also the likelihood of all the stellar properties (e.g. mass-radius relation; mode frequencies) that can be derived by assuming the matter model ℳ(X).

3.1. Informed prior from the χ-EFT band

Here, we discuss point i above in more detail. State-of-the-art χ-EFT calculations provide the energy per particle e(n)±δe(n) of pure neutron matter, where 0.02 < n < 0.2 fm⁻³ is the neutron density and δe(n) is the uncertainty associated with the specific calculation. Since different theoretical approaches yield different (overlapping) energy bands, e(n)±δe(n), we combine all the bands presented in Huth et al. (2021) into a single ‘conflated’ band, where the lower limit is given by the unitary gas approach (see Appendix A). This ensures that we do not underestimate the uncertainty associated with the theoretical calculations of e(n). Specifically, our conflated band is interpreted as a 90% confidence interval for e(n): for each n, the band is represented by a continuous probability density p(e|n) that is flat within the conflated band and has Gaussian tails accounting for the remaining 5%+5% (see Eq. (A.2)). This helps achieve a faster burn-in of the Metropolis-Hastings algorithm. Moreover, the Metropolis-Hastings procedure applied to p(e|n) allows us to to directly sample the nuclear parameters X for which e_X(n) obtained from ℳ(X) lies within our conflated band. This process starts with a flat prior⁴ for the nuclear parameters X. The resulting posterior is then used as an informed prior for filters (ii-v). This approach provides approximately 10⁹ nuclear models ℳ(X) in the informed prior, which is the most selective yet the least computationally demanding one.

3.2. Low-density filters from nuclear phenomenology

Each metamodel instance ℳ(X) can be used to calculate the mass M_N​Z(X) of a nucleus with N neutrons and Z protons. To do so, a compressible liquid drop model is used, supplemented by five extra surface and curvature parameters (Carreau et al. 2019b; Dinh Thi et al. 2021b). Therefore, to implement a filter (iii), we compare M_N​Z(X) with the measured nuclear masses $M_{NZ}^{\mathrm{AME}}$ listed in the AME2020 mass table (Wang et al. 2021). Following Dinh Thi et al. (2021a), we assign a partial likelihood of zero – that is, M(X) is discarded – if it is impossible to find values for the five curvature and surface parameters that are consistent with nuclear phenomenology. Otherwise, the partial likelihood is the goodness of the fit:

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{AME}}\propto exp\left(-{\mathit{\chi }\left(\mathbf{X}\right)}^2/2\right),\end{array}$$(3)

where the cost function χ² is

$$\begin{array}{c}\hfill {\chi }^2(\mathbf{X})=\frac{1}{N_{\mathrm{AME}}{\sigma }^2}{\displaystyle \sum _{NZ}}{(M_{NZ}(\mathbf{X})-M_{NZ}^{\mathrm{AME}})}^2.\end{array}$$(4)

Here, σ is a measure of the theoretical error on nuclear masses⁵, and the label N​Z runs over all the N_AME nuclei listed in the mass table.

The model distribution after applying filters (i-ii) yields a posterior distribution for the nuclear parameters X that is consistent with nuclear physics information up to the saturation density. At this stage, for every ℳ(X) that has not been excluded, we can extract the unified β-equilibrated EoS and all relevant outputs in Eq. (1) for all layers, including the crust.

3.3. High-density filters from astrophysics

Astrophysical constraints are applied through filters (iii-v), which are more sensitive to how ℳ(X) describes matter above the saturation density. The first check is hard, in the sense that the partial likelihood is either zero or one; the Tolman-Oppenheimer-Volkoff (TOV) equations are solved, and the maximum TOV mass m^*(X) is extracted. We assign a unit multiplicative contribution to the total likelihood in Eq. (2) for any model that satisfies causality and thermodynamic stability (i.e. 0 < v_β² < v_FR² < 1; see Camelio et al. 2023 for formal proof) and has a non-negative symmetry energy in the range of 0 < n_B < n^*(X), where n^*(X) is the central density of the star with mass m^*(X). Otherwise, the model’s likelihood ℒ(X) is set to zero (i.e. the instance ℳ(X) is discarded).

After this preliminary hard filter, we can go through the remaining filters (iii-v), which require the mass-radius relation and tidal deformability and are implemented as in Dinh Thi et al. (2021a), Scurto et al. (2024), Char et al. (2023) and Davis et al. (2024). We briefly list them below and refer to previous work for further details.

To implement filter iii, we require the maximum TOV mass m^*(X) to exceed the measured mass of PSR J0348+0432, M = 2.01 ± 0.04 M_⊙ (Antoniadis et al. 2013). The resulting contribution to the total likelihood is

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{J}0348}=\frac{1}{0.04\sqrt{2\pi }}{\displaystyle {\int }_0^{m^{\ast }(\mathbf{X})/M_{\odot }}}exp\left(-\frac{{(x-2.01)}^2}{2\times 0.{04}^2}\right)\mathrm{d}x.\end{array}$$(5)

Filter iv uses data from GW170817 and is based on the comparison between the effective dimensionless tidal deformability $\stackrel{\sim }{\mathrm{\Lambda }}$ calculated with ℳ(X) and the data of the Ligo-Virgo Collaboration (LVC). The likelihood takes the following form (see Appendix B for details):

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{LVC}}={\displaystyle {\int }_{0.73}^1}P(\stackrel{\sim }{\mathrm{\Lambda }},q)\mathrm{d}q,\end{array}$$(6)

where $\stackrel{\sim }{\mathrm{\Lambda }}$ is the effective tidal deformability, q < 1 is the mass ratio of the lighter object over the heavier one, and $P(\stackrel{\sim }{\mathrm{\Lambda }},q)$ is the observational joint posterior distribution reported in Abbott et al. (2019).

Finally, in filter v we check if the mass-radius relation R_X(m) obtained with the nuclear model ℳ(X) is consistent with the updated NICER estimates of the joint mass-radius distributions for three pulsars:

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{NICER}}={\displaystyle \prod _{i=1,2,3}{\int }_{0.7M_{\odot }}^{m^{\ast }(\mathbf{X})}}{P}_i(m,R_{\mathbf{X}}(m))\mathrm{d}m,\end{array}$$(7)

where P₁ is the joint probability distribution of mass and radius of the PSR J0030+0451 pulsar (Vinciguerra et al. 2024), P₂ refers to PSR J0437−4715 (Choudhury et al. 2024), and P₃ refers to PSR J0740+6620 (Salmi et al. 2024).

4. Frozen and equilibrated normal modes

The Bayesian procedure outlined in the previous section allows us to assign a likelihood to each model ℳ(X) based on its compatibility with nuclear physics phenomenology and astrophysical constraints. We now proceed to compute the normal mode frequencies, with the double aim of checking the impact of chemical transfusion and obtaining a posterior predictive distribution based on ℒ(X) for the mode frequencies.

4.1. Numerical scheme for the normal mode frequencies

First, we recall how to determine the frequencies of the f and p₁ normal modes of a spherically symmetric non-rotating NS in the relativistic Cowling approximation (McDermott et al. 1983). Since the space time remains unperturbed, no gravitational waves are emitted and, consequently, the radiation damping is absent. Moreover, the two equilibrated and frozen limits we considered are non-dissipative regimes (e.g. Gavassino et al. 2021, Sect. II-D), implying that there is no bulk viscosity damping due to reactions (e.g. Sawyer 1989; Haensel et al. 2002; Schmitt & Shternin 2018; Alford & Harris 2019; Alford et al. 2023). This limits our study to purely real frequencies.

Following Sotani et al. (2011), and consistently with the more complete full-GR derivation in Lindblom & Detweiler (1983), Sotani et al. (2001), the spherical space-time metric is

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=-e^{2\mathrm{\Phi }}\mathrm{d}{t}^2+e^{2\mathrm{\Lambda }}\mathrm{d}{r}^2+r^2\mathrm{d}{\theta }^2+r^{2}sin\theta \mathrm{d}{\varphi }^2,\end{array}$$(8)

while the Lagrangian fluid displacement (see Thorne & Campolattaro 1967) is a three-vector ξⁱ(t, r, θ, ϕ) defined with respect to the space-like part of the coordinate basis:

$$\begin{array}{c}\hfill {\xi }^i=r^{-2}(\stackrel{\sim }{W}{e}^{-\mathrm{\Lambda }}{Y}_{\mathit{lm}},-\stackrel{\sim }{V}{\partial }_{\theta }{Y}_{\mathit{lm}},-\stackrel{\sim }{V}{sin}^{-2}\theta {\partial }_{\varphi }{Y}_{\mathit{lm}}),\end{array}$$(9)

where $\stackrel{\sim }{W}(r,t)=W(r)e^{i\omega t}$ and $\stackrel{\sim }{V}(r,t)=V(r)e^{i\omega t}$ characterise the amplitude of the perturbation and Y_lm are the spherical harmonics, as in Lindblom & Detweiler (1983) and Sotani et al. (2001). Using these variables, the equations for the oscillation modes are

$$\begin{array}{cc}\hfill \frac{\mathrm{d}W}{\mathrm{d}r}& ={\left(\frac{\mathrm{d}P}{\mathrm{d}ϵ}\right)}^{-1}[{\omega }^2{r}^2{e}^{\mathrm{\Lambda }-2\mathrm{\Phi }}V+\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r}W]-l(l+1)e^{\mathrm{\Lambda }}V\hfill \\ \hfill \frac{\mathrm{d}V}{\mathrm{d}r}& =2\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r}V-\frac{1}{r^2}{e}^{\mathrm{\Lambda }}W,\hfill \end{array}$$(10)

where we take l = 2, since we focus on quadrupolar oscillations. The different regime of the balance between oscillation frequency and reaction rate is determined by the term dP/dϵ, the squared speed of sound, which encodes information equivalent to that in the adiabatic index (e.g. Haensel et al. 2002; Andersson 2019).

Boundary conditions at the star centre and surface are required to solve the system in Eq. (10). An inspection of the system shows that W(r) = Cr^l + 1 + 𝒪(r^l + 3) and V(r) = Cr^l + 𝒪(r^l + 2) for r → 0, with C being an arbitrary constant. The other boundary condition is obtained by demanding that the pressure perturbation vanishes at the stellar surface, which leads to

$$\begin{array}{c}\hfill {\omega }^2{r}^2{e}^{\mathrm{\Lambda }-2\mathrm{\Phi }}V+\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r}W=0.\end{array}$$(11)

With this condition, the problem becomes an eigenvalue one, which we solve using a standard shooting method. After determining the metric functions and stellar structure by solving the TOV equations, we solve the system in Eq. (10) using an initial guess for the (purely real and positive) pulsation ω. We then refine the value of ω with a bisection method, iterating the process until we find the exact ω that satisfies Eq. (11).

4.2. Frozen speed of sound and the thermodynamic stability-causality condition

For perturbations that are slow enough, matter will be always almost in β equilibrium, and dP/dϵ in Eq. (10) can be taken to be v_β², the sound speed arising from a purely barotropic EoS:

$$\begin{array}{c}\hfill \frac{\mathrm{d}P}{\mathrm{d}ϵ}=\frac{\mathrm{d}P(n_B,\delta (n_B))/\mathrm{d}{n}_B}{\mathrm{d}ϵ(n_B,\delta (n_B))/\mathrm{d}{n}_B}=v_{\beta }^2(n_B).\end{array}$$(12)

On the other hand, in a fast oscillation regime the composition of each fluid element has no time to relax and return to chemical equilibrium, and one should use the sound speed v_FR² at frozen composition, which we conveniently write as

$$\begin{array}{c}\hfill \frac{\partial P}{\partial ϵ}{|}_{\delta (n_B)}=\frac{\partial P(n_B,\delta (n_B))/\partial n_B}{\partial ϵ(n_B,\delta (n_B))/\partial n_B}=v_{\mathrm{FR}}^2(n_B).\end{array}$$(13)

As discussed in Camelio et al. (2023), the velocity v_FR coincides with the ‘maximal characteristic speed’ (the speed defining the Courant-Friedrichs-Lewy condition) for a signal propagating in a chemically reacting fluid mixture, implying that the mixture is both thermodynamically stable and causal⁶ only if

$$\begin{array}{c}\hfill 0<v_{\beta }^2(n_B)<v_{\mathrm{FR}}^2(n_B)<1.\end{array}$$(14)

From the point of view of global oscillations, particularly in g modes (e.g. Lai 1994; Jaikumar et al. 2021; Tran et al. 2023), the same criterion guarantees the local convective stability of the star (see Eq. (A.12) of Camelio et al. 2023 with Eq. (4.17) in Lai 1994). Therefore, as mentioned in Sect. 3.3, we only retain the metamodel instances ℳ(X) that satisfy the fundamental thermodynamic stability-causality condition 0 < v_β²(n_B) < v_FR²(n_B) < 1 at least up to the central density of the NS with maximum TOV mass.

5. Results and discussion

For each nuclear model, ℳ(X), we extracted the f and p₁ normal mode frequencies in the Cowling approximation, with the purpose of testing the QU relations with a large set of metamodel instances and quantifying the potential impact of assuming frozen or equilibrated composition. Finally, we used the known QU relation in full GR to estimate a more realistic posterior predictive distribution for the mode frequencies.

5.1. Differences between frozen and barotropic frequencies

We evaluated the Cowling frequencies of the f and p₁ modes in the two ideal limits of frozen and equilibrated composition, as outlined in the previous section. The results in the β-equilibrated case are shown in Fig. 1, where the prediction of the model associated with the highest likelihood are given by solid lines. We can see that, though the two modes are clearly separated, accounting for the uncertainty in the nucleonic model leads to an important dispersion of the predictions, particularly for the p₁ case. As a consequence, the discrimination between hadronic and strange stars from the measured value of the p₁ frequency might be harder than expected in first works that only considered a limited set of hadronic models (e.g. Vásquez Flores & Lugones 2014).

Fig. 1. Probability density distributions for f-mode frequencies (left) and p1-mode frequencies (right), both obtained within the Cowling approximation in the barotropic limit. The three shaded regions refer to the 68%, 95%, and 99% percentiles. The solid black line represents the model with the highest likelihood.

The frequencies of the f mode are almost unaffected by the equilibration assumption, with differences smaller than 0.5%, as shown in Fig. 2. On the other hand, the p₁ mode exhibits a more interesting behaviour, where the differences between the two cases are more evident and tend to increase with mass, as can be seen in Fig. 3. However, for the models with high likelihood, the ones in the darkest region of the plot, the frequencies calculated in the frozen limit remain close to the ones obtained by assuming the barotropic sound speed.

Fig. 2. Posterior probability density for frequencies of frozen f-mode (upper panel). The lower panel shows the relative difference between the frequencies in the frozen limit and the barotropic limit. The shaded regions represent the 68%, 95%, and 99% percentiles, respectively, while the black line indicates the model with the highest likelihood.

Fig. 3. Posterior probability density for frequencies of frozen p1 mode (upper panel). The lower panel shows the relative difference between the frequencies in the frozen limit and the barotropic limit. The three shaded regions correspond to the 68%, 95%, and 99% percentiles, while the black line represents the model with the highest likelihood.

Based on these results, only the frozen frequencies are presented in the subsequent discussion, as the differences are negligible for the f mode and less than 5% for the p₁ mode in reasonable mass ranges (not too close to the maximum TOV mass). For the same reason, the present analysis confirms – on the basis of a large set of nuclear models – that when the frozen speed of sound or the frozen adiabatic index is unavailable, the β-equilibrated speed of sound v_β can be used with minimal error. Namely, f and p₁ modes obtained with agnostic barotropic models can be trusted within 5% or better, especially for masses below ∼2 M_⊙.

5.2. Test of proposed quasi-universal relations

Andersson & Kokkotas (1998) proposed a QU relation for the for the p₁ mode, where the mode pulsation ω times the NS mass M was expressed in terms of the compactness M/R. This same scaling was later used to look for a QU relation for the f mode by Tsui & Leung (2005):

$$\begin{array}{c}\hfill \omega M(\mathrm{rad}/\mathrm{s}\mathrm{km})=a_3{\left(\frac{M}{R}\right)}^3+a_2{\left(\frac{M}{R}\right)}^2+a_1\frac{M}{R}+a_0,\end{array}$$(15)

where the coefficients a_i are obtained from a fit over a limited number of barotropic EoS models. This empirical expression was recently tested with a₂ = a₃ = 0 for the f mode in the Cowling approximation (Pradhan & Chatterjee 2021) and with a₃ = 0 for the f mode in full GR (Pradhan et al. 2022). Moreover, Sotani (2021) applied the empirical relation (15), including all coefficients, to the p₁-mode frequencies in full GR (see Table 1).

As a preliminary check, we evaluated the accuracy of the QU relation (15) using our set of metamodel instances in order to assess the quality of the proposed fits and the dispersion of the metamodel instances around them. The results are presented in Fig. 4, where we compare our findings with the fits from the aforementioned works (the coefficients of these fits are listed in Table 1, along with a fit of the p₁ mode to our numerical results). The first panel of Fig. 4 shows the density map of mode frequencies resulting from our Bayesian filtering, with the f mode in blue-green and the p₁ mode in orange-yellow, alongside the various QU relations mentioned earlier.

Fig. 4. Relationship between rescaled pulsations ωM and compactness M/R for the frozen f modes (blue-green) and frozen p1 modes (orange-yellow). The dashed lines correspond to the fits presented in Table 1, while the solid black line represents the model ℳ(X) with the highest likelihood ℒ(X). The fit lines overlapping with the distributions are based on the Cowling approximation, whereas the others, obtained in full GR, are shown for comparison. The differences between the rescaled pulsations ωM calculated in the Cowling approximation and the corresponding QU fits are illustrated in the two lower panels. Each panel also includes three shaded regions representing the 68%, 95%, and 99% quantiles of the distribution.

To quantify the dispersion of the metamodel instances around the QU fits, the two lower panels of Fig. 4 display the differences between our numerical results and the Cowling QU fitting formula. In both lower panels of Fig. 4, the dispersion around the proposed f-mode QU fit is minimal, demonstrating that our extensive set of EoSs adheres to it with the expected level of precision, with errors smaller than 2.5%. However, a structure in the residuals remains visible, which can be attributed to the choice of a linear fit. In contrast, the functional form for the fit of the p₁-mode QU seems appropriate, as there are no evident underlying structures observed in the dispersion of the residuals. Nevertheless, it is noteworthy that in this case the precision to which the QU relation is realised is lower, with errors ranging from approximately 5% to 10%.

For completeness, we also tested an alternative empirical relation, linking the f-mode frequency and the average density of the star (Andersson & Kokkotas 1998; Pradhan et al. 2022):

$$\begin{array}{c}\hfill f=a+b{\left(\frac{\overline{M}}{{\overline{R}}^3}\right)}^{1/2},\end{array}$$(16)

where $\overline{M}=M/1.4M_{\odot }$, $\overline{R}=R/10$ km and the constants a and b are obtained from a fit to the numerical results. This relation has been tested by different authors with different barotropic EoSs (not all compatible with the constraint imposed by the measured mass of PSR J0348+0432). Therefore, we verify whether Eq. (16) is a QU relation by using our filtered set of nuclear models. This is shown in Fig. 5, where the upper distribution represents our Cowling results, which is compared to the one obtained by reversing relation Eq. (15) with the coefficients provided in Pradhan et al. (2022). We also compare these distributions to the empirical relations presented in Pradhan et al. (2022), Benhar et al. (2004), and Andersson & Kokkotas (1998). Since these empirical relations are all derived from fits to frequencies extracted in full GR, they are obviously not compatible with our Cowling results. In contrast, the Cowling relation presented in Pradhan & Chatterjee (2021), obtained within the Cowling approximation, is closer to our results. It can be observed that this relation strongly depends on the selected set of EoSs, resulting in a significant spread around the relation in Eq. (16).

Fig. 5. Distributions of f-mode frequencies as a function of the average density $\overline{M}/{\overline{R}}^3$. The upper colour map refers to the frequencies we obtain within the Cowling approximation, while the lower distribution is obtained by reversing the relation in Eq. (15) with the coefficients ai from Pradhan et al. (2022). For comparison, the lines are the several fits listed in Tab-II of Pradhan et al. (2022). The three shaded regions contain the 68%, 95%, and 99% quantiles of the distribution.

5.3. Estimation of full-GR mode frequencies

Because of the excellent agreement between Eq. (15) and the metamodel result in the Cowling approximation, we can assume that the dispersion observed in Fig. 4, due to the different softness levels of the nuclear models, will equally affect the degree of validity of the QU relations in full GR. Under this assumption, the QU relation in Eq. (15) can be used to quickly estimate the frequencies for ℳ(X) in full GR directly from the R_X(M) relation, as long as the opportune parameters a_i are used. We denote these frequencies as ‘synthetic’ since they are not obtained by solving the eigenvalue problem but rather by simply unpacking the QU relation Eq. (15) via the mass-radius relation R_X(M) of each ℳ(X).

More precisely, the procedure used to recover the synthetic frequencies f_GR in full GR (i.e. beyond the Cowling approximation) is

$$\begin{array}{c}\hfill f_{\mathrm{GR}}(M)=\frac{1+\mathrm{\Delta }(M)}{2\pi M}{U}_{\mathrm{GR}}(M),\end{array}$$(17)

with Δ given by

$$\begin{array}{c}\hfill \mathrm{\Delta }(M)=\frac{M{\omega }_{\mathrm{C}}(M)-U_{\mathrm{C}}(M)}{U_{\mathrm{C}}(M)},\end{array}$$(18)

where ω_C is the mode pulsation that we found within the Cowling approximation in the frozen limit; U_C, GR is the Cowling (C) or full-GR (GR) quasi-universal relation, namely the right hand side of Eq. (15) with the appropriate coefficients a_i listed in Table 1.

The prescription Eq. (17) for the synthetic frequencies f_GR is designed so that we do not underestimate the spread of the p₁ frequencies, as discussed in Appendix C. Essentially, we unpack the QU relation in Eq. (15), with the coefficients a_i extracted from numerical results in full GR, and transporting the spread of our Cowling calculation onto the unpacked results.

Figure 6 shows the estimated probability density of the synthetic frequencies for the f mode and p₁ mode, together with the prediction of the model associated with the highest likelihood. As expected, the f-mode frequency increases more rapidly with mass than the p₁-mode one, which remains relatively flat. Consequently, extrapolating NS features from the p₁-mode frequencies is expected to be much more challenging. To quantify this further, in Fig. 7 we show the posterior distributions of frequencies for an NS with masses of M = 1 M_⊙, 1.4 M_⊙, 2 M_⊙ in the frozen limit. The three distributions for the p₁ mode are nearly indistinguishable, as they almost completely overlap. In contrast, the three distributions for the f mode show only partial overlap, suggesting that it may be possible to constrain the mass of an NS, despite uncertainties in the nuclear EoS. On the other side, the quasi-universality of the p₁ frequency in the purely hadronic hypothesis opens the compelling possibility of it being challenged in hybrid or strange stars, as proposed in Vásquez Flores & Lugones (2014) and Wilson & Ho (2024).

Fig. 6. Probability distributions of synthetic full-GR frequencies obtained from Eq. (17). The left panel shows the f mode, while the right one shows the p1 mode. The three shaded regions in each panel contain the 68%, 95%, and 99% percentiles. The solid black line represents the model with the highest likelihood.

Fig. 7. Distribution of synthetic full-GR frequencies at 1 M⊙, 1.4 M⊙, and 2 M⊙ obtained from Eq. (17). The left panel refers to the f mode, while the right one refers to the p1 mode.

6. Conclusions

With the advent of next-generation interferometers, it becomes important to evaluate how future GW detection from oscillating NSs could be used to constrain nuclear models of neutron star interiors or infer the mass of the object. To address this, we adopted the nuclear metamodel framework (Margueron et al. 2018) for cold npeμ matter, generating a large set of unified equations of state, together with their β-equilibrated composition and the two (barotropic and frozen) sound speeds. These nuclear models were then assigned likelihoods ℒ(X) through a sequence of Bayesian filters, designed to weigh each metamodel instance based on its consistency with established nuclear and astrophysical phenomenology. Given this posterior for the EoSs, we find the posterior predictive distributions for the f- and p₁-mode frequencies – shown in Fig. 7 – by inverting known full-GR quasi-universal relations. More precisely, Fig. 7 is our ‘synthetic’ full-GR prediction of the f- and p₁-mode frequencies as a function of the NS’s mass: while an f-mode detection could constrain the NS mass, this information is almost completely lost for the p₁ mode.

The generation of a large set of metamodel instances and the relative stellar structures for different masses also allowed us to check another point, which is more related to the physical assumptions underlying the computation of the modes. Proposed QU relations pertaining to mode frequencies have been found using barotropic models or, equivalently, non-barotropic nuclear models that are always at strict β equilibrium, namely using dP/dϵ = v_β² in Eq. (10). Hence, we checked the impact of the more realistic (Haensel et al. 2002) frozen limit assumption, dP/dϵ = v_FR², to see if it could introduce any deviation from the known QU relations for the f and p₁ modes. This check is a first, albeit partial, step towards a more systematic study of the impact of nuclear reactions on NS oscillation spectra (see e.g. Counsell et al. 2024a); in principle, reactions introduce mode damping, whose strength depends on the details of the nuclear model and physical conditions of temperature and density (e.g. Haensel et al. 2002; Schmitt & Shternin 2018; Alford & Harris 2019; Alford et al. 2023, 2024). However, in the two ideal limits considered here, any possible bulk-viscous effect is exactly zero (Gavassino et al. 2021; Camelio et al. 2023). This is a caveat to be kept in mind.

Our analysis shows that both the f and p₁ modes do not significantly depend on whether the sound speed used is the barotropic or frozen one. This is in contrast with what is known for g modes, where both velocities have to be used to find the frequency spectrum (e.g. Reisenegger & Goldreich 1992; Tran et al. 2023; Zhao & Lattimer 2022; Counsell et al. 2024b). Therefore, we conclude that studies assuming purely barotropic agnostic models for the EoS are accurate to within a few percent. This behaviour is reflected in the goodness of the QU relation, which can thus be used to estimate the mode frequencies without solving the perturbation equations – a crucial advantage in Bayesian studies that involve millions of agnostic EoSs.

Finally, the posterior set of metamodel EoSs obtained through the filtering procedure represents a refinement over previous similar studies (Dinh Thi et al. 2021a; Davis et al. 2024), owing to the implementation of the more stringent stability-causality condition 0 < v_β² < v_FR² < 1 derived in Camelio et al. (2023). This improved posterior set may also serve as a useful input for further studies on potential constraints on NS interiors, such as those derived from pulsar glitches (Antonelli et al. 2022).

Acknowledgments

We thank Hoa Dinh Thi, Micaela Oertel, Chiranjib Mondal and Debarati Chatterjee for useful discussion, and Philip John Davis for technical support. Partial support comes from the IN2P3 Master Project NewMAC, the ANR project “Gravitational waves from hot neutron stars and properties of ultra-dense matter” (GW-HNS, ANR-22-CE31-0001-01), and the CNRS International Research Project (IRP) “Origine des éléments lourds dans l’univers: Astres Compacts et Nucléosynthèse” (ACNu).

References

Appendix A: The chiral band of neutron matter

The χ-EFT ab-initio calculations taken into account in this work are presented in Fig. 1 of Huth et al. (2021): the i-th approach provides an estimate of e_i(n)±δe_i(n), the energy per baryon of pure neutron matter in the range 0.02 < n < 0.2fm⁻³ where all approaches are expected to provide reliable results. For each metamodel instance ℳ(X), we can easily extract e_X(n), and compare it with the theoretical microscopic results e_i(n)±δe_i(n). To this end, we have to conflate all the bands e_i(n)±δe_i(n) reported in Huth et al. (2021) into a single one, e(n)±δe(n): the lower limit e₋(n) = e(n)−δe(n) is given by the unitary gas model, while the upper bound e₊(n) = e(n)−δe(n) is

$$\begin{array}{c}\hfill e_{+}(n)={\displaystyle \underset{i}{max}}\{e_i(n)+\delta e_i(n)\}.\end{array}$$(A.1)

In order not to underestimate the theoretical systematic error, we interpret [e₋(n),e₊(n)] as the 90% confidence interval where e_X(n) should lie arising from a smooth probability distribution to be used within the metropolis-Hastings algorithm. Namely, we consider the following normalized distribution:

$$\begin{array}{c}\hfill p(e|n)=Q_n\{\begin{array}{cc}exp(-\frac{{(e-e_{-}(n))}^2}{2{\sigma }_n^2})\hfill & \mathrm{if}e\in (-\infty ,e_{-}(n)]\hfill \\ 1\hfill & \mathrm{if}e\in (e_{-}(n),e_{+}(n)]\hfill \\ exp(-\frac{{(e-e_{+}(n))}^2}{2{\sigma }_n^2})\hfill & \mathrm{if}e\in (e_{+}(n),\infty )\hfill \end{array}\end{array}$$(A.2)

where

$$\begin{array}{c}\hfill Q_n=\frac{0.9}{e_{+}(n)-e_{-}(n)}{\sigma }_n=\frac{e_{+}(n)-e_{-}(n)}{9\sqrt{2\pi }}\end{array}$$(A.3)

In this way, the central plateau of the distribution accounts for the 90% while each tail for the remaining 10%, in accordance with other prescriptions used previous studies (Dinh Thi et al. 2021a; Carreau et al. 2019b; Scurto et al. 2024). Then, the partial likelihood of ℳ(X) is given by the geometric product integral

$$\begin{array}{c}\hfill \mathcal{L}{\chi }_{\mathit{EFT}}(\mathbf{X})=exp{\displaystyle \int d}nlogp(e|n)\end{array}$$(A.4)

over the 0.02 < n < 0.2 fm⁻³ density range. In practice, the density range is divided in N equally spaced slices at densities n_j and the resulting likelihood is

$$\begin{array}{c}\hfill \mathcal{L}{\chi }_{\mathit{EFT}}(\mathbf{X})\propto {\displaystyle \prod _{j=1}^N}p\left(e_{\mathbf{X}}(n_j)|n_j\right).\end{array}$$(A.5)

Appendix B: Scheme for the LVC constraint

For completeness, we provide a schematic presentation of the LVC constraint in Eq. (6). This may help the reader sort the details and complement the sketch given in previous works that adopt the same prescription (Dinh Thi et al. 2021a; Mondal & Gulminelli 2022; Scurto et al. 2024; Char et al. 2023; Davis et al. 2024).

The analysis in Abbott et al. (2019) provides the observational joint posterior $P(\stackrel{\sim }{\mathrm{\Lambda }},q)$ for the effective tidal deformability $\stackrel{\sim }{\mathrm{\Lambda }}$ and the mass ratio q of GW170817. In principle, both quantities can be determined from the masses m_k and tidal deformabilities Λ_k of the two NSs (k = 1, 2) using known analytical expressions, $\stackrel{\sim }{\mathrm{\Lambda }}(m_k,{\mathrm{\Lambda }}_k)$ and q(m_k). The GW170817 data enabled a relatively precise determination of the chirp mass m_c (treated as a given constant in the following), which can also be expressed analytically in terms of the two masses, m_c(m_k).

To implement the constraint imposed by knowledge of $P(\stackrel{\sim }{\mathrm{\Lambda }},q)$ and m_c, the first step is to recognize that we are adopting a framework where the mass M is treated as an independent variable, and ℳ(X) can be used to obtain the relations R_X(M) and Λ_X(M). This is a natural and convenient choice, considering that R_X(M) and Λ_X(M) are genuine functions, whereas M_X(R) or Λ_X(R) can be multivalued.

Now, the observational information we have is $P(\stackrel{\sim }{\mathrm{\Lambda }},q)$ and the value of m_c, but both q and m_c depend only on the masses that, in our framework, carry no dependence on ℳ(X). Therefore, the nuclear model dependence can only enter via $\stackrel{\sim }{\mathrm{\Lambda }}$, leaving us with the possibility of marginalising over q. For any given instance ℳ(X):

1. From m_c and q we find m_k(m_c, q) for the two NSs, k = 1, 2.

2. We can use the model-specific relation Λ_X(M): the two tidal deformabilities are Λ_k = Λ_X(m_k).

3. At this point we can compute ${\stackrel{\sim }{\mathrm{\Lambda }}}_{\mathbf{X}}(q,m_c)=\stackrel{\sim }{\mathrm{\Lambda }}(m_k,{\mathrm{\Lambda }}_k)$, where the dependence on ℳ(X) enters via Λ_k. The arguments of ${\stackrel{\sim }{\mathrm{\Lambda }}}_{\mathbf{X}}$ are q and m_c because of step (i).

4. The likelihood is given by the marginalization over q, namely ${\mathcal{L}}_{LVC}(\mathbf{X})\propto {\displaystyle \int d}qP({\stackrel{\sim }{\mathrm{\Lambda }}}_{\mathbf{X}}(q,m_c),q)$ over the whole range of possible q values. This is exactly the prescription in (6).

Appendix C: Testing the prescription for the synthetic full GR frequencies

Given a QU relation for the mode frequencies, as the one in (15), it is possible to obtain the mode pulsation ω(M) simply by using the mass-radius relation R_X(M) of each nuclear model ℳ(X). However, such a method will give the exact ω(M) if and only if the QU relation is exact, that is in the limit of negligible dispersion of the model predictions around the QU line. Since this is clearly the case for the f-mode (see Fig. 4), we have followed this strategy to produce the lower density colour map of Fig. 5. The validity of the procedure is shown by the fact that the colour map overlaps with the fit, and also the spread of the frequencies is close to the one obtained by Pradhan & Chatterjee (2021).

Fig. C.1. Posterior distribution of the frequencies at 1 M⊙, 1.4 M⊙ and 2 M⊙ for the p1-mode. The left panel refers to the prescription described in Sec. 5.3 while the right is obtained by inserting the RX(M) relation of each model in (15) with the QU coefficients ai of Sotani (2021). The frequencies in the right panel have a narrower distribution, meaning that the prescription in (17) is necessary to make the spread of our synthetic full GR frequencies similar to the one found by Sotani (2021).

However, if the relation is only quasi-universal, as manifestly it is the case for the p₁-mode shown in Fig. 4, the inversion method will lead to an underestimation of the dispersion of the predictions, and the strategy proposed in Sec. 5.3 should instead be adopted. To check the validity of this statement, we have estimated the distribution of the p₁ frequencies in full GR by simply injecting the mass-radius relation R_X(M) of the models into the QU relation obtained by Sotani (2021).

The p₁-mode distribution obtained in this way – i.e., by “unpacking” the QU relation with R_X(M) – is shown in the right panel of Fig. C.1. For each mass, this distribution is narrower than the original spread between the frequencies for different EoS found by Sotani (2021). On the other hand, when we use (17) to transfer the dispersion around the QU relation obtained in Cowling to the full GR prediction (as done in the left panel of Fig. C.1), we qualitatively recover the same spread of frequencies reported in Sotani (2021), after we remove the Shen EoS (Shen et al. 1998) used therein. This EoS is particularly soft and not compatible with the 2.01 M_⊙ observation, consequently it cannot be reproduced by our data.

All Tables

All Figures

	Fig. 1. Probability density distributions for f-mode frequencies (left) and p1-mode frequencies (right), both obtained within the Cowling approximation in the barotropic limit. The three shaded regions refer to the 68%, 95%, and 99% percentiles. The solid black line represents the model with the highest likelihood.
In the text

	Fig. 2. Posterior probability density for frequencies of frozen f-mode (upper panel). The lower panel shows the relative difference between the frequencies in the frozen limit and the barotropic limit. The shaded regions represent the 68%, 95%, and 99% percentiles, respectively, while the black line indicates the model with the highest likelihood.
In the text

	Fig. 3. Posterior probability density for frequencies of frozen p1 mode (upper panel). The lower panel shows the relative difference between the frequencies in the frozen limit and the barotropic limit. The three shaded regions correspond to the 68%, 95%, and 99% percentiles, while the black line represents the model with the highest likelihood.
In the text

Fig. 4. Relationship between rescaled pulsations ωM and compactness M/R for the frozen f modes (blue-green) and frozen p1 modes (orange-yellow). The dashed lines correspond to the fits presented in Table 1, while the solid black line represents the model ℳ(X) with the highest likelihood ℒ(X). The fit lines overlapping with the distributions are based on the Cowling approximation, whereas the others, obtained in full GR, are shown for comparison. The differences between the rescaled pulsations ωM calculated in the Cowling approximation and the corresponding QU fits are illustrated in the two lower panels. Each panel also includes three shaded regions representing the 68%, 95%, and 99% quantiles of the distribution.

In the text

Fig. 5. Distributions of f-mode frequencies as a function of the average density $\overline{M}/{\overline{R}}^3$. The upper colour map refers to the frequencies we obtain within the Cowling approximation, while the lower distribution is obtained by reversing the relation in Eq. (15) with the coefficients ai from Pradhan et al. (2022). For comparison, the lines are the several fits listed in Tab-II of Pradhan et al. (2022). The three shaded regions contain the 68%, 95%, and 99% quantiles of the distribution.

In the text

	Fig. 6. Probability distributions of synthetic full-GR frequencies obtained from Eq. (17). The left panel shows the f mode, while the right one shows the p1 mode. The three shaded regions in each panel contain the 68%, 95%, and 99% percentiles. The solid black line represents the model with the highest likelihood.
In the text

	Fig. 7. Distribution of synthetic full-GR frequencies at 1 M⊙, 1.4 M⊙, and 2 M⊙ obtained from Eq. (17). The left panel refers to the f mode, while the right one refers to the p1 mode.
In the text

Fig. C.1. Posterior distribution of the frequencies at 1 M⊙, 1.4 M⊙ and 2 M⊙ for the p1-mode. The left panel refers to the prescription described in Sec. 5.3 while the right is obtained by inserting the RX(M) relation of each model in (15) with the QU coefficients ai of Sotani (2021). The frequencies in the right panel have a narrower distribution, meaning that the prescription in (17) is necessary to make the spread of our synthetic full GR frequencies similar to the one found by Sotani (2021).

In the text