The impact of nuclear reactions on the neutron-star g-mode spectrum

Mature neutron stars are expected exhibit gravity g-modes due to stratification caused by varying composition. These modes are affected by nuclear reactions, leading to complex (damped) mode frequencies and the suppression of high order g-modes, in contrast with the common non-dissipative analysis which leads to an infinite g-mode spectrum. Focusing on the transition between the fast and slow reaction regimes, we examine the impact of nuclear reactions on the g-mode spectrum. The general framework for the analysis is presented along with sample numerical results for the BSk21 equation of state and the standard Urca reactions.


INTRODUCTION
Neutron stars are highly compact objects that involve a rich and complex variety of physics.Due to this, there are (more or less) specific classes of oscillation modes associated with each aspect of the physics involved.One key feature is the varying composition of matter, introducing buoyancy as a restoring force in the equations of fluid dynamics (Reisenegger & Goldreich 1992).This, in turn, leads to the presence of low-frequency gravity g-modes.The typical frequency of the leading neutron star g-mode is of the order of a few 100 Hz (depending on the matter equation of state) and the higher overtones lie at lower frequencies.As these modes rely on stable composition stratification for their existence, they are sensitive to nuclear reactions.Reactions will strive to reinstate beta equilibrium in the neutron star matter and hence lead to each g-mode becoming damped.The standard description is really only relevant in the limit of infinitely slow nuclear reactions.While this may apply to the highest frequency g-modes (the lowest overtones) there will always be modes that have low enough frequencies that the impact of finite reactions must be accounted for.Hence it is interesting-at least from a formal point of view-to consider the impact of nuclear reactions on the g-mode spectrum.This is the issue we will address in the following.The aim is to establish a more precise understanding of how the g-modes are affected by reactions and what happens to the mode spectrum when the reaction rate becomes fast compared to the dynamics.This issue may not have immediate relevance for many astrophysical applications, but it connects with problems that involve the very high order g-modes.Two such problems immediately come to mind.First, the nonlinear saturation of modes that are driven unstable by gravitational-wave emission (Schenk et al. 2001).Second, the so-called p-g instability which has been proposed to operate in neutron star binaries (Weinberg et al. 2013) and which involves the coupling of very high order pressure p-modes and g-modes to the dynamical tide.
The importance of stellar seismology had been appreciated since the work of Cowling (1941) (see Cox 1980;Unno et al. 1989).While originally the focus was on the mathematical formulation of the problem, recent X-ray timing observations and the detection of gravitational waves from the GW170817 neutron star merger event have led to renewed focus on the problem (Abbott et al. 2017(Abbott et al. , 2018)).The vast majority of the literature on neutron star seismology has focused on oscillation modes that depend weakly on the precise nuclear physics below the crust, such as the fundamental f-mode or the pressure p-modes.Due to the complexity of the problem, different aspects of the physics have been studied one piece at a time.
A considerable body of work has also been dedicated to the neutron star g-modes (see Finn 1986;McDermott et al. 1983;Miniutti et al. 2003, to name a few).This is natural, as these modes should be present in both mature and hot young neutron stars (noting recent evidence that the g-modes may be excited during the proto-neutron star stage following a core-collapse supernova; Vartanyan et al. 2023).In general, the g-modes depend on both the internal matter composition and the state of matter (e.g. the presence of neutron superfluidity).Unlike other oscillation modes such as the f-and p-modes, the g-modes rely specifically on the conditions in the outer core of a neutron star, below the elastic crust.It is therefore of interest to ask what constraints on the nuclear physics could be made from stellar observations of g-modes.Specifically, for a mature neutron star the mode frequencies depend on the variation of the proton fraction with density, in turn subject to nuclear parameters like the nuclear symmetry energy.An observation of a specific g-mode frequency-perhaps as a tidal resonance induced during binary inspiral (Andersson & Ho 2018;Ho & Andersson 2023)-might then provide some insight into the nuclear physics and the equation of state (beyond bulk properties like the mass and radius of the star).
Building on previous proof-of-principle work by Andersson & Pnigouras (2019), this paper considers the spectrum of g-modes in a cold mature neutron star and the impact nuclear reactions have on the damping of oscillations.Instead of focusing on either the fast or slow reaction regimes, we consider how the modes behave for arbitrary (parameterised) reaction rates and examine the effect damping has on the crossover from slow to fast reactions.This provides a better understanding of the phenomenology of the problem and a clearer idea of what happens to the very high order g-modes in a realistic neutron star model.Our demonstrations are based on a specific equation of state, BSk21 (Fantina et al. 2013;Potekhin et al. 2013), which allows us to employ a realistic description of the matter stratification.
The layout of the paper is as follows: in Section 2, we outline the background equations and physics that go into the problem.Section 3 examines the fast and slow reaction regimes using a plane-wave approach.Section 4 then solves the mode-problem for a general reaction rate, presenting the dimensionless equations and results.Finally, Section 5 summarises the work and presents ideas for future continuation of this effort.

