Quasi-universal behaviour of the threshold mass in unequal-mass, spinning binary neutron-star mergers

The lifetime of the remnant produced by the merger of two neutron stars can provide a wealth of information on the equation of state of nuclear matter and on the processes leading to the electromagnetic counterpart. Hence, it is essential to determine when this lifetime is the shortest, corresponding to when the remnant has a mass equal to the threshold mass, $M_{\rm th}$, to prompt collapse to a black hole. We report on the results of more than 360 simulations of merging neutron-star binaries covering 40 different configurations differing in mass ratio and spin of the primary. Using this data, we have derived a quasi-universal relation for $M_{\rm th}$ and expressed its dependence on the mass ratio and spin of the binary. The new expression recovers the results of Koeppel et al for equal-mass, irrotational binaries and reveals that $M_{\rm th}$ can increase (decrease) by $5\%~(10\%)$ for binaries that have spins aligned (anti-aligned) with the orbital angular momentum and provides evidence for a non-monotonic dependence of $M_{\rm th}$ on the mass asymmetry in the system. Finally, we extend to unequal-masses and spinning binaries the lower limits that can be set on the stellar radii once a neutron-star binary is detected, illustrating how the merger of an unequal-mass, rapidly spinning binary can significantly constrain the allowed values of the stellar radii.


INTRODUCTION
With the birth of multi-messenger gravitational-wave astronomy and the detection of binary neutron star (BNS) mergers via the the GW170817 event The LIGO Scientific Collaboration & The Virgo Collaboration (2017), and the use of numerical simulations, it has become possible to establish a number of constraints on the equation of state (EOS) of nuclear matter (see, e.g., Margalit & Metzger 2017;Bauswein et al. 2017;Shibata et al. 2017;Rezzolla et al. 2018;Ruiz et al. 2018;Radice et al. 2018;Montaña et al. 2019;Raithel et al. 2018;Tews et al. 2018;Malik et al. 2018). Due to the degeneracy between tidal and stellar spin effects on the inspiral waveform, low-or high-spin priors are needed to infer the properties of the binaries from the gravitational-wave measurements (The LIGO Scientific Collaboration & The Virgo Collaboration 2017). Strong emphasis has been given to the low-spin priors under the expectation that the neutron stars (NSs) lose a significant fraction of their spin angular momentum through electromagnetic dipolar radiation well before the merger. While these assumptions favor irrotational tootle@itp.uni-frankfurt.de NSs in the late inspiral and merger of BNSs, the parameter estimations from the gravitational-wave detections still strongly depend on the given a-priori distribution of expected spins, with the high-spin prior leading to large uncertainties in the mass ratio and effective spin of the binary (Abbott et al. 2020).
Much of the theoretical modelling of BNS mergers has been concentrated on mass ratios q := M 2 /M 1 0.7 and on irrotational constituents (Shibata & Taniguchi 2006;Rezzolla et al. 2010;Bauswein et al. 2013;Dietrich et al. 2015;Radice et al. 2016;Lehner et al. 2016, see, e.g., ), in accordance with the (limited) sample of observed binary pulsar systems. However, the modelling of smaller ratios q 0.45 (Dietrich et al. 2017b;Most et al. 2020a;) and of higher spins has also started (Kastaun et al. 2013;Bernuzzi et al. 2014;East et al. 2016;Dietrich et al. 2017a;Most et al. 2019c). Given its impact on the electromagnetic counterpart, a particularly important prediction of the theoretical modelling has been the determination of whether the BNS underwent a prompt collapse at merger. The threshold mass discerning a prompt from a delayed collapse has been investigated thoroughly for irrotational binaries (see Baiotti & Rezzolla 2017;Burns 2020, for some reviews) using a number of EOSs. The natural expectation that M th can be parametrized in terms of the maximum mass of a nonrotating NS, M TOV (Bauswein et al. 2013), has been refined by more advanced parametrizations (Koeppel et al. 2019;Agathos et al. 2019) and by incorporating the effect of asymmetric binary systems (Bauswein et al. 2021). The latter resulted in a lowering of M th at small mass ratios, i.e., q 0.6, depending on the stiffness of the EOS (Bauswein et al. 2021). Furthermore, simulations have suggested that the lifetime of merger remnants increases with the binary spin (Kastaun et al. 2013;Bernuzzi et al. 2014;East et al. 2016;Kiuchi et al. 2019), and decreases significantly at very small mass ratios (Rezzolla et al. 2010;Dietrich et al. 2017a;Kiuchi et al. 2019;Bernuzzi et al. 2020).
Using three temperature-dependent EOSs compatible with present astronomical observations, we report on a systematic exploration the effects that spin and mass asymmetry have on M th of BNS configurations with mass ratio in the range 0.5 ≤ q ≤ 1, concentrating mostly on configurations in which the primary is spinning, while the secondary is irrotational. The large amount of data collected in this way has allowed us to derive a quasi-universal relation for M th in terms of the mass ratio and spin of the binary.

