Universal relations for gravitational-wave asteroseismology of proto-neutron stars
Alejandro Torres-Forné Pablo Cerdá-Durán Martin Obergaulinger Bernhard Müller José A. Font









Max Planck Institute for Gravitational Physics 
	(Albert Einstein Institute), D-14476 Potsdam-Golm, Germany




Departamento de
	Astronomía y Astrofísica, Universitat de València,
	Dr. Moliner 50, 46100, Burjassot, Spain




Departamento de
	Astronomía y Astrofísica, Universitat de València,
	Dr. Moliner 50, 46100, Burjassot, Spain
Institut fur Kernphysik, Theoriezentrum, Schlossgartenstr. 2, D-64289 Darmstadt, Germany



Monash Centre for Astrophysics, School of Physics and Astronomy, Monash University, VIC 3800, Australia



Departamento de
 Astronomía y Astrofísica, Universitat de València,
 Dr. Moliner 50, 46100, Burjassot, Spain
Observatori Astronòmic, Universitat de València, 
 José Beltrán 2, 46980, Paterna, Spain


State-of-the-art numerical simulations of core-collapse supernovae reveal that the main source of gravitational waves is the excitation
of proto-neutron star modes during post-bounce evolution. In this work we derive universal relations that relate the frequencies of the most common oscillation modes observed, i.e. g-modes, p-modes and the f-mode, with fundamental properties of the system, such as the surface gravity of the proto-neutron star or the mean density in the region enclosed by the shock. These relations are independent of the equation of state, the neutrino treatment, and the progenitor mass and hence can be used to build methods to infer proto-neutron star properties from gravitational-wave observations alone. We outline how these measurements could be done and the constraints that could be placed on the proto-neutron star properties. 




=


04.25.D-, 
04.40.Dg, 
95.30.Lz, 
97.60.Jd  


Universal relations for gravitational-wave asteroseismology of proto-neutron stars
    José A. Font
    December 30, 2023
==================================================================================






Core-collapse supernova (CCSN) explosions are a promising source of gravitational waves (GW) and might be one of the next
discoveries of current or future ground-based GW observatories. The most common CCSN type, a neutrino-driven explosion, is expected to be observable with Advance LIGO and Virgo within our galaxy <cit.> at a rate of about three per century <cit.>. CCSN events mark the end of the life of massive stars (8-100 M_⊙) following the formation of a heavy iron core that collapses under its own gravity. The end-result is the formation of a proto-neutron star (PNS) at densities above nuclear-matter density with a radius of ∼30 km and a stalled accretion shock at ∼100 km. This situation is maintained for 0.1-2 s, while accretion proceeds onto the PNS as it cools down through neutrino emission. The shock-PNS system may suffer instabilities, in particular convection and the standing accretion shock instability (SASI). Those instabilities are crucial for shock revival, as they allow for an enhanced energy deposition in the post-shock region by the neutrinos streaming out of the PNS. At the same time these instabilities break spherical symmetry and induce perturbations of the PNS that produce a rich, stochastic GW signal that lasts until the onset of explosion, if successful, or up to the formation of a black hole. Detailed descriptions of the CCSN explosion mechanism and GW emission can be found, e.g., in <cit.> and <cit.>, respectively. 

Numerical simulations have shown that the GW signal, albeit highly stochastic, displays some clear trends in the time-frequency plane (spectrograms) as in the example shown in the left panel of Fig <ref>. This is the result of the excitation of PNS g-modes and, in some cases, of SASI modes <cit.>. Recent work <cit.> has established that the features observed in the GW spectrograms can be very accurately matched to the l=2 PNS eigenmodes, l being the order of the spherical-harmonic decomposition. Those are obtained from the solution of an eigenvalue problem at each time-step of the simulation, using angular-averaged profiles of the simulation data as a background. Particularly important for an accurate matching has been the inclusion of general relativity and space-time perturbations <cit.>. This analysis enables a detailed study of the behaviour of CCSN GW signals in spectrograms without the need of performing computationally challenging multidimensional simulations. Instead, affordable spherically-symmetric simulations serve the purpose, supplemented with the computation of the eigenfrequencies for l=2 perturbations, responsible for the GW emission.
Those simulations, however, do not allow to compute the amplitude of the waveform but only its frequency evolution. In this Letter we show that at each instant during the PNS evolution, the characteristic frequency for the different modes does not depend on 
the exact structure of the PNS but can be estimated from the general properties of the remnant. Furthermore, the relations we derive are universal, as they do not depend on the equation of state (EoS), progenitor star, or neutrino treatment. Therefore, they provide the potential to be used as a basis for parameter inference algorithms once GW observations from CCSN become available. 