THE PERTURBATION PROBLEM
Even though we are ultimately aiming for realism, here we will explore the impact of nuclear reactions on the g-modes of a stratified neutron star in the context of Newtonian gravity.This makes the results somewhat phenomenological, but as we will have to vastly exaggerate the reaction rates to demonstrate the features we are interested in this lack of precision is not a great concern.Given this context, we consider the problem for non-rotating stars, assuming that the oscillation modes-with label n and frequency ωn-are associated with a polar perturbation displacement vector (expressed in the coordinate basis associated with the spherical polar coordinates [r, θ, φ]) with where the multipole amplitudes, W l and V l , are functions of r only and Y m l (θ, φ) are the usual spherical harmonics.Along with this, all scalar perturbations are expanded in spherical harmonics.That is, using δ to indicate an Eulerian perturbation, we have the perturbed mass density with and similar for all other scalar quantities.In the following, whenever p, ρ, Φ are used without δ they refer to the value of the pressure, density and gravitational potential, respectively, of the background equilibrium star.As we are ignoring rotation, the background configuration is spherical so all associated quantities are only functions of r.
Turning to the perturbed Euler equation, we have This leads to the radial Euler component, and then from the φ component of the Euler equation we get, We also need the perturbed continuity equation, Combining the last two equations, we get Lastly, we have the perturbed Poisson equation, Let us now add nuclear reactions to the problem.As the moving fluid is no longer in equilibrium, we need to consider additional parameters in the (perturbed) equation of state.A natural option, which helps account for nuclear reactions driven by the deviation from beta equilibrium, is to introduce the new variable β = µn − µp − µe depending on the chemical potentials for neutrons, protons and electrons (labelled n, p and e).In (cold) equilibrium, we then have β = 0.This condition allows us to solve for the matter composition, e.g. the proton fraction xp for a given density.For simplicity, we assume pure npe matter and that the star is cold enough to be transparent to neutrinos, thus the relevant reactions will be the Urca reactions.As a consequence of this addition, the general equation of state will be a two parameter function p = p(ρ, xp) where xp is the proton fraction.From Andersson & Pnigouras (2019) we have, where ∆ is the Lagrangian perturbation, and where tR is the characteristic reaction time (and the minus sign is a convention).Note that this means that ∆β = 0 in the limit of very fast reactions, when tR → 0. In the opposite limit, when reactions are slow, we have tR → ∞ and therefore A → 0 and the matter composition is frozen.Moreover, given that the unperturbed star is assumed to be in chemical equilibrium we have ∆β = δβ.Using this and combining ( 11) and ( 8) we get, As it is common to work with the perturbed pressure, we rewrite this equation as and we also need leading to At this point it make sense to use the thermodynamic relation, to get Next we define the speed of sound in equilibrium and at fixed proton fraction, respectively; These quantities are related to the commonly used adiabatic indices via We also introduce the density scale height (which is convenient as we want to avoid involving an explicit stellar model in the plane-wave analysis below) With these definitions, we have When solving the equations numerically later on, we will want to remove δρ l so we need Finally, we define the Brunt-Väisälä frequency as where the local gravitational acceleration is This leads to which is the dimensionless Brunt-Väisälä frequency, and we have Finally, rewriting ( 6) and ( 9) we get The last three equations are the main equations that will be used to determine the neutron star g-modes.