NUMERICAL METHODS AND INITIAL DATA
Numerical Methods. A significant hurdle to exploring asymmetric spinning binaries is the construction of initial data (ID) consistent with the Einstein equations. Here, we have used FUKA , the first publicly available ID solver able to reliably explore the needed parameter space. See the Appendix Materials (AM) for details on the configurations considered. The evolution of the ID is performed using the Einstein Toolkit infrastructure which includes a fixed-mesh box-in-box refinement framework Carpet Schnetter et al. (2004). The spacetime evolution was handled by Antelope (Most et al. 2019b) using a constraint damping formulation of the Z4 system (Bernuzzi & Hilditch 2010;Alic et al. 2012). Finally, the generalrelativistic magnetohydrodynamic code, FIL (Most et al. 2019b), was used to evolve the fluid quantities. FIL is a derivative of the original IllinoisGRMHD code (Etienne et al. 2015) with the addition of high-order (fourth) conservative finite-differencing methods (Del Zanna et al. 2007). Importantly, FIL handles EOSs that are dependent on temperature and electron-fraction, and includes a neutrino leakage scheme to handle neutrino cooling and weak interactions. The simulations were performed using a grid setup with an extent of ≈ 3000 km, decomposed in six levels of refinement, with the finest grid spacing being ∆x ≈ 295 m. For comparison, 36 simulations of higher-resolution were performed with a finer grid spacing of ∆x ≈ 221 m, resulting in a measure of the error budget of ∼ 1% (∆M th 0.03 M ).

Space of parameters.
Pulsar observations have provided a wealth of information on the various rotation states that NSs can have including isolated pulsars with extreme spins and binary configurations with moderate to high spins. In general, most of the spin angular momentum of the binary is associated with the most massive of the two NSs, i.e., the primary (Lazarus et al. 2016). Hence, the spin configurations considered here are mostly centered on the effects χ 1 has on M th . Here, 36 out of the 40 binaries considered have χ 1 = [−0.3, 0, 0.3] and χ 2 = 0 where χ 1,2 := J 1,2 /M 2 1,2 are, respectively, the dimensionless spins of the primary and secondary with spin angular momenta J 1,2 . However, in order to gauge the impact of the spin of the secondary we have considered also four binary configurations with the same effective dimensionless spinχ := (1 + q) −1 (χ 1 + q χ 2 ). Collapse-Time Measurement. A crucial aspect of any study wishing to determine M th is a rigorous, systematic and reproducible definition of what constitutes a collapse. Such an approach was previously not always considered, and very qualitative definitions of the threshold mass have been employed in the literature. Here, instead, we follow the prescription proposed by Koeppel et al. (2019), which tracks the minimum of the lapse to compute the collapse time and compares it to the shortest possible over which a prompt collapse can take place, that is, the free-fall timescale. This approach, whose details can be found in the AM, allows not only for a precise definition and measurements, but also for the reproducibility of the results presented here.