Our analysis is based on 25 1D simulations with different combinations of numerical codes, gravity approximations, EoS, and progenitor stars. The simulations are performed with the  code <cit.> and the  code <cit.>. Both are multidimensional Godunov-based Eulerian hydrodynamics codes for spherical polar coordinates and include neutrino-transport.Both codes include multigroup treatment for three neutrino species (electron neutrinos and antineutrinos, and heavy flavor neutrinos).  uses an algebraic Eddington-factor method with an M_1 closure <cit.>, whereas  uses a less sophisticated approach <cit.> based on a stationary transport solution with only a one-moment closure. In both cases neutrino interactions include charged-current reactions of electron neutrinos and antineutrinos with nucleons and nuclei, neutrino scattering off nucleons and nuclei, and neutrino production by nucleon bremsstrahlung.  also includes electron-positron pair processes and inelastic scattering off electrons. Regarding the gravity treatment,  can use either Newtonian gravity or a pseudo-relativistic potential (TOV-A in <cit.>) while  uses general relativity in the XCFC formulation <cit.>. 
We use 6 different finite-temperature EoS in our simulations: LS220 <cit.>, GShen-NL3 <cit.>, HShen <cit.>, SFHo<cit.>, BHB-Λ <cit.> and HShen-Λ <cit.>, the latter two including Λ hyperons. However, hyperons are not included in the neutrino interaction rates. As initial data we use a selection of pre-supernova models from <cit.> including solar-metallicity models between 11 and 75 M_⊙ and one 20 M_⊙ model with 10^-4 solar metallicity. The complete list of simulations can be found in the Supplementary Material.
For each simulation and at each time after bounce, the eigenmode frequencies of the region including the PNS and the shock are determined using the code  <cit.> which includes corrections for space-time perturbations of the lapse and conformal factor in general relativity. This eigenmode calculation does not provide a classification of the modes by itself. We use the 
procedure described in <cit.> to classify modes in f-modes (^lf), p-modes (^lp_n) and g-modes (^lg_n), where 
n indicates the number of radial nodes. We restrict ourselves to l=2 modes, which are the dominant ones for GW emission.
Compared to our mode classification presented in <cit.>, we here re-label some previously misclassified modes: one of the p-modes is reclassified as ^2g_1, for all other g-modes n is increased by one (^2g_n → ^2g_n+1), and for all p-modes above the reclassified p-mode n is decreased by one (^2p_n+1→ ^2p_n). The new and wider set of simulations used in this work (<cit.> only considered 2 simulations) shows that the lowest order g-mode had been misclassified as a p-mode. This is supported by three aspects of the simulations: the energy density distribution of the mode, much more concentrated in the PNS interior, the time evolution of the frequency, which often presents crossings with other p-modes, and the behaviour of the universal relations found in this work (see below). 











We focus on the f-mode and on the lowest-order g- and p-modes, which have been shown to be the only ones observable in the GW spectrum <cit.>. We obtain the eigenfrequency of each mode considering all evolution times and all simulations as a single dataset and try to find relations between this frequency and the properties of the system, namely the PNS mass and radius 
(M_ PNS and R_ PNS), the shock radius (R_ shock), the total mass inside the shock (M_ shock), the
central density and pressure (ρ_ c and P_ c), as well as different thermodynamical quantities 
at different radii. To find the best possible relations we systematically perform fits of the eigenfrequencies with polynomials of the form 

f=a + b x + c x^2 + d x^3,