PLANE-WAVE ANALYSIS
In order to gain intuition and help explain the numerical results later, it is useful to consider a local plane-wave analysis.First, we introduce to get the, fairly concise, equations and From these equations it is evident that the problem will change when we consider finite timescale reactions.In the limit of no reactions, when A = 0, we are dealing with an eigenvalue problem for ω 2 n , so we will always have two roots ±ωn.When A ̸ = 0 the eigenvalues become complex and the symmetry of the mode pairs is less obvious.
As we are mainly interested in the qualitative behaviour at this point, we introduce the Cowling approximation (which is expected to be reasonably accurate for the g-modes); setting δΦ l = 0.Then, using (35) to remove δρ l from the problem we get and where the Lamb frequency is defined as In order to explore the nature of the waves we are interested in, we now adopt the plane-wave approach with p = δp l , Ŵ = r 2 Wl , ∂r → ik. (39) This leads to Here we make two simplifying assumptions.First we focus on short-wavelength motion, such that k|H| ≫ 1.Secondly, we assume that N 2 ≪ 1 as appropriate for weakly stratified matter (eventually leading to the anticipated low-frequency g-modes).
It is now easy to see how the expected barotropic result emerges in the N 2 → 0 limit.For fast reactions, we get leading to the dispersion relation This solution represents sound waves-the pressure p-modes in the full mode calculation later.Higher overtone modes have shorter scales (=larger k) and therefore lie at higher frequencies.In a neutron star, we expect to find an infinite set of high-frequency p-modes.
In the opposite limit of slow reactions, we have Now we instead arrive at This equation has two sets of roots.If it is also the case that (effectively focusing of dynamics slower than the sound waves) then In the opposite limit, when it is easy to see that we retain the p-modes from the barotropic case.In essence, the introduction of the stratification has added a set of low-frequency modes to the spectrum.These are the g-modes.It is easy to see that, as the wavelength decreases (=larger k) the frequency decreases.This agrees with the results of Unno et al. (1989) where in beta equilibrium, the g-mode frequency was shown to tend to In a neutron star, with stable stratification, we expect to find an infinite set of undamped, low-frequency g-modes.This changes when we consider the nuclear reactions.For finite reaction rates, we have (focussing on the ω 2 n ≪ L 2 l case) That is, or As expected, the reactions lead to complex-frequency (damped) oscillations.This is as it should be, given that the reactions lead to bulk viscosity which damps the fluid motion Schmitt & Shternin (2018).It is, however, easy to see that we retain the previous (undamped!)results in the fast/slow reaction limits.
In the general case, solving for the frequency we have This is the key result, showing how we retain the undamped modes in the slow-reaction limit.Taylor expanding for large ω0tR we see that That is, the two modes from the non-reactive problem are both damped and symmetric with respect to the imaginary axis.
The numerically determined modes retain this symmetry.Similarly, it is easy to see that both roots become purely imaginary in the fast-reaction limit.The two roots limit to ωn = 0 and i/tR (which means that ωn → i∞ as tR → 0), respectively.In fact, it is easy to see that there are no oscillatory modes below the critical reaction timescale tR = 1/2ω0.At this point, the pair of modes from (57) merge on the imaginary axis.Above the critical reaction time, they split again and one mode moves towards the origin while the other moves towards +i∞ as tR decreases.In this regime the mode solutions represent pure diffusion.Above the critical timescale, for larger values of tR, we have damped oscillations.The implications of this behaviour are-at least formally-important.While the classic analysis suggests that the g-mode spectrum is infinite, for a realistic neutron star model this cannot be so.The (potentially very) high overtones, for which ω0tR is small, will be overdamped.This accords with the results from Andersson & Pnigouras (2019).Finally, from (56) we see that, in the regime where the modes are oscillatory, we have ω 2 n ≈ ω 2 0 .In essence, we expect the modes to move along a quarter circle, from the real axis (when tR = ∞) to the imaginary axis (at the critical reaction time).This prediction is testable with numerical solutions and we now turn to that problem.