RESULTS
Dependence on Mass Ratio and Spin. While the work of Koeppel et al. (2019) provided a first well-defined and rigorous manner of determining a quasi-universal relation for M th , it was restricted to the analysis of equal mass, i.e., with q = 1 and irrotational binaries, i.e., with χ 1 = 0 = χ 2 . In this specific scenario, M th was found to depend on the EOS rather simply, so that a quasi-universal relation M th = M th (EOS) was proposed. However, it is clear that when considering the additional influence of spin and mass asymmetry, the functional dependence of M th must account also for these additional degrees of freedom so that M th = M th (EOS, q, χ), where χ := χ 1 + χ 2 is the total dimensionless spin of the binary.
Determining the exact expression for M th (EOS, q, χ) clearly requires the exploration of the space of parameters in mass ratio and spin for any given EOS. As an example, we report in Fig. 1 the dependence on the total spin and mass ratio of M th for the three EOSs considered here. More specifically, shown with different coloured symbols are the values of M th -and the corresponding error bars, which are shown in black -as determined following the prescription discussed in Sec. 2 and by Koeppel et al. (2019) (see also the AM for an alternative approach leading to the same results), while the gray-shaded surfaces represent the best fitting function to the data. Note that the fits for the three EOSs have variable errors, but all with a small chi-squared of X 2 BHBΛΦ = 0.001, X 2 DD2 = 0.001, and X 2 TNTYST = 0.003, all of which have an average (maximum) deviation from the fit that is 1% (2%). When expressed in absolute terms, the average deviations from the fits amount to ∆M th = 0.03 M . Quasi-Universal Behaviour. As can be readily appreciated from the inspection of Fig. 1 -which correspond to 40 distinct binary configurations differing in mass ratio and spin -M th shows a behaviour that is similar for the three EOSs, but also that it leads to systematically different values for each of the EOSs considered. Determining accurately the threshold mass for each configuration has required the calculation of the inspiral and merger of about a dozen simulations with varying initial mass, the computational cost associated with Fig. 1 is of about 360 simulations. Extending this work to an arbitrarily large number of EOSs is computationally prohibitive.
However, we can exploit the existence of a quasi-universal behaviour of M th (Bauswein et al. 2013;Koeppel et al. 2019;Agathos et al. 2019;Bauswein et al. 2021). More specifically, we extend the quasi-universal relation derived by Koeppel et al. (2019) by proposing that the functional dependence of M th (EOS, q, χ) can be split into a part that is dependent on the EOS and a part dependent on q and χ. From a mathematical point of view, this essentially amounts to the separability in the functional dependence and hence in adopting the following ansatz where the dependence on the EOS is expressed via a multiplicative function following the study of Koeppel et al.
where C TOV is the compactness of the nonrotating stellar configuration with the maximum mass and R TOV its radius, i.e., C TOV := M TOV /R TOV . The coefficients in expression (2) have been reported by Koeppel et al. (2019) and are a = 2b/(2 − c), b = 1.01, c = 1.34. In practice, expression (1) proposes that f (q, χ) is a surface that models M th as a function of q and χ independently of the EOS. This surface can then be rescaled via a function describing the EOS, κ(EOS) = κ(R TOV , M TOV ), depending uniquely on the stellar compactness for the maximummass nonrotating configuration. Stated differently, the quasiuniversal expression for M th is expressed asM th := M th /κ(EOS), where all the dependence on the mass ratio and spin of the binary is contained in the function f (q, χ), whose behaviour remains to be determined.
Considering the nonlinear behaviour shown by the three fitting functions in Fig. 1, it is natural to express the ansatz for f (q, χ) via a second-order polynomial, i.e., Since we wish to recover the quasi-universal fit obtained by Koeppel et al. (2019), we set a 1 = 1 in (3).
Note that the data clearly shows a non-monotonic growth of M th as the mass asymmetry in the binary changes, which has been discussed already by East et al. (2016) and Bauswein et al. (2021) and more recently by . Although the earlier works were restricted to nonspinning binaries with larger mass ratios, the evidence for a non-monotonic dependence remained tentative and, indeed, the fitting expressions proposed by Bauswein et al. (2021)  monotonic behavior should be expected due to an increase in accretion disk mass.
On the other hand, the presence of additional angular momentum in the system stabilizes the remnant against gravitational collapse (Breu & Rezzolla 2016;Weih et al. 2018). This implies that the fitting function f should monotonically increase with spin. Since the maximum dimensionless spin of uniformly rotating NSs is well constrained to be χ max 2 × 0.65, (see, e.g., the discussion in Most et al. 2020b), we can improve (3) by requiring that M th has a maximum for the maximum allowed value of the dimensionless spin, i.e., Imposing (4) on (3) leads to thus leaving only four independent fitting coefficients, whose values are a 2 = 0.11, a 3 = 0.12, a 4 = 0.07, a 5 = −0.3. We note that we obtain essentially the same coefficients either when fitting all the values ofM th for the same values of q and χ or when averaging the coefficients obtained from the three distinct fits to the different EOSs (see Tab. 3 of the AM for the coefficients of the three fits). Figure 2 reports the fit of the quasi-universal threshold mass for the combined data, showing the accurate represen-tation of the functional behaviour. Indeed, the fit yields a chisquared of X 2 = 0.028, with an average (maximum) deviation from the fit of only ∼ 2% (∼ 6%). The very good match between the data and the fitting function provides strong evidence for the existence of the quasi-universal behaviour conjectured with the separability ansatz (1) and captures nicely the dependence of M th on the mass ratio and spin of the binary. In particular, it highlights that -when compared to the irrotational case -M th increases ∼ 6% for the alignedspin binaries considered here, and decreases ∼ 10% for the anti-aligned spin binaries. We note that while the fit has considered all of the 40 binary configurations, we have verified that predictions for M th made when considering only χ 2 = 0 match equally well the numerical results obtained when χ 2 is nonzero. Lower Limits on the Stellar Radii. Following Koeppel et al. (2019), we can use Eqs. (1) -(3), to set lower limits on the radii of possible stellar models by computing sequences of M th and M TOV for fixed radii, as shown in Fig. 3. More specifically, we recall that once a value for R TOV ,q=1 is fixed, expression (1) selects a line in the (M th , M TOV ) plane. The first intersection of this line with the measured mass of a BNS with total mass M tot that has not collapsed promptly sets a lower limit on M th (Bauswein et al. 2017). This is shown in the upper panels of Fig. 3 [cf., left panel of Fig. 4 of Koeppel et al. (2019)] -one for each of the spins considered here -and where we report with a horizontal blue-dashed line the total gravitational mass estimated for GW170817, M tot,q=1 = 2.74 +0.04 −0.01 M (The LIGO Scientific Collaboration & The Virgo Collaboration 2017). Since the mass ratio is not well-known, the constraint for GW170817 is actually given by a band (blue-shaded area) whose vertical edges depend on the mass ratio q (see arrows in the top-left panel of Fig. 3). Also reported in the top panels of Fig. 3 with a grey-shaded area is the limit set by causality and that requires M TOV /R TOV 0.354.
In essence, therefore, the blue band constrains the redshaded area from above, yielding a lower limit for the radius of the maximum-mass star, R TOV ,q=1 (red solid line). The most important difference with a similar figure presented by Koeppel et al. (2019) for q = 1 and χ = 0, is that we can now exploit the dependence of M th on the mass ratio (and spin) to report the lower limit on R TOV for different values of q (dashed lines). In this way it is possible to appreciate that in the case of anti-aligned spins (e.g., for χ = −0.3), very strong constraints can be put on R TOV ,q=1 as the threshold mass is, in this case, considerably smaller (cf., R TOV ,q=1 = 10.24 km for χ = −0.3 vs R TOV ,q=1 = 9.44 km for χ = 0.3). Overall, the top panels in Fig. 3 highlight how the knowledge of the mass ratio and spin of the binary can be extremely powerful in setting lower limits on the stellar radii.
via the same scaling function f (q, χ) in Eq. (3), namely Because relation (7) should be seen as an ansatz representative of the EOSs used here, it might need additional corrections for EOSs that include strong phase transitions ( . Note that, for a given spin of the binary, the knowledge of the mass ratio can considerably increase the lower limit on the stellar radii, especially for systems with negative total spin. At the same time, given the non-monotonic nature of the dependence on q, equal-mass systems do not necessarily yield the weakest constraint, which is instead attained for q 0.8 − 0.9.

CONCLUSIONS
We have performed the first systematic study of the impact that mass asymmetry and spin have on the threshold mass of BNS systems, M th . We have done so by measuring M th for 40 different BNS configurations encompassing three temperature dependent EOSs, four mass ratios, and a systematic sampling of spin configurations. Using this data we have derived a quasi-universal relation for M th and expressed its dependence on the mass ratio and spin of the binary. The new expression recovers the results of Koeppel et al. (2019) for equal-mass, irrotational binaries, and reveals that M th can increase (decrease) by 5% (10%) for binaries that have spins aligned (anti-aligned) with the orbital angular momentum. In addition, we find evidence for a non-monotonic dependence of M th on the mass asymmetry in the system, which can be explained by an increase in accretion disk mass . Furthermore, we have extended to unequalmasses and spinning binaries the lower limits that can be set on the stellar radii once a neutron-star binary is detected, obtaining generic and analytic expressions for the minimum radii as a function of the mass ratio and of the total spin of the binary. In this way, we have highlighted that the merger of an unequal-mass, rapidly spinning binary can significantly constrain the allowed values of the stellar radii.
With more than 360 binaries simulated, this work has represented a challenge both for the associated computational costs and for the amount of data to be analysed. Furthermore, while it has provided a first sparse but complete survey of the space of parameters to be expected in BNSs, it can be improved in a number of ways. First, by investigating more carefully the role played by spin of the secondary star. Second, by increasing the number of EOSs considered so -especially if they involve a softening from phase transitions (Most et al. 2019a;Weih et al. 2020) -as to further refine the properties of the quasi-universal behaviour. Third, by elucidating the role played by the disruption of the stars at merger. Finally, by exploring the possibility that ultra-strong magnetic fields could modify the stability properties of the merged object and hence the threshold mass. We leave all of these improvements to future works. APPENDIX A. NUMERICAL METHODS AND SETUP Numerical Simulations. As mentioned in the main text, our ID is constructed making use of the recently developed initial-data solver FUKA, which is based on the KADATH spectral solver library, solving the eXtended Conformal Thin Sandwich (XCTS) formulation of Einstein's field equations (Pfeiffer & York 2003;Pfeiffer & York 2005;). The ID is initially constructed using the force-balance equations to determine a quasi-circular orbit. The spin of each NS is achieved by modeling the velocity field of the fluid as a linear combination of a purely irrotational component and a uniformly rotating component (Tacik et al. 2015;Tsokaros et al. 2015;. To minimize the eccentricity of the inspiral, we utilize 3.5 Post-Newtonian order estimates of the expansion factor,ȧ, and the orbital velocity as discussed in . This is important as quasi-circular ID for asymmetric binaries and binaries with spinning objects result in very eccentric inspirals which is also discussed in . Additionally, it has been shown that eccentricity can have an impact on the stability of the remnant (East et al. 2016), therefore, the use of 3.5 PN estimates is important to obtain accurate measurements of M th . Finally, the initial separation of the compact objects are set to 50km, thus allowing the binary to equilibrate over the course of a few orbits prior to merger due to the approximations used in the initial data construction. Realistic, hot EOSs.
The binaries simulated here have been modelled with three different realistic and temperaturedependent (hot) EOSs. In particular, they are the rather soft TNTYST EOS Togashi et al. (2017), the rather stiff BHBΛΦ EOS Banik et al. (2014), and the intermediate (HS-)DD2 EOS (Typel et al. 2010). The most salient properties of these EOSs, namely, the maximum mass of the nonrotating configuration, the corresponding radius, the compactness, and the free-fall timescales are reported in Tab. 2. Spin Configurations. Emphasis in our analysis has been placed on low-spin priors and on the role of the spin of the primary in determining the threshold mass. Hence, the large majority of our binaries has an irrotational secondary (i.e., χ 2 = 0). However, binaries with non-negligible mixed spins cannot be ruled out and, thus, we have performed additional spot tests which include configurations of constant mass ratio and effective spin, but with a spinning secondary so as to ascertain the difference in the collapse time using the BHBΛΦ EOS. All of the 40 binaries considered here and their properties are collected in Tab. 1.
With the data collected in this way we have performed a double analysis. First, we have considered only binaries with irrotational secondaries and performed the fits as discussed in the main text. Next, using the resulting coefficients we have compared the predictions of the fit with the actual numerical results of the four binaries with spinning secondary. In this way, we have found that the analytic predictions and the numerical results differ by less than 3% and are therefore well within the average error of the fit. Overall, this result has given us confidence that the spin of the secondary has a negligible impact on M th , which is instead dynamically dominated by the spin of the primary. This is a reasonable estimate for two distinct reasons. Firstly, M th effectively measures the strength of the gravitational field of the merged object and is therefore dominated the properties of the primary. Secondly, the modifications on M th due to spin are, overall, not large. As a second approach, we have performed a fit using all of the data available and this is what is actually presented in the main text; the corresponding values of the fitting coefficients and their uncertainties are collected in Tab. 3.

B. COMPARISON OF THRESHOLD MASS MEASUREMENTS
To determine M th from our numerical simulations we have elected to utilize the method of measuring the collapse time, t col , as proposed in Koeppel et al. (2019) given the physical motivation behind it. We will reference this as the "freefall", method which is discussed below. However, similar measurements can be made with minimal impact to the result of the measured value of M th by using a simpler approach, which we refer to as the "averaging" method, and which averages the masses of the two binaries closest to the critical value leading to a prompt collapse. Free-Fall Method. The original method of measuring the collapse time, t col , as proposed in Koeppel et al. (2019), is defined by tracking the minimum of the lapse,α := min(α). While this approach is adequate for equal-mass, nonspinning binaries, it leads to potential biases when considering unequal-mass and spinning binaries, since the lapse is not only dependent on the lengthscales of the system, but it is strongly impacted by asymmetries in the system. Therefore, when performing measurements of the collapse time, t col , we analyzeα when normalized by its maximum, which we define asα :=α/max(α).
Because of this rescaling, new brackets must be applied to determine the time of merger, t merge , and the time of blackhole formation, t BH . To do so, we set as the merger time the coordinate time such that t merg :α merge = 0.9 .
(B1) Included is the determined threshold mass, M th ; the corresponding chirp mass, M th ; the ADM mass, M ADM , of each NS at infinite separation; the dimensionless spin of the each NS, χ1 , χ2; and the effective spin of the binary,χ. In all cases, the spin axis is parallel to the orbital rotation axis. M th in configurations with a * were obtained using the averaging method as discussed in the Appendix Materials Similarly, we set as the threshold for black-hole formation and the related coordinate time to be such that the remnant is gravitationally unstable across all EOSs and configurations considered and the formation of a BH is certain. Note that the merged object never collapse faster than the shortest free-fall timescale (t col >= τ ff ) where the free-fall timescale is defined as (Rezzolla & Zanotti 2013) and the shortest τ ff occurs at τ TOV := τ ff (M TOV , R TOV ) and is a physically well-motivated lower limit. We recognise that both of these markers are somewhat arbitrary, but, as done in Koeppel et al. (2019), they reproduce well the behaviour of the binaries when other markers, e.g., the time of the maximum of the gravitational-wave emission and the time of the appearance of a black-hole ringdown, are used. Using these markers, we measure the collapse time as defined as which is the same as the prescription from Koeppel et al. (2019), but here it is calculated based onα and with different bracket values. In Fig. 4 we show examples of howα vs coordinate time looks for all equations with q = 0.7 , χ 1 = 0.3 , χ 2 = 0. Additionally, one panel shows the fit for the three EOSs and the extrapolation for t/τ TOV → 1, as discussed by (Koeppel et al. 2019). Averaging method. As from its name, the averaging method simply averages the masses of the two binaries that are closest to the critical threshold mass, with M sup being the supercritical value and M sub the subcritical one, i.e., To distinguish between the two cases, we define as supercritical any binary whose evolution of the normalised lapse sharply decreases to the lower limitα = 0.1 without a characteristic post-merger local maximum. On the other hand, we define as subcritical any binary whose evolution of the normalised lapse shows a local maximum after merger indicating, therefore, that the merged object has contracted but also expanded at least in one complete oscillation. Obviously, the averaging method is far simpler (it does not require any fitting or extrapolation) but also does not provide any measure of how close M sup and M sub are from M th . Interestingly, we have measured M th using the averaging method against the free-fall method and obtained the same results with an average (maximum) difference relative to the value obtained using the free fall method of < 1% (∼ 2%). The situations in which the averaging method is to be preferred is when the difference between M sup and M sub is < 1%, in which case the error in measuring M th would be dominated by the evolution resolution. This was particularly useful in the case of the TNTYST EOS, where the separation between the supercritical and subcritical solutions was very small in some specific mass ratios and spins.
Both of the methods outlined require a definition on whether or a not a dataset is supercritical or subcritical. In this work we use the definition such thatα| tBH tmerg diverges di-  Horizontal lines mark the threshold values to define merger and black-hole formation, i.e.,αmerge = 0.9 andαBH = 0.1, respectively. Bottom Right: Example of measuring M th using the " free-fall" method for the three EOSs; dashed lines denote a linear-fit result. Table 3. Values of the fitting coefficients and their uncertainty in the functional ansatz (3) in the main text. The first three rows refer to the specific EOSs considered here and the last one to the quasi-universal expression.