with x=A^α B^β C^γ D^δ and where A, ⋯, D are all different combinations of the quantities defining the system
and α, ⋯,δ are exponents ranging in the interval [-3,3] in steps of 0.5. Universal relations for each of the modes are built
using the best-fitting combinations together with our intuition on the physical processes that should determine the frequencies. The results for the fits and the combination used for x are presented in Table <ref>. 
For most of the fits we have considered a=0, i.e. imposed that the frequency goes to zero as x vanishes. d 0 was considered only in one case since in all others this coefficient did not improve the fit significantly.

Fig. <ref> shows the fits for the universal relations. In each case we observe a scatter of the data around the best fit, that can be characterised by a standard deviation σ (given in Table <ref>). This scatter does not depend systematically on the EoS or on the progenitor used for the simulation. The only systematic behaviour found is caused by employing two different codes (indicated with different symbols in the plots) which leads to slightly different relations at the upper frequency part of the fits. This systematic difference may be due to the different gravity treatment in either code (better in ) or to the different neutrino treatment (better in ). Although the errors in the fits are fairly small, the observed systematics indicates that they could be decreased even further by performing more complete numerical simulations, which might be necessary for the purpose of inference (as described below). 
 
The lowest order g-modes, ^2g_1 and ^2g_2, depend primarily on the surface gravity of the PNS, i.e. M_ PNS/R^2_ PNS. The reason is because buoyancy inside the PNS is responsible for the excitation of these modes and the phase velocity of the associated (gravity) waves depends directly on the surface gravity. Measuring the frequency of these modes gives an idea of the compactness of the PNS. 
 
The f-mode and p-modes depend on the square root of the mean density inside the shock, i.e. √(M_ shock/R_ shock^3), as their frequency is primarily determined by the local sound speed and the size of the region containing the modes. The frequency of these modes track primarily the location of the shock. The general behaviour matches the dependence seen in previous work <cit.> and is expected according to general theory of modes in stars <cit.>. We note that the ^2g_1 mode, improperly labeled as a p-mode in <cit.>, behaves as a g-mode, since its frequency depends on the PNS surface gravity. Attempts to fit it with a p-mode dependence result in very poor fits. This is further evidence that, despite of its high frequency, this mode is indeed a g-mode.

The ^2g_3 mode presents a significantly different behaviour to other g-modes. Our analysis shows that, despite of being classifies as a g-mode, it cannot 
be fitted using the PNS surface gravity (R^2=0.76). Instead, in the best fit (see table <ref>), the dominant behaviour is related to the
mean density inside the shock (similarly to p-modes) but additional corrections including central density and pressure are necessary. Without this additional correction the quality of the 
fit degradates significantly (R=0.85).
The term p_c/ρ_c^2.5 is very close to p_c/ρ_c^Γ_1, which is related to the entropy (here Γ_1 is the adiabatic index). However, neither using the entropy itself nor Γ_1 instead of 2.5,
improves the fit. Note that the frequency of this mode decreases with increasing x, since b<0. In fact the time evolution of the frequency is monotonically decreasing in all our simulations. The puzzling behaviour of this mode may require special attention, mainly because it has been observed in the GW spectrograms of several numerical simulations <cit.>. 

 
From these results we conclude that it is possible to derive relations that connect the frequencies of the PNS modes, observable in the GW signal, with intrinsic properties of the PNS. These relations are universal, within some uncertainties, since they do not depend strongly on the EoS, neutrino treatment, or progenitor star. Therefore, they could be used for parameter estimation when GW observations
 of CCSN are accomplished. This idea was first proposed by <cit.> for cold neutron stars but with the newly proposed universal relations we show its feasibility in the CCSN context.
 The procedure that could be used for inference is outline next. First, the traces of the different modes have to be extracted from the spectrograms. The dominant emission mode, as observed in all previous work <cit.>, would correspond to the ^2g_2 mode. If present, the pattern of the ^2g_3 mode is also easy to recognise, as it monotonically decreases with time. Finally, in the presence of the SASI, the f-mode can be easily identified due to its low frequency. Other modes that may appear will likely have lower amplitudes and could be classified as higher or lower overtones of those already identified. Each identified mode follows a track in the time-frequency plane, f(t), with an error associated with the measurement process of a GW signal buried in detector noise. Each of those functions f(t) can be transformed into a x(t) using the adequate universal relations f(x). The error for x(t) will depend on both the error for f(t) and the error associated with the universal relation f(x). As a result, this method will allow to measure 
 M_ PNS/R_ PNS^2 from the dominant g-modes and √(M_ shock/R_ shock^3) from the f-mode, with their uncertainties. The development of an inference pipeline based on this proposal is one of our main immediate goals and will be discussed elsewhere.
 

 