Dimensionless Form
In order to solve the mode equations numerically, it is convenient to work with a dimensionless formulation (Unno et al. 1989).Relaxing the Cowling approximation, we define the following variables, Inserting these definitions into (29)-( 31) and the perturbed Poisson equation 10 we arrive at the coupled differential equations and The equation for y3 is obtained from the trivial relationship between y3 and y4 and we have the used the same definitions as Unno et al. (1989) where plus we define with M and R being the total mass and radius of the background star, respectively.Using ( 27) we can rewrite the equations as where is a dimensionless radial coordinate.This system of equation can be compared to that used by, for example, Unno et al. (1989) and it is easy to confirm that our set limits to the usual one when A → 0.
Next we need four boundary conditions: two regularity conditions at the centre of the star and two conditions at the surface.Near the origin, as r → 0, a Taylor expansion reveals that the we should have Meanwhile, near the star's surface, as r → R, we need to impose the vanishing of the Lagrangian pressure perturbation and the continuity of the gravitational potential and its derivative.These conditions take the form and These relations follow by first requiring that δΦ l vanishes at infinity and secondly that there exists a distinct boundary of the star where ρ, p ≈ 0, thus ∆p = 0.These are the same boundary conditions as in the usual calculation (Unno et al. 1989).We can also see from this that our equations have no dependence on m which is to be expected due to the spherical symmetry of the background star.Finally, to close the system we need an equation of state.We will discuss our chosen model below.

Results
First of all, in order to gain confidence in the numerical implementation, we reproduce (and check) some results from the literature.The natural choice is to consider a Γ = 2 (n = 1) polytrope with a constant Γ1 representing the stratification Figure 1.Plots of the dimensionless overlap integral Qn on the left and dimensionless frequencies ω on the right vs Γ 1 − Γ on log-log scales for the first 6 g-modes.The background star is polytropic with γ = 2, R = 10 km and M = 1.4M ⊙ .For the overlap integrals the slopes range from 0.95 − 0.98 and for the frequencies the slopes range from 0.46 − 0.47.(Pnigouras 2017).In particular, in work by Xu & Lai (2017), the numerical results show the following relations for the l = 2 g-modes of a non-rotating Newtonian neutron star with a polytropic equation of state where Qn is the dimensionless overlap integral defined as also commonly referred to as the mass-multipole moment of each mode in the literature.
Using the results from our code we obtain the results shown in Figure 1 for the first 6 g-modes, where Γ = 2 and Γ1 is varied from 2.01 − 2.5.This allows us to confirm the suggested scaling relations.Using a linear regression function, for the overlap integrals the slopes range from 0.95 − 0.98 and for the frequencies the slopes range from 0.46 − 0.47.This is in agreement with ( 80)-(81) allowing for small numerical errors.

Realistic stratification
In order to make the model more realistic we introduce stratification motivated by the BSk family of equations of state (Fantina et al. 2013;Potekhin et al. 2013).Specifically, mainly as proof of principle, we provide results for the BSk21 model (see, also, Gittins & Andersson 2023;Andersson & Gittins 2023).First, from BSk21 table, the relationship p(n b ) and N 2 (n b ) are calculated, where n b is the baryon number density.From this, n b is converted to mass density using ρ = m b n b , where m b is baryon mass.The pressure and density profiles as function of radius then follow (as usual) from the Newtonian stellar structure equations.Having fixed the equation of state, we still need to pick a sample neutron star.The particular model we consider has central density ρ = 5.10 × 10 15 g cm −3 , which leads to a radius of R = 13.49km and a mass of M = 1.43M⊙.Equations ( 71)-( 74) along with the boundary conditions ( 76)-( 79) were solved numerically using this background star and for different reaction rates by varying A.
From the behaviour of the radial eigenfunction, W l (r), and the mode-frequency we determine if a mode is a p-or g-mode.Gravity modes tend to have small ωn and a W l (r) that shows prominent features well beneath the surface of the star.Meanwhile, p-modes will have higher ωn and a W l (r) that shows prominent features close to the surface.The fundamental f-mode can be thought of as the p-mode with the lowest frequency.One can also easily distinguish between overtones of the p-and g-modes.
As the overtone n increases, ωn will decrease for g-modes but increase for p-modes (as per the plane-wave discussion).Also by examining W l (r) again, the number of nodes in the diagram is roughly equal to n, allowing us to identify the specific overtones.
For slow reaction rates, tr > 1s, the modes exhibit minuscule damping, Re(ω) ≫ Im(ω) ≈ 0. The expected modes are found; the f-mode and the g-and p-modes with their respective overtones.When tr < 1s, the damping begins to increase for all modes.As expected, the effect is greater for the g-modes due to their relatively low frequencies in comparison to the f-and p-modes.As can be seen in Figure 2, as |A| increases, which corresponds to tr − → 0, the real part of the mode frequency decreases while the imaginary part increases.From Figure 3 we see that the g-modes tend to move in unison with each other until a mode hits the relevant critical reaction rate, at which point the real part of the frequency quickly tends to 0 and the mode becomes "purely" imaginary, Re(ω) ≈ 0.
Beyond the critical reaction rate, a pair of pure imaginary modes are found for each specific g-mode, gn.As the reaction rate becomes faster, the pair diverge, both staying on the imaginary axis but one increasing and the other tending towards the origin.This agrees with the behaviour seen in the plane-wave analysis of the fast reaction limit.The existence of the second mode, ω2, is due to ω2 = −Re(ω1) + iIm(ω1) being a solution if ω1 is a solution with positive real and imaginary components.The presence of two solutions is down to the actual eigenvalue being ω2 n .We can think of the two solutions physically as modes propagating in opposite directions around the star, which due to the lack of rotation are equivalent.This obviously does not impact the mode damping, hence why the imaginary part of the frequency is the same for both solutions.
As expected, the lowest frequency g-modes are the first to hit their critical reaction rate and reach the imaginary axis.This is akin to the "switch off" predicted in the previous discussion.In agreement with the plane wave argument, we also see that the frequency at which the g-modes first hit the imaginary axis is approximately the imaginary counterpoint of the real-valued g-mode frequencies in the slow reaction limit.Sample numerical data is shown in Table 1, showing good agreement with the predicted behaviour in the fast reaction limit.We also see that, as the g-modes sweep through the complex plane they appear to trace a circular path, with centre at the origin and radius of the initial slow-reaction limit mode frequency.The solution with a negative real part traces a symmetric pattern in the second quadrant of the complex plane.Once the g-modes hit the imaginary axis, the symmetry breaks and the two mode solutions diverge along the imaginary axis, as previously discussed.Figure 3. Results showing how the dimensionless frequencies ωn of the first four g-modes (g 1 , g 2 , g 3 , g 4 ) move in the complex plane when the reaction rate is varied.The diamonds represent the numerical values and the dashed black lines represent the theoretical prediction of a circle with a radius equal to the numerical value of the g-mode frequency when A ≈ 0.
Table 1.In the table are the dimensionless frequencies of the first 10 g-modes, calculated numerically for A = 0 and for when the mode frequency first becomes purely imaginary which we will call the critical frequency, A Crit .The ratio is defined as

CONCLUSIONS
Expanding on the work of Andersson & Pnigouras (2019) we have shown how the presence of nuclear reactions in a neutron star leads to a damping of the composition g-modes for arbitrary reaction rates.This was achieved by setting up linear perturbation equations in a dimensionless formalism similar to that of Unno et al. (1989) and others.We showed that, as the reaction rate increases, the mode frequencies sweep through the complex plane from the real axis to the imaginary axis.As a consequence, the higher order modes are the first to be removed from the oscillation spectrum.This continues until, at a certain reaction rate, there would be no oscillatory g-modes left in the neutron star.The good agreement between our numerical results and those derived from our plane-wave analysis provides strong confidence in these conclusions.These results are, however, currently mainly relevant in a phenomenological sense due to the simplifications made in the calculations, such as using Newtonian gravity.Still, we learn that we need to be careful with problems where the higher order g-modes play a role, such as the p-g instability and mode excitation during a dynamical tide, due to the impact reactions have on the mode spectrum.
Following on from these results, the next step would be to consider more realistic background stars by extending this calculation to general relativity.One could also consider the impact temperature profiles and gradients would have on the g-modes as one of the assumptions made was that we were using a mature, cold neutron star.For example, work by Krüger et al. (2015) shows that the g-mode frequencies and damping times can be dramatically affected by high temperature.Finite temperatures also bring in a new class of g-modes, thermal g-modes, in addition to the composition modes we have considered here.This future direction would open the way to considering objects like proto-neutron stars formed after the gravitational collapse of successful core-collapse supernova explosions (Ferrari et al. 2003;Vartanyan et al. 2023), systems where nuclear reactions will also have an impact on the dynamics.

Figure 2 .
Figure 2. Plots showing how the dimensionless mode frequency of the fundamental g-mode, g 1 , varies with A. Plotted separately are the real and imaginary parts of the frequency on the left and right, respectively.
) .As we can see this agrees well with the predicted value of 2