Instead of providing an example of inference, i.e. the solution of the inverse problem, out of the scope of this work, we present here an example of the direct problem. We use the 2D CCSN simulation s20 of <cit.> and using only the values of the M_ PNS/R_ PNS^2 and √(M_ shock/R_ shock^3) extracted from the simulation, we compute the predicted frequencies
 for the dominant modes. The right panel of Fig. <ref> compares those frequencies with the GW spectrogram of the simulation and shows that both agree within the errors of the universal relations. Note that data from the s20 simulation was not included in the models used to build the universal relations, which are all 1D. This is an indication that our relations are valid even when applied to multidimensional models and therefore are truly universal. 
 



The relations presented in this work may eventually allow to measure physical properties of a PNS during its first second of life, out of purely GW information. The measurement of the PNS surface gravity gives information about the EoS of nuclear matter at finite temperature, which has not been probed yet in other astrophysical scenarios (let alone at the laboratory). The only similar scenario is the post-merger evolution of binary neutron stars but the typical frequencies are somewhat higher and thus less accessible with current GW detectors <cit.>. Measuring M_ PNS/R_ PNS^2 also allows to constraint the PNS radius. The minimum mass of the iron core at the onset of collapse is set by the Chandrasekhar mass <cit.> for a gas of degenerate electrons, M_ Ch0=5.83 Y_e^2, where Y_e ≃ 0.46 is the electron fraction for iron, which results in a minimum mass of 1.2 M_⊙. The maximum mass depends on finite temperature effects <cit.>. Stellar evolution models show that it can be as high as ∼ 2.5 M_⊙, although for solar-metallicity models it does not exceed 2 M_⊙ <cit.>. The iron core collapses into the PNS in a typical timescale of a few 100 ms. After this the mass accretion rate drops significantly due to the lower density of the outer layers. Therefore, at about 0.5 s after bounce the mass of the PNS is very close to the mass of the iron core and is in the range 1.2-2.5 M_⊙. Similar arguments hold for M_ shock, which is approximately equal to M_ PNS, because of the small amount of mass between the PNS and the shock. This allows to transform the universal relations into constraints for the shock and PNS radii. Fig. <ref> shows an example of how to constraint the PNS radius by measuring the frequency of the dominant ^2g_2 mode. For the s20 model the frequency at 0.5 s post-bounce is about 1 kHz, which would correspond to a measurement of the PNS radius in the 2σ confidence interval 19 - 37 km. The PNS radius computed from the actual simulation is 30 km, which falls within this interval. Similarly, a measurement of the f-mode frequency of about 300 Hz places a constraint on the shock radius in the 2σ interval 50-100 km, at 0.5 s post-bounce to be compared with the value of 70 km obtained from the simulation. In an actual CCSN event there could be additional constraints for the mass of the iron core, which would reduce further the errors in the estimation, e.g. constraints from the neutrino luminosity and spectrum, from the observation of the progenitor star and from the direct observation of the compact remanent, possibly as a pulsar, decades after the explosion.













We thank Marie Anne Bizouard and Sanjay Reddy for useful discussions. Work supported by the Spanish MINECO (grant AYA2015-66899-C2-1-P), by the Generalitat Valenciana (PROMETEOII-2014-069) and by the European Gravitational Observatory (EGO-DIR-51-2017). PCD acknowledges the support from the Ramon y Cajal program of the Spanish MINECO (RYC-2015-19074). JAF acknowledges support from the European Union's Horizon 2020 RISE programme H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740. 
MO acknowledges support from the European Research Council under grant EUROPIUM-667912, and from the Deutsche Forschungsgemeinschaft through Sonderforschungsbereich SFB 1245 "Nuclei: From fundamental interactions to structure and stars". BM acknowledges support by ARC grant FT160100035.


apsrev4-1







